Identification of species and materia medica within Saussurea subg. Amphilaena based on DNA barcodes

Taxon sampling

In total, 20 species were sampled in the present study, including 18 from the 38 species recognized in the latest monograph on S. subg. Amphilaena (Raab-Straube, 2017), one recently published species, S. bogedaensis (Chen & Wang, 2018), and a Jurinea species, which was selected as an outgroup according to a previous study (Wang et al., 2009). Photos of some species are presented in Fig. 1. Our sample focus on medical resources and 15 species formally recorded in the medical literature were included in the analyses (Table 1). For most of the species in the ingroup, we collected from two or more populations, with more than three individuals from each population. In total, we collected 132 individuals and their details are listed in Table 2.

DNA extraction, PCR amplification, and sequencing

Genomic DNA was extracted from dried leaves in silica gel using the CTAB method (Doyle, 1987). Five regions (rbcL, matK, trnH-psbA, trnK, and ITS) (Berends, Jones & Mullet, 1990; Ford et al., 2009; Olmstead et al., 1992; Sang, Crawford & Stuessy, 1997; White et al., 1990), were amplified and sequenced using the primers listed in Table 3. A PCR reaction mixture comprising 25 µL was prepared and amplified according to the procedure described by Wang et al. (2009). The PCR products were sent to the Beijing Genomics Institute for commercial sequencing. Sequences were aligned using CLUSTALX v.2.1 (Thompson et al., 1997) with the default settings and adjusted manually with Bioedit v.7.0.5 (Hall, 1999). All of the sequences were registered in GenBank (Table 2).

Data analysis

We constructed 31 datasets for ITS, psbA-trn H, matK, and trnK, either individually or in different combinations. For the combination of ITS and each chloroplast loci, incongruence length difference (ILD) was preferred to test the incongruence (Farris et al., 1995) using PAUP version 4b10 (Swofford, 2003). For each dataset, the inter- and intraspecific genetic divergences were calculated as described by Meyer & Paulay (2005) and used to determine whether a barcoding gap was present. For each dataset, best close match (BCM) and two tree-based methods comprising neighbor-joining (NJ) and Bayesian inference (BI) were employed to analyze the five single markers and their different combinations. BCM analysis was conducted using the SPIDER package in R (Brown et al., 2012). NJ trees were constructed using PAUP with the Kimura two-parameter model (Swofford, 2003). Support for nodes was assessed based on 100,000 bootstrap replicates. BI analysis was implemented using MrBayes on XSEDE (v3.2.6) (Ronquist et al., 2012) and the optimal models for each marker were determined according to Akaike’s information criterion with jModelTest2 in XSEDE (v2.1.6) (Darriba et al., 2012). Species were considered to be identified successfully if individual samples of a species clustered in species-specific monophyletic clades.

Proposed DNA barcodes for S. subg. Amphilaena

Among the fragments tested in the present study, ITS obtained a much higher success rate compared with the other loci. In addition, all of the combinations without ITS yielded much lower success rates, regardless of the method used (Table 7). Moreover, the rate of successful PCR (92.7%) was more or less higher for ITS than the other fragments (72.9–91.6%). It has also been reported that this fragment is highly efficient in other Asteraceae genera (Gao et al., 2010; Gong et al., 2016). However, an intrinsic problem with this fragment is that an individual may have undergone recent hybridization, thereby resulting in multiple mosaic sites (Li et al., 2011). In S. subg. Amphilaena, two species failed to form monophyletic clades in the BI and NJ trees, which could be attributed to the presence of multiple mosaic sites (Fig. 3). However, ITS performed better than the other fragments in S. subg. Amphilaena, and thus we propose that this fragment should be the first or best choice when selecting only one of the current candidates.

We found that it was difficult to identify the best second choice after ITS. TrnK performed much better than rbcL in terms of its efficiency when used individually, but its combination with ITS obtained contradictory results, i.e., ITS + trnK was inferior to ITS + rbcL in terms of efficiency. This contradictory result was unexpected and it is not common in other taxa (Cao et al., 2010; Müller & Borsch, 2005). We attributed this result to higher degree of congruence of the concatenated sequences of rbcL and ITS (P = 0.12 for ILD test), in compare to trnK and ITS (P = 0.001). But it might derive from some other mechanisms, such as the higher rate of mutation for trnK that could have caused differentiation within species, but not high enough to form distinct genetic differentiation among species, and thus a failure to cluster as a monophyletic group in line with species (Naciri, Caetano & Salamin, 2012; Petit & Excoffier, 2009). Therefore, we suggest that using trnK alone is problematic and instead we propose to use rbcL as complementary to ITS because this combination could identify all 19 of the sampled species based BCM, and 17 by NJ or BI (89%) (Table 7) (Fig. 4).

