Twenty four batches of agarwood and its products were included in this study ( Table 1 ). Four batches of dry leaves of A. sinensis, A. agallocha, A. crassna, and A. subintegra were collected from the South medicinal botanical garden. Twenty batches of commercially available agarwood products (such as wood patch, chips, blocks, powder, and Chinese patent medicines) were collected. All materials were identified by Professor Yangyang Liu, and all specimens were kept in the agarwood identification center herbarium.

Samples taken from leaves (100 mg) and dried wood or products (30 mg) were cleaned and rubbed for 5 min at a frequency of 30 r/s in a case of liquid nitrogen. DNA extraction proceeded using the modified HP plant DNA kit (OMEGA), with the addition of 1,000 μL of CPL and 10 μL of β-mercaptoethanol. The concentration and purity of the obtained sample DNA were determined, and the DNA was stored at −20°C for further use. Four batches of plant leaf DNA were interrupted using ultrasound; libraries were constructed by end repair, addition of sequencing junctions, purification, and PCR amplification; and the library fragment size and quality were examined. Sequencing was performed using Illumina’s high-throughput sequencing platform, NoVaSeq 6000.

Approximately 4 GB of raw data consisting of 150-bp paired-end reads were generated. Raw data in fastq format were first processed through fastp. In this step, clean data (clean reads) were obtained by removing reads containing the adapter, reads containing ploy-N, and low-quality reads from the raw data. At the same time, the content of guanine and cytosine (GC) was calculated. All downstream analyses were based on clean, high-quality data. Chloroplast genomes were assembled using GetOrganelle v1.7.6.1 with default settings, which filtered plastid-like reads, conducted de novo assembly, purified the assembly graph, and generated the complete chloroplast. All chloroplasts were initially annotated using PGA and GeSeq, with the annotated chloroplast of A. subintegra (NC052859) selected as the reference. For confirmation, all annotations were compared with previously published chloroplasts of A. subintegra (NC052859) and manually examined. Circular genome maps were visualized using OGDRAW v1.3.132. The sequences were aligned using MAFFT v7 with default settings (strategy of FFT-NS-2).

In this study, we compared and analyzed the chloroplast genomes of two medicinal agarwood species (A. sinensis and A. agallocha) and two non-medicinal species (A. crassna and A. subintegra) obtained by sequencing. The GenBank numbers of the four species are OR608759, OR608758, OR608757, and OR608760, respectively. Meanwhile, the chloroplast genomes of 12 published species of Aquilaria Lam. were selected in this study, including A. agallocha (NC065040 and MZ145047), A. sinensis (MN720647 and MN147870), A. beccariana (MN125347 and NC052855), A. crassna (MN125348 and NC043844), A. cumingiana (MZ145048 and NC065041), A. hirta (MN125349 and NC052856), A. malaccensis (MH286934 and NC041117), A. microcarpa (MN125350 and NC052858), A. rostrata (MN125351 and NC052858), A. rugosa (MZ145049 and NC065042), A. subintegra (MN147871 and NC052859), and A. yunnanensis (MG656407 and NC052859); full-length comparisons of these chloroplast genomes were performed using MEGA11 software; and the results were used to further screen for loci that could be used in medicinal agarwood, as well as other species-specific loci. Based on the specific loci, DNA mini-barcoding was designed for the identification of medicinal agarwood species using Oligo software. Primers that were likely to have dimers, hairpin structures, or excessive annealing temperatures were abandoned using Oligo Calc software.

PCR amplification, DNA melting, and end-point fluorescence level acquiring PCR amplifications were performed in a total volume of 20 μL on a Light Cycler^® 96 System (Roche). The reaction mixture contained 1 μL of genomic DNA, 10 μL of 2× TransStart^® Tip Green qPCR SuperMix (Trans), 1 μL of 10 mM forward and reverse primers ( Table 2 ), and nuclease-free water up to the final volume. The amplification protocol consisted of 10 min of pre-incubation at 95°C followed by 45 cycles under the following conditions: denaturing step at 95°C for 10 s, annealing at 58°C for 15 s, and extension at 72°C for 15 s. After the amplification, before the HRM step, the products were heated to 95°C for 60 s and then cooled to 40°C for 60 s. The HRM conditions were 65°C–95°C, rising to 0.07°C/s with 15 acquisitions per °C. In addition, these amplified PCR products were sequenced to verify the reliability of the results. All amplified products were sequenced by Guangzhou Ige Biotech Co., Ltd. The sequence was assembled using Codon Code Aligner 5.1.5 and aligned using ClustalW and DNAMAN.

