Six and four V. emarginatum samples from Japan and Taiwan, respectively, were analyzed (Supplementary Table S1; Figures 1A–C). Leaf samples from the Amami population were obtained from ex situ stocks collected from the native habitat on Amami Island and grown in Toyama Botanic Gardens, Kyoto Botanical Garden, and private nurseries on Amami Island. At all these facilities, all plants were originated from separate wild individuals on Amami Island; the plants have not undergone a generational change since collection, and hence retain their genetic information as it was at the time of collection. For E. platyphyllum, leaf samples were collected from plants growing at least 2 m apart from each other, in Iriomote Island and Taiwan (Figures 1A,D,E), respectively. In Iriomote Island, we took 43 samples throughout the distribution area; in Taiwan, 18 and 22 samples were collected from two growing areas (Wulai and Taiping Shan), respectively (Supplementary Table S2). The leaf samples of both species were preserved in silica gel until DNA extraction. Total genomic DNA was extracted from dried leaf material using a CTAB method (Murray and Thompson, 1980). DNA concentration was measured using the Qubit® dsDNA BR Assay Kit (Invitrogen, MA, United States) and adjusted to 15 ng/μL with TE buffer.

Elatostema platyphyllum reproduces vegetatively through the rhizome, therefore, visual identification of individuals is difficult. To identify the clonal samples of E. platyphyllum, all samples of this species were examined using MIG-seq, based on the original protocol with minor modifications (Suyama et al., 2022). A sequencing library of 350–800-bp sequences was diluted to 12 pM and sequenced on an Illumina MiSeq Sequencer (Illumina) using reagents from the MiSeq Reagent Kit v3 (150 cycles, Illumina). Sequencing reads were preprocessed using Trimmomatic version 0.32 (Bolger et al., 2014) to cut two bases off the start and end of reads which are below a threshold quality and align all read lengths to 76 bp. Then, reads with the quality below 30 were removed. Subsequent SNP detection was conducted using the Stacks software package version 2.4 (Rochette et al., 2019). We set the minimum number of identical reads required to create a stack (−m 3), the mismatches allowed between stacks within individuals (−M 2), and maximum distance allowed to align the secondary reads to the primary stacks (–N 4).

The clonal samples were identified using GenoDive version 3.05 (Meirmans, 2020). To distinguish sequence errors and variations between clones, loci data with low missing values retained by 80–100% of the whole samples were used. Additionally, the likelihood of misidentification was reduced by obtaining two independent data from the same samples and removing sequences that did not match between repetitions. Clone assignment was conducted following the manual of GenoDive version 3.05. The Infinite Allele Model was used for the calculation of genetic distance between samples, and missing data were removed from the data set. The threshold to be considered as the same genet was set to 2. The unit identified as the same clone based on the analysis was considered one genet, and 28 samples were used as representative samples for subsequent analysis.

For phylogenetic and genetic structural analysis, all samples of V. emarginatum and representative samples of E. platyphyllum were used for ddRADseq. Genomic DNA was double digested using the restriction enzymes BglII and EcoRI; the digested DNA fragments and two adapters (BglII adapter and EcoRI adapter) were ligated. Digestion and ligation were performed simultaneously at 37°C for 16 h. The reaction mixture contained 20 ng of genomic DNA, 5 units of BglII (New England BioLabs, Ipswich, MA, United States), 5 units of EcoRI-HF (NEB), 1× NEB buffer2 (NEB), 1× bovine serum albumin (New England BioLabs), 0.2 μM BglII adapter, 0.2 μM EcoRI adapter, 1 mm ATP (Takara, Shiga, Japan), and 300 units of T4 DNA ligase (Enzymatics, Beverly, MA, United States). The reaction product was then purified with AMPure®XP (Beckman Coulter, CA, United States). Purified DNA was amplified by PCR with indexed primers (Sakaguchi et al., 2015) using KOD-Plus-Neo (TOYOBO, Osaka, Japan). PCR product fragments of approximately 320 bp were selected using E-Gel SizeSelect 2% (Life technologies, Waltham, MA, United States). V. emarginatum library was uniquely barcoded and sequenced with 100 bp pair-end reads in two lanes of an Illumina HiSeq2500 (Illumina, San Diego, CA, United States), and that of E. platyphyllum was sequenced with 150 bp pair-end reads in two lanes of an Illumina HiSeqX (Illumina, San Diego, CA, United States) by Macrogen (Seoul, South Korea). Sequencing reads were preprocessed using Trimmomatic version 0.39 (Bolger et al., 2014) to remove adapter sequences and bases with the following parameters: ILLUMINACLIP TruSeq3-PE-2.fa2:30:10 and SLIDINGWINDOW:4:15. Furthermore, all read lengths were aligned to 100 bp for V. emarginatum and 150 bp for E. platyphyllum. The numbers of filtered reads are shown in Supplementary Tables S1, S2.

