Genome survey of Zanthoxylum bungeanum and development of genomic-SSR markers in congeneric species

The genus Zanthoxylum, in the family Rutaceae, consists of approximately 250 species and is distributed worldwide [1]. Zanthoxylum bungeanum (ZB), also referred to as ‘Chinese prickly ash’, ‘Sichuan pepper’ is a representative species of Zanthoxylum in China and Southeast Asia (Figure 1). The pericarps of ZB have been widely used as a traditional culinary spice and Chinese herbal medicine for thousands of years due to its flavor and medicinal characteristics [2]. As an economically important species, the cultivation area of ZB and a closely related species, Zanthoxylum armatum (ZA), occupies approximately 1.67 million hectares and produces an annual output of 350000 tons of dried pericarp. The economic value of this agricultural industry generates more than 4 billion US dollars in China. This level of economic value has generated interest in increasing the molecular data available for ZB. Although some transcriptome information has been obtained from several tissues in ZB [3,4], data about the genome structure of ZB are lacking.

Because of recent advances in DNA sequencing techniques, draft genomes have been assembled for many plant species as resources for genomic and genetic research efforts [5–7]. However, the genomes of some plant species, particularly tree species, are highly heterozygous, have a complex genetic background, and have an unknown genome size. Therefore, before large-scale sequencing is initiated, basic biological characteristics of the target material are evaluated, such as chromosome ploidy analysis or low-coverage genome sequencing (also known as a genome survey), to gauge the complexity of the genome and provide a reference for whole-genome sequencing [8].

Genome surveys, which use next-generation sequencing (NGS), yield a large amount of genomic data in a rapid, cost-effective manner. Genomic data from genome surveys not only provide useful information on genome structure, such as an estimation of genome size, heterozygosity levels, and repeat contents but also establish a genomic sequence resource from which molecular markers can be developed [9–11].

Simple sequence repeats (SSRs), also referred to as microsatellites or short tandem repeats (STRs) are tandem repeats of 1–6 nucleotides that are widely distributed in eukaryotic genomes [12]. Depending on the location of an SSR, it can be classified as either a genomic-SSR (G-SSR) or an expressed sequence tag (EST)-SSR, which indicates whether the SSR is either in a non-coding region or a translated region, respectively [13]. SSR markers are used for genetic and evolutionary analysis, germplasm resource identification, genetic map construction, and marker-assisted selective breeding because of their wide distribution, high polymorphism, and co-dominance and because the cost of developing these markers is low [14]. Therefore, increasing the number of SSR markers for ZB will provide a useful resource for genetic research.

In the present study, the main objectives were (1) to obtain information about the genome size, GC content, repeat sequence rate and heterozygosity rate of ZB by a genome survey; (2) to identify SSRs in assembled genomic sequences of ZB; and (3) to develop and evaluate G-SSR markers to assess genetic diversity in ZB and ZA (a plant that is used in the same manner as ZB in Southwest China).

Healthy leaves were collected from one ZB adult plant from the Yangling comprehensive test and demonstration station of Northwest A&F University in Shaanxi, China. The leaves were cleaned with purified water, immersed in liquid nitrogen, and stored at −80°C until use. For polymorphic marker screening, ten ZB individuals (ZB01–ZB10) and five ZA individuals (ZA01–ZA05) were collected from six provinces in China (Table 1). All samples were collected as young leaves from three individual plants that were combined for DNA isolation. Genomic DNA was extracted using a plant DNA extraction kit (DP305; Tiangen, Beijing, China), and the quality and quantity of the isolated DNA were evaluated by 1% agarose gel electrophoresis and NanoPhotometer spectrophotometer (Implen NanoPhotometer, Westlake Village, CA, U.S.A.).

