De novo assembly and characterization of the transcriptome and development of microsatellite markers in a Chinese endemic Euphorbia kansui

Abstract

Euphorbia kansui Liou is an endemic and traditional medicinal plant in China. Despite its value, the genomic and transcriptomic information of E. kansui is limited. In the present study, a total of 43,211,690 high-quality reads were generated using Illumina paired-end sequencing technology. After assembly, 58,362 unigenes were recovered with an N50 length of 1,683 bp. 36,396 (62.36%) of these unigenes had at least one blast hit against the multiple public databases. Based on KEGG pathway analyses, unigenes were discovered to encode known enzymes that catalyze the formation of the terpenoid backone and diterpenoid. In addition, the unigenes related to casbene synthase, which catalyzes the formation of the major active constituents (diterpenoids) in E. kansui, were identified. For functional marker development, 7,016 candidate simple sequence repeats (SSR, microsatellites) were identified from 6,150 unigenes, and 40 SSR loci were selected randomly and amplified in two populations of E. kansui. A total of 28 loci were amplified successfully and 23 loci exhibited polymorphisms. The number of alleles of each locus ranged from 2 to 8 alleles and the mean number of alleles per locus was 3.391. The observed heterozygosity and expected heterozygosity varied from 0.100 to 1.000 and 0.099 to 0.809, respectively. The microsatellite markers developed here will be useful for studies on population genetics of this endemic species. Moreover, the information of the transcriptome will provide a useful resource for the future research on metabolism, regulation of the terpenoid biosynthesis, genetic and proteome studies of E. kansui.

Keywords:

• Euphorbia kansui
• transcriptome
• terpenoid biosynthesis gene
• microsatellite markers

Introduction

Euphorbia kansui Liou is a perennial herb of Euphorbia (Euphorbiaceae). It is an endemic species of China, and is distributed mainly in Shaanxi, Henan and Gansu provinces. The dried root of E. kansui is a commonly used and well-known traditional Chinese medicine and is listed in the Chinese Pharmacopoeia. It is widely applied to treat with edema, ascites, schistosomiasis and asthma [1]. It was reported that E. kansui showed curative effects on anti-tumor [2], anti-viral [3], inhibition cell division [4, 5], anti-fertility [6] and pesticidal activity [7]. At present, it is often used clinically to treat pancreatitis, intestinal obstruction, hydrothorax, ascitic fluid phlegm syndrome and cirrhotic ascites. In recent years, the genome sequences of some plants in Euphorbiaceae have been completed, such as Ricinus communis [8], Jatropha carcas [9] and Hevea brasiliensis [10]. However, genomic information about E. kansui is lacking, and the current lack of genomic sequence data limits the molecular research. In our previous study, we reported the development and ultrastructure of laticifers and the latex proteome of E. kansui, but identifying proteins involved in the development of laticifers has been challenging greatly because of lack of its genome and transcriptome information [11].

Transcriptome is the set of transcripts of an organism at a defined spatial and temporal locus, which provides valuable and comprehensive information on gene expression, gene regulation, and amino acid content of proteins [12]. In recent years, next-generation sequencing (NGS) technologies, such as Illumina paired-end transcriptome analysis, provide new opportunities to acquire transcriptomic data for the purpose of gene mining and molecular marker development, especially for species without a reference genome [13]. Wang et al. (2009) reported the transcriptome data of Artemisia annua, and characterized the terpenoid biosynthesis pathway genes, sesquiterpene synthase and monoterpene synthase [14]. The candidate genes associated with carotenoid biosynthesis in Lycium chinense were discovered and further confirmed based on transcriptome data [15].

Microsatellites, also called simple sequence repeats (SSR), are found in the nuclear genomes of all eukaryotes and are often abundant and evenly dispersed [16]. SSR markers are highly informative and widely used for evolution and genetics studies [17]. Transcriptome sequencing is capable of identifying massive numbers of potential SSR [18]. Illumina paired-end sequencing was successfully used for developing SSR markers of non-model organisms, including Mucuna prurien [19], Dendrobium officinale [20]. In E. kansui, only 12 genome SSR markers were developed [21], but no EST-SSR markers are available to date. EST-SSR markers predicted and developed in this study will enrich the number of molecular markers, and facilitate gene mapping, linkage map development and genetic diversity analysis of this endemic species.

