1. Introduction

Species identification in some angiosperm taxa is challenging. Several approaches have been developed to identify species in these groups, including the standard DNA barcodes (Li et al., 2011; Bieniek et al., 2014; Jiang et al., 2020), the next generation barcodes (Kane et al., 2012; Li et al., 2015; Hollingsworth et al., 2016), capturing single-copy nuclear genes (Alvarez et al., 2008; Straub et al., 2012; Liu et al., 2021), and the recently developed methods based on all reads generated from genome skims (Yi and Jin, 2013; Ondov et al., 2016; Balaban et al., 2020). Genome skimming, defined as low-depth genome sequencing, is typically applied in organelle genome recovery and species identification (Straub et al., 2012). It is a popular approach in phylogenomics, as it can extract total DNA from a small amount of fresh or silica-dried samples, or even from long-preserved herbarium specimens for low-depth genome sequencing (Zeng et al., 2018; Nevill et al., 2020). However, most nuclear DNA data produced by genome skimming are discarded in current species identification practices, reducing its discriminatory power (Rubinoff and Holland, 2005; Straub et al., 2012; Dodsworth, 2015; Yu et al., 2018; Bohmann et al., 2020; Antil et al., 2023). Genome skimming could be optimized by using assembly-free and alignment-free methods. Although methods for assembly-free genetic distance calculation exist, such as Co-phylog (Yi and Jin, 2013), Mash (Ondov et al., 2016) and APPLES (Balaban et al., 2020), these approaches require deep sequencing depth, at least 5× (Ondov et al., 2016), which increases costs and the need for computing power.

An assembly-free and alignment-free approach called 'Skmer' is proposed to calculate genomic distances using genome skims. Data with 1× sequencing depth might be sufficient for accurate species identification in most cases when using Skmer (Sarmashghi et al., 2019). It was confirmed that using Skmer for species identification with marine mollusks was feasible and efficient (Xu et al., 2022). Furthermore, Skmer has been shown to increase species identification rate in Cymbidium Sw. (Orchidaceae) (Zhang et al., 2023) and Torreya Arn. (Taxaceae) (Mo et al., 2023). It seemed that Skmer could be proposed as a potentially credible approach for species discrimination.

Species delimitation has long been controversial in Schima Reinw. ex Blume (Gordonieae, Theaceae) (Melchior, 1925; Bloembergen, 1952; Willis and Airy Shaw, 1985; Keng, 1994; Chang and Ren, 1998). Schima is distributed in subtropical and tropical areas of Asia. It is characterized by white flowers, globose to oblate fruits containing small reniform seeds with marginal membranous wings (Fig. 1; Ming and Bartholomew, 2007). It is ecologically important as one of the dominant elements of the subtropical evergreen broadleaved forests in East Asia (Yu et al., 2017a, 2017b; Li et al., 2020; Tang et al., 2020). Schima species also used in traditional Chinese medicine and have recently been shown to have anti-cancer properties (College 1985; Liang et al., 2019; Wu et al., 2019). Furthermore, these plants are used as horticultural resources, fire-resistant trees, and furniture wood (Yang et al., 2022; Zheng et al., 2022). The current classification system of Schima relies on morphological characters, which overlap in some species, making species identification difficult. Recent research has used DNA barcodes (matK, rbcL, trnH-psbA and nrITS) and next generation barcodes (plastome, nrDNA arrays) to identify species in Schima (Yu et al., 2022). However, next generation barcodes identified only three out of 11 Schima species, whereas standard DNA barcodes identified no species. Thus, the efficacy of current DNA barcodes for species identification in Schima is limited.

Here we used genome skimming data from 47 individuals of 10 species of Schima to determine the efficacy of Skmer in identifying species in a taxonomically difficult genus. We also explored the sequencing depths required for accurate species identification.

2. Material and methods 2.1. Sample collection

In total, 47 individuals representing 10 out of 13 Chinese species of Schima were represented. Leaves of eight individuals representing Chinese species of Schima were newly collected. We also used genome skimming data from 39 individuals from our previous study, representing five additional Schima species (Yu et al., 2022). These data represent at least 2× sequencing depth according to the genome size of Schima argentea as 900 Mb (unpublished data). The numbers of individuals used for each species varied from two to nine, depending on the ranges of their geographic distribution. We also collected two species of Camellia (C. mingii and C. connata) as outgroups were collected. Fresh leaves were dried by silica gel and stored at low temperature. Voucher specimens from newly collected individual were deposited in the Herbarium of Kunming Institute of Botany, Chinese Academy of Sciences (KUN) (Table S1).

2.2. DNA extraction and data collection