The locus assembly and single nucleotide polymorphism (SNP) call was performed using denovo_map.pl. from Stacks version 2.53 (Rochette et al., 2019). We set the minimum number of identical reads required to create a stack (−m 3), the mismatches allowed between stacks within individuals (−M 3), and those between stacks among individuals (−n 2). For E. platyphyllum, the phylogenetic analysis of all samples revealed some identical clones as described below. Since having the same clones in a sample set can introduce biases during SNP calls with stacks, SNP calls were performed again using only one representative sample from the identical clone in the sample set.

For phylogenetic analysis, the populations program of the Stacks package was used to extract SNPs and create a variant call format (VCF) file. We set the minimum percentage of individuals in a population required to process a locus (−r 0.8), minimum minor allele frequency (−-min-maf 0.05), and maximum observed heterozygosity required to process a nucleotide site at a locus (−-max-obs-het 0.5). The –write-single-snp option was used to restrict data analysis to only the first SNP per locus to avoid any influence of linkage among SNPs on the following analysis. The resulting VCF files of the SNP extraction were converted to PHYLIP format using PGDSpider version 2.1.1.5 (Lischer and Excoffier, 2012), and the invariant sites were removed with the script ascbias.py.¹ Subsequently, a concatenation-based maximum likelihood tree was obtained for each species by RAxML-NG version 0.9.0 (Kozlov et al., 2019), using the best model of evolution determined by MODELTEST-NG (Darriba et al., 2020). Branch support values were obtained through 1,000 standard bootstrap replicates. The resulting phylogenetic tree was drawn with FigTree version 1.4.4 (Rambaut, 2012).

Additionally, the phylogenetic networks were constructed with SplitsTree version 4.10 (Huson and Bryant, 2006). The output VCF files from the populations program of the Stacks package were converted to NEXUS format using the vcf2phylip.py script.² Subsequently, the networks were constructed using the NEXUS files with default settings in SplitsTree.

For genetic structure analysis, the populations program of the Stacks package was used to extract SNPs and create a VCF file. We set the minimum percentage of individuals in a population required to process a locus (−r 0.8), the minimum number of populations a locus must be present (−p 4 for V. emarginatum and 3 for E. platyphyllum), minimum minor allele frequency (−-min-maf 0.05), and maximum observed heterozygosity required to process a nucleotide site at a locus (−-max-obs-het 0.5). Moreover, we used the –write-single-snp option to restrict data analysis only to the first SNP per locus to avoid any influence of linkage among SNPs on the following analysis. This program also calculated the number of private alleles in each population. The VCF files were converted to structure files using PGDSpider version 2.1.1.5.

Genetic structure analysis of the population was performed using the Bayesian clustering method STRUCTURE version 2.3.4 (Pritchard et al., 2000) with K = 1 to 9, 20 times for each K, with 100,000 burn-ins for 100,000 iterations per run. The STRUCTURE Harvester (Earl and VonHoldt, 2012) was used to infer the optimal K based on the delta K method (Evanno et al., 2005). The result of the run with the highest likelihood LnP(D) among all runs for that K was adopted.

In addition, the principal component analysis (PCA) among accessions was carried out using TASSEL version 5 (Bradbury et al., 2007).