To evaluate the repeatability of the HRM analysis method, the A. sinensis, A. agallocha, A. crassna, and A. subintegra species were randomly selected to perform the PCR reaction and HRM procedure, and the assays were performed in three repetitions. The parameters of the DNA template amount were tested using 1 ng, 5 ng, 10 ng, 25 ng, 50 ng, and 100 ng of DNA template in the reaction to evaluate the variation in the amount of DNA template. After obtaining suitable primers to test the sensitivity of the developed method, HRM was performed on standard samples. The A. sinensis, A. agallocha, and A. crassna sample powders were mixed in pairs at different ratios of 1%, 5%, 10%, 25%, 50%, 75%, 90%, 95%, and 99%. Real-time PCR amplification was performed as described earlier.

For A. sinensis, A. agallocha, A. crassna, and A. subintegra, the total length of the circular chloroplast genome was determined to be 174,828–174,911 bp, containing a long single-copy sequence (LSC) region of 87,281–87,359 bp and a short single-copy sequence (SSC) region of 3,345–3,355 bp, separated by two SSCs (IRa and IRb) of 42,102–42,109 bp ( Figure 1 ). The GC content of the four chloroplast genomes was 36.7%, with inverted repeat (IR) regions having a higher GC content (38.6%–38.8%) than LSC (34.9%–35.0%) and SSC (29.0%–29.1%). This differs from the published A. agallocha (MZ145047 and NC065040) and A. subintegra (MN147871 and NC052859) genomes by only 3 bp, whereas the GC contents of the four sequenced species are the same as those of the above published species. The difference between the sequenced species and A. sinensis (MN720647 and MN147870) is also only 3 bp or 4 bp, mainly in the LSC region, whereas the GC contents are the same. A. crassna (NC043844) exhibits differences relative to the sequenced species in the length of the IR and LSC regions ( Supplementary Table S1 ). The IR/LSC and IR/SSC barcodes of the A. sinensis, A. agallocha, A. crassna, and A. subintegra chloroplast genomes were compared ( Supplementary Figure S1 ). Among the four chloroplast genomes, the rps19 gene crossed the LSC and IRa regions with a 16-bp extension into the IRa region. The ndhF gene crossed the IRa/SSC boundary, with 28 bp, 27 bp, 26 bp, and 26 bp of extension into the IRa of A. sinensis, A. agallocha, A. crassna, and A. subintegra, respectively. The rpl2 gene was in the IRa and IRb regions, the rpl32 gene was in the SSC region, and the psbA and rpl22 genes were in the LSC region.

The four chloroplast genomes encode 141–142 genes, with the main difference being the ycf15 or ycf1 gene in the IR region. These genes include 95, 96, and 97 coding sequence, 38 transfer RNA genes, and eight ribosomal RNA genes. The first major group of genes related to transcription and translation includes ribosomal protein subunit genes, RNA polymerase genes, ribosomal RNA genes, and transfer RNA genes, totaling 77 genes, of which tRNA genes have the largest number, totaling 38; the second major group of genes related to photosynthesis includes photosystem I genes, photosystem II genes, nicotinamide adenine dinucleotide phosphate (NADPH) dehydrogenase genes, cytochrome b/f complex genes, ATP synthase, and rubisco large subunit. The third major group of genes related to other biosynthesis processes consisted of 14 genes ( Supplementary Table S2 ).

The nucleotide diversity was estimated using the four chloroplast genomes of Aquilaria Lam., and the nucleotide diversity values (Pi) ranged from 0 to 0.02056, with an average of 0.001164, indicating that there is generally mild divergence among the genomes of these four species. However, as shown in Supplementary Figure S2 , seven loci showed significantly higher Pi values (> 0.008); these included ndhC-trnV (0.02056), psbM-trnD (0.01289), rps4 (0.012), psbJ-petA (0.01144), trnk-rps16 (0.00911), psaA-ycf3 (0.00911), and ndhF-rpl32 (0.008), which can be applied to the mini-barcoding design.

