A total of 21 voucher samples derived from seven species of Senna were collected from different geographical locations of Tamil Nadu, Karnataka, and Uttar Pradesh (Supplementary Figure 1). These included three individuals of each of the following: S. alexandrina, S. italica subsp. micrantha, S. spectabilis subsp. spectabilis, S. auriculata, S. uniflora, S. italica subsp. Italica, and S. tora species. The authenticity of the samples was verified by Dr. V. Sundaresan, Senior Scientist, CSIR-Central Institute of Medicinal and Aromatic Plants (CSIR-CIMAP), Research Centre, Bangalore, using the taxonomical monographs, floras and through the herbarium vouchers from the Botanical Survey of India (BSI). The reference voucher samples were deposited in the herbarium maintained at CSIR-CIMAP, Lucknow and the obtained accession details are tabulated in Supplementary Table 1. Additionally, raw drugs samples were purchased from different drug stores and markets in Bangalore (Karnataka) and Tamil Nadu, and were randomly selected for testing. The market samples were coded as HSAs, vouchered accordingly and deposited in the herbarium of CSIR-CIMAP, Lucknow. The sample codes are HSA01, HSA02, HSA03, HSA04 from Bangalore and HSA05, HSA06, HSA07, HSA08, HSA09, and HSA10 from markets in Tamil Nadu (Supplementary Table 2).

Seven species with three individuals (n = 3) were used for study. Total genomic DNA from the field samples of reference species was isolated following the protocol described by Mishra et al. (2016b). The genomic DNA from commercial samples (HSA01–HSA10) was isolated with DNeasy^® PlantMini Kit (Qiagen, Valencia, CA, United States) using the protocols supplied with the kit. The quality of the DNA was checked by electrophoresis on a 0.8% agarose gel with standard markers and quantified by spectrophotometric analysis (NanoDrop, ND-1000, United States). The DNA was diluted to working concentration of 25–50 ng/μl for PCR amplifications. Five commonly used candidate DNA barcode markers were amplified from three individuals of each species with the established primers, which included; two coding cpDNA regions rbcL and matK; one non-coding cpDNA intergenic spacer region, psbA-trnH and the nrDNA regions, ITS1 and ITS2. Details of primers and PCR conditions are listed in Table 2. PCR reactions were set up in a final volume of 50 μl with 1X Taq DNA polymerase buffer containing 1.5 mM Mgcl₂, (Genei Bangalore, India), 200 μM of each dNTP (Genei Bangalore, India), 5–10 pmol each of forward and reverse primers, 1 Unit of Taq DNA polymerase (Genei Bangalore, India) and 25–50 ng of template DNA. Successful amplicons were analyzed through electrophoresis on a 2% agarose gel. All reactions were performed in triplicates. Subsequently products of target molecular weight were purified with a Nucleospin PCR purification kit (Macherey-Nagel–07/2014, Rev.03) according to the supplied protocol and rechecked through electrophoresis on a 2% agarose gel. The PCR amplicons were sequenced from both ends through Sanger sequencing, using the Big Dye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Inc., Foster City, CA, United States) on an ABI 3130 XL genetic analyzer (Applied Biosystems, Inc., Foster City, CA, United States).

Specimen data for each barcode region were deposited in the Barcode of Life Data Systems (BOLD)¹ (Ratnasingham and Hebert, 2007) under the project CRCBS-Barcode marker for Senna authentication (Supplementary Table 1). All the related data are publicly accessible under the dataset DS-CIMAP². Lab-generated barcode sequences were deposited to the GenBank (Benson et al., 2013) under accession numbers listed in Supplementary Table 1. The electropherograms obtained for each region were base-called using PHRED vII (Ewing and Green, 1998). The proofreading of sequencing peaks and contig assembly were done in Sequencher v5.4.6 (Gene Codes Corporation, Ann Arbor, MI, United States). Finally, the contigs were analyzed through NCBI BLASTN 2.2.1+ (Zhang et al., 2000; Morgulis et al., 2008) and loaded on to BOLD using Identification Request for accessing their identity percentage with other similar sequences in database. All the barcode sequences were greater than 300 bases in length and free from contamination. The sequences were then aligned with Muscle 3.8.31 on the EMBL-EBI website³ under default parameters and alignments were adjusted manually in BioEdit v7.1.3.0 (Hall, 1999). The sequences were trimmed at both the ends to remove the primer sequences and the variable sites were reconfirmed using the original trace files.

