1. Introduction

The Cycadales (cycads), whose origin can be dated to the Late Paleozoic (~265-290 Ma), are the most primitive living seed plants. The earliest Cycad fossil in the world was found in China and dates to Lower Permian, whereas peak abundance and diversity of cycads date to the Mesozoic (Martínez et al., 2012; Gao and Barry, 1989). Extant cycads are distributed in tropical and subtropical regions of Africa, Asia, Oceania, and America. Ten genera and ~300 species are currently accepted (Hill et al., 2004; Norstog and Nicholls, 1997). Cycads, which originated from seed ferns, have many of same characteristics as ferns. Such characteristics include pinnately compound leaves, circinate vernation, unique girdling leaf traces, no axillary buds, dichotomous branching (versus axillary branching in higher plants). However, cycads have recognizable intermediate morphological traits between angiosperms and gymnosperms, and have therefore been classified as gymnosperms. Of the four order that comprise the gymnosperms-Ginkgoales, Gnetales, Coniferales, and Cycadales-Cycadales are considered the most ancestral (Nixon et al., 1994; Soltis et al., 2002). Their pollen tubes possess multiciliate sperm, and their ovules are borne on the margins of leaflike megasporophylls (Loconte and Stevenson, 1990). Clearly, cycads represent a key node in the phylogeny of seed plants.

Previous research has also indicated the importance of understanding cycad biological nitrogen fixation (BNF). Cycads have a particular root type, referred to as coralloid roots due to their 'corallike' appearance. Nitrogen-fixing cyanobacteria invade coralloid roots, where BNF processes occur (Vessey et al., 2005; Lindblad and Costa, 2002). To date, no symbiosis-related genes in cycads have been identified (or searched) and the molecular mechanisms of symbiotic nitrogen fixation in these plants are still largely unclear (Rai et al., 2000).

In non-model species with a large genome size, EST (expressed sequence tag) sequencing and annotation is a means for gene discovery and a way to understand the transcription and expression patterns of specific genes. Normalization techniques reduce the frequency of highly expressed genes and increase the rate of rare gene discovery (Soares et al., 1994). Cycads have a large genome size (~20-30 Gbp) and are evolutionarily and horticulturally important, long-lived plants (Zonneveld, 2012). Cycads are also the only early seed plants that have evolved a specialized coralloid root to host endophytic bacteria that fix nitrogen for the plant.

Microsatellites, also called SSRs (simple sequence repeats) or STRs (short tandem repeats), are highly polymorphic genetic markers. They have been extensively used for plant population genetic studies because of their co-dominant inheritance, relative abundance, multi-allelic nature, high reproducibility, and ease of detection (Powell et al., 1996). The development of SSR markers has improved the characterization and use of genetic variation in cycads (Ju et al., 2011; Yang et al., 2008). SSRs have been adopted to evaluate genetic diversity and to reconstruct the population structure, allowing researchers to design reasonable conservation and management protocols.

Cycas is at the basal node of the Cycadales (Treutlein and Wink, 2002). Cycas debaoensis, endemic to southern China, is a rare and endangered plant species (Chen and Zhong, 1997; Ma et al., 2003). In this study, we constructed a normalized full-length cDNA library using RNA derived from the coralloid roots of C. debaoensis. We assembled ESTs into contigs and singletons, and subsequently performed comparative protein annotations. Using this database, we identified a total of 22 ESTs that encode putative receptor-like protein kinase (RLK) genes, which play a variety of important defensive and symbiotic roles in plant-microbe interactions (Shiu and Bleecker, 2003; Shiu et al., 2004). Finally we used EST data to detect SSRs.

2. Materials and methods 2.1. Plants and total RNA extraction

Coralloid roots of C. debaoensis were collected from National Cycad Germplasm Conservation Center, Shenzhen Fairylake Botanical Garden. Tissue was obtained from the topsoil, washed with sterile water, and then 70% ethanol, wrapped with aluminum foil bags, frozen in liquid nitrogen for more than 2 h, and stored at -80 ℃. Total RNA was extracted from pulverized, frozen tissue using trizol reagent (Invitrogen). The integrity of the total RNA was confirmed by examining the ratio of 28S and 18S ribosomal RNA with 1% agarose gel electrophoresis. Quality and quantity of the isolated RNA were determined using a spectrophotometer. The extracted RNA was found to be of high quality (OD260/ OD280 = 2.02).