In this study, by performing full-length comparisons of the chloroplast genomes of 12 species of the genus Aquilaria Lam. and combining the highly variable regions in the four chloroplast genomes obtained by sequencing, 59 SNP loci were identified that could distinguish two medicinal agarwood species (A. sinensis and A. agallocha) from 10 non-medicinal agarwood species, among which, 18 loci could distinguish A. sinensis from other species and 38 loci could distinguish A. agallocha. The other three loci could distinguish between medicinal agarwood and other species, but loci 130573 and 130574 needed to be used in combination thymine and adenine (TA) for medicinal agarwood and thymine and guanine (TG), guanine and guanine (GG), and guanine and adenine (GA) for other species) ( Supplementary Table S3 ).

We designed 57 primer pairs (not listed) and obtained 10 improved primer pairs by qPCR screening. Five primer pairs obtained from the psbM-trnD, psbA, rps16, petN, and ndhE-psaC regions could be used to distinguish A. sinensis from other species. Four primer pairs obtained from the rps4, atpE, ycf1, and rps15-trnN regions could be used to distinguish A. agallocha from other species. In addition, the primer pairs obtained from the matK regions could distinguish between medicinal agarwood and other species.

We performed HRM analysis on four species of the genus Aquilaria (A. sinensis, A. agallocha, A. crassna, and A. subintegra) and screened 10 pairs of primers with the ability to identify medicinal and non-medicinal agarwood by plotting three melting curves. In this study, the normalized fluorescence curve ( Figures 2A1–J1 ) and normalized fluorescence difference curve ( Figures 2A2–J2 ) produced only different graphs. The normalized fluorescence curve demonstrates consistent results in the five primer pairs of the chloroplast genome (Pt197FR, 1535FR, 5987FR, 30515FR, and 125057FR), which could clearly identify A. sinensis from A. agallocha, A. crassna, and A. subintegra ( Figures 2A1–E1, A2–E2 ). The four primers pairs (48959FR, 55861FR, 118649FR, and 148002FR) clearly identified A. agallocha from A. sinensis, A. crassna, and A. subintegra ( Figures 2F1–I1, F2–I2 ). In addition, the primer pairs (M856F2R1) could simultaneously distinguish medicinal agarwood (A. sinensis and A. agallocha) from non-medicinal agarwood (A. crassna and A. subintegra) ( Figures 2J1, J2 ). To intuitively describe the curve characteristics, we constructed the negative derivative of the fluorescence versus temperature (dF/dT) curve to show the melting temperature (Tm) of different samples ( Figures 2A3–J3 ). The melting curve further magnifies the difference, making it easier to observe the difference between the four species. In short, the differences in these curves are caused by single-base variations in the sequence.

We sequenced and analyzed the PCR products of four species (A. agallocha, A. sinensis, A. crassna, and A. subintegra) amplified by 10 primer pairs, and the results were identical to the bases of the SNP loci ( Figure 3 ), verifying that the differences in HRM solubility temperature and curve were due to differences in the bases of the SNP loci, indicating that the method was feasible.

In this study, we investigated the amount of template DNA and the melting curve. The normalized fluorescence difference curve results ( Figure 4 ) demonstrated that the correct groupings were obtained in the range of 1 ng to 100 ng of DNA for the four Aquilaria Lam. species with 10 primer pairs. The melting curve of the non-medicinal agarwood slightly deviated only when 50 ng of template DNA was used, but this did not affect the identification and classification of the species. In the applied HRM analysis, the melting curves were more stable at low template concentrations (≤ 25 ng/μL) than at high concentrations (≥ 50 ng/μL) within an appropriate dynamic range. We statistically analyzed the melting temperatures for the six DNA template concentrations (1 ng/μL to 100 ng/μL) from four agarwood species ( Table 3 ).

Different primers produced different melting temperatures for the target to distinguish the species of medicinal agarwood. Among them, the melting temperatures of five primer pairs (Pt197FR, 1535FR, 5987FR, 30515FR, and 125057FR) could distinguish A. sinensis from other species; four primer pairs (48959FR, 55861FR, 118649FR, and 148002FR) could distinguish A. agallocha from other species; and primer M856F2R1 could distinguish two medicinal agarwood species (A. sinensis and A. agallocha) from other agarwood species.

