Molecular identification of herbal species belonging to genus Piper within family Piperaceae from northern Peninsular Malaysia

2.1 Plant collections

Plant samples from genus Piper were collected randomly from natural locations that represents three states in northern Peninsular Malaysia (Fig. 1). Sampling activities were conducted for consecutively three days every month during January 2014-February 2015. The leaves of a plant were excised using scissors and then stored in a sterile polyethylene bag. All samples were labelled to avoid misidentifications before transporting to School of Biological Sciences, Universiti Sains Malaysia (USM) for storage. All plant samples were identified based on morphological and taxonomical characteristics by local taxonomist from USM. The book titled Photographic Atlas of Botany and Guide to Plant Identification authored by Castner (2005) was also used to identify all plants collected in this study (Table 1). The vernacular name of all plant samples collected were identified by local people. All the document and preserved specimens were stored in the USM’s herbarium for record and references.

Close

Fig. 1

Sampling locations of all samples collected.

Export to PPT

2.2 DNA isolation and amplification

The leaves were dried in silica gel prior to DNA extraction. Approximately 300–400 mg fresh leaves were ground to a fine paste and homogenized in a DNA extraction buffer [(50 mM Tris HCL pH 8.0, 25 mM EDTA pH 8.0, 150 mM NaCl, 40.0 ml H2O) and 1 g PVP 40 mw 40 000]. Distilled water was added to the mixture up to 150 ml and HCl was used to accustom the pH to 5.0. The paste was subsequently conveyed to another sterile tube. The modified cetyl trimethylammonium bromide (CTAB) method (Cota-Sánchez et al., 2006) has been utilized for the plants total genomic DNA extraction. The purified DNA in all samples analysed was quantified using ethidium bromide (EtBr) stained band intensities against λ DNA in order to acquire a high quality DNA. This method was also used to confirm that the extracted DNA is clear from metabolites that might impede the Polymerase Chain Reaction (PCR) amplification process. The PCR method was used to amplify the target region of the gene in the nuclear or the plastid genome of plants depending on the primer pairs used. The set of primers used in the present study are 1) rpoC1-F, 5′- GGCAAAGAGGGAAGATTTCG- 3′, rpoC1-R, 5′- CCATAAGCATATCTTGAGTTGG- 3′ (Hollingsworth et al., 2009a) and 2) rbcL-F, 5′-ATGTCACCACAAACAGAGACTAAAGC-3′ (Levin et al., 2003), rbcL-R, 5′-GTAAAATCAAGTCCACCRCG-3′ (Kress et al., 2009).

The PCR reaction consisted of 10X PCR reaction buffer, 25 mM MgCl2, 1.25 mM of each dNTPs, 1 unit of Taq polymerase, 10 μM of each primer, 20 ng genomic DNA and 30 μl Milli-Q water. Thermal cycling conditions (on a T100TM Thermal Cycler; BioRad, Singapore) were 35 × [94 °C for 45 s, 40 °C − 50 °C (depends on the primer used) for 45 s, 72 °C for 1 min] and a final incubation at 72 °C for 10 min. Additional purification was performed on all samples utilizing the Qiagen DNA Mini Kit (Germany) in order to ensure the end products free from contaminants and/or other PCR suppressor). Purified PCR products were then forwarded to NHK Laboratory Inc. (Malaysia) for sequencing purposes. At NHK Laboratory, DNA sequencing was implemented using a BigDye^TM Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) and used around 20–30 ng of purified and cleaned PCR mixtures with the set of specific primer in the analysis. The end mixtures were cleaned up prior electrophoresed on an ABI PRISM 3130xl Genetic Analyzer according to the standard protocols.

2.3 Sequence alignments and analyses

