1. Introduction

Accurate plant species identification is fundamental for biological studies. DNA barcoding using standard molecular markers has been developed for species identification (Hebert et al., 2003; Smith et al., 2005; Hollingsworth et al., 2011). Researchers have proposed the plastid regions (rbcL and matK) as well as the nuclear ribosomal internal transcribed spacer (ITS) as standard DNA barcode sequences for land plants (CBOL Plant Working Group, 2009; China Plant BOL Group, 2011). Although using these standard DNA barcodes could help identify species in many land plants, their discriminatory power is low for recently originated groups via rapid speciation (Hollingsworth et al., 2011, 2016; Li and Zeng, 2015; Coissac et al., 2016). Next-generation sequencing has expanded the plant barcodes using substantial genomic sequences to enhance the discrimination efficiency (Yang et al., 2013; Li and Zeng, 2015; Hollingsworth et al., 2016). Complete plastomes and nuclear ribosomal DNA (nrDNA) sequences, which are comparable to the standard barcodes and have more information, can improve the resolution (Kane et al., 2012; Li and Zeng, 2015; Coissac et al., 2016; Hollingsworth et al., 2016). However, only a few empirical studies have assessed their discrimination efficiency (e.g., Ji et al., 2019; Fu et al., 2022).

Genome skimming is a cost-effective and user-friendly approach that uses next-generation sequencing technology to produce substantial genomic sequences at low coverage (Kane et al., 2012; Straub et al., 2012; Coissac et al., 2016). It can recover highly repetitive genomic regions, including plastid genome and nrDNA, from plant cells and extend the plant barcodes (Kane et al., 2012; Hollingsworth et al., 2016). Meanwhile, RNA-sequencing is a high-efficiency approach that can produce substantial nuclear sequences and organelle data (Ekblom and Galindo, 2011; Smith, 2013).

Bamboos are a challenging group of plants where species identification is difficult. Moreover, only a few reproductive organ samples can be collected for bamboo identification due to the long vegetative period and flowering cycle (up to 120 years) (Janzen, 1976). Therefore, vegetative characteristics are a common choice for describing and identifying bamboo species (Soderstrom, 1981; Keng and Wen, 1989; Stapleton et al., 2009). However, these features are variable and prone to environmental factors (Guo et al., 2001), complicating bamboo classification and identification.

The temperate woody bamboo tribe Arundinarieae (Poaceae: Bambusoideae) is notorious for its intractable intergeneric phylogenetic relationships (Zeng et al., 2010; Cai et al., 2012; Yang et al., 2013; Ma et al., 2014; Guo et al., 2019, 2021). All the species of the tribe originated via an allopolyploid event through hybridization, and thus, intense conflicts exist between the nuclear (Guo and Li, 2004; Wang et al., 2017) and plastid phylogenies (Yang et al., 2013; Ma et al., 2014). Recently, Zhang et al. (2020b) proposed a subtribal classification of Arundinarieae based on analyses of double digest restriction-site associated DNA-sequencing (ddRAD-seq) data. The divergence of the four out of five major lineages (the ADH, Gaoligongshania, leptomorph, and pachymorph lineages) occurred in a short period during the Miocene (13.65–11.05 million years ago, Ma) (Guo et al., 2021). Furthermore, the relationships between the closely related genera and species within the major lineages of Arundinarieae have not been fully resolved.

The utility of standard DNA barcodes (rbcL, matK, and ITS) in identifying bamboos has been tested in a previous study (Cai et al., 2012) with 27 species of the temperate woody bamboos sampled. However, the discriminatory power was low. Three-marker combinations of the plastid barcodes (rbcL, matK and trnH-psbA) provided a low level of discrimination (< 14.8%), while ITS performed much better but the universality of PCR amplification and sequencing was relatively low. Obviously, this result is far from adequate for achieving DNA identification in temperate bamboos. Therefore, more powerful markers need to be screened to improve the discriminatory power and provide practical benefits for species identification in bamboos.