We first extracted DNA by mixing different agarwood powders in proportion to each other and then verified the adulteration sensitivity of medicinal and non-medicinal agarwood using the mini-barcoding HRM method. We mixed the samples of A. sinensis, A. agallocha, and A. crassna with each other and obtained the normalized fluorescence difference curves and the negative derivatives of the fluorescence versus temperature curves by HRM analysis to determine the results. When the primers Pt197FR, 1535FR, 5987FR, 30515FR, and 125057FR were used for the adulteration verification of A. sinensis and A. crassna ( Figure 5A ), the melting curve line shape of the mixed samples deviated from that of the homozygote when other agarwood samples were added to the domestic medicinal agarwood (A. sinensis) was more obvious in the identification of the negative derivative of the fluorescence versus temperature curve. Similarly, using primers Pt197FR, 5987FR, 30515FR, 125057FR, and M856-F2R1, the mixed adulteration verification of domestic medicinal agarwood with imported medicinal agarwood revealed that the line shape of the melting curves of the mixed samples also deviated from those of the homozygote samples ( Figure 5B ). The melting curve of the HRM mixed sample using primer 5987FR showed an extra peak ( Figures 5A a3–1, B b2–1). Additionally, the adulteration of imported medicinal agarwood (A. agallocha) was verified using primers 48959FR, 55861FR, and 118649FR, and the line shape of the melting curve changed as the amount of other agarwood increased ( Figure 5C ). The detection limit of adulteration in medicinal agarwood was examined by this method, and the melting curve line shape at ≥ 1% of adulteration was compared to that of the homozygote. However, as a control, the addition of a homozygote of medicinal agarwood is essential in HRM analysis. Therefore, this method can effectively determine whether the samples are adulterated or not.

In this study, HRM analysis and PCR sequencing of 20 samples using 10 primer pairs gave consistent results. Eight batches of samples were medicinal agarwood (including one batch of proprietary Chinese medicines containing agarwood), and 11 batches were non-medicinal agarwood (including two batches of proprietary Chinese medicines containing agarwood). One batch of proprietary Chinese medicines containing agarwood may have used medicinal agarwood mixed with non-medicinal agarwood (adulteration rate between 5% and 10%), but this requires further verification at a later stage. Of these, seven batches (C001, C002, C003, C004, C011, C016, and C020) of medicinal agarwood samples were A. sinensis, and one batch (C014) was A. agallocha. Most of the 11 batches of non-medicinal agarwood originated from Vietnam, Malaysia, Thailand, Kalimantan, Papua New Guinea, and other Southeast Asian countries ( Table 1 ; Supplementary Table S4 ). Although no further species identification of these 11 batches of non-medicinal agarwood samples was performed, this does not affect the judgment of the results. The 10 primer pairs generated in this study were effective in identifying the species of medicinal agarwood samples, among which primer M856-F2R1 could directly and rapidly determine whether the species of the samples were medicinal agarwood. Nevertheless, two specific sites were simultaneously present in samples C018 (primers 30515FR, 125057FR, and M856F2R1) and C019 (primers 1535FR and 118649FR) when amplified sequences from those samples were sequenced. Furthermore, the amplified primer 5987FR sequence from the C019 sample showed heterozygous sequencing. Future research will shed more light on these unanticipated occurrences.

It has been reported that both species of the genus Aquilaria and Gyrinops (Thymelaeaceae) can produce agarwood (Lee et al., 2022), a resin-containing wood (Liu et al., 2013, 2019), and it is difficult to differentiate the species information during the trade process, especially when it is crushed and added to preparations to improve the ease of adulteration. Currently, there are only two legal species of traditional Chinese medicinal agarwood in China. Given the high use of commercially available agarwood and the fact that approximately 138 types of medicines contain agarwood, it is meaningful to develop a method to rapidly and accurately identify the species of Chinese medicinal agarwood. Our findings will provide a reference for the identification of medicinal agarwood species in Japan and the Taiwan Province of China, which is of practical significance and application value for the development of the agarwood pharmaceutical industry.

Until now, molecular identification methods for agarwood have mainly focused on DNA barcode screening and its combinations, such as ITS2+trnL-trnF (Lee et al., 2016), ITS2+matK (Kang et al., 2018), ITS+matK, ITS+rbcL, ITS+trnL-trnF (Li et al., 2018; Lin et al., 2023), and psbJ-petA+trnT-trnL (Hishamuddin et al., 2023). However, all of these methods require agarwood plant material as the basis for the identification of species using sequencing and construction of phylogenetic trees. In addition, these barcode amplification products require long lengths, such as the ITS length of approximately 673 bp, the matK length of approximately 778 bp, and the trnL-trnF length of approximately 938 bp or even longer. Moreover, as DNA extracted from agarwood is stored for a long time and can undergo severe degradation, the amplification results using these primers are often unsuccessful. Consequently, plant leaves are often required to build DNA barcode data combinations or databases, which is difficult to apply in practice. Therefore, in this study, by comparing the chloroplast genome sequences of different species of the genus Aquilaria Lam., we screened for SNP sites and used the combination of mini-barcoding and HRM, which can either visually determine the species from the melting curve or sequence the amplification products to obtain the specific base sites. This method avoids the need to construct multiple DNA barcode combination libraries, which improves the timeliness and accuracy of traditional Chinese medicinal agarwood species screening.