In the present study, we performed large-scale transcriptome sequencing using Illumina paired-end sequencing technology to obtain comprehensive genomic information. SSR markers were developed based on the transcriptome sequences for genetic diversity studies of E. kansui. After annotation, we analyzed candidate genes involved in terpenoids biosynthesis. To our knowledge, this is the first systematic report about transcriptome analysis of E. kansui. This study will provide the basis for molecular biology, genetics and further proteomic research of E. kansui.

Materials and methods

Plant material

Euphorbia kansui plants were obtained from the field at the Botanical Garden of Northwest University (34.25 N, 108.92 E), Xi’an Botanical Garden (34.21 N, 108.95 E) and Yongshou county (34.70 N, 108.15 E) in Shaanxi Province (Shaanxi, People’s Republic of China). Tender leaves, inflorescence and fruits were collected and flash frozen in liquid nitrogen and stored at −80 °C until RNA extraction.

The leaves of 20 individuals (10 individuals of each population) were sampled from the Botanical Garden of Northwest University (34.25 N, 108.92 E) and Xi’an Botanical Garden (34.21 N, 108.95 E) respectively to develop SSR markers. These plants were transplanted from different wild populations.

RNA extraction

To globally characterize the E. kansui transcriptome and enhance sequence coverage, the total RNA was isolated respectively from young leaves, inflorescence and fruits using RNA extraction Trizol Kit (Gibco BRL, Grand Island, USA) following the manufacturer’s instructions. Then the total RNA was treated with RNase-free DNase I (Takara Bio, Japan) for 30 min at 37 °C to remove residual DNA. The quantity and quality of the isolated RNA were assessed by Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). Equal amounts of RNA isolated from different tissues were mixed for the subsequent cDNA library construction step.

cDNA library construction and illumina sequencing

Transcriptome library for sequencing was constructed according to the Illumina TruSeq RNA library protocol outlined in “TruSeq RNA Sample Preparation Guide” (Illumina 2009). Briefly, poly(A) mRNA was purified using Sera-mag Magnetic Oligo (dT) Beads (part of the mRNA-seq Preparation kit, Illumina, USA). The purified mRNA was first chopped into short fragments by adding fragmentation buffer, and the cleaved RNA fragments were reverse transcribed into first-strand cDNA with reverse transcriptase and random hexamer primers. Then second-strand cDNA was synthesized by using DNA polymerase I, dNTPs and RNase H. The short double-stranded cDNA fragments were purified using a QiaQuick PCR extraction kit (Qiagen, Hilden, Germany), and these purified fragments was then washed with EB buffer for end repair and A tailing. Then, Illumina Adapters were ligated to sequencing adapters using T4 DNA ligase. Following agarose gel electrophoresis and extraction of cDNA from gels, the cDNA fragments were purified and enriched by PCR to construct the final cDNA library. Finally, the cDNA library was sequenced on the Illumina sequencing platform (Illumina HiSeq TM 2000 system) using the paired-end technology in a single run.

De novo assembly and data analysis

The raw data from the images acquired after sequencing was transformed by base-calling into raw reads and stored in FASTQ format. Clean reads were selected by removing reads containing adaptor sequences, reads with more than 5% N base (base unknown) low quality sequences. And the clean reads were assembled using Trinity and TGICL to construct unique consensus sequences with the default setting.

Gene annotation

All obtained unigenes were annotated with the basic local alignment search tool (BLASTx) analysis with an E-value cut-off of 10⁻⁵ against the databases including non-redundant protein (Nr) database, Swiss-Prot protein database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database, and the Cluster of Orthologous Groups (COG) database to retrieve proteins with the highest sequence similarity score with the given unigenes along with their protein functional annotations. COG database was used to analyze, predict and classify unigenes with probable functions, whereas KEGG database was used to assign pathways of unigenes. The Blast2 GO program was used to obtain the annotation result of all the unigenes in Gene Ontology (GO) database which can describe the information of genes through three ontologies: molecular function, biological process and cellular component. After that, statistical analysis of GO functional classification was performed by WEGO software.

Identification of SSR loci, primer design and marker validation

