Screening and identification of seed-specific genes using digital differential display tools combined with microarray data from common wheat

DDD

DDD comparison between seed libraries and non-seed libraries was used to identify unigene sets that were seed-specific. The libraries chosen for comparison were released on July 19, 2010 (Table 1 and Additional file 1). Following the formula mentioned in the methods, four hundred and seven unigene sets differentially expressed in selected libraries were retrieved. Among them, 108 (27%) unigenes were characterized gene models (CG), and 228 (56%) unigenes had no corresponding protein of known function but show similarity (identity > 45%) to proteins of other species (SG). In addition, the other 71 (17%) unigenes represented unknown transcripts gene models (UG). According to cluster sizes (the sum of ESTs), their expression abundance was classified into different levels, each level contains various amount of unigenes (Figure 1). From figure 1, it is clear that the ratio of the CG gene models is high when the size of unigene is large. In other words, it is more likely to characterize the unigene when the cluster size is large.

	Pool A			Pool B
Seed organs	Libraries	ESTs	Non-seed organs	Libraries	ESTs
Endosperm	5	11346	Callus	1	9685
Seed embryo	2	4463	Crown	4	14529
Seed	8	67394	Leaf	25	76992
Kernel	4	13850	Shoot	12	94514
Grain	9	54579	Root	16	169382
Total	28	151632	Seedlings	2	2305
			Heads	4	5143
			Liquid cultured tissue	1	10164
			Sheath	1	1068
			Spike	2	26307
			Floral organs	12	60963
			Young spikelets	1	2093
			Total	81	473145

The details for the libraries were listed in additional file 1.

Identification of genes specifically expressed in seeds using microarray data