Total DNA was extracted from about 100 to 200 mg of leaf material and purified using the Plant Genomic DNA Extraction and Purification Kit. The quality of DNA was examined by agarose gel electrophoresis and NanoDrop® ND-1000 Spectrophotometer. High-quality DNA (> 1.5 μg) was used for subsequent library construction and sequencing. For library construction, DNA was first randomly sheared by a Covaris crusher, before the desired length of DNA fragments (300–350 bp) was recovered by electrophoresis. The NEBNext Ultra II DNA Library Prep Kit for Illumina was used to construct a library. Briefly, DNA fragments were end repaired, polyA tails and sequencing junctions added. DNA fragments were then purified prior to PCR amplification. Finally, the target fragments were subjected to paired-end sequencing using the Illumina NovaSeq 6000 platform with a unidirectional read length of 150 bp. Sequencing data for each individual was about 2 Gb.

2.3. Sequence dataset preparation

Fastp v.0.23.2 (Chen et al., 2018) was used for quality control to remove contamination and filter out low-quality reads, with parameters set to default. To explore the minimum sequencing depth require for species discrimination using Skmer in Schima, we used BBtools (Bushnell, 2014) to randomly extract sub-datasets from our dataset derived from genome skimming in 47 individuals. Two different sizes of sub-datasets were generated: 0.5 Gb (0.5×) and 1 Gb (1×).

2.4. Distance calculation

To calculate mismatch rate, each dataset was analyzed separately according to the Skmer protocol. Sequencing error and coverage were estimated by calculating k-mer frequency curves in JellyFish (Marcais and Kingsford, 2011). Shallow sequencing of each genome was collected and used to create separate k-mer sets. The similarity between any two sets (k-mer sets) was computed as the J (Jaccard index J) using the hashing technique by Mash (Ondov et al., 2016). The distance matrix is constructed with the Skmer default parameters, generating a master distance matrix and 100 subsample distance matrices. Finally, the estimated distance was corrected using the Skmer command to obtain values extremely close to the results.

2.5. Phylogenetic analysis

For building a phylogenetic tree, FastME 2.0 (Lefort et al., 2015) was used to infer the backbone tree, generating one main matrix tree file and 100 sub-repeat matrix tree files. Then the Maximum-Likelihood (ML) trees were built using RAxML-1.2.0 (Stamatakis, 2014), with 1000 rapid Bootstrap replicates to compute Bootstrap support (BS) and to search for the optimal ML tree. Species identification was evaluated as being successful if all conspecific samples clustered in a clade in the phylogenetic tree with ML bootstrap support above 70% (David and James, 1993; Pamela and Douglas, 2003).

3. Results 3.1. Species discrimination of Schima estimated by Skmer analysis

Skmer identified six of ten Schima species (i.e., Schima wallichii, S. khasiana, S. sericans, S. noronhae, S. crenata, and S. remotiserrata) with high support based on the 2× (ca. 2 Gb) genome skimming data (Table 1). Specifically, Skmer divided 47 Schima samples into two main clades (BS = 100%) (Fig. 2). Clade A was formed by 27 individuals from seven species and one subspecies collected in southwestern China, Myanmar and Nepal. Eight individuals of S. wallichii clustered into a monophyletic subclade with strong support (BS = 100%), which were distributed in southwestern China (four individuals), Myanmar (two individuals), and Nepal (two individuals). Similarly, three individuals of S. khasiana (collected from Yunnan, China) formed a subclade with 100% bootstrap support. Another monophyletic subclade consisted of three individuals of S. sericans var. sericans and two individuals of S. sericans var. paracrenata. Two individuals of S. noronhae collected from Southwestern China formed a monophyletic subclade but was nested within S. sinensis, where five individuals of S. sinensis were sub-divided into two groups respectively corresponding to Yunnan (two individuals) and Sichuan (three individuals). Three individuals of S. argentea and one individual of S. brevipedicellata clustered in a subclade together.

Clade B was composed of 20 individuals of five species mainly from southeastern China. Three individuals of Schima remotiserrata and four individuals of S. crenata (two from China and another two from Cambodia) formed two distinct monophyletic subclades, respectively, with strong support (BS = 100%). However, nine individuals of S. superba were divided into two groups with 100% bootstrap support. One group consisted of four individuals of S. superba, while the other group of five individuals of S. superba was accompanied by one individual of S. argentea. Interestingly, the other two individuals of S. brevipedicellata, which were collected from southwestern China (Yunnan) and southeastern China (Guangdong), clustered into a subclade, with the remaining one individual of S. argentea from south China (Guangxi).