The two loci comprising trnH -psbA and matK were affected by the same problem as trnK, with higher mutation rates and barcode efficiencies compared with rbcL when used individually, but lower efficiency when combined with ITS. Thus, their combination with ITS + rbcL failed to significantly increase the success rate and lower results were even obtained in some cases (Table 7). However, among the combinations without ITS, the combination with higher mutation rates was more efficient than those with lower mutation rates, e.g., trnK + trnH-psbA was better than matK + rbcL, which was proposed previously as the core DNA barcode for plants (Hollingsworth et al., 2009). Therefore, if ITS is subjected to hybridization, we propose that the priority order should be the following: trnK > trnH-psbA >  matK > rbcL. Moreover, the combination with more loci performed better than that with less loci. However, even the combination of all four loci was not sufficient to discriminate each species and new fragments should be considered.

Insights into taxonomic problems based on DNA barcodes

Most of the analyses failed to identify the species within two groups, i.e., S. luae vs. S. publifolia and S. globosa vs. S. erubescens (Figs. 3–5; Table 7). We found that these failures might have been attributable to taxonomic problems. For the first group, we found that S. luae was rather heterogeneous in terms of the ITS sequences. Some cp sequences were slightly differentiated compared with S. velutina, but the others were closer to those in S. glandulosissima or S. uniflora (Fig. 5). By contrast, the ITS sequences lacked variance and after excluding the mosaic sites, they were closely related in S. pubifolia or S. bracteata (Fig. 3). These nuclear-cytoplasmic inconsistencies suggest that hybridization may have occurred among these species.

The second group comprising S. globosa and S. erubescens was often confused in previous studies because the latter resembles a smaller form of S. globosa, which has various forms across its distribution (Raab-Straube, 2017). In agreement with the morphology, the genetic distance between the cp sequences within S. erubescens was zero whereas that within S. globosa was 0.04% (Table 5), which is even larger than that between S. erubescens and S. globosa (Table 6). The ITS sequences had a very similar pattern and the rich mosaic sites in both species also indicated differentiation accompanying substantial gene flow (Naciri, Caetano & Salamin, 2012). Both the BI and NJ methods found that S. globosa formed a clade within which S. erubescens nested as a monophyletic clade (Fig. 3). Based on these results, we propose that S. globosa might be a species with a series of differentiated populations where S. erubescens represents one of the most obvious. The current delimitation might need revision on the basis of extensive morphological as well as genetic diversity across the distribution range of both species.

Identification of the medicinal species and the potential substitutes

All of the known medically important species could be identified using our proposed DNA barcodes, i.e., ITS + rbcL or ITS alone (Table 7; Figs. 3–4). Moreover, some species such as S. bogedaensis, S. glandulosissima, S. polycolea, S. wettsteiniana, and S. orgaadayi could be identified with the cp DNA barcodes (Fig. 5). This high rate of success was unexpected because some species such as the two species in the S. obvallata complex (S. glandulosissima and S. sikkimensis) have been morphologically confused for many years and they were only separated very recently (Raab-Straube, 2017). Their distinction is indicative of difference in bioactive components. Therefore, our results caution against their indiscriminating usage in medicine.

Barcode sequences can also help to identify substitutes for medically useful species because closely related species might possibly share the same or similar secondary metabolites and bioactivities (Zhou et al., 2014). Thus, we propose that nine of the 15 medically useful species might be substituted by their close relatives according to the molecular phylogenetic context. Six of these species, which formed three groups, are also morphologically similar, i.e., S. involucrata and S. orgaadayi or S. bogedaensis, S. globosa and S. erubescens, and S. wettsteiniana and S. glandulosissima (Fig. 3) (Raab-Straube, 2017). Among the remaining three species, S. bracteata appears to be closely related to S. pubifolia whereas S. iodostegia and S. nigrescens are closely related to each other according to phylogenetic tree (Fig. 3). These affinities were not expected according to their morphology, but they are possibly due to convergent evolution or radiation in Saussurea (Wang et al., 2009). Secondary metabolomes or bioactivities are wanted to confirm their similarity.