Feasibility of nuclear ribosomal region ITS1 over ITS2 in barcoding taxonomically challenging genera of subtribe Cassiinae (Fabaceae)

Taxon sampling, DNA amplification and sequencing

A total of 54 accessions of 18 species belonging to three genera viz. Cassia, Senna, and Chamaecrista from India were examined during the study. For obtaining the sequences generated from molecular experiments in our lab, a total of 18 individuals corresponding to three different genera were collected from different geographical regions of South Western Ghats and Uttar Pradesh. The species were identified and authenticated using the morphological characters described in a monographic study on Cassiinae in India (Singh, 2001) by Dr. V. Sundaresan, Senior Scientist, Central Institute of Medicinal and Aromatic Plants, Research Centre (Bangalore). For each of the species, herbarium specimens were prepared and deposited at the Herbaria of the Central Institute of Medicinal and Aromatic Plants (Table 1, CIMAP Communication No.: CIMAP/PUB/2016/24), Lucknow.

Legumes family produce a high diversity of secondary metabolites, which causes extreme difficulty in isolation of high-quality nucleic acids. Based on literature and commercial kits available, we attempted modification of several previously reported methods to isolate high quality DNA. Ultimately, total genomic DNA from individual accessions was extracted from the leaf tissues (dried in silica-gel) using the modified cetyl trimethyl ammonium bromide (CTAB) protocol with necessary major modifications (Khanuja et al., 1999) and supplementing it with the Nucleospin Plant II Maxi prep kit using the manufacturer’s protocol (MACHEREY-NAGEL, Duren Germany). The concentration of β-mercaptoethanol and PVP (Polyvinylpyrrolidone) were increased to 2% v/v and 4% w/v, respectively. An additional chloroform-isoamyl alcohol (96:4) purification step was performed to remove proteins and potentially interfering secondary metabolites. Isolated DNA was checked for its quality and quantity by electrophoresis on a 0.8% agarose gel and spectrophotometric analysis (NanoDrop, ND-1000, USA). The nuclear internal transcribed spacer (ITS1 and ITS2) regions of all the individuals were amplified according to PCR reaction conditions (94 °C, 5 min; (30 cycles: 94 °C, 1 min; 50 °C, 1 min; 72 °C, 1.5 min); 72 °C, 7 min) following guidelines from the CBOL plant-working group and sequenced using universal primers ITS5a forward 5^′-CCTTATCATTTAGAGGAAGGAG-3^′ and ITS4 reverse 5^′-TCCTCCGCTTATTGATATGC-3^′ (Kress et al., 2005). PCR amplifications for each primer set were carried out in a 50 µl volume solution containing 1x Taq DNA polymerase buffer, 200 µM each dNTPs (dATP:dTTP:dCTP:dGTP in 1:1:1:1 parts), 10 pmol of each primer (forward and reverse), 1 unit of Taq DNA polymerase and ∼25–50 ng of template DNA. The PCR fragment lengths were determined on a 2% agarose gel. The PCR products were purified with Nucleospin PCR purification kit (MACHEREY-NAGEL, Duren, Germany) as per the manufacturer’s instructions. Presence of the specific product was confirmed by running the purified PCR products on 2% agarose gel. All the purified PCR products were subjected to double-stranded sequencing using the Applied Biosystems Big Dye Terminator Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) on an ABI 3130 XL automated sequencer (Applied Biosystems).

Apart from the lab-generated sequences, all the nucleotide sequences belonging to genera Cassia, Senna, and Chamaecrista for the regions ITS1 and ITS2 were downloaded from the NCBI based on the BLAST results. The sequences were filtered on the basis of length (less than 300 bp were omitted), lack of voucher specimens as well as verification (sequences categorised as unverified in GenBank were omitted). An effort was made to include minimum five individuals for each species, but due to unavailability of sequences for few species in the NCBI database and difficulty in obtaining the species in the field, the representatives of each species were limited to three. The GenBank accession numbers used in this study are listed in Table 1.

Data analysis

