Polygonatum rhizomes were collected from 47 different germination areas in China: Anhui (5), Chongqing (2), Fujian (1), Guangdong (2), Guangxi (4), Guizhou (4), Henan (7), Hubei (2), Hunan (1), Jiangxi (3), Shaanxi (2), Sichuan (2), Yunnan (6), and Zhejiang (6) Provinces. Table 1 shows the information on the preliminary taxonomic classification, the geographic origin, and the phenotypic traits of rhizomes. Approximately 100 rhizomes free of lesions and mildew infection were selected from each region and screened primarily for morphological features such as their thickness, number of nodes, and diameter of tubercle. Fresh tissue samples (0.5 g) were cut out from five rhizomes from each region, rinsed with distilled water, dried with filter paper, and preserved at −80 °C for DNA isolation.

Almost 80 rhizomes from each region with similar characters were selected and cultivated in a PR GAP planting base (Buchang Pharma, Hanzhong, China). The experimental plots were located in Zhongchuanba Village, Wulongdong Town, Lveyang County (106°12′30″E, 33°30′9″N; 897 m above sea level). The area of each plot was 1.2 m × 5.0 m, and the spacing between each plot was 1.2 m. The rhizomes were planted in rows; the distance between the rows and that between rhizomes was 25 cm.

Genomic DNA was isolated according to the CTAB method [38]. The quality of the extracted DNA was verified by a NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA), diluted to 20 ng·μL⁻¹ and stored at −20°C for molecular marker analysis.

A total of 100 ISSR primers, designed at the University of British Columbia, were synthesized by Shanghai Sangon Biological Engineering Technology & Services (Shanghai, China). All primers were screened for their amplification efficiency and 15 primers (Table 2) with high amplification efficiency and reproducibility were selected.

PCR amplification was performed in a 96-well Veriti Thermal Cycler (Applied Biosystems, Foster City, CA, USA) in 20-μL reaction volumes comprised of 10 μL of 2 × Es Taq MasterMix (CWBIO, Shanghai, China), 1.2 μL of 1 μM primer, 1 μL of 20 ng·μL⁻¹ template DNA, and 7.8 μL of ddH₂O (CWBIO, Shanghai, China). The PCR amplification procedure was set as follows: initial denaturation at 94°C for 5 min, followed by 35 cycles of denaturation at 94°C for 45 s, annealing at Ta (screened in advance using gradient temperature PCR, Table 2) for 40 s, and an extension at 72°C for 90 s, and a final extension at 72°C for 7 min. All amplification products were separated by electrophoresis on a 1.5% agarose gel in 1 × TBE buffer, stained with ethidium bromide, and visualized under a Molecular Imager Gel Doc system enclosed in Universal Hood II (Bio-Rad Laboratories, Inc., Hercules, CA, USA).

In total, 85 SCoT primers were tested for their amplification efficiency and stability. Finally, 10 primers (Table 2) with high amplification efficiency and reproducibility were selected.

Reactions were prepared following the same protocol as that used for ISSR markers (Section 2.1). The SCoT PCR consisted of an initial denaturation at 94°C for 5 min, 35 cycles of denaturation at 94°C for 50 s, annealing at 50°C for 1 min, and extension at 72°C for 90 s, and a final extension at 72°C for 7 min. All the amplified products were separated and visualized as described in Section 2.1.

The specific phenotypic traits listed in Table 3 were directly investigated using the method we named VG (visual assessment by a single observation of a group of plants or parts of plants). If more than 95% of plants in each plot showed a specific phenotypic state, we confirmed the trait's state and recorded the data according to the stipulations in Table 3. Common quantitative traits such as plant height (cm), stem diameter (mm), leaf length (cm), and leaf width (mm), which are remarkable and convenient for statistics and analyses, were chosen for direct observations conducted nearly every 20 days after seedling emergence and until the end of the growth period when the aerial parts in each plot withered. Five plants with similar growth were randomly selected through a five-point sampling method and labeled for the observation of morphological indexes. Where appropriate, the plant height was measured with a 5-m-long measuring tape, and Vernier calipers were used for surveying the stem diameter, the leaf length, and the leaf width. For uniformly determining leaf length and leaf width, we measured two leaves located in the middle of the stem of the sampled plants.

2.6.1 Statistics and analysis of ISSR and SCoT data

2.6.2 Statistics and analysis of morphological data

