uencer. The original image data obtained by sequencing basecalling were the original sequence reads. Each read in the Solexa paired-end (PE) sequencing was 101 bp in length. There were 112.10 million reads and a 11.32 Gb original data sets produced during sequencing. After the raw data were trimmed, 57,382,380 clean reads for NJCMS1A sample and 45,599,106 for NJCMS1B sample were obtained (Table 1). All clean reads were matched to the soybean reference genome by Tophat software, allowing two base mismatches [15, 22]. As a result, 53,376,483 mapped reads for NJCMS1A sample and 42,066,351 for NJCMS1B sample were obtained, with an average matching rate of 92.64% (Table 1).To estimate whether the sequencing depth was sufficient for the transcriptome coverage, the sequencing saturation and coverage in the two cDNA libraries were analyzed. Saturation analysis showed that most genes with moderate contents of expression (genes with greater than 
3.5 FPKM) became saturated when more than 40% of sequencing reads were aligned (vertical axis numerical approached 1), which indicated that the overall quality of sequencing saturation in the two cDNA libraries was high, and sequencing amount covered the vast majority of expressed genes (Fig 1). Coverage analysis showed that two ends of the sequencing coverage in the two cDNA libraries had no significant peaks, which indicated that sequencing data among the two cDNA libraries was normally distributed (Fig 2).
A total of 88,643 transcripts were produced in the Illumina sequencing. Then 56,044 genes were obtained matching the soybean reference genome by Cufflinks 1628838-42-5 software [16] (S1 Table). “FDR 0.05 and |Log2FC| ! 1” were used as the threshold to screen 10205015 the DEGs between NJCMS1A and NJCMS1B. It was found that there were 365 DEGs between NJCMS1A and NJCMS1B (S2 Table), among which, 339 down-regulated and 26 up-regulated in NJCMS1A compared to in NJCMS1B. Furthermore, 93 down-regulated DEGs were only expressed in NJCMS1B and 9 up-regulated DEGs were uniquely expressed in NJCMS1A. Results showed that the number of the down-regulated DEGs was obviously larger than that of the up-regulated DEGs in NJCMS1A compared to in NJCMS1B. All of these RNA-Seq reads were deposited in Sequence Read Archive database (http://www.ncbi.nlm.nih.gov/Traces/sra/) under the Accession number SRP052011. Saturation analysis of sequencing data of NJCMS1A and NJCMS1B. X-axis represented the percentage of mapped reads to soybean genome (%); Y-axis represented the fraction of genes within 15% of quantitative deviation. Each color line represented the saturation curve of different gene expression level, and the gene number within different FPKM interval was displayed in the lower right corner.
Gene ontology (GO) is an internationally standardized gene function classification system used to describe properties of genes and their products in any organism, containing three ontologies: biological process, cellular component and molecular function [8]. In this study, plant GO Slim annotation was conducted by Blast2GO software (Version 2.3.5) [23]. Based on sequence homology, 242 DEGs were annotated to 19 functional categories, including 9 biological processes, 3 cellular components and 7 molecular functions (Fig 3, S3 Table). Among the biological process categories, “embryo development” was the main functional groups, followed by “cellular component organization” and “carbohydrate metabolic process”. Among the cellular component categories, “cellu