SSR markers were detected using MIcroSAtellite tool (http://pgrc.ipk-gatersleben.de/misa/). SSRs were considered to contain motifs with two to six nucleotides in size and a minimum of 4 contiguous repeat units. Then primer pairs were designed by primer 3.0 (http://www.broadinstitute.org/genome_software/other/primer3.html) in the flanking regions of SSRs. Primers were designed based on the following criteria: GC content ranging from 40% to 60%, primer length ranging from 18 to 27 bp, melting temperature ranging from 55 to 63 °C and the expected PCR product sizes ranging from 100 to 280 bp.

A total of 40 primer pairs were randomly selected and synthesized by Sangon biotechnology Co., Ltd (Shanghai, China). For validation of SSR primers, total genomic DNA was extracted from young leaves of 20 individuals from two populations described previously (ten for each one) using plant genomic DNA kit (Tiangen, Beijing, China). PCR amplifications were performed in a total volume of 10 μL containing 1 ng/μL DNA, 5 μL 2 × Taq PCR master mix (0.1 U Taq polymerase/μL; 500 µmol/L each dNTP; 20 mmol/L Tris–HCl, pH8.3; 100 mmol/L KCl; 3.0 mmol/L MgCl₂) and 0.2 μmol/L of each forward and reverse primer and 3.6 μL ddH₂O. The PCR amplification was performed under the following conditions: initial denaturation at 94 °C for 5 min followed by 30 cycles of 94 °C for 50 s, 30 s annealing at 55–64 °C, 45 s extension at 72 °C, and 5 min at 72 °C for final extension. The PCR products were separated by 10% polyacrylamide gel electrophoresis; then the gels were stained with silver nitrate. Fragment sizes of each locus were estimated using Quantity One Software (Bio-Rad Laboratories, Drive Hercules, CA, USA) with pBR322 DNA/MspI as the molecular weight standard (Tiangen, Beijing, China).

Microsatellites data analysis

The number of alleles (N_A), expected heterozygosity (H_E), observed heterozygosity (H_O), polymorphic information content (PIC) and Hardy-Weinberg (HW) equilibrium were calculated using CERVUS version 3.0.

Results and discussion

Sequence assembly and data analysis

The next-generation sequencing technology has provided a powerful and cost-efficient tool to obtain large amounts of transcriptome data from many organisms and tissue types without a reference genome [22, 23]. In the present study, RNA was extracted from fresh leaves, inflorescence and fruits. The extracted RNA from these tissues was mixed equally and sequenced with Illumina paired-end sequencing technology, so as to find more unigenes of E. kansui. A similar method has been used in transcriptome sequencing of Benincasa hispida [18] and Camelia sativa [24]. After cleaning the adaptor sequences and quality checks, about 43,211,690 clean reads were obtained amounting to 4,321,169,000 base pairs (bp), and 97.80% of Q20 bases, and these reads had an average GC content of 44.09% (Table 1). The raw data were deposited in the NCBI Sequence Read Archive (SRA) database (accession number: SRP067381).

Table 1. Summary of transcriptome sequencing dataset.

	Number	Total nucleotide bases (bp)	Average length (bp)	N50 (bp)
Clean reads	43,211,690	4,321,169,000	100	–
Contigs	61,712	56,176,504	910	1,530
Unigenes	58,362	51,550,633	1,683	1,683

After length and quality-based filtering, a total of 61,712 contigs were generated, containing 5,61,76,504 bp total nucleotides, with an average length of 910 bp, and a N50 of 1,530 bp (Table 1). The contigs with the length more than 500 bp accounted for about 53.18% (32,817), and 2,012 (3.26%) contigs were longer than 3,000 bp (Figure 1; Supplemental Figure S1). Then, de novo assembly of these contigs generated 58,362 unigenes which did not contain gap regions. Both the N50 and average length of the unigenes were 1,683 bp (Table 1). The mean length was longer than the results shown in previous studies, such as Erigeron breviscapus (916 bp) [25] and Litchi chinesis (710 bp) [26]. N50 is a value to evaluate the assembly quality, and the N50 length was longer than those in other species [18, 22]. These results suggested that the transcriptome sequencing data were effectively assembled and it is possible to obtain reliable unigene sequences. Meanwhile, these data also indicated that Illumina paired-end sequencing technology is useful for the de novo sequencing and assembly of transcriptomes of non-model plants. The length of assembled unigenes ranged from 201 bp to 12,178 bp. The length distribution of the unigenes showed that 28,224 (48.36%) were 200–500 bp long, 30,138 (51.64%) unigenes were longer than 500 bp, and those with a length of over 3,000 bp were 1,786 (3.06%) (Figure 1; Supplemental Figure S2). A total of 23,805 unigenes (40.79%) were realigned with 11–100 reads, which comprised the largest group; 14,473 (24.80%) and 4949 (8.48%) unigenes were remapped by more than 100 and 1000 reads, respectively; whereas only 273 (0.47%) unigenes were remapped with more than 10,000 reads (Supplemental Figure S3).

Figure 1. Assembly length distribution of contigs and unigenes in Euphorbia kansui.

Figure 2. Number of unigenes annonated in protein databases by BLASTx with E-value threshold of 10⁻⁵. The numbers in the circles indicate the number of unigenes annotated by single or multiple databases.

Figure 3. Gene ontology classification of assembled unigenes of Euphorbia kansui. In total 15,506 unigenes were categorized into three main categories: biological process, cellular component and molecular function.

Functional annotation of unigenes

All the unigenes were searched using BLASTx against Nr, Swiss-Prot, COG and KEGG databases with an E-value threshold of 10⁻⁵. The results suggested that among the 58,362 unigenes, a total of 36,396 (62.36%) unigenes were annotated in the public databases. 36,318 (62.23%) unigenes and 26,640 (45.65%) unigenes were homologous with the Nr database and Swiss-Prot database respectively (Figure 2). Among them, 6195 unigenes were annotated in all four databases and 8771 unigenes were annotated uniquely in Nr database, 37 uniquely in Swiss-Prot database, 2 unigenes in COG database and 30 unigenes uniquely in KEGG database (Figure 2). This implied that the transcriptome sequencing yielded a great number of unique genes, while 21,966 (37.64%) unigenes were not annotated in part because the whole genome sequence data was not available for E. kansui. The unigenes that had no hits to any sequences in the public databases may represent novel genes [27]. We also analyzed the species distribution of BLAST hits for each unigene in the Nr database. The species distribution showed that 9990 (27.51%) of the mapped unigenes were annotated by Theobroma cacao homologs, 5252 (14.46%) from Vitis vinifera, and 2531 (6.97%) from the Fragaria vesca subsp. vesca (Supplemental Figure S4).

Figure 4. Cluster of Orthologous Groups (COG) function classification of Euphorbia kansui unigenes. A total of 13,528 unigenes were annotated and classified into 25 different classifications.

Gene ontology classification and COG classification

In this study, 15,506 (42.70%) unigenes were assigned with at least one GO term under the three main categories: biological process, cellular component and molecular function (Figure 2). The functionally assigned unigenes covered 43 GO categories. Under biological process, 8635 (55.69%) unigenes were assigned to metabolic process, followed by cellular process (8341, 53.79%) and response to stimulus (3484, 22.47%), whereas only a few unigenes were assigned to locomotion (8, 0.05%) and viral reproduction (13, 0.08%) (Figure 3). For the cellular component class, the cell (9037, 58.28%) and cell part (9037, 58.28%) were the dominant terms, whereas the least abundant term was extracellular region part (11, 0.07%) (Figure 3). For the molecular function category, the major terms were catalytic activity (8,159, 52.62%) and binding (7,069, 45.59%), whereas the least abundant terms were electron carrier activity (15, 0.10%) and transcription regulator activity (17, 0.11%) (Figure 3).

All unigenes obtained in our study were searched against the COG database for functional prediction and classification. Overall, 13,528 unigenes were annotated and assigned to 25 COG categories (Figures 2 and 4). The cluster of general function prediction represented the largest one (4,552, 33.65%), followed by replication, recombination and repair (2,388, 17.65%), transcription (2340, 17.30%), posttranslational modification, protein turnover, chaperones (1,853, 13.70%), signal transduction mechanisms (1,776, 13.13%), translation, ribosomal structure and biogenesis (1,625, 12.01%), carbohydrate transport and metabolism (1,516, 11.21%), function unknown (1,439, 10.64%) and cell wall/membrane/envelope biogenesis (1,119, 8.27%), whereas only a few unigenes were assigned to extracellular structures (9, 0.07%) and nuclear structure (6, 0.04%) (Figure 4). A large number of unigenes were assigned to 43 GO categories and 25 COG classifications, which indicated that our transcriptome data represented a broad diversity of transcripts in E. kansui genome. These annotations provide comprehensive information on the unigene functions.

Function classification by KEGG

All the unigenes were analyzed by the KEGG pathway database to further identify the active biochemical pathways. In total, 9,562 unigenes were mapped to 269 KEGG pathways, and some unigenes were assigned to more than one KEGG pathway. All the KEGG pathways are listed in Supplemental Table S1. Among the 9,562 unigenes, 420 (4.39%) were assigned to spliceosome, which was the largest group among the KEGG categories. The other large pathway groups were ribosome (365, 3.82%), protein processing in endoplasmic reticulum (289, 3.02%) and RNA transport (281, 2.94%). In particular, unigenes involved in terpenoid backone biosynthesis (68 unigenes), diterpenoid biosynthesis (22 unigenes), fatty acid metabolism (54 unigenes) and flavonoid biosynthesis (27 unigenes) were all identified in the database (Supplemental Table S1).

Identification of candidate genes involved in terpenoid biosynthesis

Diterpenoids are major constituents of E. kansui [28]. The key genes involved in terpenoid backbone (isopentenyl diphosphate) and diterpenoid biosynthesis pathways were discovered, such as the genes encoding thiolase, 3-hydroxy-3-methylglutaryl coenzyme A synthase, mevalonate disphosphate decarboxylase and so on (Figure 5, Table 2). 3-Hydroxy-3-methylglutaryl-CoA reductase (HMG-CoA reductase) is a key rate-limiting enzyme in the mevalonate (MVA) pathway, which catalyzes the conversion of HMG-CoA into mevalonate. Genes that encode HMG-CoA reductase have been cloned and characterized from E. pekinensis and Hevea brasiliensis [29, 30]. Mevalonate disphosphate decarboxylase was identified in latex of E. helioscopia and E. kansui. Whereafter, mevalonate disphosphate decarboxylase gene was cloned and it was specifically expressed in different tissues in E. helioscopia [31, 32]. In addition, three downstream enzymes (geranyl-diphosphate synthase, geranylgeranyl pyrophosphate synthase and farnesyl pyrophosphate synthase) involved in terpenoid biosynthesis were also discovered in the transcriptome data.

Figure 5. Terpenoid backbone and diterpene biosynthesis pathways. Terpenoid backbone (isopentenyl diphosphate) is synthesized from two distinct pathways, methylerythritol phosphate (MEP/DXP) pathway and mevalonate (MVA) pathway. The enzymes participating in each step of these pathways are listed in Table 2.

Table 2. Identified unigenes of terpenoid backbone and diterpene biosynthesis.

Steps of the pathway^a	Enzyme name	Number of unigenes
MVA (mevalonate) pathway
1	Thiolase	1
2	3-hydroxy-3-methylglutaryl coenzyme A synthase	7
3	3-hydroxy-3-methylglutaryl-CoA reductase	4
4	Mevalonate kinase	1
5	Phosphomevalonate kinase	4
6	Mevalonate disphosphate decarboxylase	4
7	Isopentenyl-diphosphate Delta-isomerase	3
MEP/DXP (methylerythritol phosphate) pathway
8	1-deoxy-D-xylulose 5-phosphate synthase	8
9	1-deoxy-D-xylulose 5-phosphate reductoisomerase	3
10	2-C-methyl-D-erythritol 4-phosphate cytidylyltransferase	1
11	4-(cytidine 5′-diphospho)-2-C-methyl-D-erythritol kinase	2
12	2-C-methyl-D-erythritol 2,4-cyclodiphosphate synthase	1
13	4-hydroxy-3-methylbut-2-en-1-yl diphosphate synthase	1
14	4-hydroxy-3-methylbut-2-enyl diphosphate reductase	1
15	Geranyl-diphosphate synthase	2
16	Geranylgeranyl pyrophosphate synthase	6
17	Casbene synthase	1
18	Farnesyl pyrophosphate synthase	5
19	Isoprene synthase	2

Note:

^a The steps of the terpenoid backbone and diterpene biosynthesis pathways which are indicated in Figure 5.

We also identified the candidate genes which encoded casbene synthase. In Euphorbiaceae, casbene synthase (diterpene synthase) is a key enzyme in the diterpenens biosynthesis pathway. It catalyzes the synthesis of macrocyclic diterpene olefins casbene [33]. Kirby et al. (2010) discovered genes encoding putative casbene synthases in Homalanthus nutans, Sapium sebiferum, E. resinifera, E. esula and R. communis, and they considered that casbene is likely the precursor to a large number of Euphorbiaceae diterpenes [34]. Casbene synthase has been isolated and characterized, and its cDNA was cloned from R. communis [35, 36] and roots of E. fischeriana [37]. The information about the unigene sequences in the transcriptome data provides a vital resource for cloning the cDNA of enzymes, and will facilitate the future research of terpenoid biosynthesis and regulation in E. kansui.

Est-SSR discovery and validation

Microsatellite markers have numerous uses, including linkage map development, quantitative trait loci mapping, marker-assisted selection, parentage analysis, cultivar fingerprinting, genetic diversity studies, gene flow and evolutionary studies [38]. Recently, transcriptome sequencing has become a viable method for microsatellite mining and SSR marker development, and has been utilized in many species [39]. Traditional strategies for SSR marker development are labour-intensive and costly. In E. kansui, 12 genome SSR markers were developed using paired-end Illumina shotgun sequencing [21]. However, no EST-SSR markers were reported. The lack of molecular markers has impeded the genetic analysis of E. kansui, an important medicinal plant. the In the present study, all the 58,362 unigenes generated were used to mine potential simple sequence repeats, which were defined as di- to hexa-nucleotide SSR with a minimum of four repeats for all motifs (except for di-nucleotide with a minimum of six repeats, and tri-nucleotide with a minimum of five repeats). A total of 7,016 SSR loci were detected in 6,150 unigenes, out of which, 747 unigenes contained more than one SSR (Table 3). The SSRs were highly distributed over the sequences, providing a useful resource for further research of these SSR loci as genetic markers. The most abundant motif of EST-SSR varied with species. Detailed analysis showed that tri-nucleotide repeat motifs were the most abundant types with a frequency of 58.87%, and this finding was consistent with previous reports of Camelina sativa (38.21%) [24] and Chrysanthemum indicum (52.32%) [40], followed by di-nucleotide (25.06%), tetra-nucleotide (9.78%), penta-nucleotide (4.06%) and hexa-nucleotide repeat motifs (2.24%) (Figure 6a). The frequencies of SSR distributed in different repeat numbers are shown in Figure 6b. SSRs with five tandem repeats were the most abundant (2384), followed by those with six tandem repeats (1900), seven tandem repeats (1072) and four tandem repeats (886). Motifs that showed more than 12 repeats (di-nucleotide, 9; tri-nucleotide, 1; hexa-nucleotide, 1) were very rare. The most abundant motifs were tri-nucleotide repeat units AAG/CTT, ATC/ATG and AAT/ATT with frequencies of 24.10%, 10.46% and 6.74%, respectively. Among the di-nucleotide repeat units, AG/CT (18.87%) and AT/AT (4.36%) were the most abundant motifs (Figure 7). This showed that AAG/CTT and AG/CT were the most common motifs. Similar findings have also been reported in Manihot esculenta [41] and Mentha piperita [42].

Figure 6. Distribution of simple sequence repeats (SSR) in Euphorbia kansui: distribution of different repeat types (a); number of SSRs per different tandem repeat motifs (b).

Figure 7. Frequency distribution of main SSR repeat motifs.

Table 3. Summary of SSR search results.

Total number of sequences examined	58,362
Total size of examined sequences (bp)	51,550,633
Total number of identified SSRs	7,016
Number of SSR containing sequences	6,150
Number of sequences containing more than one SSR	747
Number of SSRs present in compound formation^a	259

^a Maximal number of bases interrupting two SSRs in a compound microsatellite: 100 bp.

A total of 4,146 SSR loci were designed by primer 3.0, and detailed information of the SSRs discovered in this study is described in Supplemental Table S2. A total of 40 loci (di-, tri-, tetra-, penta- and hex-anucleotide repeats) were randomly selected for primer synthesis and validation among two populations with 20 E. kansui individuals. Among the 40 primer pairs, 28 (70.00%) successfully yielded amplification products, whereas the other 12 SSR loci had no amplified products at the expected size (Supplemental Table S3). Among them, 23 primers could be easily and consistently amplified and were polymorphic in two populations. The number of alleles (N_A) of each locus ranged from 2 to 8 alleles with an average of 3.391. The polymorphic information content (PIC) ranged from 0.094 to 0.758 with an average of 0.468. The observed heterozygosity (H_O) and expected heterozygosity (H_E) varied from 0.100 to 1.000 and 0.099 to 0.809, with an average of 0.735 and 0.550, respectively (Table 4). In 23 polymorphic EST-SSRs, E27704, E4746, E57686 and E58997 showed significant departures from the Hardy-Weinberg (HW) equilibrium (Table 4). Five loci were monomorphic (Supplemental Table S4), and those SSRs had no detectable polymorphism within two populations of E. kansui. We speculated that it is possibly because of the small sample size and population (two populations). Possibly, the SSRs exhibited polymorphism among other populations which were distributed in other regions. The other 12 loci did not amplify products at the expected size, perhaps because of the presence of introns in the unigenes, or because there may have been some errors in the base calls or assembly that resulted in incorrect flanking sequences from which primers were designed [43, 44].

Table 4. Characterization of 23 polymonophic SSR markers in Euphorbia kansui.

Locus	Blast match to Nr database	Primer sequence (5'–3')	Repeat motif(s)	Tm (°C)	Product size (bp)	NA	PIC	Ho	HE	HW
E27704	PREDICTED: ubiquitin-fold modifier 1-like [Solanum lycopersicum] 720 bp	F: TGCCAAGCCCGTATCTTTAC R: ACCGGTAGCAAATTACAGCG	(AG)11	64	240–266	3	0.555	0.500	0.641	***
E25447	aux/IAA protein [Populus tremula x Populus tremuloides] 1346 bp	F: GATTGAGCCTTGGTTGTGGT R: CACCATTAGGGGGAGAGACA	(CAA)9	61	238–265	8	0.587	0.900	0.663	ND
E58869	PREDICTED: protein NIM1-INTERACTING 1-like [Cicer arietinum] 588 bp	F: TCCATCTTTTCAACCCGAAG R: CCCAAAAAGTGCTTGGAAAG	(AGC)9	55	263-272	4	0.601	0.900	0.683	ND
E43796	arabinogalactan protein 13 [Arabidopsis thaliana] 689 bp	F: CAGGAATCTCCAGCTCCAAG R: GAAAAATCCCAACCAAGCAA	(GAT)9	61	260–272	4	0.607	0.900	0.691	ND
E5019	UDP-D-glucuronic acid 4-epimerase [Gossypium hirsutum]	F: GATGGTGCACTTGCTTTTCA R: ACCGGGAGTAGAAGGGAAGA	(GAAA)6	50	223–231	2	0.342	0.650	0.450	ND
E182	PREDICTED: ribosome biogenesis protein NSA2 homolog [Vitis vinifera] 1207 bp	F: TTTGTAATCCAGTCCTCGCC R: TGTGACACCTGCTGGAAAAG	(CTAGA)7	64	241–271	4	0.672	0.800	0.738	ND
E59065	Cationic amino acid transporter [Populus trichocarpa] 2235 bp	F: TCTAGTAAACCCGGCTCGAA R: CCATGACCCCTCCCTATTTT	(AGACAT)7	62	178–184	2	0.305	0.400	0.385	ND
E4746	ARM repeat superfamily protein [Theobroma cacao] 813 bp	F: TGAAGAAGGGGTTGACCTTG R: TGCCGGTGCTTTAGGTTTAC	(AGAGGA)5	59	283–289	2	0.375	1.000	0.513	**
E8796	PREDICTED: probable xyloglucan endotransglucosylase/hydrolase protein 23-like [Glycine max] 344 bp	F: GAGGTTTCCTCAAGGGCTTC R: GATCGAAATACATACAGGGGGA	(CCTAGT)5	61	165–183	4	0.628	0.600	0.708	ND
E29656	Mediator subunit 8 isoform 1 [Theobroma cacao] 2437 bp	F: TAGCTGGAAATGCAGCAAGA R: TGTTCCAAATTGAGCTGCTG	(GCA)8	56	239–245	5	0.687	0.900	0.753	ND
E60960	RNA-binding KH domain-containing protein [Theobroma cacao] 2141 bp	F: AACCATAACTTCCAGCGTGC R: ACAGAATTTCATGGCCGAAG	(TGC)8	58	170–179	3	0.532	0.950	0.619	ND
E54173	ATP-dependent Clp protease proteolytic subunit 4 [Theobroma cacao] 1561 bp	F: CAGTGTGTGACGTTTTTGGG R: TGGGATAGGTTTTCTGGCTG	(ATT)8	56	265–286	2	0.269	0.400	0.328	ND
E42984	F21B7.22 [Theobroma cacao] 1335 bp	F: GGAAAGATCCACACTTCGGA R: GAGGTGGGATGCATGTTCTT	(GTT)7	58	215–248	4	0.562	0.900	0.631	ND
E11630	Quinolinate synthase [Theobroma cacao] 3702 bp	F: CAGCAGCCATTGTTGATGAT R: TAAAGTGTGGGCAAGTGCAA	(TGATC)5	56	215–230	3	0.345	0.450	0.422	ND
E29774	HCP-like superfamily protein with MYND-type zinc finger [Theobroma cacao] 1634 bp	F: CAGTAGACCCCAGGACCTCA R: GAGTTTAAACGGCGGTGAGA	(ATTTA)5	62	208–218	4	0.560	0.900	0.629	ND
E925	Metal ion binding protein [Arabidopsis thaliana] 1158 bp	F: AAGAGCCGGAGAAGCCTAAC R: CCGTATCCGCCATAACATTC	(TCCGGC)4	60	226–256	7	0.758	0.900	0.809	ND
E60501	PREDICTED: trichohyalin-like [Cicer arietinum] 2213 bp	F: TATGCCTCCCATTTCCATTC R: GGTTTGGGTCGAATCAAGAA	(CCCATT)5	56	279–291	2	0.374	0.950	0.512	ND
E5857	None 675 bp	F: TCTCTAAGCAACACCGCAGA R: CTCCAAGCAGCAATTTAGCC	(CAATC)5	58	243–258	3	0.094	0.100	0.099	ND
E23570	15-cis-zeta-carotene isomerase [Theobroma cacao] 1996 bp	F: CCCTCAGAGATGGCGATAAA R: ATAGGGCTTGGATGAAGGGT	(CAG)7	58	213–231	3	0.580	1.000	0.672	ND
E4339	PREDICTED: LOW QUALITY PROTEIN: 1-phosphatidylinositol-3-phosphate 5-kinase-like [Vitis vinifera] 6159 bp	F: AAGGAAAGAGACGGTGGGAC R: CGAATTGCAACCCAGTGTAA	(AAAT)6	56	173–193	3	0.303	0.400	0.344	ND
E57686	Dentin sialophosphoprotein-like protein [Arabidopsis thaliana] 4774 bp	F: TGGGCTGTCAACACATTCAT R: GCCTTGCCATAAGAAACCAA	(CTGGTT)5	62	164–170	2	0.375	1.000	0.513	**
E58997	Ubiquitin interaction motif-containing protein [Theobroma cacao] 2518 bp	F: ACATGCCCAGGAAAGTTGAA R: TTGTCAACATGCTCCTTTGG	(ATCTC)7	60	240–245	2	0.375	1.000	0.513	**
E10981	Adenine nucleotide alpha hydrolases-like superfamily protein isoform 1 [Theobroma cacao] 1354 bp	F: ATCAATGTGCCTGTGGATGA R: GGAAAGGGGAACATTCCTTC	(GGTAT)6	55	290–315	2	0.269	0.400	0.328	ND

Note: Tm, annealing temperature; N_A, number of alleles; PIC, polymorphism information content; Ho, observed heterozygosity; H_E, expected heterozygosity; HW, Hardy-Weinberg equilibrium.

EST-SSRs derived from transcribed regions are highly conserved [45]. EST-SSRs have an enhanced cross-species transferability and efficiency of amplification in comparison with SSRs located in non-transcribed regions [46]. Transferability of EST-SSRs has been reported in Euphorbiaceae plants including Hevea brasiliensis and Jatropha curcas. The cross-species transferability analysis has become a very popular strategy to the development of SSR markers in closely related species [47, 48]. As a result, the developed 23 polymorphic EST-SSR markers and potential SSR loci will likely be valuable for genetic analyses of E. kansui and other close Euphorbia species.

Conclusions

In this study, we used high-throughput sequencing data to characterize the transcriptome of Euphorbia kansui, a species for which no genomic data are available. The candidate genes involved in terpenoid biosynthesis were identified. A total of 7,016 SSR loci were detected based on the transcriptome sequences. Forty ones were validated and 23 ones were polymorphic. Transcriptome data and the developed SSR markers provide a valuable resource for medicinal ingredient biosynthesis and genetic studies in E. kansui and other closely related Euphorbia species.

Disclosure statement

The authors declare no conflict of interest.

Additional information

Funding

This work was supported by the National Science Foundation of China [Grant number 31870173], Science and Technology Program of Shaanxi Academy of Sciences [Grant number 2020 K-28], Key Research and Development Program of Shaanxi Province [Grant number 2018NY-88], the science and technology plan project of Xi'an [Grant number 2019112913CXSF007SF017] and the Education Department Special Research Project of Shaanxi Province [Grant number 18JK0953].