Development of EST-SSR markers in the relict tree Davidia involucrata (Davidiaceae) using transcriptome sequencing

RNA extraction

Four accessions of D. involucrata were selected for RNA extraction. Total RNA was extracted from young leaf tissues and was immediately frozen in liquid nitrogen. For every 100 mg tissue, 1 mL TRIzol reagent (Invitrogen, USA) was added and treated with RNase-free DNase I (TaKaRa Bio, Shandong, China) for 1 h at 37°C. The extracted RNA was then diluted in RNase-free water (Ambion, USA). From each sample, 1 µL was used to check the quality and concentration by both NanoDrop (Thermo Scientific, USA) and Agilent Bioanalyzer 2100 (Agilent, USA). The four samples were pooled by mixing the total RNA at equal volumes. This step was done to ensure the quality of the transcriptional units and to enhance the downstream processes.

cDNA synthesis and sequencing

The working concentration for the cDNA synthesis was reduced to be in the order of 50-100 ng/μL. mRNA was isolated from the pooled total RNA and purified using Micropoly (A) Purist^TM mRNA purification kit (Ambion™), following the manufacturer instructions. This step was followed by cDNA synthesis employing a slightly modified protocol based on Ng et al. (2005). To synthesize the first-strand cDNA, 10 µg mRNA template was utilized, using GsuI-oligo as reverse transcriptome primers and 1000 U Superscript II reverse transcriptase (Invitrogen). The mixture was incubated at 42°C for 1 h. The cDNA was treated with sodium periodate (Sigma, USA) and the mRNA 5' oxide caps were biotinylated using Dynal M280 beads^TM (Invitrogen). The biotin-linked mRNA/cDNA was further treated with alkaline lysis to release the first-strand cDNA. The complementary strand cDNA was synthesized using Ex Taq^TM polymerase (TaKaRa Bio Inc., Shiga, Japan) and GsuI restriction enzyme was used to cut-off poly-A tails. The double-stranded cDNA was then purified using the QIAquick PCR extraction kit (Qiagen, Hilden, Germany), followed by ligation of sequencing adaptors onto the fragments. In order to maximize the quality and to enhance the accuracy of the sequencing process, uniformity of the fragments was ensured by selecting a range of 300-500 bp. These fragments were then purified using Ampure^® beads (Agencourt Canada, USA). The final step of polymerase chain reaction (PCR) was performed to enrich the fragments and to construct a library of transcripts for sequencing. Finally, the library was sequenced using Illumina HiSeq^TM 2000 platform (Illumina Inc., CA, USA).

De novo assembly and functional annotation of unigenes

A stringent filtering process was implemented and the clean reads were assembled using the Trinity software (Grabherr et al., 2011). The EMBOSS software (Rice et al., 2000) was used to predict and identify any putative open-reading frames as well as untranslated regions within the transcripts. The predicted protein-coding sequences (unigenes) were aligned against the Swiss-Prot and TrEMBL protein databases using BLASTp, with the E-value set at <1e^-5, for authentication and annotation.

We used GoPipe (Chen et al., 2005) to assign gene ontology (GO) annotations. The metabolic pathways were constructed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2010). The conditions for bi-directional BLAST were also set at E-value <1e^-5. The unigenes were assigned KEGG orthology (KO) numbers and their involvement in the metabolic pathways was defined.

SSR mining, primer designing, and marker validation