The cluster analysis of the ISSR data classified the plants with close genetic similarities into two main clusters (Fig. 1a) with the similarity coefficient of 74% (I and II). Clade I, which included the samples from the regions of Henan and Shaanxi Provinces, with a similarity coefficient of 80.4% among them, was comprised of three subclades: Ia–Ic (S33). A relatively high similarity coefficient of 94.7% was calculated between S21 and S22. The remaining samples, all of which originated from other regions, were classified into clade II. At a similarity coefficient of 79.8% among them, this clade comprised seven subclades, IIa–IIg. In clade II, the highest similarity coefficient of 95%, was calculated between S36 and S37. According to the preliminary species identification (Table 1), the samples of Polygonatum sibiricum were classified into a group with a similarity coefficient of 80.4%, and samples belonging to P. kingianum and P. cyrtonema were grouped together with a similarity coefficient of 74.6%. Samples from the same province were clustered into subgroups with similarity coefficients ranging from 74.5% to 91%. Cluster analysis (Fig. 1a) indicated that the samples collected from Anhui, Chongqing, Guangdong, Guangxi, Guizhou, Henan, Hubei, Jiangxi, Shaanxi, Sichuan, Yunnan, and Zhejiang were classified into a subgroup with similarity coefficients of 84%, 84%, 91%, 74.5%, 74.5%, 80.2%, 78.8%, 80%, 80.3%, 85.6%, 78.6%, and 77.2% among them, respectively.

The MDS analysis based on genetic distances has been used for plant classification, and to some extent, to validate the results of molecular analyses [42]. The samples from all regions were differentiated into six groups according to their genetic distances (I–VI) (Fig. 1b). The samples with close genetic distances were grouped together, in a fashion slightly different than their geographical distribution. This fact may be because a different distribution makes the genetic information of samples of the same species vary. The samples belonging to Polygonatum sibiricum (I) were spatially isolated in the two-dimensional MDS plots from the other samples. The samples of P. kingianum and P. cyrtonema were divided into five small groups (II–VI). Group VI (S14 and S18) showed a spatial distribution with in −0.716 to −0.612 in dimension 1 and −2.162 to −1.825 in dimension 2. The samples in group V were spatially isolated in–0.632 to 0.459 in dimension 1 and −1.547 to–0.797 in dimension 2. Samples of Group IV indicated a spatial distribution within–0.624 to 0.235 in dimension 1 and −0.274 to 0.442 in dimension 2. Samples in group III were spatially isolated within −1.760 to −1.034 in dimension 1 and −0.268 to 0.618 in dimension 2. The remanent samples (II) revealed a spatial distribution within −1.282 to 0.528 in dimension 1 and 1.000 to 1.697 in dimension 2. These results confirmed that ISSR markers could distinguish PR from other species of Polygonatum, thus corroborating previous reports [28,29].

The similarity coefficient between the samples based on the SCoT data ranged between 67% and 93%. The cluster analysis also classified the samples into two main branches with a similarity coefficient of 67% (I and II; Fig. 1c). Cluster I comprised S6, S7, and S17 with the similarity coefficient of 75.1%; this branch was further divided into two subclades, Ia (S6 and S7) and Ib (S17). The similarity coefficient between S6 and S7 was relatively high, at 89.8%. The samples from other regions were classified into cluster II with a similarity coefficient of 74.5%. Cluster II was divided into nine subclades: IIa–IIi. In cluster II, the highest similarity coefficient of 93% was observed between S27 and S45. According to the preliminary species identification (Table 1), the samples of P. sibiricum were clustered into a group with a similarity coefficient of 78.8%, whereas the samples belonging to P. kingianum were grouped together with a similarity coefficient of 77.8%, and the samples from P. cyrtonema were classified into a group with a similarity coefficient of 67%. The samples from the same province were classified into a subgroup with a similarity coefficient of 67%–90.4%. Cluster analysis (Fig. 1c) indicated that the samples from Anhui, Chongqing, Guangdong, Guangxi, Guizhou, Henan, Hubei, Jiangxi, Shaanxi, Sichuan, Yunnan, and Zhejiang were grouped as sharing a similarity coefficient of 82.8%, 90.4%, 88.6%, 73%, 67%, 79.2%, 71%, 80%, 81.4%, 68.8%, 76.4%, and 79.6%, respectively.

MDS analysis based on the SCoT resolved the samples into six groups (Fig. 1d). Thus, the spatial representation of the data using MDS corroborated the results based on genetic distance among the plants inferred from the SCoT markers.

The two DNA-based markers showed similar genetic capacity for the identification of the PR species. The cluster analyses based on the two types of markers resolved the samples from P. sibiricum into groups with comparable coefficients of similarity (80.4% ISSR and 78.8% SCoT). Similarly, the samples of P. kingianum were resolved in a single clade with a similarity coefficient of 74.6% (ISSR) and 77.8% (SCoT), and P. cyrtonema samples were grouped in a clade with similarity coefficients of 74.6% (ISSR) and 67% (SCoT). Therefore, both the ISSR and SCoT markers are suitable for the identification of PR species. However, the results were not completely unanimous. The ISSR and SCoT-based topologies were incongruent in that the samples of P. sibiricum were resolved within cluster I on the ISSR-based dendrogram or within subclade IIb on the SCoT-based dendrogram. In addition, the similarity coefficients among the samples were slightly different between the two types of markers. These differences may due to the fact the two schemes and binding sites of DNA were different [45].