3.2. Sequencing depths for Schima

To test suitable sequencing depths of Skmer in Schima, two fixed datasets (0.5 Gb and 1 Gb) were extracted randomly from the original dataset (about 2 Gb, ca. 2×). Two phylogenetic trees constructed by 0.5 Gb (0.5×) or 1 Gb (1×) dataset exhibited similar topologies to that constructed by the original dataset (Figs. 2-4). Skmer successfully discriminated six species of Schima (Schima wallichii, S. khasiana, S. sericans, S. noronhae, S. crenata, and S. remotiserrata) using each dataset (0.5×, 1× and 2× depth). For datasets based on genome skimming at the original (2×) and 0.5× sequencing depths, all 47 individuals of 10 Schima species were divided into two main clades (Figs. 2 and 4), which corresponded well with their geographic distribution. In contract, fixed 1× sequencing depth yielded three main clades, with S. wallichii clade not clustering with other species found in clade A (Fig. 3). As expected, bootstrap support values for some clades for most clades or subclades were higher in the original 2× dataset than for the 0.5× and 1× datasets. However, all three phylogenetic trees formed six monophyletic subclades. Surprisingly, these six monophyletic subclades obtained complete 100% bootstrap support, even in the 0.5× dataset (Fig. 4).

4. Discussion 4.1. Skmer approach was efficient in species identification of Schima

Next generation DNA barcodes have been used well in species identification of some taxa, such as Taxus (Taxaceae) (Fu et al., 2019), Panax (Araliaceae) (Ji et al., 2019), and Calligonum (Polygonaceae) (Song et al., 2020). However, these barcodes have not been effective for species identification in taxonomically challenging taxa, such as Rhododendron (Ericaceae) (Fu et al., 2022) and Fargesia (Poaceae) (Lv et al., 2023), or our previous study in Schima, where only 3 out 11 Schima species were identified (Yu et al., 2022). Here, Skmer approach, based on the unassembled data, improved species identification rate in Schima, identifying six out of ten species (60%) (Table 1). Our results indicate that Skmer is more efficient than the next generation barcodes at species discrimination in Schima. This finding is similar to a previous study that showed that the Skmer approach improved upon the species identification rate of next generation barcodes (72% vs. 68%) in Cymbidium (Zhang et al., 2023). Our finding is also consistent with previous studies that showed that the discrimination rate of Skmer approach (77.7%) is the same as that of plastid genomes, but relatively higher than that of nrDNA cistrons (62.5%) (Mo et al., 2023).

The improved identification rate of Schima species might be attributed to an extended exploration of the unassembled nuclear reads from genome skims. Nuclear genes have been used to infer of angiosperm phylogeny (Hughest et al., 2006; Zimmer and Wen, 2012; Dong et al., 2021). However, nuclear genome sequencing is expensive and requires significant computational power. For example, the minimum sequencing depth required that falls in the range of at least 50× for a relatively straightforward diploid organism (Sohn and Nam, 2018). As a result, a fully assembled nuclear genome may not be the most suitable option for widespread species identification purposes. Genome skimming generates a large amount of nuclear genomic sequence data that are not used for species discrimination and plant phylogenetic studies, but instead discarded. In fact, it was ever proposed that ca. 1× sequencing depth or 2 Gb genome skimming data might be sufficient for species-level identification in most cases when Skmer was developed as a new approach (Sarmashghi et al., 2019; Bohmann et al., 2020). The Skmer approach allows for species identification based on the genome skimming data ranging from 0.1 Gb (ca. 0.50×–1.68×) to 4 Gb (ca. 0.93×–2.78×) with an increasing depth corresponding to a strong support (Xu et al., 2022). In our study, the phylogenetic trees based on data ranging from 0.5 Gb (ca. 0.5×), 1 Gb (ca. 1×), and 2 Gb (ca. 2×) were similar in topology, with each dataset identifying six Schima species with strong support (BS = 100%). Therefore, Skmer might be an effective, credible, and economic approach for exploring the potential of nuclear genomic data.