SSR motifs were detected using MiSa (http://pgrc.ipk-gatersleben.de/misa/), a Perl language-based program. Mononucleotides were excluded in our search. The search criteria were set at 6, 5, 3, 3, and 3 for di-, tri-, tetra-, penta-, and hexanucleotide repeats, respectively. We used the Primer 3.0 software (Rozen and Skaletsky, 2000) to design SSR primers under the following conditions; primer length was set between 18-22 bp, 40-60% G-C composition, 45°-60°C annealing temperature, and an expected product size of 100-300 bp.

Firstly, we randomly selected 18 SSR primer pairs and tested them for amplification and specificity using the four accessions of D. involucrata. Genomic DNA was extracted using MagicMag Genomic DNA Micro Kit (Sangon Biotech Co., Shanghai, China) following the manufacturer protocol. The 13 best primers were selected and the sense sequences were labeled with 6-FAM fluorescent dye on the 5'-end. These were used to genotype 24 accessions of D. involucrata, collected from three distant natural populations in Sichuan and Shanxi Provinces in China. The total 20 µL PCR mixture was composed of 50-80 ng/µL genomic DNA, 2.5 µL 10X Taq Buffer (with Mg²⁺), 0.8 µL 10 mM dNTPs, 0.8 µL 10 mM each primer, 0.2 U Taq Polymerase enzyme, and double-distilled water.

The PCR protocol was set as follows: an initial denaturation step of 3 min at 95°C, followed by 30 cycles of 30 s denaturation at 95°C, 30 s annealing at 50°-55°C, and 30 s extension at 72°C. A final extension step of 10 min at 72°C was added at the end of the program. Quality was checked on a 1.5% agarose gel. The PCR products were separated using an ABI 3730 XL automated sequencer (TsingKe Biotech, Beijing, China) and visualized using the GeneScan system (Applied Biosystems, Foster, CA, USA). The number of alleles (N_A), observed heterozygosity (H_O), expected heterozygosity (H_E), for each of the EST-SSR markers were calculated using GenAlex 6.5 (Peakall and Smouse, 2012).

Illumina de novo sequencing and functional annotations and classification

A total of 17.4 million paired-end raw reads (3.4 billion bp) were obtained for D. involucrata. The reads were assembled into 209,238 contigs that were reduced to 142,950 after removal of low-quality redundant contigs. The contig length ranged from 201 to 45,506 bp, with an average of 491 bp. In total, 142,432 total protein coding sequences were predicted. Of these, 25,220 (17.7%) genes had a clear biological function. Of these, 9123 (6.4%) transcripts significantly matched 3373 GO terms. These were organized into three major categories; biological processes (28,730), molecular function (25,267), and cellular component (22,870) (Figure 1). Of the predicted proteins, 2496 had significant matches in the KEGG database and could be assigned to 25 major pathways organized into five main classifications. These included metabolism, genetic information processing, environmental information processing, cellular processes, and organismal systems (Figure 2). In general, metabolism (2090) was highly represented, followed by genetic information processes (1061). A number of unigenes were found to be active in more than one biological process. The specific pathways involving the majority of the unigenes included signal transduction mechanisms (6634) and general function (5618). Cell motility (69) was the least represented pathway.

Figure 1

Classifications of assembled unigenes on gene ontology (GO). Of the 9123 predicted proteins, 3373 significantly matched GO annotations.

Figure 2

Kyoto Encyclopedia of Genes and Genome (KEGG) pathway assignment of 2496 unigenes to metabolic pathways, based on cellular process, genetic information processing, environmental information processing, metabolism, and organismal system.

EST-SSR marker development, characterization, and validation

From 142,950 assembled contigs and singletons, a total of 68,325 SSRs were identified in 52,630 sequences. Of these, 11,471 sequences contained more than one SSR motif. Among the SSRs, 25,806 (37.8%) were dinucleotide repeats, 5600 (8.2%) were trinucleotide repeats, the rest (0.8%) were tetra-, penta-, and hexanucleotides. The A/T (50.4%) was the most abundant repeat motif, followed by AG/CT (23.7%) and AAG/CTT (2.1%). The frequency analyses of the developed SSRs are shown in Table 1. A total of 12,208 primers were designed from the 30,411 identified SSRs. Of these, 7887 primers were designed from dinucleotide repeats, 3472 were designed from trinucleotides, whereas the rest were designed from tetra-, penta-, and hexanucleotides. The N_A per locus ranged from 3 (Locus1) to 12 (Locus13), whereas H_O and H_E varied from 0.083 to 1.0 and 0.102 to 0.69, respectively (Table 2).