The five candidate DNA barcode markers and their 20 possible combinations in multigene and tiered barcoding approaches were evaluated based on the methods proposed by the consortium for the barcode of life (CBOL). Candidate diagnostic nucleotides to classify sequences from specimens to species using a set of classification rules were identified using character-based machine learning approach in the program BLOG2.0 (Bertolazzi et al., 2009; Weitschek et al., 2013). The different barcode datasets used in this study were subjected to 90% slicing within species-level. The maximum iterations was set to 500 (GRASPITER = 500) with maximum given time of 5 min for complete analysis (GRASPSECS = 300). Among the derived logic formula sets, the one with the lowest false positive rate against the reference dataset was used as identification basis in Senna.

Bayesian-inference (BI) analysis was performed on ITS marker using MrBayes v.3.2.2 (Ronquist et al., 2012) at the CIPRES Science Gateway⁴. The best fit substitution model GTR+G was selected for analysis using the jModelTest v2.1.7 (Posada, 2008). The Bayesian analysis with the metropolis-coupled Markov Chain Monte Carlo (MCMC) was run for 10,000,000 generations, saving at every 1000th generation (the first 25% of trees were discarded as burn-in). The generations were checked until the average deviation of split frequencies reached under 0.01 and the potential scale reduction factor (PSRF) for all parameters approached 1.0. Convergence of runs was assessed using Tracer v. 1.6 (Rambaut et al., 2014) to generate a consensus tree with Bayesian posterior probabilities (PP) values. The values ≥0.95 were considered and included for each marker and concatenated topologies. To verify the results of the HRM analysis, the neighbour joining (NJ) tree-based method was used for species identification analyses. The NJ tree was constructed using MEGA6.0 with Kimura-2 parameter (K2P) model. The reliability of each node was assessed by performing a bootstrap analysis set to 1000 pseudo-replicates (Felsenstein, 1988).

Based on a previous study (Mishra et al., 2016b) and the results from the present dataset, ITS1 barcode reflected significant sequence divergence among the seven Senna species and were thus selected for subsequent HRM analysis. Conserved regions flanking the variable sites were identified visually and HRM-suited potential primers were designed with AlleleID (v.7.7, Premier Biosoft International, Palo Alto, CA, United States) using its SYBR Green Design function. The target amplification was set between 100 and 250 bases and primer length was set to 15 and 30 bases with an estimated melting temperature of 55.0 ± 5.0°C. The internal primers were screened for optimal coverage of polymorphic sites and conserved sequence at primer sites. The designed primers were analyzed using the IDT OligoAnalyzerv3.1 tools⁵ for the primer properties, hairpins and self/cross hybridization. The primers were verified using the Primer-BLAST NCBI to ensure specificity and were synthesized through Eurofins Analytical Services India Pvt. Ltd. (Bangalore, India).