SNP molecular marker technology can be used to identify abundant loci and has been applied to the analysis of germplasm resource evolution and species identification (Li et al., 2023; Palasciano et al., 2023). Chloroplasts are maternally inherited, and the genetic information of each species can be stably inherited. By analyzing the differences in SNP and designing the corresponding primers, it is possible to accurately identify different species. In the present study, 59 SNP loci were obtained from the chloroplast genomes of 12 published species of Aquilaria Lam. and could be used to distinguish traditional Chinese medicinal agarwood (A. sinensis and A. agallocha) from non-medicinal agarwood species. The SNP loci were obtained correctly along with the four agarwood chloroplast genomes obtained by sequencing in this study. We further screened to obtain the 10 best primer pairs, of which five pairs could distinguish A. sinensis from other species, four pairs could distinguish A. agallocha from other species, and another one pair could distinguish two medicinal agarwood species from non-medicinal agarwood. These primers were from the psbM-trnD, psbA, rps16, petN, ndhE-psaC, rps4, atpE, ycf1, rps15-trnN, and matK regions, and their amplification products were < 200 bp with high amplification success.

HRM analysis is a simple, rapid, and cost-effective method for identifying genetic variants. The dissociation of double-stranded DNA is measured in real time by HRM, and the polymorphism of the sample DNA can be easily monitored according to the shape and position of the melting curve (Gundry et al., 2003; Vossen et al., 2009; Suesatpanit et al., 2017; Kim et al., 2023). In this study, we successfully screened the chloroplast genome from the genus Aquilaria Lam. for species identification of SNP loci that can be used for medicinal agarwood and validated them by applying HRM and sequencing. We conducted methodological investigations using these 10 primer pairs, and the melting curve line shapes and melting temperatures were correctly grouped for agarwood sample DNA concentrations in the range of 1 ng to 100 ng. In the applied HRM analysis, the melting curves were more stable at low template concentrations (≤ 25 ng/μL). We also examined the detection limits of adulteration using different ratios of agarwood powder. The melting curves and temperatures deviated from those of the homozygote sample after the addition of non-medicinal agarwood to the medicinal agarwood samples, and the detection limit of adulteration was 1%. The method was validated in 20 batches of commercially available samples, among which eight batches were identified as medicinal agarwood (seven batches were traditional Chinese medicinal agarwood and one batch was imported medicinal agarwood), and 11 batches were non-medicinal agarwood. Additionally, we unexpectedly found that one batch of proprietary Chinese medicine containing agarwood used a mixture of medicinal and non-medicinal agarwood, with an adulteration rate between 5% and 10%, although this needs to be further verified at a later stage.

However, it is usually during HRM analysis that false-positive results are most likely to occur. The main practices that are effective in avoiding false positives are the use of multiple characterized SNP loci assessment, amplified sequence sequencing validation, and the use of target species controls. In addition, the occurrence of false-positive results is also related to PCR amplification conditions and target DNA concentration and purity. In this study, the above situation was evaluated, and the 10 pairs of primers obtained in this study can satisfy the identification of two medicinal agarwood species, and, according to the purpose, one or more of the primers can be selected to carry out species identification. If necessary, the reliability of the method can be verified by direct sequencing of the amplification products, which, in this study, was shown to be reliable by the sequencing results. In addition, methodological validation was performed in this study, and because of the differences in DNA concentration and purity between different samples and the reaction reagents themselves, we suggest that, when using this method, it will be necessary to add positive homozygote samples for analysis.

Therefore, compared to the traditional DNA barcode analysis method, the method established in this study is time-sensitive and can be used to quickly, intuitively, and accurately determine the medicinal agarwood species. In the future, in addition to the traditional Chinese medicine agarwood species identified in this study, there are many other agarwood species that need to be identified, and the mini-barcoding of each species can be further developed along the lines of this study according to the usage or industry needs.