The elite purple-flesh tubers and the conventional white-flesh tubers of yam (D. alata) were cultivated in a yam producing region (Wenzhou city, Zhejiang province, China; 121°09′48.82″ E and 28°27′53.62 N). Both were planted at the same time and cultivated in similar conditions. Tubers were harvested 10 days after new tuber emergence (DAM) and used for transcriptome analysis. Tubers of each cultivar were collected from five different plants, with a total of 15 tubers per cultivar. The tubers were washed and their skin was peeled off. The samples were labeled as DP (purple-flesh tuber) and DW (white-flesh tuber), then immediately frozen in liquid nitrogen, and stored at −80°C prior to use.
Total RNA from the DP and DW samples was extracted using the RNAiso kit for polysaccharide-rich plant tissue (Takara Biotechnology (Dalian) Co., Ltd.) and purified using RNeasy plant mini kit (Qiagen, Valencia, CA) to avoid DNA contamination. The RNA quality was analysed by measuring the absorbance at 260 nm/280 nm (A260/A280) using a ND-1000 spectrophotometer (Nano-Drop Technologies, Wilmington, DE, USA). Further, RNA Integrity Number (RIN) values were determined using a Bioanalyzer 2100 (Aligent Technologies, Santa Clara, CA) to make sure all samples had a RIN greater than 8.5. Two separate RNA pools for the DP and DW cultivars were prepared for cDNA library construction, each comprising 15 RNA samples from 15 tubers of five plants per cultivar.
Two sequencing libraries were constructed using a cDNA Synthesis kit (Illumina Inc., San Digo, CA, USA) following the manufacturer’s instructions. Paired-end (2 × 150 bp) sequencing of the cDNA libraries was performed on the Illumina HiSeq 2000 (Illumina Inc., San Diego, CA, USA). Libraries from both the cultivars yielded more than 4 GB of clean data. Sequencing was completed by the Hangzhou Woosen Bio-technology Co. Ltd. (Hangzhou, China).
The clean reads were obtained by read trimming of raw data by removing adaptors, reads in which unknown bases were more than 10%, low-quality reads with quality scores less than Q30, and low-quality bases less than (Q30) at the 3′ end. Next, the high-quality filtered reads were further assembled using a de novo assembly program Trinity (released 2011-05-19, http://trinityrnaseq.sourceforge.net/) with the main parameters “K-mer = 25, group_pairs_distance = 500, min_glue = 2, min_kmer_cov = 1” [42]. .Briefly, for each library (DP and DW), short reads were first assembled into longer contiguous sequences (contigs) according to their overlap regions, and then these reads were mapped back to the contigs based on their paired-end information. With paired-end reads it is possible to detect the contigs from the same transcript as well as the distances between these contigs. Afterwards, the contigs were further assembled, and the assembled sequences that could not be extended on either end and were defined as unigenes. Finally, the potential unigenes from DP and DW library were clustered using the TGICL clustering tool [43] to acquire a single set of non-redundant unigenes. In addition, to obtain assembly statistics profile about reads that could be mapped back to the assembled unigenes, TopHat (version 2.0.8) (released 2013-02-26, http://tophat.cbcb.umd.edu/) [44] with the parameter “mate-inner-dist = 250”, was used to align short reads to the constructed transcripts by de novo assembling.
All assembled unigenes were annotated by matching against the NCBI non-redundant protein (NR), the Arabidopsis thaliana protein dataset of NR (ATNR), Gene Ontology (GO), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) using the BLASTX analysis with a cut-off E-value of 10−5. Based on NR annotation, the Blast2GO software (version 2.3.5) was used to obtain GO annotations according to molecular function, biological process and cellular component ontologies (http://www.geneontology.org) [45]. The unigene sequences were subsequently matched against the COG database to predict and classify possible functions. The KEGG pathway annotation was also performed by comparison against the KEGG database using the online KEGG Automatic Server (KAAS) (http://www.genome.jp/tools/kaas/) [46,47].
In order to assess the differential expression between the two investigated yam cultivars, TopHat (version 2.0.8) was first used to match against the assembled unigenes, which was followed by estimation of total mapped reads [44]. After the alignment, cufflinks (version 2.1.1) (released 2013-04-11, http://cufflinks.cbcb.umd.edu) was used to estimate the abundances of unigenes as Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) [33]; and cuffdiff was carried out to perform pairwise comparisons between different investigated cultivars. Differentially expressed genes (DEGs) were further characterized and estimated using the R software module edge R (R v2.14; edgeRv 2.3.52) in term of the results from cufflinks [48]. False discovery rate (FDR) <0.05 and an estimated absolute log2 fold-change (log2 FC) ≥ 1 were used as threshold for determining significant difference in gene expression between the purple- and white flesh tubers of yam. Moreover, all DEGs were mapped to terms in the KEGG database and searched for KEGG terms to identify pathways related to purple-flesh trait in yam tubers.
Total RNA was extracted from the white and purple flesh-tubers of yam as described above. Approximately 2 mg of total RNA per sample was treated with DNaseI (Takara), and reverse transcribed into cDNA using Promega A3500 reverse transcription system. Eight DEGs were selected for Quantitative real-time PCR (qRT-PCR) analysis to verify the expression patterns revealed by the RNA-seq analysis. Gene specific qRT-PCR primers (Additional file 7) were designed using Premier 5.0 software (Premier Biosoft International, Palo Alto, CA). qRT-PCR was performed using SybrGreen qRT-PCR Master Mix (Ruian Biotechnologies, Shanghai, China) in an ABI 7500 FAST Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Amplification program comprised an initial denaturation step at 95°C for 2 min, followed by 40 cycles of denaturation at 95°C for 10 s and annealing at 60°C for 30 s. Three replicates were performed, and the amplicons were subject to melting curve analysis to determine amplification specificity. The relative expression level of the selected unigenes were normalized to UBC gene and calculated using the 2-ΔΔCt method [49].
Simple sequence repeats (SSRs) in unigene sequences were identified using MIcroSAtellite package (MISA, http://pgrc.ipk-gatersleben.de/misa). The SSR search parameters were defined to identify di-, tri-, tetra-, penta- and hexa-nucleotide motifs with a minimum of 6, 5, 4, 5 and 5 repeats, respectively. Subsequently, primer pairs were designed for genes with SSRs using the Primer3 (version 2.23) (http://sourceforge.net/projects/primer3/) with default settings [50], and the PCR product size ranging from 100 to 280 bp.
All clean reads generated by Illumina sequencing have been deposited in the Sequence Read Archive (SRA) database (http://trace.ncbi.nlm.nih.gov/Traces/sra/) under the accession ID SRX652481 for DP, and SRX652483 for DW.