Fargesia is the largest genus within the Arundinarieae (Li et al., 2006; Bamboo Phylogeny Group, 2012) and belongs to the pachymorph lineage (Zhang et al., 2020b), with more than 90 species. At least 78 species of Fargesia, of which 77 are endemic, are distributed in China (Yi, 1996; Li et al., 2006). The Fargesia species are shrubby or sub-arborescent and mainly found in the subalpine regions (1400–3800 m altitude) in Southwest China, and therefore adapted to the high mountain environment (Yi, 1996; Li et al., 2006). As the dominant species in the shrub layer of the forests, these bamboos play an indispensable role in ecosystem services, providing food and shelter for the giant panda and many other endangered animals (Yi, 1985a, 1985b; McNeely, 1996). The generic boundary and infrageneric classification in Fargesia had faced an extended period of controversy (Soderstrom, 1979; Chao et al., 1980; Hsueh and Li, 1987; Soderstrom and Ellis, 1987; Yi, 1996; Chao and Renvoize, 1989; Li et al., 2006). Previous molecular studies based on plastid loci and few nuclear genes and with relatively limited taxon sampling did not resolve the infrageneric phylogeny of Fargesia (Guo and Li, 2004; Triplett and Clark, 2010; Yang et al., 2013; Zhang et al., 2019a, 2019b; Zhou et al., 2020). In a recent study, Ye et al. (2019) used ddRAD-seq data to divide Fargesia into four well-supported clades and two isolated species (V-Fargesia1, V-Fargesia2, V-Fargesia3, and V-Fargesia4; Fargesia angustissima and Fargesia lincangensis, respectively). The four main clades of Fargesia originated between 2.94 Ma and 3.88 Ma with rapid radiation (Ye et al., 2019). A prior phylogenetic framework has been established for Fargesia; however, it is still a big challenge for species identification in such a complex genus.

In this study, we would like to take the genus Fargesia as a pilot test for the species discriminatory ability of the complete plastomes and nrDNA sequences in woody bamboos. Specifically, we aimed to address the following questions: (1) How about the discriminatory power of the complete plastomes and nrDNA sequences in Fargesia? (2) What are the differences among the standard DNA barcodes, the complete plastomes and nrDNA sequences in species identification capability for the genus? (3) To what extent the complete plastomes and nrDNA sequences can improve the phylogenetic resolution of Fargesia with more taxon sampling? The findings of this study will provide an empirical case for the use of complete plastomes and nrDNA sequences as organelle barcodes in bamboos and other taxonomically difficult plant groups.

2. Materials and methods 2.1. Taxon sampling

We sampled 167 individuals from 62 species of Fargesia, which covered approximately 70% of the total described species in the genus. Species were mainly sampled from the mountains of Southwest China and adjacent areas. Two to nine individuals per species were sampled for 32 species, and a single individual was sampled for the remaining 30 species. The 167 samples included 129 silica-dried leaf materials used for total genomic DNA extraction and 38 fresh leaf materials, immediately frozen and stored in liquid nitrogen, used for RNA extraction. All vouchers were deposited in the Herbarium of Kunming Institute of Botany, Chinese Academy of Sciences (KUN). Another seven individuals with whole plastid genome data and 22 individuals with ITS data were downloaded from the GenBank. In addition, two individuals of Phyllostachys edulis were sampled as outgroup taxa. Detailed information on sampling, species, vouchers, and data sources is provided in Table S1.

2.2. Sequencing, assembly, and annotation

Total genomic DNA was extracted from silica-dried leaf materials using the modified CTAB method following Yang et al. (2014). Purified DNA samples were sheared into 500 bp fragments to construct short-insert libraries, and 2 × 150 bp paired-ends sequencing was conducted on an Illumina HiSeq 2000 platform, generating at least 2 GB of data for each sample. High-quality reads were assembled into whole plastid genomes and nrDNA sequences using the GetOrganelle pipeline (Jin et al., 2020). Plastomes were annotated with the Plastid Genome Annotator (PGA; Qu et al., 2019) and manually corrected using Geneious v.9.0.2 (Kearse et al., 2012). We also reannotated plastomes downloaded from the GenBank database using PGA. The ribosomal DNA sequences of Oryza sativa (KM036285) were used as the reference to assemble and annotate the nrDNA sequences.

RNA extracted from the fresh leaf materials was used for library preparation and Illumina sequencing. High-quality Illumina reads were de novo assembled into contigs, transcripts, and unigenes with Trinity v.2.10.0 (Grabherr et al., 2011). All transcriptomes were mapped to the reference ribosomal DNA sequences, and the nrDNA sequences were extracted using Geneious v.9.0.2. For incomplete nrDNA sequences, only ITS regions were extracted.