All the samples were classified into two main populations based on leaf arrangement — those with alternate and those with verticillate leaves (Supplementary Table). The samples collected from 13 regions (S19–S25, S32, S33, S36, S37, S39, and S41) were with verticillate or mostly verticillate leaves. Four of them (S36, S37, S39, and S41) with red leaf sheaths were identified as P. kingianum (Supplementary Fig. S1a). Meanwhile, S24, which was characterized by a curved leaf apex, was identified as P. cirrhifolium (Supplementary Fig. S1b). The remaining samples (S19–S23, S25, S32 and S33) were identified as P. sibiricum (Supplementary Fig. S1c). The leaves of the remaining samples (S1–S18, S26–S31, S34, S35, S38, S40, S42–S47) were alternate or mostly alternate. The leaf apexes in plants from six regions (S6, S7, S15, S17, S34 and S35) were, similarly to P. sibiricum, strongly cirrhous or curved. Those samples were identified as P. alternicirrhosum Hand-Mazz. (Supplementary Fig. S1d). Finally, samples from the remaining regions (S1–S5, S8–S14, S16, S18, S26–S31, S38, S40, and S42–S47) were identified as P. cyrtonema (Supplementary Fig. S1e).

HCA based on morphological indexes further verified the results of morphological observations and divided all sampled plants into 14 groups (Fig. 2). The samples from all the regions were classified into two clusters with a squared Euclidean distance of 25; clade II comprised 19 samples and clade I included the rest of the samples (I). Clade I included samples with alternate leaves, whereas clade II comprised mainly plants with verticillate leaf arrangement. The exceptions were S6, S7, S15, S17, S34, and S35, which have alternate leaves, but because of their cirrhous or curved leaf apexes were classified into the cluster with verticillate leaves by HCA. Cluster I was further subdivided into two subclusters, of which one included samples from a single region (S43, Id), while the samples from the remaining 27 regions were further clustered into two small groups of nine and 19 regions. Moreover, the 1st group (Ic), with nine regions, was further subdivided into two subgroups containing one (S13) sample and eight other subgroups at a squared Euclidean distance of 2. The second group was also subdivided into two subgroups, the first one with eight regions (Ib), and the other one with the remaining regions (Ia) with a squared Euclidean distance of 3. The smaller cluster (II) contained two small subclusters of four regions (IIa), and the other group of 15 regions sharing a node with them at a squared Euclidean distance of 9. The first subcluster with four regions was subdivided into three groups sharing a Euclidean square distance of 2. The second subcluster was subdivided into two small groups, one containing five regions and another 10 regions (IIc), with a squared Euclidean distance of 6. The first group (IIb) was subdivided into two subgroups, and the second group was subdivided into three subgroups, sharing a node with a squared Euclidean distance of 2. Therefore, the analysis of the samples using selected morphological data successfully identified the species of PR.

Compared with morphological identification, the ISSR marker data assigned the samples (S19, S20, S21, S22, S23, S25, S32 and S33) of cluster I (Fig. 1a) collected from Henan and Shaanxi Province to P. sibiricum, except for the sample S24, which was identified as P. cirrhifolium. In cluster II, S36 and S37 were clustered with S39 and S41, which were all identified as P. kingianum. The samples S6, S7, S15, S17, S34, and S35, identified as P. alternicirrhosum, did not form a single clade, although S34and S35 were grouped together. The analysis of the SCoT markers produced similar results. Clade I included three samples, S6, S7, and S17, which were identified as P. alternicirrhosum, whereas Clade II included S19, S20, S21, S22, S23, S25, S32, and S33, which were all identified as P. sibiricum. The samples S36, S37, S39 and S41, which were identified as P. kingianum, formed a clade.

Consequently, the HCA of morphological data and the cluster analysis of ISSR and SCoT markers using the UPGMA are suitable for the identification of the germplasm resources of PR. The samples of P. kingianum were grouped with those of P. sibiricum, P. alternicirrhosum, and P. cirrhifolium (S24) based on morphological characters. However, the molecular marker-based identification grouped P. kingianum with P. cyrtonema. These results showed that the genetic distance between P. kingianum and P. cyrtonema was smaller than that between P. kingianum and P. sibiricum. This is possibly caused by the geographical distribution of these species. According to the Flora of China [46], P. cyrtonema is mainly located in southern China, P. kingianum is mainly distributed in Sichuan, Guizhou, and Yunnan province in western China, and P. sibiricum is mainly found in northern China. Hence, although many morphological differences exist between P. kingianum and P. cyrtonema, their genotypes showed greater similarity and smaller genetic distance.

The species identification of PR that relies on rhizome traits is uncertain and prone to errors. Furthermore, accurate identification of PR based on morphological characters can be completed only after cultivating the rhizomes for 3–6 months. In contrast, identification using molecular markers is faster. Both ISSR and SCoT markers recovered similar values for genetic distance among the species. Specific qualitative and quantitative characters that were selected for morphological identification in combination with molecular marker identification were efficient in determining the PR species from other species of Polygonatum.