Electropherograms corresponding to raw sequences of individual accessions from both the forward and reverse primers were assembled and edited using CodonCode Aligner v.3.0.1 (CodonCode Corporation). Sequences were clipped at the end to avoid the presence of variable sites introduced by the sequencing artefacts. Due to its well-conserved nature, the 5.8S gene region was removed from any sequence so that the ITS1 and ITS2 regions could be analyzed separately and concatenated. The edited sequences were then aligned with MUSCLE 3.8.31 on the EMBLEBI website (http://www.ebi.ac.uk) with default parameter and adjusted manually in BioEdit v7.1.3.0 (Hall, 1999). All the variable sites were rechecked on the original trace files. To evaluate the effectiveness of ITS1, ITS2 and their combination (ITS1 + 2) as barcodes in the concerned genera, three widely used methods viz. distance-based (PWG-distance), similarity-based and tree-based were applied.

Genetic distance-based method

To evaluate the measure of effective barcode locus, DNA barcoding gap was calculated using TaxonDNA software with a ‘pairwise summary’ function under K2P nucleotide substitution model (Meier et al., 2006). The pairwise genetic distance were calculated at the observed levels of intra- and inter-specific divergence for each barcode. To test the accurate species assignments, the distributions of the pairwise intra- and inter-specific distances with 0.005 distance intervals were generated. The histogram of distances vs. abundance were plotted to estimate the presence of any barcoding gaps. For the PWG-distance method, the genetic pairwise distance was estimated by MEGA version 6 (Tamura et al., 2013) using the Kimura two-parameter distance model (K2P) with pairwise deletion of missing sites (Kimura, 1980). Average inter-specific distance was used to characterize inter-specific divergence (Meyer & Paulay, 2005; Meier, Zhang & Ali, 2008) and ‘all’ intra-specific distance, mean ‘theta’ and coalescent depth were used to characterize intra-specific distances. Finally, the obtained inter- and intra-specific distances were plotted with frequency distribution in bin interval of 0.05 to illustrate the existing DNA barcoding gap (Meyer & Paulay, 2005; Lahaye et al., 2008).

DNA sequence similarity-based method

To test the potentiality of ITS regions to identify species accurately based on sequence similarity, the proportion of correct identifications were calculated using SpeciesIdentifier program from the TAXONDNA software package with ‘Best match’ (BM), ‘Best close match’ (BCM) and ‘All species barcodes’ functions (Meyer & Paulay, 2005). The tool examines all the sequences present in aligned data set and compares each successive sequence with all the other sequences to determine the closest match. The ‘Best match’ modules than classifies the sequences as correct and incorrect based on the indicated pair from the similar species or different species respectively. While the various equally best matches from different species are referred to be as ambiguous. The ‘Best close match’ module works on the intra-species variability criterion and considered to be the more rigorous method in TaxonDNA. The sequences classified as ‘no match’ are the results above the calculated threshold value (Meier et al., 2006).

Tree-based method

To evaluate the ability of candidate barcode to delimit the species into discrete clades or monophyletic groups, three different optimality criteria (tree-building method) viz Neighbour-joining with minimum evolution (NJ), maximum likelihood (ML) and Bayesian inference (BI) were employed. To test the reliability of the result, NJ and ML trees were constructed and compared with two different softwares: (i) In MEGA using the K2P distance as model of substitution (Tamura et al., 2013) and (ii) In PAUP 4.0 with the HKY-gamma substitution model (Swofford, 2003). The reliability of the node was assessed by a bootstrap test with 1,000 pseudo-replicates with the K2P distance options (Felsenstein, 1988). Bayesian sampling was performed in BEAST1.8 using the operators: HKY substitution model with four gamma categories, a constant-rate Yule tree prior and 10,000 chain lengths and all other priors and operators with the default settings. Coalescent tree priors were used for population-level analysis and speciation prior were applied to estimate relationships and divergence times of inter-species data. Trees were sampled for every 5,000 generations resulting in a total of 10,000 trees, and a burn-in of 50,00,000. Beast file was created using the BEAUti program v1.8.2 within Beast and performance of each run was further analysed with the program Tracer (Rambaut et al., 2014). The resulting Beast tree files were annotated through TreeAnnotator v1.8.2 and visualized and edited with FigTree v1.4.2. (Rambaut, 2014, http://tree.bio.ed.ac.uk/software/figtree). Visualization and analysis of all the resulting trees through PAUP 4.0 was done in Dendroscope3 (Huson & Scornavacca, 2012). Gaps were treated as missing data for all the phylogenetic analysis.

PCR amplification and sequence characteristics

The sequence characteristics of ITS regions evaluated in this study showed good success rates (90%) for PCR amplification (ranging from 571-1153 bp with mean size ∼707 bp; gel images can bé provided on request) and sequencing in both the direction using a single primer pair ITS5a forward and ITS4 reverse. The presence of large amount of secondary metabolites, polysaccharides and polyphenolic compounds in the plants of sub-family Caesalpinioidae, hindered the isolation of pure nucleic acids. Therefore few samples had to be excluded from the study after 3–4 initial amplification attempts that failed due to the presence of inhibitory components. The present study generated 15 new sequences belonging to 15 different species of Cassia, Senna, and Chamaecrista. The sequences were submitted to NCBI (www.ncbi.nlm.nih.gov/genbank/) and corresponding GenBank accession numbers were obtained for each species. A total of 64 sequences corresponding to 18 different species of Cassia, Senna, and Chamaecrista for ITS regions (ITS1 and ITS2) were obtained from NCBI and included in the study (Table 1). The ITS1 region had an aligned length of 315 bp (Alignment S1) which was greater than that of ITS2 with 258 bp (Table 2; Alignment S2). The combined region ITS1+2 showed an align length of 573 bp (Alignment S3) with 80.1 % of pairwise identity (Table 2). The aligned ITS1 matrix consisted of 315 bp with 206 parsimony sites. The number of variable sites was 210. The maximum intra-specific divergence was observed among the individuals of Senna siamea with 0.023 PWG-distance while minimum inter-specific distances were recorded between Senna hirsuta and Senna occidentalis with 0.039 PWG-distance. The species of genus Chamaecrista showed lowest K2P distances (Table 3). Overall the summary statistics for DNA alignments and DNA sequences for the ITS dataset evaluated in this study are summarized in Tables 2 and 3 respectively.

Genetic divergence and Barcoding gap

The presence of DNA barcoding gap based on the concept of an inter-specific distance being larger than the intra-specific distance for a species, directly reveals the species discrimination ability of candidate barcodes. In this study, the relative distribution of frequencies of K2P distances for three ITS datasets using TaxonDNA software showed a significant pattern with the inter-specific distance being higher and did not fully overlap with the intra-specific distance resulting in the presence of an identified barcoding gap in the genera. The observed pattern of ITS1, ITS2 and ITS1 + 2 results are presented in Fig. 1. The mean intra- and inter-specific genetic divergence based on PWG distances through MEGA, for ITS1 varied in the range from 0.023 to 0.000 and 0.033–1.185 respectively (Table 3).

Species discrimination based on different analytical methods

In accordance with the CBOL PWG-distance method, a favourable barcode should possess a high inter-specific divergence to distinguish different species. The result obtained through the different datasets showed significant pattern of inter-specific divergence, whereby ITS1 was concluded to be the best among the candidates. The mean pairwise inter-specific distances were found to be higher in comparison to intra-specific distances in all the barcodes, resulting in the presence of a clear barcode gap. The distance distribution range of all inter- and intra-specific distances for all markers are shown in Fig. 2.

Compared with the PWG- distance method, the BM and BCM functions of TaxonDNA showed the better discrimination success. All the three datasets presented same success rate of species identification when BM was selected in comparison to BCM. The highest and same rate of discriminatory power (81.6%) was observed for ITS1 on both BM and BCM functions. The other two datasets; ITS2 and ITS1 + 2 datasets recovered 75.0% and 77.4% BM respectively (Table 4).

The tree building methods for the evaluation of barcode sequences were estimated based on the correct assignment of individuals forming a monophyletic clade (Fig. 3 and Fig. S1). Among the different phylogenetic methods, BI recovered the highest value for species monophyly in all the datasets. While in the combination of ITS1 + 2, all the three methods viz. NJ, ML and BI provided near similar topology, concluding 77.41% of individuals identified correctly (Fig. 4). The resulting bootstrap value lends support to our findings. Comparing the potentiality of the ITS datasets and the phylogenetic algorithms employed, the highest discriminatory power was observed when ITS1 was used alone, which successfully maintained the genera (Cassia, Senna, and Chamaecrista) monophyly with few exceptions (Fig. 5). The coalescent and speciation tree priors intrinsically correlated the rate of evolution and time in inferring genetic differences between species. It is interesting to conclude that all the species from genera Senna and Cassia framed in two different clusters viz. Cluster I and II according to traditional morphology. The phylogenetic tree presented a slight divergence in the clustering of Chamaecrista absus accession obtained from GenBank which might be due to the mis-identification of samples. Referring to the species relationships within genera; to some extent, the phylogenetic relationships obtained were in consistent with the result obtained from the traditional morphological classification method. The clustering pattern of three different genera Cassia, Senna, and Chamaecrista within the subtribe Cassinae based on the nuclear ribosomal region ITS1, proved to be successful in comparison to the infrageneric clustering of taxa. The clustering of Senna tora, Senna uniflora and Senna obtusifolia accessions based on molecular algorithm of ITS1 complies with the morphological similarity occurs among them, while in ITS2, Senna uniflora showed little divergence (Fig. 3). Also, we were not able to find out the clear pattern of lineage of respective species within the genus at a molecular level, as according to traditional taxonomy. Worthy to note here, that the resulting pattern within the individuals of same species and high reliability value obtained for their nodes concludes the existence of genetic similarity among them. Framing of Senna occidentalis and Senna hirsuta into the individual cluster through ITS1, were in consistent with the key classification (Fig. 3).

Besides, all the tree species belonging to genus Cassia, undertaken in this study framed an individual cluster (Cluster II) according to their diversity there by concluding the importance of molecular characterization in corroboration with morphological methods in biosystematics study. The analysis conducted in subtribe Cassiinae with the tree based, similarity based and distance based methods showed that BI phylogenetic method and BM similarity methods outperformed the PWG- distance method when using these barcode loci (Fig. 4).

Discrimination success

Hitherto several different analytical methods were framed for the assessment of the species discrimination ability, which includes tree-based (NJ, MP, Bayesian), distance-based (PWG-distance, p-distance, K2P-distance) and sequence similarity-based methods (BLAST and TaxonDNA), etc., and all of them show different discrimination power on the same data set (Little & Stevenson, 2007; Austerlitz et al., 2009; China Plant BOL Group, 2011; Sandionigi et al., 2012). In this study, sequence analysis of ITS datasets using Bayesian inference (BI) tree-based method gave the highest species resolution based on the topology with the highest product of posterior clade probabilities across all nodes followed by BM and BCM model of TaxonDNA, which too presented equally efficient results either in single or combination of barcodes. Similarly, patterned results have been obtained in different DNA barcoding studies in various plant groups (Yan et al., 2014; Giudicelli, Mäder & De, 2015; Xu et al., 2015; Yan et al., 2015). The clustering algorithm of the Bayesian framework provides a flexible way to model rate variation and obtain reliable estimates of speciation times, provided the assumptions of the models be adequate (Drummond et al., 2012).

The PWG-distance method based on simple pairwise matching recommended by CBOL Plant Working Group as a universal and robust method for the assessment of clear barcoding gap indicated the significance of ITS1, thereby highest number of variable and informative sites (210 and 206, respectively) were obtained. Moreover, the rate of species discrimination is equally efficient when ITS1 and ITS2 are concatenated. These results were expected, considering the complexity of the genera and directly reflected on the performance of ITS1 and ITS2 as barcode markers in Cassia, Senna, and Chamaecrista. The possible reason behind the results might be the inter-specific sharing of identical sequences or failure of conspecific individuals to group together. Besides, many other aspects have also been reported for unclear barcoding gap such as imperfect taxonomy, inter-specific hybridization, paralogy and incomplete lineage sorting (Yan et al., 2015). However, ITS region has proved to be a suitable marker in authentication of Cassia species in the commercial herbal market (Seethapathy et al., 2014). The strong identification ability of nuclear region ITS have been verified in many complex groups (Baldwin et al., 1995; Alves et al., 2014; Wang et al., 2014; Giudicelli, Mäder & De, 2015). Therefore, we suggest that ITS1 itself could be the first option for DNA barcoding in subtribe Cassiinae, though ITS2 should not be discarded.

Moreover, the differences among the three methods compared here, have their possible cause in the theories behind their algorithms and the matter of comprehensive sampling. Thus the comparison of species resolution between studies without consideration of the methods should be avoided for one or the other reasons discussed, as species resolution is an important criterion for assessment of robust barcodes.

Biological implications of ITS based signalling in Cassiinae

The corroboration of morphological, ecological, geographical, reproductive biology and DNA sequence information paved the successful path for constructing robust taxonomy for diverged plant taxa (DeSalle, Egan & Siddall, 2005; Fazekas et al., 2009; Hollingsworth, Graham & Little, 2011). The ITS region appears to evolve more rapidly than coding regions in interpreting phylogenetic relationships at lower taxonomic levels (Inter-generic and Inter-specific). Species discrimination for the genera Cassia, Senna and Chamaecrista sampled in this study was high with the strong identification ability of nuclear region ITS. All the three genera maintained the monophyly of the clade either alone or in combination of barcoded loci. The resulting bootstrap value lends support to our findings. To some extent, the divergence of species within the genus did not outperformed as designated according to key taxonomy. The possible reasons behind the findings could be the complexity of the genus with large number of highly polymorphic species which has been found to devise greater interspecific variation (Mohanty et al., 2010). Sometimes interspecific hybridization and gene introgression had accounted for the limited barcoding event at genus level. Moreover genera Cassia and Senna accounts for high morphological complexity based on species polymorphism, which have been reported in few studies in the past. Successful PCR amplifications, sequencing strategy and alignment matrix obtained from the present study provided further evidence to support the separation of species and genera. The robust phylogenetic signalling of ITS region seems obvious in Cassinae. Although an earlier study (excluding ITS) did not report any single novel region to differentiate the existing Cassia species (Purushothaman et al., 2014), our findings provide the potentiality of the ITS region with data support. The delineation of genera based on ITS regions provided a basic framework to have an authentication prospect of correct species at the industrial level.