2.3. Data analysis

For each annotated plastid genome, we extracted the large single copy (LSC), small single copy (SSC), inverted repeat (IR) regions, and the plastid DNA barcode regions rbcL, matK, and trnH-psbA using PhyloSuite v.1.2.2 (Zhang et al., 2020a); the ITS regions were subsampled from the nrDNA sequences. We compared species discriminatory power of the following seven data sets: (1) the complete plastome; (2) the large single copy (LSC) region; (3) the small single copy (SSC) region; (4) one inverted repeat (IR) region; (5) a combination of plastid barcode regions (matK, rbcL, and trnH-psbA); (6) the ITS region consisting of ITS1-5.8S-ITS2; and (7) the nrDNA sequences consisting of 18S-ITS-26S. All sequences from these seven data sets were aligned using MAFFT v.7.22 (Katoh and Standley, 2013) with default parameters and manually adjusted in Geneious.

Further, we evaluated the species discrimination success by using a tree-based method and distance-based analysis. Phylogenetic analysis with the seven data sets was performed using the maximum likelihood (ML) method in RAxML v.8.2.12 (Stamatakis, 2006) with the GTR+Γ model and 1000 rapid bootstrap (BS) replicates. To compare "barcoding gaps" among data sets (for species with at least two individuals per species in each data set), we calculated pairwise distance using Kimura's two-parameter (K2P) substitution model with a pairwise deletion in MEGA X (Kumar et al., 2018).

3. Results 3.1. Characteristics of data sets

We obtained 136 plastomes for Fargesia, including 78 novel and 58 newly published ones (Ye et al., 2021). These plastomes are circular and highly conserved in size and gene content (Fig. S1). There was little variation in the length of the whole plastomes (139, 404 bp to 139, 718 bp), or the LSC (82, 999 bp to 83, 302 bp), SSC (12, 786 bp to 12, 861 bp), and IR (21, 793 bp to 21, 804 bp) regions (Table S2). Each plastid genome encoded 131 identically arranged genes, including 84 protein-coding genes, 39 transfer RNAs, and 8 ribosomal RNAs (Fig. S1).

Data sets produced from the plastid sequences included 138 samples representing 63 species. The matrix of plastomes yielded 141, 988 positions, with only 1426 (1.00%) variable sites and 577 (0.41%) parsimony-informative (PI) sites. The SSC region had the highest variation rate (1.61%) with 212 variable sites and 87 (0.66%) PI sites, followed by the LSC region with 1134 (1.33%) variable sites and 468 (0.55%) PI sites. The IR region had the lowest variation rate, with 39 (0.18%) variable sites and 12 (0.05%) PI sites. The data set based on the combined plastid barcode regions matK+rbcL+trnH-psbA was 3514 bp long, with 39 (1.11%) variable sites and 16 (0.46%) PI sites (Table 1).

We newly obtained 110 nrDNA sequences and generated two nuclear data sets, ITS and nrDNA. The nrDNA data set recovered from 110 samples included 90 assembled from genome skimming and 20 extracted from transcriptomes; the length of the sequences ranged from 5790 bp to 5803 bp (Table S2). A total of 93 (1.60%) variable sites and 69 (1.19%) PI sites were identified in the 5813 bp alignment matrix (Table 1). The ITS data set consisted of 150 taxa, including 110 subsampled from nrDNA, 18 recovered from transcriptomes, and 22 downloaded from the GenBank database (Table S2). Alignment of ITS showed the highest percentage of variation, with 67 (10.60%) variable sites and 50 (7.91%) PI sites (Table 1).

3.2. Phylogenetic reconstruction

Phylogenetic analysis using the plastome and LSC sequences produced nearly identical topologies, with the tree based on plastomes recovering higher branch support values (Figs. S2 and S3). However, the phylogenetic resolution was low in both data sets. Here, V-Fargesia1, V-Fargesia3, V-Fargesia4 plus F. angustissima, and F. lincangensis formed a large group, with species from these clades nested within each other. V-Fargesia2 was paraphyletic with Fargesia incrassata, forming a monotypic clade, and the remaining V-Fargesia2 species plus Fargesia funiushanensis formed another clade (Figs. S2 and S3). Ye et al. (2019) did not include F. funiushanensis, which clustered into the V-Fargesia2 clade in the present study. Due to insufficient phylogenetically informative sites, the SSC region, IR region, and the combined plastid barcode regions (matK+rbcL+trnH-psbA) were not used in phylogenetic analysis.

