Considering that L. macranthoides in Flos Lonicerae has many cultivated varieties, a wide planting area, and more stable and representative genetic characters. In contrast, L. hypoglauca, L. confuse, and L. fulvotomentosa currently have no cultivated varieties, and there are few wild varieties with low representativeness and other factors. In this experiment, L. macranthoides was selected as the representative of Flos Lonicerae. The experiment was conducted with 39 wild and cultivated varieties of L. japonica and L. macranthoides existing in China. There were 35 varieties of L. japonica and four varieties of L. macranthoides ( Figure 1 ). The specific variety information is shown in Table 1 . For each germplasm accession, seven vigorous plants were randomly selected at the seedling stage, and young leaf samples totaling 5 g were collected on tinfoil, frozen immediately in liquid nitrogen, and transported to the laboratory, where they were stored at −80°C until they were used for DNA extraction.

After grinding the leaf samples, genomic DNA was extracted using a Plant Genomic DNA Kit (TIANGEN, China), and the quality and concentration of DNA were measured using a NanoDrop2000 UV spectrophotometer (Thermo Fisher, Waltham, MA, USA). Then, a paired-end library with a length range of 300–500 bp was constructed via double-digested (ddRAD) library construction of qualified sample DNA. First, 500 ng of genomic DNA was reacted with 0.6U EcoRI (NEB), T4 DNA ligase (NEB), ATP (NEB), and EcoRI connectors (including index sequences of differentiated samples) at 37°C for 3 h and annealed at 65°C for 1 h. The restriction enzyme NlaIII (NEB) and the NlaIII connector were then added and allowed to react for 3 h at 37 °C. After the reaction, the endonuclease was inactivated at 65°C for 30 min in a polymerase chain reaction (PCR) amplifier. Digested products of 400–600 bp were recovered via agarose gel electrophoresis and quantified using Qubit 3.0 (Life Technology). After mixing 39 samples in equal quantities, an Illumina TruSeq kit was used to construct a DNA library of the mixed products (Hanania et al., 2004).