Real-time PCR amplification followed by DNA melting and fluorescence measurements was performed on the StepOnePlus^TM Real-Time PCR System (Applied Biosystems, Inc., Foster City, CA, United States). The real-time PCR assays consisted of 5 μl of 2x MeltDoctor HRM Master Mix (Applied Biosystems, Inc., Foster City, CA, United States), 0.2 μl of 10 mM forward and reverse primers, 1 μl (1–50 ng) of genomic DNA and DNAse-free water to make up the final volume to 10 μl. Positive (containing a known amount of genomic DNA from each species) and negative controls were included. PCR runs were conducted using an initial denaturing step at 95°C for 5 min followed by 35 cycles of 95°C for 30 s, 57–60°C for 30 s and 72°C for 20 s, then a final extension step of 72°C for 2 min with collection of fluorescence signal at the end of each cycle. For HRM analysis, the PCR products were denatured at 94°C for 10 s and then annealed at 50°C for 15 s to randomly form DNA duplexes. The melting analysis was performed with the temperature increasing from 60 to 95°C at a ramp rate of 0.1°C/s. Fluorescence data were acquired at the end of each melting phase and processed using High Resolution Melt Software v3.0 (Applied Biosystems, Inc., Foster City, CA, United States). All reactions were done in triplicates and reference samples for HRM profiles analysis were included. The samples providing cycle threshold (C_t) values below 30 were considered suitable for HRM analysis. The -dF/dT (negative derivative of fluorescence F over temperature T) curve were plotted to derive the characteristic T_m for every species; the normalized raw curve depicted the decrease in fluorescence with increasing temperature. To normalize the raw melting curves, pre- and post-melt normalization regions were set and adjusted to define the temperature boundaries of the used plots. The characteristic melting temperatures (T_m) were recorded for each species of the Senna clade. S. alexandrina was set as a reference species. The developed Bar-HRM method was authenticated on the commercial Senna samples for detection of constituent species. DNA extracted from the raw drug powder of S. alexandrina was pooled with DNA of other six Senna species, viz. S. italica subsp. micrantha, S. spectabilis, S. auriculata, S. uniflora, S. italica subsp. Italica, and S. tora in percentage of 4, 8, 12, 25, and 50% w/w and tested for their specificity.

Seven potential species of the genus Senna representing a total of 21 individuals were successfully amplified and sequenced using five DNA barcodes, viz. rbcL, matK, psbA-trnH, ITS1, and ITS2 with 100% PCR and sequencing success rate (Table 1). The present study generated 105 new sequences, which were submitted to BOLD database and GenBank (Supplementary Table 1) and their diagnostic characteristics have been tabulated in Table 1. The PCR amplicons of all the five barcodes ranged according to the average size of the respective marker (Table 2). NCBI BLAST hits of all the seven species shared maximum similarity ∼98–100% with other species and to other genera of Cassiinae group. ITS sequences ranged from 607 to 721 bases with 738 aligned sites comprising of 304 variable and 303 parsimony informative sites. The sequences of complete ITS regions were annotated and trimmed to the regions of ITS1, 5.8S and ITS2. The sequence length of ITS2 ranged from 457 to 462 bases with 474 aligned sites. The numbers of variable and parsimony informative sites were 154. Both the markers showed several indels in the range of 1–4 bases within the aligned region. The coding regions rbcL and matK were highly conserved with 682/705 sites 726/785 sites, respectively, and were without indels. The intergenic spacer psbA-trnH showed high sequence length variation with 341–384 bases and the number of aligned sites were 437. The aligned region comprised several indels of 1–4 bases and 359 conserved sites.

BOLD-based genetic distance analysis of the four barcode regions matK, ITS1, ITS2 and psbA-trnH revealed that mean intra-specific distances were less than the distance to nearest neighbour (NN) of each Senna species. Individuals of S. italica subsp. italica exhibited maximum intra-specific divergence of 0.77–16.03%. The highest divergence range recorded is also due to inclusion of the individuals of subsp. micrantha. S. auriculata and S. tora share a maximum identity of 0.14% in terms of nearest neighbor (Table 3). Among the single barcode markers, only the ITS region resulted in the existence of a clear barcode gap, which is ideal for species identification (Figure 1). Besides, the two-barcode combination increased the resolution percentage in most of the tested combinations. S. alexandrina and S. italica subsp. micrantha+italica recorded the lowest NN distance of 5.74% among them, which reflects the high genomic similarity among them. The nuclear region ITS1 and ITS2 showed maximum inter-specific distances (4.97–19.34%) among all the seven species, thereby qualifying as potent marker for discriminating the species of genus Senna. However, the coding region rbcL showed only 0.14–0.43% divergence with any of the nearest neighbor among all seven species (Table 3). Figure 1 depicts the scatter plot of the maximum intra-specific distances against the NN distances to confirm the existence and magnitude of the barcode gap for all five candidate barcodes. Maximum intra-specific distances were less than 2% in all the species, except S. spectabilis (2.73%) and S. italica subsp. italica (16.03%). Based on the utility of individual markers, ITS1 and psbA-trnH were the favorable choices in the genus Senna and the regions were combined with other markers to assess their resolution rate (Table 4).

