1. Introduction

Celtis L. is a genus of 60–70 species of trees, or rarely shrubs, commonly known as hackberries or nettle trees, widespread in warm temperate regions of the Northern Hemisphere and parts of central Africa and South America (Keerler, 1902; Hwang et al., 2003; Martins et al., 2015). Celtis species are grown as ornamental plants, used in traditional medicines (Kim et al., 2005), provide important wildlife habitat and winter food to birds and mammals (Correll and Johnston, 1970), and serve as pollen sources (Krajicek and Williams, 1990). Despite their high economic and ecological importance, the genomic resources necessary for examining the phylogenetic relationships of Celtis are currently unavailable.

Plants in Celtis show a great diversity of growth forms, along with significant divergence in physiological characters caused by strong ecological differences (Whittemore, 2005). Furthermore, there are apparent intergradations between species, which many scholars have deduced are likely a result of interspecific hybridization (Fernald, 1950; Whittemore and Townsend, 2007). Consequently, the genus Celtis has long been known for its taxonomic difficulty and is one of the most poorly understood plant groups (Hooker, 1885; Whittemore, 2005; Sattarian, 2006; Lee et al., 2011). The genus was previously put in either the elm family (Ulmaceae) or a separate family (Celtidaceae). However, the Angiosperm Phylogeny Group Ⅲ (APG Ⅲ) system transferred it into an expanded hemp family (Cannabaceae) (Angiosperm Phylogeny Group Ⅲ, 2009). Subsequent phylogenetic studies have confirmed that Celtis is a member of Cannabaceae and closely related to Trema Lour., Parasponia Miq., Chaetachme Planch., Pteroceltis Maxim., Humulus L. and Cannabis L. (Yang et al., 2013; Zhang et al., 2018). To date, few studies have analyzed the phylogenetic relationships within Celtis. A molecular phylogeny of Cannabaceae based on four plastid loci included six Celtis accessions, which formed two sister clades, with C. madagascariensis Sattarian, C. iguanaea (Jacq.) Sarg. and C. ehrenbergiana (Klotzsch) Liebm. 2 grouped into one clade, and C. biondii Pamp., C. tournefortii Lam., C. sinensis Pers. and C. ehrenbergiana 1 clustered into the other (Yang et al., 2013). Chen et al. (2017) reported a new species (C. neglecta Zi L. Chen & X. F. Jin) from Zhejiang province, and revealed a topology of ((C. neglecta, C. julianae C.K.Schneid.), C. vandervoetiana C.K.Schneid.) with moderate support based on the atpB-rbcL and trnL-F sequences. Clearly, further phylogenetic analysis is required to improve our understanding of the phylogeny and evolution of Celtis.

Traditionally plant lineages have been inferred using plastid introns and inter-genic spacer regions with highly variable loci (Baldauf et al., 2000; Moncalvo et al., 2002). More recently, whole plastomes have increasingly been used for phylogenetic analyses (Williams et al., 2016; Liu et al., 2017; Gu et al., 2019). Studies have shown that whole plastomes are effective when resolving recalcitrant phylogenetic relationships at shallow taxonomic levels (Givnish et al., 2015; Duvall et al., 2016) and at deeper levels (Jansen et al., 2007; Moore et al., 2010). However, the uniparental inheritance of plastomes means that they are not ideal for identifying phylogenetic relationships of species that have undergone interspecific hybridization (Gruenstaeudl et al., 2013; Johnson et al., 2019).

Additional molecular resources can be obtained by genome skimming, which is one of the genome-level sequencing methods first identified by Straub et al. (2012). The high-copy fraction of a genome (such as plastome, mitogenome and nuclear repetitive elements) is sequenced through sampling of low-coverage genome data, and partial sequences with low-copy nuclear loci recovered from genome skimming are provided to design PCR primers or probes for amplification (Cronn et al., 2012; Straub et al., 2012). Moreover, the internal and external transcribed spacers (ITS and ETS), as well as polymorphic nuclear simple sequence repeats (nSSRs) and plastid genomic hotspots are commonly obtained using this technique (Liu et al., 2018, 2020). Genome skimming has been widely used at different taxonomic levels, for intraspecific 'ultra-barcoding' (Kane et al., 2012), intergeneric (Mcpherson et al., 2013; Li et al., 2017; Liu et al., 2018) or phylogenomic studies (Besnard et al., 2013; Malé et al., 2015; Xu et al., 2018; Yu et al., 2018).