2.2. Construction of full-length normalized cDNA library

The cDNA library was constructed using Creator SMART cDNA library construction kit (Clontech, USA). Double stranded cDNA was synthesized according to the manufacturer's protocol. DSN normalization was applied according to the instructions of Trimmer-Director kit (Evrogen, Cat. No. NK002). The cDNA inserts were directionally cloned in pDNR-LIB vector and transformed using DH10B electrocompetent cells of Escherichia coli. The recombinants were selected in LB agar plates supplemented with ampicillin. The quality and quantity of the isolated plasmid DNA was confirmed on 0.8% agarose gels before sequencing.

2.3. EST sequencing and assembly

Plasmid DNA was sequenced by 5'end single-pass sequencing on automated DNA capillary sequencer ABI 3730XL (Applied Biosystems) according to the manufacturer's instructions using T7 primer (50TAATACGACTCACTATAGG3'). All EST sequences were deposited in the GenBank database under the accession numbers dbEST JZ917721eJZ922731. Sequence analysis was performed at Beijing Luhe Huada Gene Sci-Tech Company. ESTs were scanned and trimmed for vector sequences using NCBI's VecScreen tool. Low quality and short (< 100 bp) sequences were also removed. The processed sequences were assembled into contigs and singletons using Contig Assembly Program CAP3. Then EST and contig redundancy was calculated.

2.4. Annotation and classification of singlets

All contigs and singletons were compared against the NCBI nonredundant protein (Nr) and Swiss-Prot database (Boeckmann et al., 2003). Based on the Nr annotations, GO annotations of the unigenes were obtained using the Blast2GO program. GO terms were assigned to each unigene and classified into three functional annotation categories: biological process, cellular component, and molecular function. All unigenes were also compared to the proteins in the Cluster of Orthologous Groups (COG) databases (Tatusov et al., 2000) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008).

2.5. EST-SSR search and primer design

