The genus Pistacia L., belonging to the family Anacardiaceae, is consisting of more than 13 species (Rauf et al. 2017). P. vera is a well-known species from the genus Pistacia and has been studied in various aspects due to its agronomical importance. P. chinensis is distributed in China, Taiwan, and Pakistan as well as North America. This small tree has been used for landscape and shade tree and its seed oil possesses not only biodiesel properties, but also pest repelling properties (Rashed et al. 2016). Tender burgeon of P. chinensis has been used as vegetables, and seeds can be used to make confectionery or vegetable oil. Also, P. chinensis has been used as a rootstock for P. vera because it has a strong adaptability and resistance to adverse environment (Tang et al. 2012).

Several Pistacia species, including P. chinensis, have been used as folk medicines to treat various illnesses (Akhtar et al. 2013; Bozorgi et al. 2013; Rauf et al. 2017). P. vera has been used to treat abdominal ailments and rheumatism (Bozorgi et al. 2013; Rauf et al. 2017). P. khinjuk and P. integerrima have been used to treat hepatitis and liver disorder (Rauf et al. 2017). P. lentiscus has been used to relieve various illnesses such as coughs, sore throats, eczema, and kidney stones (Rauf et al. 2017). Similarly, all parts of P. chinensis can be used for medicinal purpose to relieve dysentery, inflammatory swelling, psoriasis and rheumatism (Tang et al. 2012) and also to treat jaundice and liver diseases (Akhtar et al. 2013). Other Pistacia species such as P. terebinthus, P. palaestina, P. eurycarpa, P. weinmannifolia, and P. atlantica have also been used to relieve illnesses (Rauf et al. 2017).

Various studies have shown that most Pistacia species contain diverse and valuable secondary metabolites (Bozorgi et al. 2013; Rauf et al. 2017). Like other Pistacia species, various metabolites were identified from P. chinensis. Identified metabolites include new N-phenyl-pyrrolidone derivative (Liu et al. 2008), gallic acid and 6-0-galloyl arbutin-quercitrin (Shi and Zuo 1992), β-sitosterol, luepol, myricetin 3-O-α-rhamnoside (Rashed et al. 2016), 4-aryldihydrocoumarins (Nishimura et al. 2000). These chemicals have been used or have potentials for pharmaceutical drugs. Even though several Pistacia species are known to have valuable pharmaceutical and industrial properties, genomic and genetic studies of Pistacia species have not been vigorously studied except for P. vera. Therefore, it is essential to enrich genomic and genetic resources to delve into genetic insight of the genus Pistacia and also to identify not only key genes but also biosynthetic pathways involved in useful metabolite biosynthesis.

Advanced RNA sequencing transcriptome analysis tool is in preference to get transcriptome information of the target organism because it is cost- and time-effective (Lee et al. 2015). Because RNA sequencing method generates transcriptome information in a short time, this tool is frequently utilized to analyze transcriptome information of non-model organisms including medicinal plants (Bae et al. 2018; Eum et al. 2019). An increasing number of previously unexplored medicinal plants have been sequenced through this advanced sequencing technology, providing genomic resources for unravelling genes and biosynthetic pathways involved in metabolite biosynthesis in various medicinal plants (Bae et al. 2018; Eum et al. 2019; Kotwal et al. 2016; Loke et al. 2016; Rai et al. 2016). However, a vast majority of medicinal plants are yet to be studied.

In this study, RNA-seq transcriptome analysis was performed to characterize genomic features of P. chinensis. Generated RNA-seq data were used de novo assembly, and resulting unigenes were used for gene ontology (GO) analysis, KEGG metabolic pathway search, and SSR mining. In addition, the reference transcriptome of P. chinensis would be very useful resources for enriching and facilitating genetic/genomic studies, molecular marker discovery, and various genetic/biological studies of not only P. chinensis but also the genus Pistacia and close relatives.

Plant materials

Fresh leaf tissues of a fully grown P. chinensis Bunge were harvested June 2015 at Wunnan, China and submerged into liquid N₂. For long-distance transportation, the leaf samples were transferred into RNAlater solution (Ambion Ins, USA), which was stored in -20°C freezer prior to mRNA extraction for RNA-seq analysis.

Preparation of RNA-seq library

Leaf samples of P. chinensis were used for total RNA extraction using TRIzol reagent by following the manufacturer’s instructions. The remaining procedures were performed according to the methods described by Bae et al. (2018) and Eum et al. (2019).

De novo assembly and unigene annotation, and