All five candidate barcodes and their 20 possible combinations were compared with TaxonDNA and BLOG to conclude on their discrimination percentage. All the barcoding datasets represented equal number of individuals corresponding to respective species and the rates of correctly identified, misidentified and not identified percentage were recorded for each datasets. Averaged over both the methods, ITS1 and psbA-trnH produced the highest success rate (100%) among the tested single barcodes (Table 4). The coding regions rbcL and matK resulted in very poor discrimination success for single as well the two-combination barcodes. However, supplementing the combination of rbcL+matK with the non-coding ITS1region, itself being at the first position, afforded 100% correct identification (100/0/0) through both TaxonDNA and BLOG based methods. Thus the tiered approach of barcoding proved a promising way to barcode the species of the genus Senna.

For estimating the evolutionary divergences among the species of the genus Senna, we employed character-based methods on all the barcode regions carried out using the BI model in MrBayes. The consideration of barcode marker based on the computational phylogenetics depicted similar hypothesis in agreement with the method based on the presence of barcoding gap. As shown in the phylogenetic tree the most favorable barcode dataset ITS1 presented highest level of discrimination at the species level. Bootstrap values clustered at the species level ranged from 64 to 100%. In particular the individuals belonging to same species maintained the species monophyly reflecting no intra-specific divergence among them. S. alexandrina and S. italica subsp. micrantha framed to be 100% similar (Supplementary Figure 2), which is consistent with the results from the BLAST analysis. Also the market samples of raw drug of S. alexandrina showed the presence of S. italica subsp. micrantha, when analyzed at the molecular level (Figure 2). Based on the indel polymorphism, we obtained marker nucleotides that can be used for species discrimination. Species-specific SNP positions and indel mutations detected in the ITS regions (total 143 sites) allowed the simultaneous discrimination of S. italica subsp. italica and S. alexandrina with transition mutations C ↔T and G ↔ A at 448 bp and 584 bp positions in the aligned length of antisense primer region. At the large polytomy, all the sister species of S. alexandrina shared maximum identity with 77–100% posterior probabilities. Individuals of S. italica subsp. italica framed the nodal cluster at the base of the tree (Supplementary Figure 2).

Based on the ability and performance of ITS1 DNA barcodes, the developed ITS-HRM primers (Table 2 and Supplementary Figure 3) that allowed the successful discrimination of the S. alexandrina from its allied species S. italica subsp. micrantha, S. spectabilis, S. auriculata, S. uniflora, S. italica subsp. Italica, and S. tora (Figure 3A and Table 5). Negative controls did not indicate any presence of plant species template. All the replicate specimens per species resulted in nearly similar C_t values and shapes of melting curves, which were further confirmed by the sequencing of the end products. S. alexandrina and S. italica subsp. micrantha depicted similar C_t values resulting in similar amplification efficiencies among species. Further the identifications derived by the HRM assays were confirmed by sequencing of the region of all the analyzed DNA extracts. The corresponding sequences were submitted to GenBank under the accession numbers listed in Supplementary Table 1, and the interpretations of both the sequencing and HRM profiles are found in Table 5. Both the techniques concluded 80–100% match concerning the targeted species. For samples SIM026 and SII031, which are difficult to discriminate visually by morphological classifications, the developed melting profile assigned them as the correct plant species consistent with the sequencing results (Figure 3).

Aligned sequences of seven ITS markers were used to screen the suitable primer sets for HRM analysis. As a result we were able to find four primer sets specific to species S. italica subsp. micrantha (SIM026HRMF/R), S. uniflora (SU029HRMF/R), S. tora (ST030HRMF/R) and S. italica subsp. italica (SII031HRMF/R), which could amplify 79 bp, 78 bp, 112 bp (Supplementary Figure 3), and 97 bp amplicons, respectively. The three best primer sets SIM026HRMF/R, ST030HRMF/R, and SII031HRMF/R were able to amplify with all the targeted species, while the primer set SU029HRMF/R was found to work for S. uniflora, S. auriculata, and S. tora, only. Thus the previous three primer sets were found to be suitable for the HRM analysis in Senna due to their universality among the species. The amplified PCR product of these primer sets were found plausible with further evaluation, which contains the variations of nucleotide sequences when compared within the target species.