In this study, our aim was threefold. First, we aimed to reveal patterns of plastome evolution in Celtis. For this purpose, we used genome skimming sequencing to reconstruct the plastomes of five Celtis species widely distributed in East Asia. We then examined whether plastome data effectively resolves phylogenetic relationships within Celtis. Finally, we developed universal genomic resources (e.g., plastid hotspot regions, polymorphic nuclear SSRs, single or low copy genes) for the genus. Our results provide powerful tools for future studies on the phylogeny, plastome evolution, and population genetics of Celtis.

2. Materials and methods 2.1. Plant materials

Five Celtis species (C. biondii, C. cerasifera C.K.Schneid., C. julianae, C. koraiensis Nakai and C. sinensis) were sampled at Jinming Campus of Henan University (114°18′27.96″E, 34°49′17.59″N), Kaifeng, Henan, China. Fresh leaves of each species were collected for genomic DNA extraction using modified CTAB reagents (Plant DNAzol, Shanghai, China) according to the manufacturer's protocol. Voucher specimens for the five species were deposited at the Herbarium of Henan University (Table 1). For comparative analysis of plastomes, genomic resources development and phylogenetic inference, we downloaded publicly available plastomes of two additional Celtis species (C. biondii and C. tetrandra) (Zhang et al., 2018; Wang et al., 2019) from GenBank database.

2.2. Illumina sequencing, genome assembly and annotation

High quality DNA was sheared to yield fragments less than or equal to 800 bp, and fragment quality was checked using Agilent Bioanalyzer 2100 (Agilent Technologies). The 500 bp short-insert length paired-end library was prepared and sequenced on an Illumina HiSeq X10 to obtain reads of 150 bp at Beijing Genomics Institute (BGI, Wuhan, China).

Raw reads were first screened for Phred score < 30 to remove low quality sequences. All remaining reads were assembled into contigs implemented in the CLC Genomics Workbench (CLC Inc. Aarhus, Denmark). The parameters set in CLC were as follows: 200 bp for minimum contig length, 3 for deletion and insertion costs, bubble size of 98, 0.9 for length fraction and similarity fraction, 2 for mismatch cost. Principal contigs representing the plastome were then screened from the total contigs using a BLAST (NCBI BLAST v.2.2.31) search, with the published plastome of Celtis biondii (GenBank accession number: MH118119) as a reference. Mapped contigs were oriented and ordered according to the reference genome to yield the draft genomes, and plastomes were finalized by re-mapping cleaned reads to the draft genomes.

We used Geneious R11 (Kearse et al., 2012) to annotate the plastomes, and putative starts, stops, and intron positions were identified by comparison with homologous genes of the reference genome. The tRNA genes were verified with tRNAscan-SE v.1.21 (Schattner et al., 2005) using the default setting. We drew the graphical map of the annotated circular plastomes using the OrganellarGenomeDRAW program (OGDRAW, Lohse et al., 2013).

2.3. Comparative analysis of plastomes in Celtis

To analyze sequence variation within the genus, multiple sequence alignments of the seven Celtis plastomes were performed in MAFFT v.7.017 (Katoh and Standley, 2013) with standard parameters. Alignments were visually inspected and manually adjusted in Geneious R11. Sequence identities of the seven Celtis plastomes were plotted using the mVISTA program in LAGAN mode (Frazer et al., 2004). Plastid DNA rearrangement analyses were performed via whole genome alignment in Mauve (Darling et al., 2004).

2.4. Genomic resources development for Celtis

To identify plastid hotspot regions within the genus Celtis, multiple alignments of the seven plastomes were performed using MAFFT v.7.017. Nucleotide diversity (Pi) was determined by calculating the total number of mutations (Eta) and average number of nucleotide differences (K) using DnaSP v.5.0 (Librado and Rozas, 2009).