Genomic DNA samples were randomly sheared into two collections of average fragment size (250 or 350 bp) using ultrasound (Covaris, U.S.A.) and used to construct four libraries (two libraries from each fragment collection) for sequencing. Library construction and sequencing were performed by Beijing Novogene Biological Information Technology, Beijing, China (http://www.novogene.com/) using an Illumina HiSeq 2000 platform. According to the genome size of Zanthoxylum oxyphyllum (3.7Gb) [15], we estimated that the genome size of ZB was approximately 4 Gb. Consequently, to ensure at least 50× coverage of the ZB genome, 238.5 Gb raw data were generated. After filtering to remove adaptors, poly(N) sequences, and low-quality reads, the remaining clean reads were used for K-mer analysis. Based on the results of the K-mer analysis, 17-mers (K = 17) were used to estimate the genomic characteristics, including genome size, repeat sequence rate, and heterozygosity rate. Furthermore, clean data were assembled (K = 41) into contigs and scaffolds using the De Bruijn graph-based assembler, SOAPdenovo (Version 1.05, BGI, Beijing, China) [16]. The GC content was calculated with contigs longer than 500 bp. More details regarding the analysis procedures employed in this study have been described by Bi et al. [17].

SSR loci were detected from scaffolds longer than 1000 bp using MISA software (Microsatellite, http://pgrc.ipk-gatersleben.de/misa/). For mono-, di-, tri-, tetra-, penta-, and hexanucleotide SSRs, the search parameters were set with minimum repeat numbers of 10, 6, 5, 5, 5, and 5, respectively. Primers were designed using Primer 3.0 software (https://primer3.sourceforge.net/webif.php), following the following parameters: 100–300 bp for amplification length, 18–27 bp for primer size, 40–70% for primer GC content, and 57–63°C for melting temperature. The PCR reaction volume was 20 μl and consisted of 10 μl 2× Taq Master Mix (Vazyme, Nanjing, China), 5 μl genomic DNA (20 ng/μl), 1.5 μl (2 μmol/l) each forward and reverse primers and 2 μl ddH₂O. The PCR program was 5 min at 95°C; followed by 35 cycles of 95°C for 30 s, 55°C for 40 s, and 70°C for 40 s; and a final extension of 72°C for 8 min. The PCR products were separated by electrophoresis on an 8% denaturing polyacrylamide gel, visualized by silver staining, and fragment sizes were estimated using pBR322 DNA/MspI markers (Tiangen, Beijing, China).

Genetic diversity parameters, such as the number of alleles (N_a), observed heterozygosity (Ho), and expected heterozygosity (He) were calculated using POPGENE1.32 [18], and the polymorphism information content (PIC) values were calculated using the formula $PIC=1−Σfi2$⁠, where f_i is the frequency of the i-th allele [19]. The dendrogram of 15 individuals was constructed by UPGMA clustering using NTSYSpc2.10 software to reveal genetic relationships [20].

In the present study, 238.5 Gb of raw sequence data were generated by four small-insert libraries. After removing low-quality reads, 197.16 Gb of clean data were used for the K-mer analysis. In the four small-insert libraries, the Q20, Q30, and GC contents were 97.28–98.10%, 93.88–95.66%, and 38.38–38.68%, respectively. Moreover, the error rate was 0.02% for each library (Table 2). With an Illumina platform, the overall accuracy of the sequencing is indicated by having Q20 and Q30 values of at least 90% and 85%, respectively [21]. Therefore, the sequencing accuracy of the ZB genome survey in the present study was high.

In the 17-mer frequency distribution, the K-mer number was 176134142868, and the K-mer depth was 44. Based on the empirical formula (genome size = K-mer number/K-mer depth), the initial estimate of genome size of ZB is 4003.05 Mb. After excluding the effects of erroneous K-mers, the revised genome size is 3971.92 Mb (Table 3). Furthermore, based on the K-mer map, a high peak (22) was observed at half the K-mer depth (44), which indicates that the ZB genome has high heterozygosity, and the heterozygosity rate is estimated to be 1.73%. In addition, a fat tail was observed in the K-mer analysis, and the repeat sequence rate was calculated to be 86.04% (Figure 2A). The heterozygosity and repeat sequence rates of the ZB genome are much higher than those reported from the genome survey data of other woody plants, such as Acer truncatum (1.06%; 48.80%) [8], Xanthoceras sorbifolium (0.89%; 62.00%) [17], and Betula platyphylla (1.22%; 62.20%) [22]. Because of the high values for heterozygosity and repeat sequence rates combined with the chromosome number of ZB (2n = 136) [23], we speculate that ZB has a very complex genome.

Preliminary genome assembly was performed using clean reads. The software SOAPdenovo generated 11185221 contigs and 11086925 scaffolds using a K-mer length of 41. The maximum length and N50 length of contigs are 13151 and 183 bp, respectively, and for scaffolds, these values are 13628 and 186 bp, respectively (Table 4). Although a large amount of clean data (197.16 Gb) were used for assembly, the assembly results were unsatisfactory. The N50 lengths of contigs and scaffolds are notably shorter than those calculated in other similar studies [8,17,22]. A likely reason for these findings is that the ZB genome contains 68 chromosomes (2n = 136), has a high heterozygosity rate, and has a large number of repeated sequences; furthermore, the insert sizes for the sequencing libraries are relatively short (250 and 350 bp). Collectively, these factors likely contributed to the unsatisfactory assembly results.

GC content analysis was performed with contigs longer than 500 bp. Figure 2B shows the relationship between GC content and sequencing depth. The data suggest that the GC content of the ZB genome is approximately 37.21%, which is higher than the GC content in Citrus plants (32.0–35.0%) within the same family (Rutaceae) as ZB [24,25]. A scatter plot of GC content shows that the data segregate into two layers, a result that is likely due to the high heterozygosity rate (1.73%) [11]. GC content can influence the quality of genome sequencing. Because different densities of GC content can reduce the sequencing coverage in certain genomic regions, sequencing bias can occur on the Illumina sequencing platform and affect the genome assembly [26,27]. However, when GC content is between 30% and 50%, there is no significant influence on genome sequence quality [28]. Consequently, the 37.21% GC content of the ZB genome is not likely to have influenced the assembly results in the present study.

Based on the complex characteristics of the ZB genome, we recommend using second-generation sequencing (Illumina) combined with third-generation sequencing (PacBio) and using the Hi-C technique and BioNano Genomics for supplement in future whole-genome sequencing studies.

Assembled scaffolds (57288) with a minimum length and total length of 1000 bp and 86.51 Mb, respectively, were selected for SSR searching by MISA software, and 27153 putative SSRs were identified on 17593 scaffolds. The mean distance (3.19 kb) between G-SSR loci in the ZB genome was longer than mean distances measured in Ziziphus jujuba (0.93 kb) [29], Dioscorea zingiberensis (1.94 kb) [30], and Saccharina japonica (2.2 kb) [31]. One possible reason is that the scaffolds we analyzed only account for 2.18% of the total size of the ZB genome.

Mononucleotides are the most abundant type of SSR and account for 72.57% (19706) of the total. Dinucleotides are the second-most abundant type of SSR (18.98%, 5154), followed by trinucleotides (6.85%, 1860), and tetra-nucleotides (1.14%, 309). Pentanucleotides and hexanucleotides SSRs are found at 60 and 64 loci, which account for 0.22% and 0.24% of the total, respectively. As the repeat motif length increases, the number of SSR loci decreases. Among mononucleotides, A/T motifs are the predominant type (98.41%). Among dinucleotides, the most frequent motifs are AT/AT (62.36%), followed by AG/CT (19.95%) and AC/GT (17.56%); however, CG/CG accounts for just 0.13%. Among trinucleotides, the most frequent motif is AAT/ATT (56.53%), followed by AAG/CTT (20.82%), and ATC/GAT (7.92%). ACG/CGT motifs are the least frequent, with only 16 loci (0.89%). Moreover, the most abundant motifs among tetra-, penta- and hexanucleotides are AAAT/ATTT (51.46%), AAAAT/ATTTT (45.00%), and AAAAAT/ATTTTT (15.63), respectively. According to these results, motifs containing only A and T residues are more common than those containing at least one C or G base, especially among di- and trinucleotides (Table 5).

There is a certain relationship between the diversity of SSRs and motif sequence types. Because the energy required to destabilize the two hydrogen bonds between A and T is lower than that required to destabilize the three hydrogen bonds between G and C, the slippage rate of A/T is higher than G/C between the two DNA strands. The elevated slippage rate causes A/T to be more frequently observed in SSR motifs [32]. Other possibilities include the insertion of 3′-terminal poly(A) sequences into the genome or the conversion of methylated C residues to T residues [33].

To identify SSR loci that could have utility as potential markers, 21243 non-redundant primer pairs for 23475 G-SSR loci were designed using Primer 3 software (Supplementary Table S1); the remaining loci may have had flanking sequences that were too short or inappropriate for primer design.

To assess the degree of polymorphism for potential SSR markers, 100 primer pairs were selected randomly and validated across 15 individuals (ten ZB and five ZA individuals). Mononucleotide SSRs were not considered. PCR amplification showed that 85 primer pairs produced fragments that were clear and stable. Of these, 36 SSR loci were polymorphic after amplification products were separated (Table 6 and Figure 3). The polymorphism rate (polymorphic markers/number of markers used for polymorphic screening; 36/85) of the G-SSR markers tested is higher than that of EST-SSR markers developed in previous studies (18/55 and 15/44) [34,35]. One reason may be that exon sequences are more conserved than intron or intergenic sequences [36]; this suggests that in conserved regions, the relatively low frequency of polymorphisms may limit the utility of EST-SSR markers; consequently, the development of G-SSR markers is necessary [37,38].

In total, 126 alleles with a range of 2 to 7 per loci (mean = 3.5) were obtained from the 36 polymorphic SSR loci. The PIC, Ho, and He values per locus ranges from 0.16 to 0.75 (mean = 0.48), 0.00 to 0.86 (mean = 0.28), and 0.19 to 0.81 (mean = 0.56), respectively. According to the classification criteria of Bostein et al. [39], loci polymorphisms can be divided into three degrees: low (PIC < 0.25), moderate (0.25 < PIC < 0.5), and high (PIC > 0.5). Among the 36 polymorphic G-SSR markers, 1 (2.78%), 22 (61.11%), and 13 (36.11%) were demonstrated to have low, moderate, and high polymorphism, respectively, in 15 individuals. Because these SSR markers were developed based on a single ZB genome sequence, three markers (ZBg17, ZBg35 and ZBg45) did not amplify any fragment, and nine markers (ZBg02, ZBg03, ZBg11, ZBg31, ZBg34, ZBg36, ZBg40, ZBg74 and ZBg76) amplified only one fragment (i.e., were not polymorphic) in the five ZA individuals. Interestingly, four markers (ZBg46, ZBg83, ZBg86, and ZBg97) that are not polymorphic in ZB are polymorphic in ZA, and three markers (ZBg04, ZBg07 and ZBg98) produce more alleles in ZA than ZB, suggesting that some loci might be more likely to mutate in ZA relative to ZB. Similar results have been reported for Vernicia fordii [40], Taxus wallichiana [41], and Saxifraga sinomontana [42].

Based on the 36 polymorphic SSR markers identified, the genetic relationships among the 10 ZB individuals and 5 ZA individuals were investigated using UPGMA clustering. The dendrogram shows that the genetic similarity coefficient (GSC) ranges from 0.55 to 0.93 and the 15 individuals are distributed into two major clusters by species (Cluster I: ZA; Cluster II: ZB) (Figure 4). When the genetic similarity coefficient (GSC) value is approximately 0.72, Cluster I and Cluster II each divide into two subclusters: I-A, I-B, II-A, and II-B. Cluster I-A includes four individuals (ZA01-ZA04) that are from adjacent provinces (Sichuan and Chongqing), and Cluster I-B includes ZA05 (Goujiao). In Cluster II-A, the three individuals (ZB01-ZB03) are from adjacent provinces (Sichuan and Guizhou). However, Cluster II-B consists of seven individuals (ZB04-ZB10) that are from three provinces (Shaanxi, Shanxi, and Shandong). Hancheng of Shaanxi is one of the main producing areas of ZB. Therefore, we speculate that ZB05 and ZB07 were probably introduced from Hancheng. In addition, although 36 G-SSR markers were used, ZA08 and ZA10, which came from the same area with different name, are indistinguishable from each other, suggesting that they came from the same individual. The clusters support the expected classification of ZA and ZB and demonstrate the efficacy of the G-SSR markers developed in the present study.

In the present study, genomic characteristics of ZB were obtained by a genome survey, and G-SSRs were identified and developed simultaneously from the sequence data. The results showed that ZB has a notably complex genome. Its genome size is 3971.92 Mb, with a heterozygosity rate, repeat sequence rate, and GC content of 1.73%, 86.04%, and 37.21%, respectively. For future whole-genome sequencing, we propose that using Illumina and PacBio sequencing technologies combined with Hi-C and BioNano for assist will yield better genome assembly results. A total of 27153 G-SSRs were identified, and 21243 non-redundant primers pairs were designed. Thirty-six of one hundred randomly selected primer pairs showed polymorphism among ten ZB individuals and five ZA individuals. An UPGMA-derived dendrogram showed that the clustering of these 15 individuals was consistent with their species of origin. These findings will be useful for future genomic and genetic studies in ZB.

The authors declare that there are no competing interests associated with the manuscript.

This research was supported by the Doctor Faculty Inaugurating Project of Northwest A&F University [grant number 2452015296]; National Key Research and the Development Program Project Funding [grant number 2018YFD1000605]; and Innovative experiment program for College Students of Northwest A&F University [grant number 201710712023].

S.Q.L., J.M.L., L.J.K., L.H.W., A.Z.W., and Y.L.L. performed the experiments. J.M.L., S.Q.L., and Y.L.L. analyzed the data, prepared figures and tables, wrote the paper and reviewed drafts of the paper. All authors read and approved the final manuscript.

EST

expressed sequence tag

GSC

genetic similarity coefficient

He

expected heterozygosity

Hi-C

high-through chromosome conformation capture

Ho

observed heterozygosity

K-mer

a sequence of k characters in a string

Na

number of alleles

NGS

next-generation sequencing

PIC

polymorphism information content

SSR

simple sequence repeat

STR

short tandem repeats

UPGMA

unweighted pair group method with arithmatic mean

ZA

Zanthoxylum armatum

ZB

Zanthoxylum bungeanum