Our HRM analysis using the primer pair ST030HRMF/R distinguished between the species of Senna clade with characteristic T_m values recorded for each species (Figures 3A,B). S. italica subsp. italica differed in the PCR product size from the other six species due to the presence of indels (at positions 345 and 346) among them, resulting in a T_m value (79.45°C) different from that of S. alexandrina (T_m 78.80°C). Also the change of nucleotides from T to C between them (Figure 3C) resulting in their differences in T_m. Amplicon from S. spectabilis and S. uniflora showed the lowest T_m value difference (79.62°C and 79.02°C, respectively) due to almost identically sized PCR products. The change of nucleotides C to A at position 355 in the forward primer region in the species S. alexandrina and S. italica subsp. micrantha from the rest of the species, resulted in almost identical T_m values (78.80°C and 78.77°C, respectively) for both the species. To better visualize the small difference between the individual melting curves, HRM software was used to calculate a difference plot for each species (Figure 3B). S. alexandrina was used as a reference species for genotyping with its melting curve as the baseline. The difference was obtained by subtracting the difference graph area of the reference species from the rest of the species derived melting curve. The genotype confidence level was measured on the cut-off value of 90% to assign a specific genotype for each barcode region. Thus the ITS1 barcode coupled to HRM primer pair amplified all the seven species with sufficient discrimination and confidence level (Figure 3B).

Herbal plant identification becomes more challenging when the plants are in processed or dry form, which often is due to the high morphological similarity of the drug part with respect to its adulterants (Supplementary Figure 4). To test the blend of other species in S. alexandrina samples, HRM coupled with the ITS1 barcode method was employed on the pooled DNA of S. alexandrina contaminated with the other six Senna species. The contamination percentage was measured as 1–50%. Consequently, when the sample of S. alexandrina was mixed with the S. italica subsp. micrantha species in range of 1–10%, the T_m values of PCR product started deviating gradually toward S. italica subsp. micrantha species (Figure 4A). The identification of the contaminant was not easy when mixed in range of 1–5%. However, the differentiation slightly increased above 6%, to detect the admixture of other species in the market samples of S. alexandrina. The limit of detection was recorded between 8 and 12%. Besides, 10 commercial samples (HSA01–HSA10) from local markets were also tested with the method to confirm their specificity. T_m values of PCR product from samples HSA01 and HSA06 were found to be in same range (Figure 4B). In contrast, both the samples presented similar melting curve profiles and cannot be differentiated further. Both the Senna samples were found to be contaminated with plant species S. italica subsp. micrantha. Sample HSA08 and HSA09 had T_m values almost identical to the T_m value recorded in the field samples of S. alexandrina. Thus both the samples were putatively found to be S. alexandrina. Sample HSA08 was found to be contaminated with S. tora, with T_m values of 79.30°C. Thus the present investigation resulted in the development of a robust and rapid method to detect the admixture of S. alexandrina raw drugs contaminated with other Senna species, based on their unique melting curves (Figure 4A).

The findings from the HRM results were further validated with bidirectional sequencing of the HRM-PCR products. All the obtained sequences were checked for their similarity with the BOLD database. The BLAST analysis of the HSA sample sequences were found to be in congruence with the results of the HRM analysis. A phylogenetic tree was constructed on the sequence dataset employing the NJ cluster algorithm (Figure 2). HSA04 and HSA06 were found to be highly contaminated with S. italica subsp. micrantha. HSA03 was mixed up with S. tora and S. uniflora with 60% contamination. HSA01, 02, 08, and 09 were the putative S. alexandrina species, while HSA05, 07, and 10 were found to be clustered with S. italica subsp. micrantha. The corresponding sequences have been submitted to DNA barcoding database BOLD for species authentication. These results demonstrate the usefulness of ITS1 barcode in species authentication in the genus Senna.