Plastid and mitochondrial sequences of the Celtis were removed from the assembled contigs using NCBI BLAST with the plastome of C. biondii and mitochondrial genome of Cannabis sativa L (GenBank accession number: KU310670) as references. Candidate polymorphic nSSRs within the genus Celtis were screened based on multiple assembled sequences using the CandiSSR v20170602 software (Xia et al., 2016) set at default parameters. Primers for each target SSR were automatically designed based on the Primer3 package built-in installation to the pipeline (Untergasser et al., 2012).

For targeted sequence capture, we used probes from 353 putative single-copy protein-coding genes from 42 angiosperms (Johnson et al., 2018), which have been shown to be universal in 283 species (the information on sequences are published at github. com/mossmatters/Angiosperms353). Our cleaned sequence data for each Celtis species was searched using pipeline HybPiper v.1.2 (Johnson et al., 2016) with default parameters. Under ideal conditions, searches identify only single contigs for each Celtis gene sequence.

2.5. Phylogenetic analysis

Phylogenetic inferences were conducted using Bayesian inference (BI) and Maximum-Likelihood (ML) methods for three different data sets: (1) a combined matrix of four plastid regions (psbA-trnH, matK, rbcL and trnL-trnF; mostly downloaded from GenBank, see Table S1) for 32 Celtis species (to reduce missing data, we only used species with at least two regions); (2) whole plastome sequences of seven Celtis accessions; (3) 78 protein-coding gene sequences commonly shared among the seven accessions. Two species from Cannabaceae, Trema orientalis (L.) Blume and Pteroceltis tatarinowii Maxim., were used as outgroups according to Zhang et al. (2018).

The best-fit nucleotide substitution model (GTR + G for three datasets) was determined by jModelTest v.2.1.4 (Posada, 2008). ML analysis was performed with RAxML-HPC v.8.1.11 on the CIPRES cluster (Miller et al., 2010), with 1000 bootstrap replicates. BI analysis was performed using MrBayes v.3.2.3 (Ronquist and Huelsenbeck, 2003). The Markov chain Monte Carlo (MCMC) algorithm was run under the parameters set as follows: five million generations, every 1000 generations for tree sampling, and convergence was determined when the average standard deviation of split frequencies reached 0.01 or less. The first 25% of the trees were discarded as a burn-in to generate the consensus tree.

3. Results and discussion 3.1. Illumina sequencing and genome features

More than 10 Gbp of high-quality data with clean reads (Q30 > 95%) were generated from the Illumina sequencing platform for each species. Between 419, 052 (Celtis koraiensis) to 527, 009 (C. cerasifera) contigs were generated after de novo assembly. All plastomes were reconstructed from three initial contigs corresponding to the large single-copy region (LSC), small single-copy region (SSC), and inverted repeat regions (IR_B/IR_A), with no Ns or gaps detected. All five plastome sequences have been submitted to GenBank (accession numbers shown in Table 1).

Plastomes of all five Celtis species possessed a typical quadripartite structure similar to the majority of land plant plastid genomes (Li et al., 2017; Liu et al., 2017). The genome size ranged from 158, 989 bp in C. biondii to 159, 082 bp in C. koraiensis (Fig. 1, Table 1), consisting of two copies of IR regions (26, 861–26, 862 bp) separated by a LSC region (86, 075–86, 171 bp) and an SSC region (19, 129–19, 190 bp). All five plastomes contained 112 unique genes, including 78 protein-coding genes, 30 tRNA genes, and four rRNA genes with 20 genes duplicated in IR regions (Table 2). Of the 112 genes, nine protein-coding genes and six tRNA genes contained one intron, and three protein-coding genes (rps12, clpP and ycf3) contained two introns.

3.2. Comparisons of the plastomes in Celtis

Multiple sequence alignment of seven Celtis plastomes showed no evidence of genome structure variation (e.g., gene loss, translocation, inversions, or genome rearrangement) (Fig. 2). Plastome sequences were relatively conserved especially the IR regions were more conserved than both single copy regions (Fig. 3).