To evaluate genetic diversity at the individual level, we calculated SNP frequency within a single plant, by extracting and outputting SNPs in VCF format using the populations program of the Stacks package with the same options for genetic structure analysis. Then, the Genotype Summary in TASSEL version 5 was used to calculate the percentage of heterozygous SNPs in each individual.

Additionally, ADZE version 1.0 (Szpiech et al., 2008) was used to estimate private allelic richness with a standardized sample size of two, which is the smallest sample size from a single population. In ADZE, the generalized number of private alleles expected in a population is estimated using the rarefaction approach when sample sizes differ across populations (Kalinowski, 2005).

A total of 37,772,894 paired-end reads were obtained from 10 samples, averaging 3,777,289 reads per sample (Supplementary Table S1). We extracted 5,887 SNPs for phylogenetic analysis with RAxML-NG from this dataset. Based on the model select of MODELTEST-NG, the maximum likelihood phylogenetic tree was created using TVM + ASC_LEWIS model. The phylogenetic tree showed that the Japanese and Taiwanese populations formed distinct clades (Figure 2). The Japanese population showed short branches between individuals, indicating small genetic differences within the Japanese population. The phylogenetic network constructed by SplitsTree based on 12,878 SNPs also suggesting the same phylogenetic structure (Supplementary Figure S1).

For genetic structure analysis, we extracted 4,341 SNPs. The STRUCTURE Harvester analysis showed that delta K peaked at K = 2 (Supplementary Figures S2, S3). Both clustering (Figure 3) and the phylogenetic tree showed a clear genetic differentiation between the Japanese and Taiwanese populations. The result of PCA (Supplementary Figure S4) indicated the differentiation within Taiwanese populations in addition to the genetic structure of Japanese and Taiwanese populations.

The individual heterozygosity of the Taiwanese population was 0.202 ± 0.036 on average, while that of the Japanese population was much lower, 0.014 ± 0.003 (t-test, p = 0.000) (Figure 3; Supplementary Table S3). Although the Japanese population had lower heterozygosity, numbers of private allele and private allele richness were larger than Taiwanese samples. The numbers of private alleles in the 4,341 SNPs extracted from Amami Island, Dahang Shan, Mingchih, and Yuanzui Shan samples were 1,188, 727, 569, and 220, respectively. The private allelic richness of each was 0.245 ± 0.006, 0.184 ± 0.005, 0.089 ± 0.004, and 0.083 ± 0.004, respectively.

An average of 397,838 reads per samples was obtained by MIG-seq and clonal identification was conducted using 144 SNPs. The results showed that the Iriomote population was considered to be composed of one genet and the samples from Taiwan consisted of 27 genets (Supplementary Table S2). One sample from each genet was selected as representative and examined using ddRADseq for the following analysis.

A total of 100,705,390 paired-end reads, averaging 3,596,621 reads per sample, were obtained from the representative 28 samples by ddRADseq (Supplementary Table S2). From the dataset, 3,539 SNPs were extracted for phylogenetic analysis with RAxML-NG. Based on the model select of MODELTEST-NG, the maximum likelihood phylogenetic tree was created using SYM + G4 + ASC_LEWIS model. Within the resulting phylogenetic tree (Figure 4), the individual from Iriomote and Taiping Shan were in the same clade, while the Wulai population and others were clearly separated. The divergence between clades was smaller to that within populations, and the differentiation between the Japan and Taiwan populations was not clear. The phylogenetic network constructed by SplitsTree based on 6,867 SNPs also showed the same situation (Supplementary Figure S5).

The estimated optimal K value calculated on the extracted 3,404 SNPs was two (Supplementary Figures S6, S7). This result indicated the differentiation between the two populations in Taiwan, but no clear genetic structure over Japanese and Taiwanese populations (Figure 5). This feature is consistent with the architecture in the phylogenetic tree and the result of PCA (Supplementary Figure S8), where differences between samples were relatively larger than differences between clades.