Based on extensive EST collections, Affymetrix GeneChip platforms have now been developed for wheat [18]. The microarray data for the genes screened above with DDD were searched against PLEXdb (http://www.plexdb.org) [19]. A total of 322 probesets corresponding to 256 (63%) unigenes were analyzed here. To better assess the data, a heat map was generated through the GenePattern program (http://www.broadinstitute.org/cancer/software/genepattern/index.html) [20] (Figure 2). Figure 2 indicated that not all probesets showed specific expression in seeds. The probesets with expression levels significantly higher in seeds than any other organs (more than 4-fold) were chosen for further analysis. Finally, 236 probesets corresponding to 201 unigenes were selected (Additional file 2). Among them, 47 (23%) were CG gene models, 119 (59%) were SG gene models and the other 35 (17%) were UG gene models. Furthermore, among the 201 unigenes, most (173 unigenes) showed specificity in the 22 DPA caryopsis (embryo and endosperm), while 14 unigenes showed specificity in the 3-5 DPA caryopsis.

Microarray-alone analysis was also performed to screen probsets corresponding to genes with seed-specific expression. The 236 probsets screened above were used as "seeds" to find seed-specific genes in wheat microarrays. The searching tool was "Profile Neighbors" in PLEXdb (http://www.plexdb.org//modules/glSuite/gl_main.php) [19]. One hundred and seventy-four additional probsets were retrieved by removing the redundant information and probsets with expression profiles that were not seed-specific. The EST profiles were also screened for the unigenes corresponding to the 174 probesets. However, 95 of the corresponding unigenes were not for seed-specific ESTs.

GO enrichment analysis

Further GO enrichment analysis was performed for the 201 seed-specific unigenes through GOEAST (Figure 3). GOEAST is a web based software providing easy to use, visualizable, comprehensive and unbiased Gene Ontology (GO) analysis (http://omicslab.genetics.ac.cn/GOEAST/index.php) [21]. The details of the enrichment are listed in table (Additional file 3). Most of the unigenes were assigned into several categories, including stress response, defense response, multi-organism process, pathogenesis, extracellular region, nutrient reservoir activity, enzyme inhibitor activity, antioxidant activity and oxidoreductase activity (Figure 3).

According to the expression specificity analysis results, 14 of 201 unigenes were specially expressed in the 3-5 DPA caryopsis. At 3-5 DPA, the caryopsis is at the coenocytic stage during which free nuclear division occurs within the primary endosperm cell to give a coenocytes. In wheat, over 2000 nuclei will be generated within 72 hours after fertilization [22]. During this stage the seed is far from mature and has not begun accumulating grain storage molecules. This explains why the 14 unigenes whose expressions were specific in the 3-5 DPA caryopsis didn't contain any storage protein transcripts. Meanwhile, 173 (86%) unigenes were found to be specifically expressed in the 22 DPA caryopsis. During the grain filling stage (14-24 DPA) the major metabolic activity is the synthesis and accumulation of storage molecules such as starch and protein. Therefore, the genes belonging to the "nutrient reservoir" category all represented in this stage. Ten genes, which have specificity in the 22 DPA caryopsis, were found in the category of "response to stress". This is consistent with the result of Becerra et al. [23] who found 6 genes specially expressed in the Arabidopsis seed that were involved in response to abiotic stress, which clearly indicated the importance of genes contributing to stress-tolerance during seed development. Interestingly, there are also 9 genes with function annotation as "serine-type endopeptidase inhibitor activity" showed specific expression in the 22 DPA caryopsis. Three genes of this category encode serpins, which are inhibitors of exogenous proteinases (capable of breaking down seed storage proteins); play important roles in the defence of specific cell types of vegetative tissues [24, 25]. Another three genes encode alpha-amylase inhibitor which plays a similar role as protease inhibitors and protect cereal seed from the attack by endogenous hydrolase [26]. Thus, during the filling stage, in addition to transcript accumulation of genes encoding storage protein, genes involved in protecting the main nutritional resource of the developing seeds from endogenous as well as exogenous attack are also seed-specific and active.

The expression patterns of homologous gene groups in rice

Both rice and wheat belong to the plant family Gramineae and rice is the best-characterized experimental model for monocot plants. Further more, these two species are evolutionarily related. To investigate whether the expression patterns of homologous gene groups involved in seed development are conserved between wheat and rice, the 201 unigenes which showed seed-specific in wheat were compared with homologous genes in rice. For some unigenes, no significant homology was found in the rice genome, while for others, more than one unigene corresponded to the same homologous gene in rice. In the end, 121 (60%) unigenes corresponding to 106 putative homologous genes in rice were inferred based on sequence similarity. To facilitate the comparison of 106 homologous genes between wheat and rice, they were divided into three categories as follows: 13 gene models with high similarity (HS: >90% in the aligned region), 52 gene models with moderate similarity (MS: 70-90% in the aligned region) and 41 gene models with weak similarity (WS: 45-70% in the aligned region) to their counterpart genes.

To examine whether the three gene models have a seed-specific expression pattern in rice, the IDs of probesets representing the homologous genes and the microarray data presenting on the Affymetrix rice genome array were retrieved using an online tool of internet based database - Rice Multi-Platform Microarray Search tool (http://www.ricearray.org/matrix.search.shtml) [27]. A total 98 homologous sequences were detected, while 8 of the 106 homologous sequences had no corresponding probe on the rice arrays. Hierarchical clustering analysis based on array data had also been carried out (Figure 4). The results indicate that the expression patterns of homologous gene groups are conserved at different degrees between wheat and rice. Though not all the analyzed genes showed seed-specificity in rice, a large percentage of corresponding genes in each gene models were specially expressed in the seed of rice. Genes are defined as seed-specific only if their expression levels in the seed is significant higher (more than 4-fold) than that in any other organs. Finally, the rice counterparts of 9 out of 13 (69%) HS gene models, 31 out of 48 (65%) MS gene models, and 23 out of 37 (62%) WS gene models were all found to be seed-specific. The results also indicated that the percentage of seed-specific genes among HS, MS and WS gene models both differed and were lower where less sequence similarity occurs.

In order to assess the quality of the list, we also performed a sensitivity analysis. Seventy genes were identified as seed-specific in rice by the same method (DDD combined the microarray) which was used in wheat. Meanwhile, microarray-alone analysis in rice was also performed using the same method applied in wheat and 358 probsets were screened as seed-specific genes. The libraries chosen for DDD comparison were released on Apr 5, 2011 (Additional file 1). The 70 seed-specific genes in rice were compared with the genome of wheat and 43 (61%) homologous genes were found. Among the 43 homologous genes in wheat, 27 were seed-specific and detected in the list of 201 wheat seed-specific unigenes screened here.

GO enrichment analysis was also performed for the 70 seed-specific genes in rice through the GOEAST [21] (Additional file 4). The result showed that most rice seed-specific genes were also associated with oxidation-reduction process, response to stress, defense response, envelope, extracellular region, nutrient reservoir activity and enzyme inhibitor activity. These results suggest that most of the seed-special genes in both wheat and rice have shared functions.

Experimental validation of the expression patterns of 10 selected unigenes by qRT-PCR

To confirm the results of DDD and microarray screening, ten unignes were selected for validation of their expression specificity with quantitative real-time PCR (qRT-PCR). Two transcripts (Ta.54227 Ta.2291), which were previously shown to have consistent expression across samples, were chosen as internal controls for all the data normalization [28]. Among the 10 genes, 6 were randomly selected from the CG gene models generated by both DDD and microarray selection (Figure 5A), two were randomly selected from the group generated by DDD based selection but without microarray data (Figure 5B). All of the eight genes tested by qRT-PCR showed the same seed specificity as in silico analysis. Meanwhile, the other two genes which didn't follow the formula used in the DDD selection, that is the rations in poolB were indicated as <0.00005 while the sum of ESTs in poolA is less than 4 times than that in poolB, were also selected for qRT-PCR analysis (Figure 6). The results showed that none of the two genes showed seed specificity.

Plant materials and growth conditions

Triticum aestivum L. cv. Jingshuang 16 (winter wheat) was cultivated at the experimental field of the Chinese Academy of Agricultural Sciences from October, 2009 to June, 2010 under natural growing conditions.

The following tissues were collected for RNA extractions: (i) flag leaves and stems; (ii) young roots and leaves of the 20-day-old seedlings, 20 days after planting; (iii) single floral organs (glumes, palea, lemma, stamens, pistil) from fully emerged spikes; (iv) fully emerged spikes before flowering; (v) developing seeds of 5, 10, 15, 20, and 30 days post-anthesis (DPA). Collected samples were immediately frozen in liquid nitrogen and stored at -80°C for future use. The main spikes of the plants were tagged at anthesis and only seeds in the middle of each spike were harvested.

RNA extraction

Total RNA was extracted from frozen samples by using TRIZOL Reagent (Invitrogen). Genomic DNA was removed by digesting each sample (20 μg of total RNA) with DNaseI (Takara) according to the manufacture's instruction. Eight microlitres of RNA treated with DNaseI was reverse transcribed according the protocol of Omniscript Reverse Transcriptase Kit (Takara).

qRT-PCR

qRT-PCR was performed in a 20 μl volume containing 10 μl 2 × SYBR^® Premix Ex Taq™ (TaKaRa), 2 μl 50-fold diluted cDNA, 0.15 μl of each gene-specific primer and 8 μl ddH₂O [33]. The PCR conditions were as follows: 95°C for 3 min, 40 cycles of 95°C for 15 s, 59°C for 15 s and 72°C for 20 s. Three replicates were used for each sample. Reaction was conducted on a CFX96 Real-Time PCR Detection System (Bio-Rad). All data were analyzed using the CFX Manager Software (Bio-Rad).

DDD

The DDD tool (http://www.ncbi.nlm.nih.gov/UniGene/ddd.cgi) was used to screen the seed-specific genes in both wheat and rice. Libraries chosen for DDD comparison are listed (Additional file 1). The output provided a numerical value in each pool denoting the fraction of sequences within the pool that mapped to the unigene cluster. In order to control the expression level of unigenes in poolB, only unigenes whose expression ratio in poolB was <0.00005 and the quantity of ESTs in poolA was >4-fold higher expression in seeds compared with poolB, were chosen in wheat. In addition, a 3-fold criteria was set for the DDD selection of wheat (data not shown), which resulted in an increase in the false positive rate. Since there are fewer seed libraries available for rice compared with wheat, the 3-fold criteria did not result in an increase in false positive rate, and thus a threshold of > 3-fold higher expression in seeds of rice was selected.

Extraction of wheat and rice microarray data

The expression data for the genes studied in different developmental contexts in wheat were obtained from PLEXdb (http://www.plexdb.org) [19] experiment TA3 [18]. The expression values from the following tissues: root, leaf and crown from seedling, immature inflorescence, floral bracts, anthers and pistil before anthesis, 3-5 DAP caryopsis, 22DPA embryo and endosperm were retrieved. Eighty-two of the 407 selected unigenes were not used in this analysis because their ESTs have no corresponding probeset in the Wheat Affymetrix 61 K GeneChip^®. Among the remaining, one or more probesets were found to be corresponding to the individual unigene. However, a number of probesets whose name are suffixed by "_s_at", "_x_at" or "_a_at" which is known to cross-hybridize in one way or another were disregarded, and those probesets whose name are suffixed by ".A1" were also disregarded, because they are predominantly of the wrong orientation [18]. Finally a total of 322 probesets corresponding to 256 (63%) unigenes were found.

The IDs of probesets representing the homologous genes and the microarray data presenting on the Affymetrix rice genome array were retrieved. Rice Multi-Platform Microarray Search (http://www.ricearray.org/matrix.search.shtml) [27] tool available at National Science Foundation Rice Oligonucleotide Array Project was used. The experiment GSE6893, which was used for analyzing the spatial and temporal gene expression in various tissues and various stages of reproductive development of rice [34], was also chosen in our study. The expression values from the following tissues and development stages: seedling, seedling root, mature leaf, Y leaf (leaf subtending the shoot apical meristem [SAM]), SAM, and various stages of panicle (P1-P6) and seed (S1-S5) development were retrieved. For accuracy of the expression analysis, only the expression data from those probesets with uniquely matched rice gene model in the rice genome were considered.

Microarray-alone analysis was also performed as such: the probsets corresponding to the 201 seed-specific genes in wheat and 70 seed-specific genes in rice were used as "seed" probsets to retrieve all the seed-specific genes in the microarray. The web based tool was called "Profile Neighbors" and found in PLEXdb (http://www.plexdb.org//modules/glSuite/gl_main.php) [19]. The neighbours were calculated using Pearson Correlation with a minimum correlation value of 0.8.

To assimilate the data generated from wheat and rice, log2 digital data (the average of three replicates per sample were taken) were normalized, and heat maps were generated through GenePattern program (http://www.broadinstitute.org/cancer/software/genepattern/index.html) [20]. Hierarchical clustering method was used to discriminate and visualize the patterns of gene expression in various organs. Hierarchical clustering [35] allows for the analysis of the relationship of expression patterns among different genes on the array.

Homology search between rice and wheat

For the 47 CG gene models, the protein sequences were used to blast against rice protein database. Of these 47 models, the 35 best matches were extracted based on sequence identity (identity > 45%). For the remaining unknown function clusters, the ProtEST section of each unigene allows the user to explore pre-computed protein similarities for the cDNA sequences found in a cluster. Possible protein products for the gene are suggested by providing protein similarities between one representative sequence from the cluster and protein sequences from elected model organism [12]. Among the remaining unknown function clusters, 86 (56%) can provide a similar protein from rice. Finally a total of 121 (60%) unigenes could find best-matched proteins in rice, and then the gene accessions corresponding to the proteins were retrieved. In rice, the proteins of 70 seed-specific genes were used to blast against wheat genome database and 43 (61%) homologous genes (identity > 45%) were found.

GO enrichment

The probsets of 201 and 70 seed-specific genes in wheat and rice, respectively, were used for GO enrichment analysis by the GOEAST (http://omicslab.genetics.ac.cn/GOEAST/index.php) [21]. In this study, hypergeometric distribution was used to calculate the p-value of GOID enrichment, and a p < 0.05 cut off value was applied. The smaller the p-value is, the more significant the GO term is enriched in the dataset. And the graph size was reduced by condensing non-significant nodes to points in figure 3.