Most angiosperm plastomes range between 135 and 160 kb (Palmer, 1985). The most typical cause of variation in the length of plastomes is the contraction or expansion of the IR into or out of adjacent single copy regions (Plunkett and Downie, 1999). Therefore, we compared the exact IR border positions of seven Celtis plastomes and the positions of adjacent genes (Fig. 4). The genes rps3-rpl22-trnH and ycf1-ndhF were located in the boundaries of LSC/IR and SSC/IR regions. The ycf1 gene spanned the SSC/IR_A region. In C. biondii (MH118119) and C. tetrandra, the pseudogene fragment of ^ψycf1 was 1093 bp, whereas in the remaining Celtis plastomes it was 1065 bp. In all seven Celtis plastomes, the ndhF gene shared a 28 bp sequence with the ^ψycf1 gene. The rps3 gene spanned the LSC/IR_B region, resulting in 124 bp to 139 bp nucleotides located in IR_B region. In C. julianae, the rpl22 gene was divided from the LSC/IR_B boundary by a 93 bp spacer and this spacer was 87 bp in the remaining species. For all species except C. tetrandra, the trnH-GUG gene was separated from the IR_A/LSC boundary by a spacer of 3 bp. Compared to plastomes of other genera within Cannabaceae (Zhang et al., 2018), the genome size and genome boundary structure of Celtis plastomes are relatively conserved.

3.3. Development of genetic resources for Celtis

Although plastomes are usually an effective tool to investigate species divergence, identify species and trace demographic history (Thomson et al., 2010), traditional screening approaches are time consuming and inefficient (Alkan et al., 2011). Here, we screened divergence hotspots in the plastomes of seven Celtis accessions by comparing coding genes, non-coding regions, and intron regions. A total of 117 loci (56 inter-genic spacers, 45 coding genes, and 16 intron regions) with a length more than 200 bp were generated (Fig. 5). Nucleotide diversity (Pi) values for each locus ranged from 0.00022 (ycf2) to 0.02121 (psbC-trnS). Six of these loci with the Pi ≥ 0.01, including trnH-psbA, psbA-trnK, trnG-trnR, psbC-trnS, cemA-petA and rps8-rpl14, showed relative high nucleotide diversity values and can be used as highly informative phylogenetic markers for the genus Celtis or higher taxonomic levels.

Genomic SSR (gSSR) markers have higher levels of polymorphism than EST-SSRs, and have thus garnered more attention (Bae et al., 2015). In this study, a total of 11, 559 candidate polymorphic nSSRs were identified for the genus Celtis. After discarding loci with either no available primers or sequence similarity < 99%, we obtained 4097 polymorphic nSSRs with a standard deviation between 0.4 and 3.7 (Table S2, Table S3). Of these polymorphic nSSRs, di- (2940), tri- (997), tetra- (140), penta- (13) and hexanucleotides (7) accounted for 71.76%, 24.33%, 3.42%, 0.32% and 0.17%, respectively (Fig. 6). Dinucleotide repeats included AT/TA (64.56%), CT/TC (15.56%) and AG/GA (15.07%) motifs. SSRs in genus Celtis are strongly biased towards AT-rich repeat motifs, not only in dinucleotides but also in tri-, tetra-, penta- and hexanucleotide repeats.

Several plastid markers are frequently used for phylogenetic studies, including single genes (Soltis et al., 1993; Cameron et al., 1999), whole plastid exomes (Gitzendanner et al., 2018; Medina et al., 2018) or whole plastome sequences (Bernhardt et al., 2017; Liu et al., 2017). Importantly, phylogenetic relationships based on plastome data alone may not reflect the evolutionary divergence of plants in Celtis, which have likely undergone extensive interspecific hybridization (Fernald, 1950). We used the HybPiper pipeline to identify putative single-copy genes for phylogenetic analysis of genus Celtis. In total, 353 target genes were successfully mapped from the raw clean reads, and contigs were formed for C. biondii (11), C. cerasifera (14), C. julianae (26), C. koraiensis (11) and C. sinensis (27) (Fig. 7; Table S4). We recovered coding sequences with a length of at least 50% of the target length for C. biondii (2), C. cerasifera (4), C. julianae (8), C. koraiensis (3) and C. sinensis (6). In addition, we recovered three coding sequences (TJLC, THHD and YKQR) in C. sinensis longer than the corresponding target genes (Table S4). Although the gene recovery rate appears relatively low and no recovered genes were shared in the five sequenced Celtis species, we were still able to obtain homologous genes within the genus through PCR amplification using universal primers. In combination with plastome data, these genes will provide useful molecular tools for further studies on the history of species divergence in Celtis.