The individual heterozygosity of the Japanese sample was 0.214 (standard deviation is not shown because all samples were an identical clone and the representative sample was only one for the Japanese samples), and that of the Taiwanese samples averaged 0.247 ± 0.012 (Figure 5; Supplementary Table S3). The values for the Japanese sample were only slightly lower than those for the Taiwanese sample. Of the extracted 3,404 SNPs, the private alleles were not found in Iriomote Island, while 241 and 349 were found in Wulai and Taiping Shan populations, respectively. The private allelic richness of each population was 0.115 ± 0.004 in Iriomote, 0.159 ± 0.003 in Wulai and 0.141 ± 0.003 in Taiping Shan.

Phylogenetic and genetic structure analysis revealed much smaller differentiation within the Japanese population of V. emarginatum than in Taiwan (Figures 2, 3). Furthermore, the frequency of heterozygous sites at the individual level in the Japanese population was extremely lower compared to that of the Taiwanese population (Figure 3). Since small populations with low genetic diversity may reduce the population’s survival abilities (Willi et al., 2006), the Japanese population is in a critical situation because of its low genetic diversity at the population and individual levels.

Although Tsutsumi (2011) revealed similarities between samples from Japan and Taiwan in nuclear Internal transcribed spacer and chloroplast sequences, the present analysis based on ddRADseq indicates clear genetic differentiation between Japanese and Taiwanese populations (Figures 2, 3). Furthermore, of the 4,341 loci extracted, we found that the Amami population had 1,188 private alleles. The private allelic richness was also the highest in the Amami population. This indicates that the differences between the Amami and Taiwanese populations in phylogenetic analysis were not only due to lower heterozygosity of Amami population, but also due to the own genetic variations in the Amami population, suggesting the genetic uniqueness. From a conservation point of view, genetically different intraspecific populations require separate genetic management (Moritz, 1995). Therefore, the Japanese and Taiwanese populations should be managed as separate conservation units.

The heterozygosity of populations and individuals decreases more rapidly over generations in small populations (Frankham et al., 2002); the much smaller heterozygosity in Japanese individuals suggests that the Amami population remained small during repeated generational changes. The species was discovered in 1959 in Amami Island, and the low individual heterozygosity of the sample collected in 1959 (YI-0001) indicates that genetic diversity was already low at the time of discovery. It is unlikely that the species was introduced to Amami Island by humans, since it grows on the tops of the trunks of tall trees growing in primary forests that have not been affected by human activities. The habitat of V. emarginatum on Amami Island is highly pristine, with few naturalized plants. Because of its long generation time, it is considered that the Japanese individuals had lost genetic diversity over a long time before the discovery in 1959. The lower heterozygosity of the Amami population can lead to expression of recessive mutations that reduce fitness (Keller and Waller, 2002; Charlesworth and Willis, 2009), which can have a negative impact on the population’s viability.

Although the Amami Island has been registered a UNESCO World Heritage site since 2021, and its nature is being protected more than ever before. Continuous monitoring of the species’ wild habitat is necessary to prevent further reduction of the number of trees on which the species grows.

Elatostema platyphyllum has characteristics that contrast with V. emarginatum. First, the population of E. platyphyllum on Iriomote Island consisted of a single clone. Second, there were no significant phylogenetic differences between the Japanese and Taiwanese samples (Figure 4) and the STRUCTURE analysis also showed no clear genetic structure between them (Figure 5). Therefore, the Japanese population of E. platyphyllum had lower phylogenetic uniqueness than that of V. emarginatum. Third, only slight differences in heterozygosity were observed between Japanese and Taiwanese individuals. Although sexual reproduction and generational changes among a small number of individuals reduce genetic diversity, the reduction was not observed in the Japanese population despite the small number of genet (a single clone). Therefore, it is considered that the Japanese population originated from a single or a small number of individuals immigrated from Taiwan to Iriomote Island in the relatively near past, and has been maintained by vegetative reproduction or with few generational changes. Silvertown (2008) suggested high frequency of clonality in populations of rare or endangered species, of alien plants, or at the edges of the species’ geographical ranges. The Japanese population meets all these conditions.

The Iriomote population is located in the central part of the island, where indigenous vegetation has not been influenced by human activity, and the population seems derived from natural distribution. However, its history is likely as short as a few generations, which is almost equivalent to transient naturalized plants.