For assembly, Trimmomatic tool (Bolger et al. 2014) was used to remove low quality reads (< Q20) and the read with a length < 50 bp. The following methods were performed according to the methods described by Bae et al. (2018). Briefly, three different assemblers were used for de novo assembly of the trimmed raw reads. Identified unigenes were annotated by comparing similarity to known proteins deposited to the NCBI non-redundant (NR) protein database.

GO analysis and KEGG pathway search

The Blast2GO analysis tool was used for GO annotation of unigenes, of which annotation information was used for functional classification of unigenes using WEGO software. For KEGG pathway search (http://www.genome.jp/kegg), all unigene sequences were run on the annotation by BlastKOALA (KEGG Orthology And Links Annotation), from which assigned K numbers were used to construct KEGG pathway search.

Identification of SSR and repetitive sequence in P. chinensis unigenes

To examine SSR accumulation in the unigenes from P. chinensis, SSRs were identified using MISA, MIcroSAtellite identification tool (pgrc.ipk-gatersleben.de/misa). Criteria of MISA were set to a minimum of four motif repeats and ≥ 12 bp. To examine a content of repetitive sequence in the unigenes, RepeatMasker (v. 4.0.7) was run on the unigenes with a default mode using the reference library, RepBaseRepeatMaskerEdition-20170127 (www.girinst.org). To compare SSR distributions of P. chinensis to other plant species, transcriptome data from several plant species (Mangifera indica, Sterculia lanceolata, Clausena excavata, Arabidopsis thaliana, and Oryza sativa) were used. Transcripts of Mangifera indica (Mango) were retrieved from NCBI nucleotide database. Transcript sequences of both Sterculia lanceolata and Clausena excavata were obtained from NCBI transcriptome shotgun assembly data. Transcripts of Oryza sativa (version 7.0) and Arabidopsis thaliana (ATH_cDNA_sequences_20101108) were retrieved from the Rice Genome database (http://rice.plantbiology.msu.edu) and from the TAIR database (https://www.arabidopsis.org), respectively.

De novo assembly of RNA-seq data and transcriptome annotation

RNA-seq whole transcriptome sequencing generated 20.8 million of raw reads (~2.6 Gb) from P. chinensis. De novo assembly using trimmed raw reads by Trimmomatic trimmer (Bolger et al. 2014) resulted in 19 million clear reads with a total of read length of ~2.36 Gb. From assembly using clear reads, a total of 18,524 unigenes with a length of 16,174,683 bp were generated (Table 1), of which GC content was 40.7%. N50 was 1,104 bp, and average length of unigenes was 873 bp (Table 1). The length of the unigenes ranged from 300 to 9,942 bp. A total of 6,353 (34.3%) unigenes was found between 297~500 bp, followed by 3,456 (18.7%) between 401~500 bp, and 2,348 (12.7%) between 701 ~ 900 bp (Fig. 1). A total of 515 unigenes (2.8%) were longer than 2,501 bp (Fig. 1).

Table 1 . Summary of sequencing and assembly data.

Data description	Data summary
Total number of raw reads	10,621,059
Total length of raw reads (bp)	2,676,506,868
Number of filtered reads used for assembly	9,625,676
Total length of filtered reads (bp)	2,358,007,393
Number of assembed contigs (Unigenes)	18,524
Total length of assembed contigs (bp)	16,174,683
Average length (bp)	873
Length of largest contig (bp)	9,942
N50 (bp)	1,104
GC content (%)	40.7

Figure 1.

Length distribution of unigenes from transcriptome of Pistacia chinensis

To annotate the unigenes, protein sequences of the unigenes were searched for similarity against NCBI non-redundant (NR) protein database. Among 18,524 unigenes, 17,814 unigenes (96.2%) were aligned to protein sequences from other organisms, whereas 710 unigenes (3.8%) did not show similarity to other known proteins (Fig. 2). Top-five plant species with most hits with annotated unigenes were Citrus sinensis with 8,977 unigenes (48.5%), Theobroma cacao with 2,655 (14.3%), Vitis vinifera 930 (5.0%), Populus euphratica 879 (4.7%), and Ricinus communis 677 (3.7%), respectively. As shown, a half of annotated unigenes showed higher similarity to Citrus sinensis (Fig. 2).

Figure 2.

Top five plant species with higher homologous genes

Accumulation of repetitive sequences and simple sequence repeats

To examine the content of repetitive sequences, RepeatMasker (http://www.repeatmasker.org) was run on unigenes, resulting in that 149,205 bp (0.92%) were occupied by repetitive and low complexity sequences (Supplementary Table S1). Simple sequence repeats were the most abundant repetitive sequences among identified repetitive sequences: 2,670 elements occupying 116,432 bp (0.72%) of the P. chinensis transcriptome (Supplementary Table S1). A total of 18 interspersed repeat elements were identified from P. chinensis unigenes.

SSR search using MISA SSR search tool identified a total of 2,629 perfect SSRs from 2,041 unigenes (Supplementary Table S2). Among SSR-containing unigenes, 393 had more than one SSR. The frequency of all identified SSRs was 162.5 per one million base pairs (Mbp). Tri-nucleotide SSRs were the most abundant SSRs with 2,343 (89.1%) occurrences, followed by di-nucleotide SSRs with 142 (5.4%) occurrences. Frequency of tri-nucleotide SSRs were 144.9 per Mbp (Fig. 3). The highest SSR frequency by motif type was AAG/CTT motif with 43.5 occurrences per Mbp, followed by ACC/GGT motif with 22.8 per Mbp (Fig. 4; Supplementary Table S2). Among di-nucleotide repeats, AG/CT motif showed the highest 7.6 occurrences per Mbp, while AAAG/CTTT motif showed the highest frequency with 0.6 per Mbp among tetra-nucleotide motifs (Fig. 4; Supplementary Table S2).

Figure 3.

SSR distribution and frequencies by repeat unit size

Figure 4.

SSR distribution by motif types

Functional classification of unigenes by GO analysis

To examine the functional classification of annotated genes, unigenes were assigned by the GO functional term using Blast2GO software (Conesa et al. 2005) and classified using WEGO tool. By functional categorization of annotated unigenes, 9,020 were classified into molecular functions at level one, 6,133 into biological processes, and 2,622 unigenes classified into cellular component, respectively. Most unigenes belonging to molecular function at level one were classified into two major categories, binding with 5,999 and catalytic activity with 4,493 unigenes (Fig. 5). As for unigenes in biological process, most genes were classified into three sub-categories, metabolic process with 4,728 genes, cellular process with 4,139, and single-organism process with 2,706 (Fig. 5). In cellular component, a majority of genes were classified into four sub-categories, cell with 1,659 genes, membrane with 1,225, organelle with 1,109, and macromolecular complex with 810 (Fig. 5; Supplementary Table S3).

Figure 5.

Gene Ontology functional categorization of annotated unigenes

Analysis of KEGG metabolic pathway of unigenes

To identify unigenes involved in the KEGG metabolic pathway, the KEGG BlastKOALA online tool was used to assign KEGG Orthology (KO) number, by which a total of 6,553 unigenes of P. chinensis were assigned. Among them, 4,061 unigenes were mapped into 391 different KEGG metabolic pathways (Table 2), in which 2,561 unigenes were shown to be involved in more than one pathway. Major portions of unigenes (6,715) and pathways (138) were found to be included into the pathways of metabolism category (Table 2).

Table 2 . Categories of KEGG metabolic pathways and their associated entries and unigenes.

Category	No. of sub categories	No. of pathways	No. of entry	No. of associated genes
Metabolism	12	138	3324	6715
Genetic Information Processing	4	22	859	1520
Environmental Information Processing	3	35	424	1242
Cellular Processes	5	31	549	1280
Organismal Systems	10	84	607	1608
Human Diseases	12	81	935	2131

Among all pathways, top four pathways with most entry enzyme hits were metabolic pathways with 821 entry enzymes, biosynthesis of secondary metabolites with 381, biosynthesis of antibiotics with 194, microbial metabolism in diverse environments with 144 entries, respectively (Fig. 6).

Figure 6.

Top 10 KEGG pathways with the most entry enzymes identified from P. chinensis transcriptome. Number for each pathway denotes total entry enzymes identified and the number in parenthesis indicates total unigenes encoding entry enzymes

Biosynthesis of many secondary metabolites is tightly correlated to some of metabolic pathways including the pathways involved in the metabolism of terpenoid, the biosynthesis of phenylpropanoid, or flavonoid. Twenty-two unigenes were found for 16 entries of carotenoid biosynthesis pathway (Table 3). For phenylpropanoid biosynthesis pathway, 63 unigenes were found to encode enzymes for 17 entries, while 24 unigenes were found for 14 entries of the flavonoid biosynthetic pathway (Table 3; Supplementary Table S4).

Table 3 . KEGG metabolic pathways related to the biosynthesis of various medicinal metabolites.

Metabolic pathways	KEGG map ID	Number of entry	No of unigenes*
Metabolism of terpenoids
Carotenoid biosynthesis	00906	16	22
Sesquiterpenoid and triterpenoid biosynthesis	00909	5	8
Diterpenoid biosynthesis	00904	5	6
Monoterpenoid biosynthesis	00902	2	5
Biosynthesis of phenylpropanoid and flavonoid
Phenylpropanoid biosynthesis	00940	17	63
Flavonoid biosynthesis	00941	14	24
Isoflavonoid biosynthesis	00943	1	2
Flavone and flavonol biosynthesis	00944	1	1

^*Number of unigenes indicates that they can encode enzymes for corresponding entry.

Medicinal plants are normally rich in traditional knowledge about medicinal usage, but there is very limited genetic information available for most traditional medicinal plants except for well-known medicinal plants. Medicinal plants are getting more interest to identify new metabolic compounds that possess important medicinal properties. In the absence of genomic information, however, it is very difficult to identify new lead molecules for pharmaceutical drug development from medicinal plants and to delve into how those molecules are synthesized in those plants. Therefore, enrichment of genomic resources as well as genetic information is crucial for studying medicinal properties and for identifying potential lead molecules from unexplored medicinal plant species. As medicinal plants are getting more interest, increasing number of those species are getting sequenced by advanced sequencing technologies.

Advanced RNA-seq technology has spurred transcriptome analysis of increasing number of medicinal plants that have not been of interest previously. In this study, transcriptome P. chinensis was analyzed by RNA-seq technology. De novo assembly generated 18,524 unigenes (Table 1), of which 17,814 unigenes were annotated. The number of identified unigenes in the present study is significantly lower than that from previous study in which 127,545 unigenes were identified from P. chinensis (Dong et al. 2016). However, compared to reports from other plant species, the number of unigenes reported by Dong et al. (2016) is very high, so further verification may be needed. On the other hand, the number of unigenes in our study is about 60% of the number of predicted genes from P. vera genome, in which 31,784 genes were reported by Zeng et al. (2019).

Many unigenes from the present study showed higher similarity to those of Citrus sinensis. Roughly, half of the annotated unigenes, 8,977 unigenes (48.5%), showed similarity to proteins from Citrus sinensis (Fig. 2), suggesting that P. chinensis is closely related to C. sinensis. Genome evolution analysis using gene families from P. vera showed a similar results among nine species from different genus, in which P. chinensis and C. sinensis showed that they diverged more recently among tested nine species (Zeng et al. 2019).

Many SSR markers have been developed mostly from P. vera (Zaloglu et al. 2015; Ziya Motalebipour et al. 2016), many of which were successful for cross-amplification among various Pistacia species (Zaloglu et al. 2015; Ziya Motalebipour et al. 2016). In this study, a total of 2,629 perfect SSRs were identified from P. chinensis, of which frequency was 162.5 per one Mbp (Supplementary Table S2). Frequency of tri-nucleotide SSRs was 144.9 per one Mbp (Fig. 3) and was the most abundant SSRs as reported in various other studies (Zhang et al. 2019; Bae et al. 2018; Eum et al. 2019; Kotwal et al. 2016). Repetitive sequence analysis showed that repetitive and low complexity sequences occupied 149,205 bp (0.92%) of the P. chinensis transcriptome, and the most abundant repetitive sequences were SSRs (Supplementary Table S1).

Various secondary metabolites possess important medicinal effect, and many metabolites are known to be derived from pathways involved in terpenoid metabolism or secondary metabolite biosynthesis pathways. A total of 4,061 unigenes were assigned into 391 different metabolic pathways through KEGG pathway analysis. Among the assigned unigenes, 131 genes were found to be involved in some of metabolite biosynthesis pathways including terpenoid, phenylpropanoid, and flavonoid (Table 3; Supplementary Table S4). Various genes involved in these pathways were identified from the present study as shown in Supplementary (Supplementary Table S4). Reconstruction of flavonoid and phenylpropanoid biosynthesis using the unigene information in Supplementary Table S4 revealed that most key genes involved in those pathways are well conserved in P. chinensis (Supplementary Fig. S1 and S2). Such genes include those that encode enzymes required for biosynthesis of kaempferol and quercetin known to have anti-cancer activity (Batra and Sharma 2013). Also, the biosynthetic pathway of luteoforol (Deoxyleucocyanidin) is well conserved, which is known to have an antimicrobial activity (Spinelli et al. 2005). Nevertheless, since P. chinensis contains various medicinal compounds having potentially medicinal effects (Liu et al. 2008; Nishimura et al. 2000; Rashed et al. 2016; Shi and Zuo 1992), our data will be invaluable to study how those metabolic compounds are synthesized.

The authors declare that there is no conflict of interests.

This study was supported by KRIBB initiative program of Republic of Korea and 2017 research grant to JKN from Kangwon National University.

This transcriptome shotgun sequencing data was deposited at BioProject: PRJNA566127 in NCBI GenBank. SRA accession number is SRR10136265.