3.4. Phylogenetic inference

ML and BI analyses based on a combined matrix of four plastid regions (3940 bp) for all sampled species generated identical trees (Fig. 8). The 32 Celtis species were divided into five major clades, but the ML bootstrap support (BS) and BI posterior probability (PP) for most clades were weak. Nine of the eleven species distributed in China were assigned to two clades (Clade A and Clade D), and the six species with available plastomes were assigned to one subclade of Clade A.

We also used a plastome data set and a concatenated matrix of 78 shared protein-coding genes for ML and BI analyses. The length of aligned whole plastome sequences (nine individuals) was 166, 417 bp, and contained 9826 polymorphic sites, which can be further divided into 6548 singleton variable sites and 3278 parsimony informative sites. The concatenated matrix of 78 shared protein-coding genes was 68, 849 bp in length, comprising 2434 singleton variable sites and 1325 parsimony informative sites. The tree topologies recovered from both the whole plastome and the concatenated data sets were completely coincident (Fig. 9). The monophyly of Celtis had maximal support (MLBS = 100%, BIPP = 1). Within Celtis, C. biondii was first to diverge from the six representative species, and the remaining five species formed the topology ((((C. koraiensis, C. sinensis), C. tetrandra), C. julianae), C. cerasifera), with strong support.

To eliminate sampling interference, we removed all Celtis species except for six species (seven accessions), which were retained to reconstruct the phylogenetic tree based on four plastid regions (Fig. S1). The interspecific relationships of the six Celtis species resolved by plastomes and 78 shared protein-coding genes were consistent with the results recovered by the four combined plastid regions but with higher BS/PP support values.

In this study, we used plastome data to reconstruct the phylogenetic relationships of six Celtis species and compared these phylogenies to those generated using four common barcoding regions. The tree topologies inferred from the whole plastomes and 78 shared protein-coding genes are identical, and every branch of these trees is supported by full BS/PP values (Fig. 9). We also found that the phylogenetic relationships of the six species were consistent with that revealed by the combined matrix of four plastid regions (Fig. S1), although these data generated trees with much lower BS/PP values. Thus, our findings indicate that whole plastome or CDS genes have more informative sites and are more effective at resolving the phylogeny of Celtis.

4. Conclusions

In this study, five Celtis species were sampled for genome skimming. Combined with two published plastomes, six Celtis species (seven accessions) were used to generate abundant genetic resources, including plastid hotspots, polymorphic nuclear SSRs and low or single copy gene fragments. All the sequenced plastomes had a typical quadripartite structure with similar size and organization to other sequenced angiosperms. In addition, phylogenetic relationships of Celtis species were well-resolved based on whole plastomes or commonly shared protein-coding genes. Our phylogenetic tree of the genus Celtis and the genetic resources developed herein will benefit further studies of phylogeny, population genetics, and conservation biology for this important woody genus.

Author contributions

PL and LXL designed the study. LXL conducted the sampling. LXL and YHZ produced and analyzed the data. LXL and PL wrote the manuscript. YHZ revised the manuscript. All authors approved the final manuscript.

Declaration of competing interest

The authors declare that they have no competing interests.

Acknowledgments

We sincerely thank Dr. James R.P. Worth and Ryan A. Folk for revising the manuscript. This research was supported by the National Natural Science Foundation of China (Grant Nos. 31900188, 31970225), Natural Science Foundation of Zhejiang Province (Grant No. LY19C030007).

Appendix A. Supplementary data

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