SSR loci of unigenes were identified and analyzed by MISA software (http://pgrc.ipk-gatersleben.de/misa/). Premier 5.0 software (PREMIER Biosoft International, Palo Alto, CA) was employed to design PCR primers for the conserved flanking regions of the SSRs. The minimum repeats of sequences with di-, tri-, tetra-, penta- and hexanucleotide motifs was set as 6, 5, 5, 5 and 5, respectively. The range of PCR product size was set between 100 and 300 bp.

2.6. Analysis of sequences

Databases used to perform BLASTN analyses included Cycas rumphii PUT (10, 901), Zamia vazquezii PUT (7, 657), Gnetum gnemon PUT (6, 193), and Ginkgo biloba PUT (10, 210). All databases as described above were downloaded from plantGDB (http://www.plantgdb.org). All BLAST searches were subject to an expect value < 1e-5.

3. Results and discussion 3.1. EST sequence quality, contigs and singlets

A total of 5393 ESTs were generated from this cDNA library; removing vector and low quality sequences resulted in 5011 highquality ESTs. The EST length was distributed from 100 to 500 bp (19.04%), 501 to 800 bp (80.42%) and 801 to 840 bp (0.54%). The cluster analyses of ESTs generated 2366 singlets and 618 contigs. The average lengths of singlets and contigs were respectively 538.87 bp and 720.66 bp. The redundancy of the library was calculated as 40.45% ((1 - Number of Unigenes/Number of ESTs) × 100%) (Tatusov et al., 2000). This is lower than that reported in other gymnosperms studies like C. rumphii (Brenner et al., 2003), Ginkgo biloba (Brenner et al., 2005), and Picea glauca (Birol et al., 2013). The number of ESTs in a contig ranged from 2 to 150. About 96.28% of the contigs had less than 10 ESTs with 338 having only 2 sequences. EST sequencing and assembly results are shown in Table 1 and the distribution of ESTs in contigs is shown (Fig. 1).

3.2. Functional gene annotation and gene discovery in unique sequences

All unique sequences (including 2366 singlets and 618 contigs) were subjected to BLASTX searches against the NCBI Nr (nonredundant) protein database. A total of 2333 unigenes had significant hits. Among the 2333 unique sequences, 577 matched proteins of known function and 1756 matched predicted proteins of unknown function. These unique sequences were also compared to proteins in Swiss-port, KEGG, and COG databases, which revealed there were 1, 901, 2255 and 951 unigenes that showed high similarity, respectively. Of these, 938 sequences had annotations in all four databases, 2343 sequences in at least one of the four database, and 641 sequences had no annotation in any database (Fig. 2).

3.3. GO and COG categories

Gene Ontology (GO) is a standard system defining gene classes (Ashburner et al., 2000). Fifty percent of 2984 unique sequences were successfully annotated with 4082 GO terms (Fig. 3). The unigenes were thus functionally classified with one or more ontologies: 954 sequences (63.8%) were assigned GO terms associated with biological processes, 1329 sequences (88.9%) were involved in molecular functions, and 415 sequences (27.8%) were cellular components. Under the category "biological process", 51.2% sequences were associated with "physiological process", 39.4% with "cellular process", 5.7% with "regulation of biological process", and 3.5% with "response to stimulus". Under the category "molecular functions", sequences associated with "binding" (805), "catalytic activity" (739), and "transporter activity" (111) were respectively 35.8%, 32.8%, and 4.9%. Under the category "cellular component", 57.3% sequences were associated with "cell", 23.7% with "organelle", and 17.8% with "protein complex".

The COG protein database is an attempt to classify orthologous gene products (Tatusov et al., 2000). All unigenes were compared to proteins in the COG database. The results showed that a total of 847 sequences were assigned to 22 COG categories (Fig. 4). The cluster for "general function prediction only" (154, 18.18%) was the largest group, followed by "posttranslational modification, protein turnover, chaperones" (132, 15.58%), "translation, ribosomal structure and biogenesis" (90, 10.63%), "carbohydrate transport and metabolism" (73, 8.62%), "energy production and conversion" (61, 7.2%), and "amino acid transport and metabolism" (56, 6.61%). In addition, the clusters for "signal transduction" (29, 3.42%) and "defense mechanisms" (5, 0.59%) were found.

3.4. Gene expression and highly abundant transcripts

The expression levels of corresponding genes can be estimated preliminarily according to the distribution frequency of transcribed ESTs in the cDNA library (Mann et al., 2013). Among the annotated sequences, 22 contained at least 10 supporting ESTs (Table 2). A family of metallothionein-like proteins was the most abundant transcript, with 470 ESTs detected. These are small metal ion binding proteins. By binding heavy metals, metallothionein-like proteins can effectively reduce the toxicity of heavy metals to the body (Zhou et al., 2005; Liu et al., 2002). Other highly abundant transcripts were genes associated with DJ-1 protein, antimicrobial proteins, germin-like proteins, aquaporin proteins, ERD (early responsive to dehydration) proteins, WAT1-related proteins, and membrane steroid-binding proteins. High abundance of genes involved in stress and defence responses is quite expected because the cDNA library was constructed from coralloid roots of C. debaoensis under natural conditions. The genes identified in our study provide a valuable transcriptomic resource for structural and functional genomics studies in C. debaoensis.

3.5. Sequence similarity and evolutionary relations

The BLASTX comparisons of unique C. debaoensis sequences were conducted against the published PUT (putative unique transcripts) of four major gymnosperms (C. rumphii, Z. vazquezii, Ginkgo biloba, and Gnetum gnemon). Phylogenetially, C. debaoensis is closest to C. rumphii, next to Z. vazquezii, Ginkgo biloba, and last to Gnetum gnemon.

Accordingly, C. rumphii and C. debaoensis had the highest matches in 1354 sequences (45.4%). The similarity results indicate that 45.4% of C. debaoensis genes discovered are homologues (orthologs) of C. rumphii genes and may have originated from a common ancestor. Ginkgo biloba followed, with 436 (14.6%) unique sequence similarities. The number of unique C. debaoensis sequences that were similar to Z. vazquezii PUTs was 429 (14.4%). Ginkgo biloba had higher similarity than Z. vazquezii since the EST resource generated for Ginkgo biloba (10, 210 PUTs) is much larger than that of Z. vazquezii (7, 657). Only 50 (1.7%) sequences had statistically significant similarity between C. rumphii and Gnetum gnemon (Fig. 5).

3.6. Simple sequence repeats

A total of 94 different 2-6 nucleotide repeats were developed from unigenes obtained in this experiment (Table 3). Trinucleotide repeats were the most frequent (46, 49%), followed by dinucleotide repeats (38, 40%), tetranucleotide repeats (8, 8.5%), pentanucleotide repeats (1, 1.1%), and hexanucleotide repeats (1, 1.1%). Among the various SSRs, AG/CT repeats were the most abundant (23, 24.5%), followed by AAG/CTT (10, 10.6%), and ATC/ATG (9, 9.6%). Our results are in agreement with previous studies that have shown that angiosperms and gymnosperms have a higher abundance of AG/CT and AAG/CTT motifs. These sequences are also known methylation targets in plants (Ranade et al., 2014). In gymnosperms, some studies report the AG/CT motif as the most abundant in Cycas (von Stackelberg et al., 2006) and Gnetum (Victoria et al., 2011). In other studies, the AAG/CTT motif has been shown to be the most abundant trimer in Picea (Rungis et al., 2004) and Cycas (von Stackelberg et al., 2006). Finally, 60 pairs of primers were designed for PCR amplification. Of these primer pairs, 20 successfully amplified the target regions using genomic DNA of C. debaoensis (Table 4). As newly developed molecular markers, they provide valuable resources for the population and conservation genetics of C. debaoensis and other cycads.

3.7. Potential candidate RLK genes involved in signaling in symbiosis and defense

The receptor-like protein kinase (RLKs) is the main receptor for plant extracellular signals (Shiu et al., 2004; Shiu and Bleecker, 2003). Plant RLKs function in diverse signaling pathways, including the responses to microbial signals in symbiosis and defense (Antolínllovera et al., 2012). Because the cDNA library was constructed from coralloid root tissues, genes involved in plantemicrobe interactions were expected to be expressed. In this study, we focus on RLKs involved in plantemicrobe interactions. Three types of the RLKs highlighted here include the leucine-rich repeat (LRR) type, lysin-motif (LysM) type, and lectin (Lec) domain type. The LRR-RLKs contain a tandemly repeated (9-26) Leu-rich motif, which plays an important role in plant development, defense symbiosis, and other biological processes (Liu et al., 2017). RLKs with lysin-motif (LysM) ectodomains confer recognitional specificity toward N-acetylglucosamine-containing signaling molecules, such as chitin, peptidoglycan (PGN) and rhizobial nodulation factor (NF), which induce immune or symbiotic responses (Antolínllovera et al., 2012). The lectin receptor-like kinases (LecRLKs) possess a characteristic extracellular carbohydratebinding lectin domain. LecRLKs play important roles in plant development and innate immunity (Vaid et al., 2013; Prashant and Laurent, 2013). Using this database, we identified a total of 22 ESTs that encoded putative RLK genes (LRR-RLKs, LysM-RLKs, and LecRLKs), which play a variety of important defensive and symbiotic roles in plant-microbe interaction (Table 5). Further study of these genes will help us to understand the signaling pathways leading to symbiosis and defense.

4. Conclusions

This study is the first to successfully construct high quality cDNA library using RNA derived from coralloid roots of C. debaoensis. We have obtained 2984 unigenes, including 641 sequences with no nucleotide similarity with public databases. These sequences are an important addition to existing databases. We have identified highly expressed genes (mainly stress-responsive genes and antimicrobial genes). A total number of 94 SSR loci were detected, of which 20 loci were successfully amplified in C. debaoensis. These SSR markers can be used for various population and conservation genetics studies of C. debaoensis and other cycads. The cDNA library also provides an excellent resource for discovering genes involved in signaling in symbiosis and defense. A total of 22 ESTs that encoded putative receptor-like protein kinase (RLK) genes were identified.

Author contributions

Conceived and designed the experiments: Yunhua Wang, Nan Li, Ting Chen; Performed the experiments: Yunhua Wang, Ting Chen; Analyzed the data: Yunhua Wang, Yiqing Gong; Wrote the paper: Yunhua Wang.

Acknowledgments

This work was supported by the Grant (201522) from Shenzhen Urban Management.