Interestingly, our findings based on Skmer analyses conflict with those of our previous study in which we used next generation barcodes. For example, next generation barcodes grouped all individuals of Schima wallichii into two subclades (Yu et al., 2022), however, Skmer approach grouped these same individuals into one highly supported monophyletic subclade (BS = 100%) (Fig. 2). Similarly, our plastomic tree failed to group S. khasiana into a subclade (Yu et al., 2022), whereas Skmer tree placed all three individuals within a monophyletic subclade (BS = 100%) (Fig. 2). The same discrepancies occurred when we used these approaches to identify S. noronhae, S. sericans, and S. remotiserrata. These results indicate that previous failure of plastome for species discrimination of Schima might be caused by hybridization or introgression among species, and thus the maternally inherited plastome could not trace the underlying species boundary (Small et al., 2004; Greiner et al., 2015). However, in this study, the genome skims based on nuclear reads in Skmer analyses could provide additional information for species identification and improve the species discrimination rate, avoiding limitations coming from chloroplast capture (Percy et al., 2014; Yi et al., 2015; Ogishima et al., 2019; Liu et al., 2020). Our previous study also found ancient introgression between the common ancestor of Gordonia and that of Schima (Zhang et al., 2022). It remains unclear whether potential hybridization has played a role in Schima speciation, although we suspect this might be an important factor that limited the efficacy of plastome data in species identification in the genus (Twyford and Ennos, 2012; Richard and Erica, 2014). Even though Skmer analysis could be an alternative of next generation barcodes, using more single/low copy genes as next generation nuclear barcodes would provide more robust framework of species discrimination, as numerous studies have suggested (Ruhsam et al., 2015; Ji et al., 2019; Fu et al., 2022; Yu et al., 2022; Zhang et al., 2023).

4.2. Skmer approach is reliable in species discrimination of Schima

Species discrimination of Schima by Skmer corresponded well with morphological and geographic evidence. The six Schima species identified by Skmer included two species endemic to China (i.e., S. sericans and S. remotiserrata) and four species widely distributed in Yunnan, China and Southeast Asia (Figs. 2-4). Previous studies have divided S. sericans, which is distributed in Xizang and Yunnan, China, into two varieties only based on morphological data, including the presence or absence of trichome on current year branchlets, petioles, pedicels and outside of sepals (Ming and Bartholomew, 2007). Furthermore, S. sericans var. sericans and S. sericans var. paracrenata are distributed in overlapping ranges, with the latter narrowly distribution in the Gongshan mountain of Yunnan Province in China, suggesting that relying solely on morphology may be insufficient for species delimitation. Skmer analysis showed that two individuals of S. sericans var. paracrenata were nested with three individuals of S. sericans var. sericans and strongly formed a subclade. Thus, our findings indicate that the two varieties of S. sericans should be merged as a single species. We also found that three individuals of S. remotiserrata also form a monophyletic subclade (Figs. 2-4). These individuals, narrowly distributed in the southeastern China, share morphological traits (e.g., leaf adaxially glabrous, margin apically sparsely obtusely serrate).

The identification of Schima crenata was strongly supported, as two individuals from China and two from Cambodia clustered into a subclade (Figs. 2-4). This finding is consistent with our previous analysis based on plastomic next generation barcodes (Yu et al., 2022). In China, S. crenata is only found on Hainan Island, although it is common throughout the Mainland Southeast Asia, Malaysia and Indonesia. It was distinguished morphologically from other Schima species by its long pedicel (4–6 cm), leaf margin apically undulate crenate, and ovary tomentose (Ming and Bartholomew, 2007).

Schima wallichii is widely distributed in southwestern China, Bhutan, N India, Laos, Myanmar, Nepal, Thailand, and Vietnam (Ming and Bartholomew, 2007). It was previously treated as the only one species in Schima, but regarded as a 'complex-polymorphic' species, which could be divided into nine geographically separated subspecies and three varieties (Bloembergen, 1952). The chromosome number (2n = 36) of S. wallichii was recorded (Bezbaruah, 1971), which was congruent with that of the whole genus of Schima (Yang et al., 2004). Previous studies used next generation barcode to divide ten individuals of S. wallichii from southwestern China, Myanmar, and Nepal into two subclades (Yu et al., 2022). Here, Skmer analysis grouped eight of these individuals into one subclade with strong support (Figs. 2-4). Our results support prior taxonomic treatment that regard S. wallichii as a good species, considering its distinct features, including the current year branchlets pubescent, leaf blade with entire margin, capsule subglobose with white lenticellate, etc. (Ming and Bartholomew, 2007).

Schima khasiana is mainly distributed in Yunnan of China, extending to Myanmar, India, as well as to Vietnam. It was distinguished by the largest flower (diameter > 6 cm) and the largest fruit (diameter ca. 2.5–4 cm) within Schima. Skmer analysis grouped three individuals into a monophyletic subclade (Figs. 2-4). In addition, two individuals of S. noronhae, which mainly occurred in Southeastern Asia, clustered as a subclade by Skmer (Figs. 2-Fig. 4). This identification is consistent with morphological trait, i.e., pale pink flowers with a long pedicle (> 5 cm) (Ming and Bartholomew, 2007).