Phylogenies based on nuclear ribosomal data (ITS and nrDNA) resolved Fargesia angustissima as sister to the remaining Fargesia species. V-Fargesia1 was resolved as a highly-supported monophyletic clade. V-Fargesia2 plus F. funiushanensis formed a monophyletic group. Within V-Fargesia3, one monophyletic clade was formed by Fargesia erecta, F. mairei, F. muliensis and F. yajiangensis, whereas another monophyletic clade was formed by F. yulongshanensis and F. xianggelilaensis. The remaining V-Fargesia3 species were nested within the polyphyletic V-Fargesia4 clade (Figs. S4 and S5).

3.3. Species discrimination 3.3.1. Species discrimination based on ML tree analysis

Species identification by the tree-based method was considered successful when a species with at least two individuals was resolved as monophyletic with a branch support value of over 50%. The plastid genome sequences that showed the highest discriminatory power were the complete plastome and the LSC region sequences (28.6% of 35 species discriminated; Fig. 1; Table 2). Sequence data from the SSC (20.0%) and IR regions (8.6%) showed a lower discriminatory power. The sequences with the lowest discriminatory power were matK+rbcL+trnH-psbA (5.7%) (Table 2). Out of the 35 species with two or more than two individuals sampled, only two species, Fargesia glabrifolia (BS = 85%) and F. yajiangensis (BS = 93%) were resolved as monophyletic with > 50% bootstrap support; a third species, F. hygrophila, had 30% bootstrap support (Table 3). By contrast, in addition to F. glabrifolia (BS = 69%), F. hygrophila (BS = 96%), and F. yajiangensis (BS = 98%), another seven species were resolved as monophyletic (100% bootstrap support for F. angustissima, F. hsuehiana, F. incrassata, F. pubifolia, F. qinlingensis, and F. spathacea, and 83% for F. canaliculata) in the complete plastid genome (Fig. 1; Table 3). The LSC region showed a discrimination rate similar to that of complete plastomes with only slight differences in support values for identical monophyletic groups. The SSC region sequences resolved F. yajiangensis, F. hygrophila, and five additional species (F. angustissima, F. hsuehiana, F. incrassata, F. pubifolia, F. qinlingensis) as monophyletic, but the three F. glabrifolia individuals failed to cluster together despite the combined plastid barcode regions supporting their monophyly. The IR sequences showed a low discrimination rate, with only F. angustissima, F. hsuehiana, and F. pubifolia resolved as monophyletic groups with 100% bootstrap support; furthermore, these species were different from those resolved by the combined plastid barcodes (Fig. 1; Table 3).

The ITS data set resolved 17 out of 36 species as monophyletic, but only 13 (36.1%) species had a branch support value higher than 50% (Table 2 and S3). The nrDNA data set showed the highest discriminatory power, with 57.7% of species distinguished (15/26; Fig. 2; Table 2).

3.3.2. Species discrimination based on genetic distance

In the distance-based analysis, a species was considered successfully identified when the minimum interspecific distance was greater than its maximum intraspecific distance. The discriminatory power in the distance-based analysis was lower than the tree-based analysis using the complete plastome (22.9% of 35 species discriminated vs. 28.6% in tree-based method), the LSC region (22.9% vs. 28.6%), the SSC region (17.1% vs. 20.0%), and matK+rbcL+trnH-psbA (2.9% vs. 5.7%). In contrast, the discriminatory power of the distance-based analysis was higher than that of the tree-based analysis for ITS (47.2% of 36 species discriminated) and nrDNA (65.4% of 26 species discriminated) data. A similar discrimination rate was obtained for the IR (8.6%) regions between the distance-based and tree-based analyses (Tables 2 and 3).