An Illumina NovaSeq 6000 PE150 was used for sample sequencing after library construction, and data quality control was performed on the original sequenced reads (Dellaporta et al., 2016). Fastp software (version: 0.20.0) was used to remove reads with an unknown base number N <5, reads with a length of bases <50%, quality value <5, connector sequences, and other low-quality sequences to obtain clean data. The detailed parameters were set to -q5 -n5. Burrows–Wheeler Aligner (BWA, 0.7.17-R1188) was then used to compare the sequenced reads with the reference genome (Byers et al., 2012). The reference genome used was the GWHAAZE000,000,000 genome Fasta (L. Lonicera) (download website: https://ngdc.cncb.ac.cn/search/?dbId=gwh&q=SAMC097356), and the parameters were set as -M -R. The insert size and coverage depth of each sample were counted, and variation was detected by comparing the positions of clean reads on the reference genome (Chao et al., 2009). Then, the same files generated by the comparison were converted to bam format using Samtools software (version: 1.9). Finally, Picard MarkDuplicates (version 2.21.2) was used to detect duplicate tags, and high-quality reads were retained for subsequent analysis.

The Genome Analysis Toolkit software (GATK) (https://gatk.broadinstitute.org/hc/en-us; version 4.1.4.1) HaplotypeCaller function was used to call variable outliers. The FilterVCF function was used to filter low-quality mutation sites to obtain the final SNP site set and obtain SNP statistics (Huang et al., 2013). The genome file and the structure annotation gff file were used to annotate the variable outliers using the eff mode in the snpEff software (version 4.3t) (Pariasca-Tanaka et al., 2015). SNPs and INDELs were screened according to the following criteria (Elbasyoni et al., 2018): average sequencing depth ≥5× (Davey et al., 2011), minor allele frequency (MAF) ≥0.05, information integrity ≥0.70, SNP quality value Q ≥30, QD <2.0, MQ <40.0, FS >60.0, SOR >6.0, MQRankSum <−12.5, and ReadPosRankSum <−8.0 (Angel et al., 2012).

Based on the high-quality SNPs obtained, GCTA software (version 1.92.1, http://cnsgenomics.com/software/gcta/#Overview) was used to perform principal component analysis (PCA) (Abdelhafez et al., 2019) with the following parameter settings: – make grm – autosome (Gong et al., 2016). The maximum likelihood (ML) method in FastTree software (version 2.1.9) was used to construct an evolutionary tree of 39 Lonicera samples (Saxena et al., 2012). Finally, admixture software (Sorkheh et al., 2007) (version 1.3.0) was used to analyze the population genetic structure (Jaganathan et al., 2015).

DNA fingerprint construction was performed based on the high-quality SNPs obtained (Xu et al., 2017). As few markers as possible were used to identify as many varieties as possible to achieve the purposes of simplicity, efficiency, and economy (Wang et al., 2009). Fifty core markers with a high detection rate and significant amounts of polymorphism that could distinguish all varieties were screened out based on the size of the polymorphic information content (PIC) value and distribution frequency (Pavana et al., 2012), and then a DNA fingerprint was constructed. Furthermore, the 15 leanest SNP combinations were further screened out (Jones and Mackay, 2015) to identify Lonicera varieties with lower costs and faster speeds.

Primers for the 15 SNP loci were designed using Primer Premier 5 software ( Table 2 ), and the DNA of the validation test sample was extracted to conduct PCR. The parameters were set as follows: length, 18–30 bp; Tm value, 55–65 degrees; and GC content, 40%–70%. The primers were synthesized by Sangon Bioengineering Co., Ltd. (Shanghai, China). Four samples were randomly selected from 39 known samples, and three unknown samples were added ( Table 1 ). A total of seven verification DNA samples were selected as templates for PCR amplification. The total volume of the PCR mixture was 25 μl, containing 1 μl DNA template (100 ng/μl), 1 μl each of Primer F and Primer R (10 pmol/L), 0.2 μl Taq Plus DNA polymerase (5 U/μl), 2.5 μl 10× PCR buffer with Mg²⁺, and 1 μl dNTP (10 mM). The PCR reaction program involved pre-denaturation at 95°C for 5 min, denaturation at 94°C for 30 s, annealing for 30 s to 63°C, and extension at 72°C for a total of 30 s; this was repeated for 30 cycles, after which there was a final extension at 72°C for 10 min before storage at 4°C. After PCR amplification, 5 μl of PCR products were taken and subjected to 1% agarose gel electrophoresis at 150 V and 100 mA. After 10–20 min of observation, the target PCR band was cut from the gel, recovered, and then sequenced using 3730XL. The sequencing data were searched in the results group, and sequence analysis software SnapGene was used for analysis. Finally, SeqMan software was used to observe the peak map to test the reproducibility, genetic stability, and specificity of 15 SNPs.

The sequencing of 39 Lonicera samples on the Illumina NovaSeq 6000 PE150 yielded a total of 84.88 Gb of clean data, and Q30 reached more than 93.25% per sample ( Table 3 ). The obtained clean reads were mapped to the reference genome, and the mapping efficiency reached 99.38% ( Table 4 ).

After sequencing, GATK software detected many SNP variants in 39 Lonicera samples; sequencing generated a total of 13,143,268 SNPs ( Figure 2A ). A total of 3,374,929 filtered SNPs were obtained, distributed on nine chromosomes ( Figure 2B ). Finally, the SNP distribution map on each chromosome was drawn according to the number and density of SNPs ( Figure 2C ). There were about 400,000 SNPs on chromosome GWHAAZE000000000 1 and 600,000 SNPs on chromosome GWHAAZE000000000 2. For chromosomes GWHAAZE000000000 3 to 9, the numbers of SNPs on these six chromosomes were similar, and there were about 300,000 SNPs on each chromosome.

The PCA of high-quality SNPs screened from 39 Lonicera samples was conducted using GCTA software. L. macranthoides samples from Hunan were clearly clustered on one side ( Figure 3A and Table 1 ), while L. japonica samples from the same species were clustered on the other side ( Figure 3B and Table 1 ). For the two clusters classified by species, detailed clustering classification was further conducted according to a variety of characteristics of different Lonicera samples. For example, honeysuckle varieties with red flower buds were clustered together. Crossbred samples were completely clustered together. Admixture software was used to analyze the population structure of the 39 Lonicera samples. The cross-validation error rate was the lowest when K = 4, a result that was consistent with the variety characteristics and source classification of the 39 Lonicera samples ( Figure 3C and Table 1 ). The 39 Lonicera germplasm resources were divided into four subgroups ( Figures 3D, E ).

Finally, FastTree software was used to construct an evolutionary tree, and the 39 Lonicera samples were clustered into four groups ( Figure 4 ), consistent with the results of the PCA and the population structure analysis. In addition, through further analysis, we found that the four clusters also had high similarity in phenotype, species, and origin. Each group was supported by a high bootstrap value.

Based on the results of simplified sequencing, SNP loci with a PIC value greater than 0.30 that were uniformly distributed on nine chromosomes were screened as the core loci for the construction of the Lonicera variety fingerprint. Finally, 50 SNP loci were selected for constructing the fingerprint, and all varieties were distinguished ( Figure 5A ). Among the 50 core SNPs, loci Markers 8, 17, 23, and 32 had the highest PIC value of 0.38. Marker 4 had the lowest PIC value of 0.30. The average PIC value was 0.35, indicating moderate polymorphism. Using this set of core SNP loci combinations, 39 samples of Lonicera were compared in pairs. The statistical results for the number of different loci between the samples showed many different loci between each pair. To quickly and cheaply distinguish Lonicera varieties, in this study we screened out the combination with the minimum number of SNPs that could distinguish 39 Lonicera varieties from the obtained 50 core loci; we obtained the most simplified SNP combination that contained 15 SNPs distributed on eight chromosomes ( Figure 5B ). When using the combination of 15 SNPs to distinguish samples, there was at least one different SNP between each of the two samples that could effectively distinguish the 39 Lonicera varieties.

Some Lonicera samples could only be distinguished based on one SNP. For example, “Longhua” (SPL022), “Yincuilei” (SPL019), and “Jincuilei” (SPL020) could only be distinguished by Marker 7; “Yujin 2” (SPL010) and “Yate 1” (SPL004) could only be distinguished by Marker 6. “Yujin 5” (SPL014) and “Yujin 5 1-2” (SPL028) could only be distinguished by Marker 11. These results also illustrate that it is difficult to develop fingerprints against the background of complex and chaotic genetic relationships in L. japonica germplasm.

We randomly selected four samples, “Yujin 1,” “Yujin 2,” “Mixian Damaohua,” and “Telei 1,” from a set of 39 samples to verify the reproducibility of SNPs. Fifteen SNPs had the same site specificity as the previous test results ( Table 5 ), with a clear and clutter-free locus peak map, 100% reproducibility, and ideal results. Next, we tested the genetic stability of 2-, 5-, 7-, 8-, 9-, and 10-year-old plants of “Yujin No. 1.” The results showed that the genetic stability of honeysuckle of the same variety in different years was as high as 97.33% ( Table 6 ). Finally, three new varieties, “Yateliben,” “Bainong 2,” and “Weizi,” were added to verify the specificity of 15 SNP loci ( Table 7 ). The results showed that 15 SNP loci could effectively distinguish the genotypes of the newly added varieties ( Supplementary Figure 2 ).

As the original and major producer of Lonicera, China has three main production areas: Fengqiu in Henan Province, Julu in Hebei Province, and Pingyi in Shandong Province (Song, 2020). The planting area of L. japonica in 2022 reached 1,066.7 km² (Lv, 2020; Ma et al., 2022), with an annual output value of nearly 10 billion CNY (Liu, 2021; Liu, 2022). At present, many Lonicera varieties have been widely planted, but there has been a lack of standardized and decentralized management. Therefore, it was necessary to clarify the genetic relationship and population structure of existing Lonicera varieties.

The 2015 Edition of the Chinese Pharmacopoeia details the differences between L. japonica and Flos Lonicerae in medicinal history, sources, characteristics, chemical components, and other aspects, showing that Flos Lonicerae is classified as L. macranthoides, L. hypoglauca, L. confusa, and L. fulvotomentosa. L. macranthoides is more recognized because of its easy cultivation, large number of plants, and wide planting area. Lu et al. (2017) established the fingerprints of five batches of L. japonica and three batches of L. macranthoides using proton nuclear magnetic resonance (1H-NMR) that provided a basis for the quality control analysis of L. japonica and L. macranthoides. Ren et al. (2022) used microscopy and thin-layer chromatography (TLC) to evaluate the quality of L. japonica varieties in Henan, Hebei, Shandong, and other areas and concluded that the index components of Shandong L. japonica were the highest. This work provided a scientific basis for the selection of raw materials for honeysuckle decoctions. Deng et al. (2021) conducted a comparative study on tree-shaped and wild-type L. japonica using HPLC fingerprints and found that wild L. japonica contained more chemical components, while tree-shaped L. japonica was more consistent. Li et al. (2020) conducted a comparative analysis of the characteristics and pharmacological effects of L. japonica and L. macranthoides, providing data support for their application in traditional Chinese medicine. The related studies focused on Lonicera in the early stages; most of the studies on the germplasm focused on morphological analysis, index component differences, and quality evaluation among varieties, while there was almost no research on the genetic relationships and population structure of Lonicera species and varieties.

In this study, 39 existing Lonicera germplasm resources in China were collected and divided into four subgroups via simplified sequencing, PCA, population structure analysis, and evolutionary tree construction. Subgroup I contains the red variety L. acuminata of L. japonica. Subgroup II contained L. macranthoides. Subgroup III comprises the main, wild, and hybrid varieties of L. japonica from Shandong and Henan. Subgroup IV contained cultivars of L. japonica that comprised varieties with high yields, strong resistance to pests and diseases, and good characteristics that have been recognized on the current market, for example, SPL009, SPL002, SPL003, and other varieties. Based on these results, it should be clear that the lineages of wild Lonicera species were extremely complex, while the lineages of systematically bred L. japonica were relatively simple. For example, “Mixian xianhua” (SPL006), “Mixian wild” (SPL007), and “Mixian Damaohua” (SPL039) are wild varieties of L. japonica with extremely complex genetic relationships. In the population structure analysis, they all had ancestral lines from the four subgroups. However, the genetic relationship between “Yujin 1” (SPL009) and “Yate liangzhong” (SPL035) was simple and derived from one pedigree. After further analysis, this study also clarified the controversial genetic relationship of “Longyao” (SPL023), a variety produced in Hunan, as belonging to L. macranthoides. However, our results confirmed that “Longyao” (SPL023) was a cross-species hybrid of a female parent L. japonica and a male parent L. macranthoides, which explained why it simultaneously had the pedigree of two subgroups in the PCA and population structure analysis. The genetic relationship of “Yujin 6” (SPL027) was further clarified, and the results showed that it was crossbred from the female parent “Telei 1” (SPL011) and the male parent “Yujin 1” (SPL009). This finding well explained its position in Group IV of the evolutionary tree. In addition, four other hybrids were used in this study. Population structure analysis showed that all hybrids have genetic genes from their parents, which was largely consistent with expectations. At present, this research group is conducting a more in-depth analysis of hybrid characteristics.

Different from most previous studies focusing on honeysuckle morphology and quality evaluation (Li et al., 2022), the present study focused on aspects including analyzing the population structure of L. japonica, excavating the differences among varieties of L. japonica, and using SNP to distinguish varieties of L. japonica. This work can not only explain the population relationships of existing Lonicera germplasm resources but can also further clarify the genetic relationships between Lonicera varieties. These results provide a favorable reference and a basis for the management of Lonicera germplasm resources and have important significance for the breeding of new Lonicera varieties in the future.

In this study, SNP molecular marker technology and DNA fingerprinting technology were used; 3,374,929 SNPs were obtained from 39 samples of Lonicera; 50 core SNPs with a strong ability to identify Lonicera germplasm were screened; and DNA fingerprints were constructed. A fingerprint map constructed with 50 core SNP markers effectively identified and distinguished differences among different species and varieties of Lonicera. In addition, this study further filtered out 15 SNPs as the most compact site combination. The DNA fingerprints constructed based on the 15 SNP loci were used to detect different Lonicera samples, and ≥1 different locus was found in each Lonicera variety that could be used to distinguish different species and varieties rapidly, accurately, and at a low cost.

Considering the extremely complex genetic background of Lonicera germplasm, this study conducted further research on varieties that had several different loci ≤1 in the fingerprints of 15 SNP loci. After further analysis, it was found that these varieties not only had high repeatability in DNA fingerprints but also had high similarity in origin, species, and characteristics. For example, “Mihua 1” (SPL017) and “Fengjin 1” (SPL012) were introduced and domesticated from “Mixian Damaohua” (SPL008). The varieties “Longhua” (SPL022), “Jincuilei” (SPL020), and “Yincuilei” (SPL019) originated from Hunan Province and were L. macranthoides “Yujin 2” (SPL010) and “Yate 1” (SPL004) had purplish red flower buds and were red varieties of L. japonica. “Yujin 5” (SPL014) and “Yujin 5 1-2” (SPL028) were both hybrids, and their parents were the same varieties ( Table 1 ). Therefore, it was concluded that when the difference in the DNA fingerprint between two varieties is ≤1 locus, it means that the two varieties have certain similarity in species, origin, parents, flower color, etc., which is very likely to be two identical or similar varieties.

To verify the reliability of the DNA fingerprints comprising 15 SNP loci, four samples were randomly selected from the original 39 samples, and 15 SNP loci were analyzed via PCR and Sanger sequencing. The 15 SNP verification results were completely consistent with the previous test results. Subsequently, genetic analysis was conducted on 10 honeysuckle samples of the same variety from different years. The results showed that the genetic stability of these honeysuckle samples was extremely strong, reaching 97.33%, a level that was consistent with the characteristics of honeysuckle as a vegetatively propagated plant. However, it is worth noting that the genotypes of the two-year-old samples “J01,” “YJ02,” and “YJ03” and the five-year-old sample “YJ07” changed from A/A to A/G in Marker 10, a result that may have been related to the second allele of the SNP. However, as a vegetatively propagated plant, the genetic traits of honeysuckle are very stable, and the probability of genetic mutation between the same varieties is very low; thus, this phenomenon was most likely caused by the duality of SNP sites. Finally, three new varieties verified the specificity of SNPs, further demonstrating the reliability of the SNPs and DNA fingerprints developed in this study.

In general, the verification results of the 15 SNP loci and the genetic stability of honeysuckle were in line with expectations; the analysis not only confirmed that the 15 SNP loci constituted the authenticity and reliability of DNA fingerprints but also verified the high genetic stability of honeysuckle, providing a data reference and theoretical support for the identification of honeysuckle germplasm and the SNP fingerprint identification method.