Skmer did not successfully identify 22 individuals sampled from other four Schima species (i.e., S. superba, S. argentea, S. brevipedicellata and S. sinensis; Figs. 2-4). The delimitation of these four species has been controversial due to their wide range distribution and morphological overlap (Ming and Bartholomew, 2007). For example, S. superba is distributed widely from the southern and southeastern China and Taiwan Island in China to Japan. It could not be easily identified by morphology because its leaves were characterized by margin undulately obtusely crenate from basal 1/2 apically but varied heavily. Furthermore, S. superba shares features with S. remotiserrata and S. crenata, e.g., variously toothed leaf margins and glabrous abaxial surface. Skmer grouped and mixed four individuals of S. argentea and three S. brevipedicellata in two distant subclades, with one individual of S. argentea grouping within the S. superba subclade (Figs. 2-4). The geographical distribution of S. argentea and S. brevipedicellata overlap in Jiangxi and Yunnan of China as well as in Vietnam. These two species also share common features, such as leaf margin entire and abaxial surface glaucous, but are distinguished by the evident length of pedicle. The pedicle of S. argentea (1–1.5 cm) is longer than that of S. brevipedicellata (< 1 cm). S. sinensis is endemic to the southwestern China, and is characterized by solitary flowers, a trait shared with S. khasiana (Ming and Bartholomew, 2007). Skmer also failed to successfully identify individuals of S. sinensis, which clustered into a subclade but with two individuals of S. noronhae nested within (Figs. 2-4). S. sinensis is morphologically similar to S. noronhae, with both having leave that are glabrous on the abaxial surface and long pedicel (> 3 cm).

The inability of Skmer to discriminate individual of four Schima species might be attributed to their low genetic variation which is reflected in their overlapping morphological characters. In addition, Schima species may have undergone interspecific hybridization or rapid radiation. Despite our progress, taxonomy and species identification of Schima remains problematic. Future taxonomic studies should integrate morphological and molecular data as well as sample populations that cover the entire geographic distribution of the species (Dayrat, 2005; Cleusa et al., 2018; He et al., 2022).

5. Conclusions

In the study, Skmer identified six out of ten Schima species. Furthermore, species identification was consistent with the unassembled genome skims ranging from 0.5× (ca. 0.5 Gb), 1× (ca. 1 Gb), to 2× (ca. 2 Gb). The species identified by Skmer were also consistent with morphological data and corresponded to their distribution patterns. Our results effectively improve the species identification rate of Schima to 60% compared with the previous next generation barcoding research (27.27%). Additionally, our comparative results showed that 0.5× sequencing depth data successfully captured sufficient nuclear genomic data for downstream species identification analysis in the case study of Schima. Therefore, Skmer can make full use of the scattered nuclear genomic data from shallow sequencing depth and thus improve species identification rates. However, comprehensive taxonomic revisions for the objective taxa would serve as significant benchmark for Skmer performance.

Acknowledgements

This work was supported by National Natural Science Foundation of China (No. 32070369), the Youth Innovation Promotion Association CAS of China (No. 2021393), the Yunnan Revitalization Talent Support Program "Young Talent" Project, the Applied Fundamental Research Foundation of Yunnan Province (202301AT070308) and the Fund of Yunnan Key Laboratory of Plant Reproductive Adaptation and Evolutionary Ecology (YNPRAEC-2023006). The authors are grateful to Prof. Liang Fang (Jiujiang University) and Prof. Xue-Jun Ge (South China Botanical Garden, CAS), to Prof. Hua Peng and colleagues En-De Liu, Ting Zhang, Ji-Dong Ya, Yong-Jie Guo, Cheng Liu and Li Chen (Kunming Institute of Botany, CAS) for assistance with sample collection. We thank the staff at KUN for providing access to study specimens of Schima, and the iFlora High Performance Computing Center of the Germplasm Bank of Wild Species (iFlora HPC Center of GBOWS, KIB, CAS) for computing.

Data availability statement

The data set supporting the conclusions of this article is available in the FigShare (https://doi.org/10.6084/m9.figshare.24848427).

CRediT authorship contribution statement

Han-Ning Duan: Writing – original draft, Formal analysis. Yin-Zi Jiang: Writing – review & editing. Jun-Bo Yang: Writing – review & editing. Jie Cai: Writing – review & editing. Jian-Li Zhao: Writing – review & editing. Lu Li: Writing – review & editing. Xiang-Qin Yu: Writing – review & editing, Conceptualization.

Declaration of competing interest

The authors have no conflict of interest.

Appendix A. Supplementary data

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