The minimum interspecific genetic distance for the plastome data set ranged from 0.000007 to 0.000366. For the LSC data set, a single species pair Fargesia nitida versus F. robusta showed the minimum interspecific genetic distance of zero. In contrast, the minimum interspecific genetic distance of zero was observed more frequently in other data sets (the SSC region, 14 out of 35 species; the IR region, 30 out of 35 species; matK+rbcL+trnH-psbA, 32 out of 35 species; ITS, 13 out of 36 species; nrDNA, 8 out of 26 species) (Table S4).

4. Discussion 4.1. DNA barcoding in Fargesia

Our results show that complete plastomes and nrDNA sequences discriminate species in Fargesia at higher rates than do standard barcodes; however, both complete plastome and nrDNA sequences were unable to discriminate all Fargesia species.

The organelle markers from the plastid genome had low power for species discrimination, ranging from 5.7% (2/35, combined plastid barcodes) to 28.7% (10/35, complete plastome) for Fargesia. The complete plastome provided substantial information and notably higher species discriminatory power than the combined plastid DNA barcodes. This increase is mainly due to the additional sequence characters. Sequences from three regions of the plastome (LSC, SSC, and IR), which have different levels of sequence variation, demonstrated lower discriminatory power that differed among them. The discriminatory power does not appear to be always consistent with substitution rates of the three plastid regions. Our results showed that the IR region had the least sequence variation (Table 1). Zhu et al. (2016) also reported the lowest molecular evolutionary rate in the IR region compared to the LSC and, SSC, due to the duplication nature of this region. And the discrimination power of data from the IR region was indeed lowest among the three main regions of the plastome (8.6% for IR, 20.0% for SSC, and 22.9% or 28.6% for LSC; Table 2). However, the SSC region had the highest evolutionary rate and the shortest sequence length among the three regions (Zhu et al., 2016; Table 1), but the discrimination rate of data from SSC was lower than that of the LSC region (17.1% vs. 22.9% and 28.6%; Table 2). The LSC region had the longest sequence length and the highest number of informative characters, and had nearly the same discrimination rate as the whole plastome. This indicates that the total number of informative sites, rather than the substitutional rate, determines the level of discrimination rate. However, the entire plastome sequence was not able to distinguish most of the species (25 out of 35) in the genus. These observations indicate that plastid markers as barcodes have very limited power in species identification in Fargesia.

In contrast, our data showed that the nrDNA markers were twice or three times as efficient as the plastid genome sequences in species discrimination in Fargesia (Table 2). Compared to ITS, the nrDNA sequences had a higher discrimination rate and significantly increased support values for the monophyletic species (Table 3), likely due to a greater number of informative sites (Table 1). Therefore, these nuclear markers are more promising than plastid markers for species identification in Fargesia. Even the ITS data alone had higher discriminatory power than complete plastome sequences. Neither ITS nor nrDNA data was able to discriminate all species, although nrDNA sequences had the highest discrimination rate among the data sets compared (up to 65.4%; Table 2).

Several factors potentially led to the low resolution of species discrimination in Fargesia. One important factor may be the long generation time associated with the low rate of DNA sequence variation in the Arundinarieae (Ma et al., 2017). Many of the species may have recently diverged, and thus lack sufficient sequence differences for discrimination. Furthermore, hybridization is prevalent among bamboo species (Triplett and Clark, 2010; Yang et al., 2013; Triplett et al., 2014), which presents another challenge for discriminating species. Hybridization can lead to sharing plastid genomes and nuclear genes between different species (Du et al., 2009), blurring the species boundary (Petit and Excoffier, 2009; Hollingsworth et al., 2016). Recent studies have suggested that Fargesia recently experienced rapid radiation and extensive hybridization, resulting in the divergence of the four main clades between 2.94 Ma and 3.88 Ma (Ye et al., 2019, 2021). In addition, Fargesia is a taxonomically complex group with narrowly defined species, and these species are mainly distributed in remote montane areas, making sampling difficult. Moreover, the delimitation of some species in Fargesia may be unreasonable (Zhang et al., 2014).