The amplified sequence of rbcL and rpoC1 were blast in GenBank utilizing BLASTn analysis and implemented in the GenBank web interface (http://blast.ncbi.nlm.nih.gov). The previously published sequence of the genus Piper in GenBank was also used as a reference in this study. The sequences were selected based on the availability of voucher specimens deposited in GenBank and the length of the sequences (>300 bp). All sequences were collapsed in haplotypes implementing the Collapse software version 1.2 (Provan et al., 2005). All haplotypes were then aligned with MUSCLE (Edgar, 2004) with default parameters and all aligned sequences was then manually synchronized and trimmed in Mega X version 10.1 (Kumar et al., 2018). Alignments were then manually revised in an attempt to minimize the positional dissimilarity. All missing data and gaps within the sequences were removed.

2.4 Species partitioning and analysis of genetic divergence

Each sample was determined for correct identification and employed in the Automatic Barcode Gap Discovery method (ABGD, Puillandre et al., 2012) in an effort to analyse the species accurately. ABGD is an automatic procedure that is based on the genetic distance method to determine a pause/gap for barcode that will partitions a putative species in the dataset by confirming that intra and inter specific genetic distance do not overlap. This tool quantifies all pairwise distances and arranges the values in rank. ABGD engages a two-stage approaches which primarily fraction DNA segment into Operational Taxonomic Units (OTUs) according to a statistically deduced barcode gap (e.g. primary partitioning), and afterwards employs a second cycle of partitioning (for e.g. recursive partitioning). There are three crucial parameters in ABGD namely (1) X , relative gap width estimate, (2) minimum and (3) maximum values of prior intraspecific divergence, P that are important to determine barcode gap. In this research, the values for the prior P (prior maximum divergence of intraspecific diversity) has been setup ranging from 0.001 to 0.1 and X  = 1.0. The whole data set will be interpreted as one species if the P value was set too high (Puillandre et al., 2012). The genetic distance examination was quantified based on Jukes-Cantor 69 (Jukes and Cantor, 1969) analysis and performed at the web interface (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html).

To retrieve the constructiveness of marker discriminatory capability, the emergence of monophyletic succession and inheritance was determined utilizing a phylogenetic tree analysis. For likelihood based analysis, the best model of substitution was resolved following the Akaike Information Criterion (AIC) and employed in MrModeltest 2.4 (Nylander, 2004). For both genes (rbcL and rpoC1), MrModelTest found the best model is the Kimura 2-parameter (Kimura, 1980). On the basis of the model, the Maximum Likelihood (ML) tree was adopted utilizing the MEGA X version 10.1 (Kumar et al., 2018) software with a Maximum Composite Likelihood (MCL) approach. The MCL is a total of related log-likelihoods in which this method was used to accurately determine the pairwise distance and related substitution parameters (Tamura et al., 2004). The evolutionary relationships of all samples were also performed employing the Neighbour Joining (NJ) approach with the same model as in ML analysis. In order to make the phylogenetic tree more comprehensive and informative, other sequences from different genera namely Peperomia tetraphylla, Zippelia begoniifolia and Verhuellia lunaria were included in the analysis for rbcL gene. However, only a single sequence (Peperomia pellucida) was included in the rpoC1 phylogenetic analysis as no reference sequences found for the genus Zippelia and Verhuellia in GenBank. Asarum caudigerum (Genebank accession no.:JF940914.1) and Asarum forbesii (Genebank accession no.:GQ436063.1) are low-growing herbs distributed mainly in Asia and were used to root the NJ and MP tree regenerations in this study. The estimate of evolutionary divergence (genetic distance, D_s) between and within samples was computed based on Maximum Composite Likelihood and implemented in MEGA X version 10.1 (Kumar et al., 2018).

In order to test the validity of phylogenetic deduction and examinations of evolutionary postulate, we conduct a Disparity Index (DI) analysis. This was implemented for all samples collected (not includes sequence retrieved from GenBank). A greater variances in base structure bias than expected will be indicated by a value more than 0 and this will be quantified on the basis of the evolutionary deviation between sequences and by chance alone. The analysis was employed in MEGA X version 10.1 (Kumar et al., 2018).

3.1 Sampling data and sequence information

A total of 11 fresh plants from genus Piper were successfully collected from six different locations in northern Peninsular Malaysia (Fig. 1). Samples were collected from Alor Setar (n = 2) and Baling (n = 2) in Kedah, Padang Besar (n = 1) and Kangar (n = 2) in Perlis, Balik Pulau (n = 2) and Bayan Lepas (n = 2) from Penang (Table 1). Additionally, based on the selection criteria mentioned earlier, a total of 22 (rbcL) and 10 (rpoC1) sequences of the genus Piper (including outgroup sequences) were also retrieved from GenBank. The aligned sequence length varied from 488 bp (rbcL) to 418 bp for rpoC1. The rbcL gene performed a full outcome (100%) in amplification reactions, but rpoC1 gene shows only 91% successful amplification and exhibited lower reaction efficiency compared to rbcL gene. This is as expected because one sample needs approximately two trials to confirm the volume of PCR reagents (i.e., DNA template, dNTPs, Taq DNA Polymerase), which exhibited a notable variability in rpoC1 amplification success among different plant samples. All the amplified sequences obtained from rbcL gene show high efficiency and quality sequences (100% sequencing success) while rpoC1 showed only 80% success.

The chloroplast region of rpoC1 was notably utilized as a genetic marker to discriminate many plants (Rydberg, 2010). However, in the present study, the rpoC1 region showed lower sequencing success when compared to rbcL region. An identical outcome was also reported by other researchers that recorded a low sequencing success for rpoC1 (see e.g., Sass et al., 2007; Rydberg, 2010; Hollingsworth, 2011; Tripathi et al., 2013). The interspecific variation for this locus is too low in some plant groups and the difficulty in amplifying and sequence has been identified as the main limitation for rpoC1 locus (Sass et al., 2007). This pitfall consequently has led to a disputation for contemplating rpoC1 as not an effective and favourable standard gene marker for plants (Sass et al., 2007).

3.2 Species identification

BLASTn analysis for all sequences showed that all samples have been correctly identified up to a species level, demonstrating that all plant samples preliminary identified based on their morphological characteristic matched with the scientific names retrieved from the conspecific sequences deposited in GenBank (Table 2). Likewise, rbcL gene consistently shows the highest rate of success for species identification (99.8%) based on BLASTn analysis. Conversely, the identification for rpoC1 was notably low at both levels (genus and species) which is only 66.7% of all the collected samples. Based on these results, the present study reveals that rbcL is the reliable DNA marker for the Piper sp. collected from northern Peninsular Malaysia (Table 2).

The variation of names used for the same plants throughout different races, ethnics and geographic range is likely the most important reason for the mismatch in species identification by BLASTn analysis. (Mankga et al., 2013). For example, in Malaysia, the medicinal plant named ‘sirih’ or ‘sireh’ is mostly referred as Piper betle and/or Piper sarmentosum. The plant was also called ‘sirih’ in Indonesia while in the Philippines it was called ‘Ikmo’. Nagori et al. (2011) showed that P. betle exhibit high antioxidant activity and possessed antifungal, antiseptic and anthelmintic properties to serve as a contraceptive for humans. Moreover, consumption of food products contained P. betle could contribute to the additional antioxidant needed in the body to enhance defence system, as well as an additional nutraceutical supplement in patients diagnosed with rheumatism and leucorrhoea (Nagori et al., 2011). However, inaccurate identification of plants will compromise the therapeutic value of medicinal plants, thus would endanger human health. For example, adulterant Bunium cylindricum that has been homogenated with B. persicum and trade in the retail outlet in Iran have resulted in the deterioration of the superiority and effectuality of the products (Joharchi and Amiri, 2012). Many cases of toxicity have been reported and mostly due to species misidentification (Viljoen, 2013). The situation, where a proliferation of herbal remedies that have been adulterated or substituted with other plant materials, has stressed the need for a quality control (Raterta et al., 2014).

3.3 Species partitioning and discrimination

The interspecific genetic distances (D_s) were significantly higher (0.00–0.095) for samples amplified by rpoC1 gene as compare to rbcL (D_s = 0.000–0.023), manifesting that all samples analysed in this study were genetically unrelated (Table 3). In contrast, intraspecific genetic distances (D_s) were significantly low, ranging from 0.000 to 0.033 (rbcL) and 0.000–0.095 (rpoC1), signifying the exclusion of a cryptic species (Table 3). The validity of rbcL gene as the reliable barcode marker for Piper sp. was further supported by the ABGD analysis. The number of species defined by ABGD ranged with the change prior thresholds from 0.001 to 0.050 prior intraspecific divergences (P) for rbcL marker (Fig. 2A) and within the range of 0.000 to 0.470 for rpoC1 (Fig. 2B). For rbcL, the lowest threshold values varied from 0.00 to 0.03 and grouped all haplotypes as a different species (11 species). However, this is not the case for rpoC1 where the lowest threshold values (0.00–0.91) clustered all haplotypes into eight groups only. This has resulted in the same grouping as shown by the Maximum Likelihood (ML) and Neighbour-Joining (NJ) phylogenetic tree model of all species. Thus, the 11 and eight species clusters that were demarcated by the ABGD approach via the high prior threshold for rbcL and rpoC1 genes, respectively are corroborated by the outcome of ML and NJ analysis. The NJ tree shows a same topology as in ML analysis, thus only the ML phylogenetic tree was represented for both genes (Figs. 3 and 4).

Close

Fig. 2

Distribution of rbcL (A) and rpoC1 (B) pair-wise JC69 distance values performed by ABGD, showing the barcode gap between the putative maximum co-specific divergence and the minimum congeneric divergence.

Export to PPT

Close

Fig. 3

The evolutionary history of Piper sp. collected from northern Peninsular Malaysia (including sequences retrieved from GenBank) inferred by using the Maximum Likelihood method and Kimura 2-parameter model for rbcL gene. The tree with the highest log likelihood (-974.05) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches.

Export to PPT

Close

Fig. 4

The evolutionary history of Piper sp. collected from northern Peninsular Malaysia (including sequences retrieved from GenBank) inferred by using the Neighbour-Joining method and Kimura 2-parameter model for rpoC1 gene. The tree with the highest log likelihood (-1840.11) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. shows sample with missing corresponds sequence from GenBank. indicate samples with ambigous identification.

Export to PPT

The ML and NJ phylogenetic trees for each barcode marker used in this study were diverged each species by delineating a monophyletic clade, where all clades revealed evidence separate from other clades with a high bootstrap value (Figs. 3 and 4). For example, the clade containing P. auritum and P. umbellatum as displayed by rbcL (Fig. 3) and rpoC1 (Fig. 4) NJ phylogenetic trees were well supported. The rbcL genes demonstrated successful identification both at all levels. Nevertheless, the verification rate of rpoC1 barcode is only 66.7% (Table 2, Fig. 4). In this case, some of the species’ positions are paradoxical in which they were grouped within a different species and/or genus. For example, two samples namely AS2 and KA2 that were confirmed as Piper sarmentosum (based on BLASTn and rbcL gene) was largely diverged from its ancestor in which it was positioned with Piper protrusum (sequence retrieved from GenBank) (Fig. 4). Several samples, namely BA2, AS1, KA1, PB and BA1 were successfully amplified, but the identity of the samples was not resolved as no correspond sequence available from GenBank.

In phylogenetic reconstruction, one of the most notable difficulties that can generate unreliable in regards to the accurate evolutionary history of the organism or genes analysed is the extensive existence of inconsistency among approaches. This was also corroborated by the incompatibility of genes or genomic regions used to infer the phylogenetic output. Accordingly, such drawback hindered the principal aims of evolutionary research (Som, 2014). This is further exacerbated by the advancement of genomic and molecular approaches that has escalated the variation of categorizations instead of lessening the issues or complications. Consequently, many important nodes remain unresolved. Thus, the fact that some samples were misidentified in this study may not be due to a misidentification or phylogenetic approaches used, but probably due to unidentified node.

The disparity index (DI) analysis of both markers shows in Table 4 was further corroborated the above analysis in which it indicates that the rbcL gene has lower range of disparity index than rpoC1. The lowest DI value was found between the PB (Perlis) and AS1 (Kedah) samples with the value of 0.0082 for rbcL gene. However, the highest value was between PB (Perlis) and BP2 (Penang) for rpoC1 gene, demonstrating the greatest variations in the composition of bases in each sequence. This might be the consequences of a high heterogeneous substitution structures constructed within the samples as amplified by rpoC1 gene.

A perfect barcode of DNA should possess sufficient conserved segment to design a universal primer and should have high variations to be utilized for discrimination of a species and must able to distinguish between complex species. This can be attained if a species has remarkably high genetic distance compared to the other intraspecific individual within the group (Hebert et al., 2004; Mankga et al., 2013). The phylogenetic tree-based approaches have been substantially used in DNA barcode studies in order to assign a species to its respective taxa and the widely utilized phylogenetic tree is NJ, where the essence of evaluation was relied on morphological distance and the species’ documentation of evolution (Liu et al., 2014).

Traditionally, approach on molecular systematic have depended on contrasting a restricted number of the orthologous sequence in order to acquire evaluation of species relationships across the tree of life (Edger et al., 2014). Additionally, a supposition of tree-based analysis of nucleic acid sequences is that each position is unconnected with the other positions was always made in analysing phylogenetic relationship of a species (Alvarez and Wendel, 2003). The outcomes of this research demonstrated that the majority of sequence sites for rpoC1 are not evolving independently, but more precisely are surviving in another position/site in order to conserve the secondary composition of a molecule. Thus, in future studies, effort should be made to determine all coevolving positions and relevant adaptations are needed before utilizing rpoC1 region as a barcode marker.

Nowadays, there are various DNA marker accessible for researchers in the plant sciences field. However, many researchers have disputed that it is very strenuous to find a universal barcode for the identification of all plant species due to morphological and geographical variations as well as reticulate evolution ( Mosa et al., 2019). A review on the role of DNA barcoding as a powerful tool for plant biodiversity analysis by Mosa et al. (2019) revealed that the ITS and rbcl gene have been recognized as core barcode markers. For example, Asterella could be distinguished by the rbcL gene alone with 90% successful rate (Table 5), while approximately 70% of angiosperm and only 32% of gymnosperm species (especially Araucaria) were successfully identified by the association of several chloroplast DNA regions (Hollingsworth et al., 2009a,b; Ran et al., 2010). This further support our conclusion on the reliability of rbcL as a DNA barcode marker. However, there is no universal barcode method for plants and researchers are utilizing several different gene markers for improving success in species identification.

Another important point has raised in this study is that some species (for e.g. P. umbellatum and P. auritum) were region-specific, in which they can be found specifically in Penang and absent in other sampling locations (Fig. 3). This is probably due to the habitat conditions that are suitable for both of the species. They usually occurs in the undergrowth of evergreen rainforest, but also in clearings and on river banks in which they always occurs in damp localities. That some of the species were not region-specific (for e.g. P. sarmentosum) probably due to the fact that they are easy to grow and therefore are easily available. Additionally, P. sarmentosum is an invasive species and easily colonized areas where they were found.