Complete plastome and nrDNA sequences have been proposed as tools to address the resolution limitations of standard barcodes in plants (Kane et al., 2012; Li and Zeng, 2015; Coissac et al., 2016; Hollingsworth et al., 2016); however, few studies have assessed their actual discrimination rate. Moreover, the discriminatory power revealed in studies using multiple congeneric species varies. For example, either complete plastome or nrDNA is used to distinguish species in Taxus (Fu et al., 2019) and Theobroma (Kane et al., 2012) successfully. However, in recently diverged and more complex plant groups, such as Panax (Ji et al., 2019) and New Caledonian Araucaria (Ruhsam et al., 2015), both plastome and nrDNA failed to improve the discriminatory power significantly. Furthermore, in Rhododendron (Fu et al., 2022) and Fargesia, the discrimination rate was substantially higher with plastome and nrDNA sequence data than with standard DNA barcodes; however, many species remain unresolved. These case studies suggest that complete plastomes and nrDNA sequences can increase the discriminatory power and practical benefits for species identification beyond the standard barcodes. However, factors such as slow sequence evolutionary rate, interspecific hybridization, recent origins, and rapid radiation can significantly hamper the power of DNA barcodes for identification of closely related species (Ruhsam et al., 2015; Hollingsworth et al., 2016). In such groups, a substantial number of nuclear markers are needed to improve the resolving power.

4.2. Phylogenetic implications

In our study, plastome and nrDNA sequences provided poor phylogenetic resolution of Fargesia, indicating that plastome and nrDNA sequences both lack sufficient genetic differentiation. This is likely due to the low rate of molecular evolution and recent divergence of many species in the genus (Ye et al., 2019; Zhou et al., 2019). Additional sampling alone may not help increase the resolution of Fargesia phylogeny. Nonetheless, the topology of the phylogenetic tree based on nrDNA sequences was more similar to that of the well-resolved ddRAD-seq phylogeny (Ye et al., 2019) than that based on plastomes. Moreover, the nrDNA data also identified two clades, i.e., V-Fargesia1 and V-Fargesia2 as in the ddRAD-seq phylogeny. For the plastome, its topology differed sharply from the ddRAD-seq phylogeny, with only a small number of species in consistent phylogenetic position. This phenomenon of cytonuclear discordance was also observed by Zhou et al. (2020) and Ye et al. (2021), which may reflect the hybridization and introgression events and the evolutionary history of rapid radiation in the genus. Notably, in all data sets analyzed F. funiushanensis was clustered in the V-Fargesia2 clade. These findings indicate that F. funiushanensis, which was not sampled by Ye et al. (2019), should be a member of the V-Fargesia2 clade (Zhou et al., 2020).

5. Conclusions

Fargesia represents a taxonomically difficult group with DNA barcoding difficulties due to long generation time, interspecific hybridization, recent origin, and rapid radiation. Screening of more markers is needed to increase species discriminatory power and practical benefits of specimen identification. Here, we comprehensively evaluated the discriminative ability of complete plastome and nrDNA sequences compared with the standard DNA barcodes. They provided more informative characters, substantially increasing the discriminatory power in Fargesia. However, they were not powerful enough to distinguish all species with the highest discrimination rate (65.4%), as expected due to the complex evolutionary history of Fargesia. More nuclear markers are desirable to gain further improvement for species identification in Fargesia and other similar plant groups.

Acknowledgements

We thank Ms. Li-Ying Luo, Mr. Zu-Chang Xu, Ms. Hong Wu, and Ms. Shuang-Xiu Xu for assistance with fieldwork, and Prof. Jenny Xiang and three anonymous reviewers for critical suggestions to improve the presentation of the manuscript. We also thank the Laboratory of Molecular Biology and iFlora High Performance Computing Centre of the Germplasm Bank of Wild Species (GBOWS) for molecular experiments and data analysis support, respectively. This work was supported by grants from CAS' Large-scale Scientific Facilities (Grant No. 2017-LSF-GBOWS-02), the Key R & D Program of Yunnan Province, China (Grant No. 202103AC100003), and Ten Thousand Talent Program of Yunnan Province (Grant No. YNWR-QNBJ-2020-297).

Author contributions

D.-Z. L. and P.-F. M. designed and conceived this study; S.-Y. L. and X.-Y. Y. collected samples; S.-Y. L. and P.-F. M. analyzed data; S.-Y. L. drafted the manuscript; P.-F. M., Z.-H. L. and D.-Z. L. modified and improved the manuscript.

Data availability

The data presented in this study are available in the supplementary material.

Declaration of competing interest

The authors declare that they have no conflicts of interest.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2022.04.002.