Transcriptome Analysis for Abnormal Spike Development of the Wheat Mutant dms

Background Wheat (Triticum aestivum L.) spike development is the foundation for grain yield. We obtained a novel wheat mutant, dms, characterized as dwarf, multi-pistil and sterility. Although the genetic changes are not clear, the heredity of traits suggests that a recessive gene locus controls the two traits of multi-pistil and sterility in self-pollinating populations of the medium plants (M), such that the dwarf genotype (D) and tall genotype (T) in the progeny of the mutant are ideal lines for studies regarding wheat spike development. The objective of this study was to explore the molecular basis for spike abnormalities of dwarf genotype. Results Four unigene libraries were assembled by sequencing the mRNAs of the super-bulked differentiating spikes and stem tips of the D and T plants. Using integrative analysis, we identified 419 genes highly expressed in spikes, including nine typical homeotic genes of the MADS-box family and the genes TaAP2, TaFL and TaDL. We also identified 143 genes that were significantly different between young spikes of T and D, and 26 genes that were putatively involved in spike differentiation. The result showed that the expression levels of TaAP1-2, TaAP2, and other genes involved in the majority of biological processes such as transcription, translation, cell division, photosynthesis, carbohydrate transport and metabolism, and energy production and conversion were significantly lower in D than in T. Conclusions We identified a set of genes related to wheat floral organ differentiation, including typical homeotic genes. Our results showed that the major causal factors resulting in the spike abnormalities of dms were the lower expression homeotic genes, hormonal imbalance, repressed biological processes, and deficiency of construction materials and energy. We performed a series of studies on the homeotic genes, however the other three causal factors for spike abnormal phenotype of dms need further study.


Results
Four unigene libraries were assembled by sequencing the mRNAs of the super-bulked differentiating spikes and stem tips of the D and T plants. Using integrative analysis, we identified 419 genes highly expressed in spikes, including nine typical homeotic genes of the MADS-box family and the genes TaAP2, TaFL and TaDL. We also identified 143 genes that were significantly different between young spikes of T and D, and 26 genes that were putatively involved in spike differentiation. The result showed that the expression levels of TaAP1-2, TaAP2, and other genes involved in the majority of biological processes such as transcription, translation, cell division, photosynthesis, carbohydrate transport and metabolism, and energy production and conversion were significantly lower in D than in T.

Conclusions
We identified a set of genes related to wheat floral organ differentiation, including typical homeotic genes. Our results showed that the major causal factors resulting in the spike abnormalities of dms were the lower expression homeotic genes, hormonal imbalance, repressed biological processes, and deficiency of construction materials and energy. We performed a series of studies on the homeotic genes, however the other three causal factors for spike abnormal phenotype of dms need further study.

Introduction
University, Zhengzhou, Henan, China and conserved in our laboratory (Fig 1). The other traits such as tiller number, internode number, floret structure and fertility of both the semi-dwarf and tall plants were almost the same as that of 'Zhoumai 18', However, the dwarf plants (D) were significantly different from 'Zhoumai 18' and T plants, The average plant height, tiller number, and spike number per plant of D were only 0.29 m, 3.9 and 2.4, respectively, they decreased by 62.4, 80.5 and 84.6% compared with those of 'Zhoumai 18' and their spikes were sterile (Table 1) [25]. Individual plant lines of mutant dms were sown from 2009-2011. Because D plants do not produce seeds, the mutant gene was kept by sowing seeds from the M plants every year. The field experiments were carried out in a completely randomized design. Forty populations (eleven and 29 derived from the T and M plants of dms respectively) were sown in Henan Agriculture University experimental field in 2012-2013. The segregating populations were sown in plots of 3.0 m in length and 2 m in width, the distance between rows was 21 cm, 30 seeds were planted in each row. A hundred and twenty seeds of each population (the progeny of M, that segregated into T, M, and D at ratio of 1:2:1; the progeny of T that did not segregate any more) were sown. Fertilizer and weed management were similar to wheat breeding [25].

Observation of spikes and stem tips
The young spikes and young stem tips were observed with an inverted microscope (Olympus 3111286), and images were captured by a camera (Nikon Coolpix 4500). The young spikes at stages from glume to stamen / pistil primordium differentiation and young stem tips (small fragment just connected with young spikes) of D and T genotypes were taken from the segregating populations (progeny of M) (Fig 2). Each of the four super bulks was consisted of 20 individuals to average gene expression background and biological variability. Young spikes of D, M and T at different developmental stages [26] were taken for real time quantitative RT-PCR.

Transcriptome mRNA sequencing
The transcriptome mRNA sequencing and analysis were carried out in BioMarker Company (Beijing, China). Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Four super bulk samples: T1, young stem tips of D at Z31 stage (from 14 to 20 March 2013); T2, young spikes of D in the phase from glume to stamen / pistil primordium differentiation (from 15 to 22 March 2013); T3, young stem tips of T at Z31 stage (from 9 to 15 March 2013) and T4, young spikes of D in the phase from glume to stamen / pistil primordium differentiation (from 9 to 15 March 2013). Pair-end sequencing was carried out on Illumina HiSeq TM 2000 using the standard protocol provided by Illumina company (San diego, CA, USA).

Sequencing data analysis
The sequencing data were assembled and analyzed as no reference genome sequences. Transcriptomes (including contigs, transcripts and unigenes) of the four bulk samples were assembled using software Trinity [27] with parameters of K-mer = 25, group pairs distance = 500 (maximum length expected between fragment pair). Unigenes were clustered using Blat (The BLAST-Like Alignment Tool; http://genome.ucsc.edu/cgi-bin/hgBlat) with parameters of tile size = 8, step size = 5. Open reading frames (ORFs) of the unigenes in the assembled unigene library were analyzed using a program Getorf in a software pack EMBOSS (http://emboss. sourceforge.net/apps/cvs/emboss/apps/getorf.html). Sequence data of each sample were compared with gene sequences in unigene library, gene expression enrichment was analyzed based on the different numbers of reads. Quantitative gene expression was calculated with the method of RPKM [28]. The formula is: Differentially expressed genes (DEGs) were searched based on gene expression enrichments in different samples using software IDEG6 (http://telethon.bio.unipd.it/bioinfo/IDEG6/),

Real time QRT-PCR
To study the gene expression profiles in the whole spike developmental stage, total 64 samples of each genotype were prepared (S1 Table). Total RNA was extracted as described for sequencing. cDNA was synthesized from 1.

Results
The expressions of most genes in differentiating spikes of D were lower Through whole transcriptome mRNA sequencing, a total 27.59 G data, 136,596,533 pair reads were obtained from the four super bulk samples: young stem tips of D, young spikes of D, young stem tips of T and young spikes of T (sample number: T1, T2, T3 and T4). Total 128,619 unigenes were de novo assembled, 20,296 of them were over 1 kb in length. 90,424 unigenes were annotated for biological functions. The transcriptome data of the young spikes on T and D were similar, including reads number (35,703,321; 35,329,399), base number (7,211,384,222; 7,135,855,487) and de novo assembled unigenes (65,316; 63,800). Generally, the gene expression levels of differentiating spikes of T were relatively higher than that of D (Fig 3).

The genes involved in wheat floral organ specification
To identify the genes involved in wheat floral organ specification, the transcriptomes of young spikes on D and T were compared with those of their stem tips. Total 419 genes highly expressed in spikes had been screened out. Functional classification referred to GO showed that the following aspects were significantly different between wheat differentiating spikes and their stem tips: (1) cellular component: organelle part, macromolecular complex, membraneenclosed lumen, cell junction and nucleoid; (2) molecular function: nucleic acid binding transcription factor, molecular transducer, structural molecule, receptor, enzyme regulator, protein binding transcription factor and transporter activity; (3) biological process: biological regulation, metabolic process, cellular component organization or biogenesis, developmental process, locomotion, multicellular organismal process, establishment of localization, reproductive process, reproduction, death and cell proliferation (Fig 4). Except for molecular transducer, receptor and antioxidant activity, all the others were more active in differentiating spikes than in their stem tips. These suggested that the most gene expressions, cell structure and major biological processes were activated when transition from vegetative to reproductive growth. The active biological processes of biological regulation, developmental process, reproductive process and cell proliferation exactly reflected the typical characters of the differentiating spikes, which were at the early stage of wheat spike specification. The KEGG analysis suggested that the purine and pyrimidine metabolism, DNA replication (S1 Fig (Table 2).
According to gene functional annotation, 40 genes among the 419 genes highly expressed in differentiating spikes were putatively involved in wheat spike development ( Table 3). The genes of DNA replication, transcription, translation factors (gene expression), ABC transporter B family, protein kinase, carbohydrate metabolism, cytochrome, signal transduction and general function prediction only were involved in the development of various components of wheat floral organ (Table 3). Both 'gene expression' and 'carbohydrate metabolism' had nine unigenes each (total 18 genes) that implied their important roles in spike differentiating. The genes of ABC transporter B family might also played an important role because three significant DEGs had been found, that suggested transmembrane production transport was key for floral organ differentiation.
Noticeably, nine of the 40 floral organ differentiation related genes belonged to the typical floral organ determining genes, homeotic genes, including MADS-box transcription factor, AP2-like ethylene-responsive transcription factor, FLORICAULA/LEAFY (FL) and DROOP-ING LEAF (DL) genes. To further explore the homeotic genes in the unigene libraries, we searched them in the DEGs of T1 vs T2 and T3 vs T4, four more homeotic genes were found   (Table 3, marked with ' Ã '). The homolog of Oryza sativa MADS47 (OsMADS47) and homolog of OsMADS55 were significantly highly expressed in stem tips than in differentiating spikes of both D and T; homolog of OsMADS6 and homolog of OsMADS34 almost only expressed in differentiating spikes. A homolog (T4-56463) of OsMADS18 was highly expressed in differentiating spikes than in stem tips of T, contrarily, no copy was detected in D genotype. One more AP2-like ethylene-responsive transcription factor gene (T2-51524) was found in T, but none in D. The significantly high expression or exclusive expression in spikes suggested the homeotic genes played a determining role in wheat floral organ differentiation ( Fig 5).

The significant DEGs between spikes of T and D
The spikes of D are male sterile with 1-6 pistils in each floret, the other genotypes are normal and fertile. To explore the genes whose expressions were significantly different in spikes of T and D, the differentially expressed genes (DEGs) of differentiating spikes were carefully screened. A total 143 DEGs were obtained. Regarding D as control, 94 genes were significantly up regulated and 49 genes were down regulated in differentiating spikes of T (Fig 6; Table 4; S3 Table).
To screen out spike specifically expressed genes in the 143 DEGs (S3 Table), they were compared with the 419 genes significantly highly expressed in both spikes of T and D (Fig 9). Seven specifically expressed DEGs in spikes were obtained ( Table 5). The homologic gene search in NCBI indicated that the unigene T2-63525 was the gene of dbj|AK334518.1| from Chinese Spring, which was a gene specifically expressed in inflorescence (http://www.ncbi.nlm.nih.gov/ UniGene/ESTProfileViewer.cgi?uglist=Ta.67718). The dynein light chain LC6 like gene family was mainly expressed in wheat spike and stem (http://www.ncbi.nlm.nih.gov/UniGene/ ESTProfileViewer.cgi?uglist=Ta.21814). Noticeably, five of the seven genes belonged to biological process of cell cycle control, cell division and chromosome partitioning. Three of the five genes were wheat dynein light chain LC6 (TaLC6). Furthermore, the expressions of all the three TaLC6 genes were significantly lower. It should be an important gene family that played a fundamental role in differentiating spikes.
The functions of the DEGs were annotated referred to seven databases (Table 6). Twentysix out of the 143 DEGs were involved in wheat differentiation and development according to the functional annotation referred to GO (Table 4; S3 Table).
Functional classification of the DEGs in differentiating spikes of T and D were carried out referred to databases COG and GO (Fig 10). The results indicated that, (1) cellular components of extracellular region part, extracellular matrix and extracellular matrix part were significantly changed; (2) molecular functions of transporter, nucleic acid binding transcription factor, molecular transducer, structural molecule, receptor, enzyme regulator, antioxidant and nutrient reservoir activity were significantly changed; (3) biological processes of response to stimulus, biological regulation, developmental process, multicellular organismal, multi-organism,  growth, immune system, death, locomotion, pigmentation, cell proliferation, and biological adhesion were significantly changed (Fig 10).
Regarding young spikes of D as control, KEEG pathway enrichment analysis indicated that 23 pathways were up regulated (many DEGs in the pathway were up regulated), 7 pathways were down regulated, 3 pathways had up and down regulated DEGs in differentiating spikes of T (Tables 6 and 7). Based on the computer analyses of the DEGs referred to databases, manual analysis was also carried out (Table 4; S3 Table). These data demonstrated that mutant significantly repressed some major metabolism pathways, that especially were photosynthesisantenna proteins, photosynthesis, carbon fixation in photosynthetic organisms, glycolysis / gluconeogenesis, galactose metabolism, fructose and mannose metabolism, these were close associated with production of 'energy' and 'construction materials' for plant differentiation and growth.

Expression profiles of twelve genes in differentiating spikes
To validate the reliability of the transcriptome sequencing results, and further explore the expression profiles of the important DEGs (Table 4), twelve DEGs considered as important genes involved in wheat spike differentiation were selected first to do expression profiling in spikes of D, M, T and 'Zhoumai 18' (Table 4, Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table.).  Table 3. doi:10.1371/journal.pone.0149287.g005 Mutant significantly suppressed the expression of a putative alpha-L-fucosidase 1-like gene (T4-52638, K), however it was highly expressed in T, M and 'Zhoumai 18', reached to several ten thouthand (control gene grdph) or several hundred thousand times (control gene actin) compared with the internal control genes, that suggested its transcript was abundance and it   The expression of a heat shock protein 81-1-like gene was increasing with spike development, its expression level was relatively lower in D (T1-38135, A). TaLC6 was suppressed at early stage in D, but highly expressed at later stage (T1-47968, B). A gene broadly similar to an uncharacterized 35.5 kDa protein in transposon Tn4556 was an unstable expressed gene (T1-49256, C). A probable alpha, alpha-trehalose-phosphate synthase gene was lower expressed at early stage and highly expressed at later stage in D (T1-71202, E). An UDP-glucose 4-epimerase 3-like gene was highly expressed at early middle stage in T and later in D (T4-20296, J). The expression profiles of this five genes were different between D and other three genotypes, but not sharp as that of the above seven genes. Their expressions were also significantly effected by mutation, so they might also play important roles during spike development. T2: gene nubmer of spikes on D, T4 gene number of spikes on T. * Genes whose spatio-temporal expression profile was analyzed by real time QTR-PCR (Figs 7 and 8).
** The accession number with 'sp' is that of database SwissProt, with 'tr' is that of TrEMBL doi:10.1371/journal.pone.0149287.t004  Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table. doi:10.1371/journal.pone.0149287.g007

Super-bulks at appropriate developmental stages guaranteed a reliable result
Compared with cDNA microarray [30], next generation sequencing (NGS) technologies are high throughput and more direct approaches for fine-mapping and gene isolation in many species [31], but sequencing technology can not eliminate biological variability [32]. An elaborately designed RNA-sequencing experiment can eliminate biological variability and produce reproducible results. In this paper, mRNA sequencing combined with super bulked segregant analysis (S-BSA) successfully profiled the whole transcriptomes of differentiating spikes on D and T plants of a wheat mutant dms. The genes specifying the original meristem of stamens and pistils should express at stage of pistil / stamen primordium, thus we pooled 20 individual spikes in this phase for each superbulked sample, which was more than that in most researches with only one or several individuals [32]. The transcriptomes of spikes and their stem tips of both T and D were parallel experiments in fact (Fig 6). Thus the approach not only guaranteed the key genes for floral organ differentiation were in the pools, but also mostly eliminated the confounding of the background genes and biological variability. A total 143 significant DEGs were screened out between D and T genotypes, which was smaller than that of the similar studies [30]. However, the expression profiling of the twelve genes completely confirmed the results of the RNA sequencing without exception (Fig 7A-7L Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table). Furthermore, most expression profiles of Transcriptome Analysis for Abnormal Spike Development of dms the homeotic genes in D and T were highly consistent (Fig 5), these demonstrated the sequencing was successful. The global expression profiles in this study truly reflected the characteristics of the mutant, and provided a set data of transcriptomes for further studies.

Elements specify wheat floral organ and influence spike development
We identified 419 genes highly expressed in wheat differentiating spikes, which highlight the elements specifying floral organ and influencing spike development.
Homeotic genes specify wheat floral organ. We have identified eight MADS-box genes in this study, six were significantly highly or exclusively expressed in differentiating spikes, they were homologs of OsMADS6 [33], OsMADS15 [33], OsMADS18 [33], OsMADS32 [34,35] and OsMADS34 [34,35] (Table 3, Fig 5). Transition from the vegetative to the reproductive growth phase in winter wheat needs vernalization, that is determined by allelic differences in the MADS-box transcription factor VERNALIZATION1 [32], which is closely related to the three paralogous Arabidopsis floral meristem identity genes APETALA1 (AP1), CAULI-FLOWER (CAL) and FRUITFULL (FUL). VRN1 is also homologous to FUL2 (= OsMADS15) and FUL3 (= OsMADS18), VRN3 is homolog of FLOWERING LOCUS T (FT) in A. thaniala [36]. In rice, a MADS-box gene PANICLE PHYTOMER2 (PAP2 = OsMADS34) combined with the AP1/FUL2/FUL3 homologs were needed for the transition of the shoot apical meristem to the reproductive stage [37]. Coexpression meta-network of OsMADS34 in rice and knock-out phenotypes characterized in surrounding network for OsMADS34 indicate its key role in spike development. The OsMADS6 is AGAMOUS-LIKE6 (AGL6) gene, which regulates floral organ identity and floral meristem determinacy [38]. The WFUL2 (wheat FUL-like genes) of wheat has class A functions in specifying the identities of floral meristems and outer floral organs (lemma and palea) through collaboration with class B and class E genes [39]. The stem tips used in this study were very short fragments connected with the differntiating spikes (Fig 2), and the homologs of OsMADS34 (TaSEP-5), OsMADS6 (TaAGL6) and OsMADS15 (TaAP1-3) almost exclusively expressed in differentiating spikes (Fig 5). The spike specific expression patterns demonstrated the six MADS-box genes played determining roles in wheat floral organ differentiation.
Two homologs of O. sativa AP2-like ethylene-responsive transcription factor gene (= AIN-TEGUMENTA, ANT), one homolog of O. sativa transcription factor FLORICAULA/LEAFY (FL), and one homolog of O. sativa protein DROOPING LEAF (DL) were identified in this study (Table 3), they all were significantly highly expressed in differentiating spikes (Fig 5). The AP2/ETHYLENE RESPONSE FACTOR (ERF) transcription factor genes are typically expressed in the organ primordia that they help to specify floral organ [40]. In Arabidopsis, integuments depend on the expression of the AP2 /EREBP transcription factor ANT. The FL genes encode plant-specific transcription factors and play major roles in the reproductive transition as well as floral organ development [41][42][43]. In grasses the homeotic function for carpel differentiation is attributed to a YABBY domain gene named DROOPING LEAF (DL) [44][45][46][47]. The expression patterns suggested that the homologs of ANT, FL and DL identified in this study were also key factors for wheat floral organ differentiation (Fig 5).
Floral organ development studies indicate that most genes involved in the differentiation of floral meristem and organ identity and in the transition to flowering in some model plants such as Arabidopsis, Antirrhinum, and Petunia belong to conserved and widespread gene families encoding transcription factors, such as the MADS-box genes, FLO-like, and AP2-like genes [48]. Currently, the ABCDE model for flower development proposes that floral organ identity is defined by five classes of homeotic genes, named A, B, C, D and E, class A genes specify the identity of sepals, class A and B genes specify petal identity, class B and C genes determine stamen identity, and class C genes determine carpel identity, class D genes are crucial for ovule Energy production and conversion 7 0 7 Replication, recombination and repair 3 2 1
Similarly, using a multi-pistil and male sterility mutant dms, we demonstrated wheat homologs of MADS-box MFO1 (E class), FUL2 (A class), FUL3 (A class), PAP2 (E class) etc, AP2 (A class), FL and DL (CRC) genes played determining roles at early stages of floral organ differentiation in wheat. This was confirmed in some extent by temporal expression profiles of many homeotic genes during spike development of D, M and T plants (unpublished data). Combined with other researcher's results [29,39], the representative genes of all A, B, C, D, E classes have been identified in wheat. Although diversity may exist, it is clear that wheat floral organ differentiation also be determined following ABCDE model.
Other gene expression regulators play important roles in wheat floral organ development. In addition to transcription factors of typical homeotic genes described above, we identified some other DNA replication, transcription and translation factor genes significantly highly expressed in differentiating spikes (Table 3). According to GO annotation, they are involved in plant floral organ development. For example, E3 ubiquitin-protein ligase, mostly expressed in inflorescence. AtPUB4, a U-box/ARM repeat-containing E3 ubiquitin ligase, is a Transcriptome Analysis for Abnormal Spike Development of dms novel player in male fertility in Arabidopsis. Loss of AtPUB4 function causes hypertrophic growth of the tapetum layer. The Atpub4 mutation also leads to incomplete degeneration of the tapetal cells and strikingly abnormal exine structures of pollen grains [58]. E3 ubiquitin ligase participates in CpG methylation-dependent transcriptional regulation and epigenetic transcriptional silencing [59]. BTB/POZ (Bric-a-brac, Tramtrack, and Broad Complex/Pox virus and Zinc finger) domain-containing protein NPY1 (= ENP, MAB4) may act as a substrate-specific adapter of an E3 ubiquitin-protein ligase complex (CUL3-RBX1-BTB) which mediates the ubiquitination and subsequent proteasomal degradation of target proteins. It expresses in the inflorescence meristem and at various stages of floral development. Disruption of NPY1 does not cause obvious defects in organogenesis, but npy1 npy3 npy5 triple mutants failed to make flower primordia, NPY genes define one key step in a pathway that controls YUCCA1 (YUC)-mediated organogenesis in Arabidopsis  Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table). The homologs of these genes (E3 ubiquitin ligase, BTB domain-containing NPY1, FVE etc) identified in this study suggested that except homeotic genes, there were many other 'gene expression' regulators and protein interaction genes involved in wheat floral organ development.
Phytohormone and other signaling pathways influence wheat spike differentiation. Phytohormone is a kind of very important regulator for plant floral organ development. For example, GIBBERELLIC ACID REQUIRING 1 (GA1) was associated with the length of petals, stamens, and to a lesser extent style-stigma length [64]. AP2 is a typical homeotic gene. We identified two AP2-like genes, one auxin response factor 14-like gene highly expressed in differentiating spikes. The relationship between stem development and GA metabolism of the wheat mutant dms was studied by measuring GA concentrations in the upper internodes of dms using High Performance Liquid Chromatography (HPLC), and investigating the effects of exogenous GA 3 on plant height in field and coleoptile length in lab (the data will be published elsewhere). It confirmed the GA signal pathway in D plants was different from that of T plants. These data suggested that phytohormone signaling was important when spike development in wheat also.
We identified two LRR receptor-like serine / threonine-protein kinase (RLK) genes highly expressed in differentiating spikes (Table 3). LRR-RLK is a large group important signalling molecules involved in various biological processes. In Arabidopsis, LRR-RLK ERECTA (= ER, QRP1, QRS1, TE1) strongly expressed in organ primordia and immature organs but weakly in mature organs. It was observed in shoot apical meristems (SAM) at low levels during the vegetative growth with an increase at the transition to the reproductive growth phase. At the reproductive stage, localized in the young developing flowers, expressed in inflorescence meristem and is up-regulated during flower initiation and formation of flower organs. Receptor kinase that, together with ERL1 and ERL2 (two functional paralogs of ERECTA), regulates aerial architecture, including inflorescence. Loss of the entire ERECTA family genes led to striking dwarfism, reduced lateral organ size and abnormal flower development, including defects in petal polar expansion, carpel elongation, anther and ovule differentiation. These defects are due to severely reduced cell proliferation [65]. The LRR-RLK genes in ERECTA family, as well as the zinc finger protein genes identified in this study may be involved in the GRNs of wheat floral organ differentiation.

Molecular basis for abnormal spike differentiation of mutant dms
In this study, we found mutation significantly affected the expressions of homeotic genes, many transcription factors, signaling related genes, phytohormone response genes, cell division related genes, carbohydrate transport and metabolism, protein synthesis and photosynthesis related genes (Table 4; S3 Table). These were the molecular basis for abnormal spike differentiation of the mutant dms.
Two class A homeotic genes were suppressed in D. Homeotic genes play important roles in wheat transition from vegetative to floral organ growth [66]. Noticeably, the TaAP1-2 was highly expressed in spikes of T, but no detected in D. The TaAP2 was significantly highly expressed in spikes of T than in D (Fig 5). Both of AP1 and AP2 were key floral organ identity genes of class A [36,67]. Murai [29] postulated that the mutual suppression between class A and C genes existed in wheat. Class A genes determine palea and lodicule development, class C genes determine the steman and pistils. The expressions of class A genes TaAP1-2 and TaAP2 were inhibited that possibly activated the class C genes, and resulted in forming additional pistils.
Other transcription factors. Most biological processes in a eukaryotic cell or organism are finely controlled at the transcriptional level by transcription factors [68]. Genome-wide studies have revealed that many of the genes that are regulated by the floral organ identity factors encode transcriptional regulators. In Arabidopsis the percentage of transcription factorcoding genes identified in microarray experiments on floral tissues is 5.5%~26% of DEGs [40]. In this study, transcription factor-coding genes consisted of a big group of DEGs in mutant dms (23 unigenes, 16.08% of DEGs; S3 Table). The transcription factors are complex. The following were some of them.
BTB/POZ domain scaffold proteins are present in many organisms. BTB-containing proteins are numerous and control cellular processes that range from actin dynamics to cell-cycle regulation. In spite of their high sequence divergence and apparently unrelated functions, many BTB-containing proteins have at least one shared role: the recruitment of degradation targets to E3 ubiquitin ligase complexes [69]. In Arabidopsis they form a family of eighty proteins that have been implied in a variety of signaling pathways. The BTB and TAZ (transcriptional adapter zinc finger) domain protein 1 (BT1) is interactor and regulator of the PINOID (PID) kinase. BT1 has a typical domain structure that is only observed in land plants. The Arabidopsis BT family, which comprises five members, named BT1 to BT5, of which at least four are encoded by auxin responsive genes. There is considerable functional redundancy among the BT members, and that BT function is crucial for early steps in both male and female gametophyte development. BT proteins, as regulators of PID activity, are likely part of a feed-back loop that affects auxin distribution [70]. The expresstion of a BTB/POZ and TAZ domain-containing protein 2-like gene (T2-47003) in spikes of D was significantly suppressed at early stage, but highly expressed at later stage ( Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table). It may be interacted with auxin signal pathway as most BT genes in Arabidopsis.
The ASYMMETRIC LEAVES2 (AS2) gene is required for the generation of the flat and symmetrical shape of the leaf lamina in Arabidopsis. AS2 encodes a plant-specific protein with an AS2/LATERAL ORGAN BOUNDARIES (AS2/LOB) domain that includes a cysteine repeat, a conserved single glycine residue and a leucine-zipper-like sequence in its amino-terminal half. Some members of the AS2/LOB family might have redundant function(s). The AS2/LOB domain of AS2 cannot be functionally replaced by those of other members of the family, and that only a few dissimilarities among respective amino acid residues of the AS2/LOB domain of AS2 and those of other members are important for the specific functions of AS2 [71]. The expression of a LOB domain-containing protein 37-like gene was significantly suppressed in spikes of D. The functional annotation implied it could regulate transcription.
Four LOC284861-like genes (T2-16245, T1-64318, T4-61132 and T2-31732) were significantly suppressed in D, unigene T2-16245 almost only expressed in T, suggested a very important role in wheat floral organ differentiation. Uncharacterized protein LOC284861 is a new protein of Homo sapiens. Blastn indicated that the unigene T2-16245 was moderately similar to a sequence (emb/FN564432.1) in wheat 3B-specific BAC library by 70% identity over 354 bp region. Such, unigene T2-16245 was a novel gene in wheat. The functional annotation implied it could regulate transcription, and its function needs to be explored in future.
Phytohormone signalling. Abscisic acid (ABA), auxin, ethylene, and gibberellin signaling related genes were affected in mutant, the later two might be more important for floral organ development, because AP2 [40,48,72], and GA1 [60] were certainly function in plant floral organ development. The expressions of two auxin-repressed 12.5 kDa protein genes were upregulated in D, indicated that the concentration of auxin in the spikes of D might be lower, and auxin pathway was repressed. The dioxygenase for auxin oxidation (DAO) gene, encoding a putative 2-oxoglutarate-dependent-Fe (II) dioxygenase, is essential for anther dehiscence, pollen fertility, and seed initiation in rice [73]. This indicated the important role of auxin signaling in floral organ development. Repressed auxin signal pathway would decrease cell division and proliferation, thereby affected growth. An ABA inducible protein PHV A (cold acclimation protein WCOR615) was suppressed, one gibberellin receptor GID1, two gibberellin stimulated transcript (GAST1; = Gibberellin Stimulated-Like proteins 1, GSL1, = snakin-1) genes were upregulated in D. In tomato (Lycopersicon esculentum Mill.), gibberellic acid (GA3) increases the abundance of GAST1 RNA while treatment with ABA blocks this effect [74]. In Arabidopsis, GASA5 suppresses GA responses, GA induces or suppresses the transcription of different GAST1-like genes and the encoded proteins regulate the redox status of specific components to promote or suppress specific GA responses [75]. Potato (Solanum tuberosum L.) GSL1 and GSL2 genes are very highly conserved. Potato transformation with antisense knock-down expression cassettes failed to recover viable plants (lethality; similar to D genotype in this study in some extent). GSL protein family plays a role in several aspects of plant development [76]. Gibberellin pathway plays an important role in wheat transition from vegetative to floral growth [66]. The relationship among gibberellin, VRN1, FT, SUPPRESSOR OF OVEREXPRES-SION OF CONSTANS1-1 and FL in determining wheat flowering has been studied. GA and FL play critical roles during early stages of wheat flower development [77]. These cases indicate that the balance of hormone signalings is pivotal for floral organ differentiation. In this study, the expressions of GID1 and GAST1-like genes were significantly activated, indicated that the gibberellin signaling pathway was significantly altered. Generally, the signal pathways of ABA, ethylene and auxin were inhibited in D, but gibberellin was activated, thereby resulted in hormonal imbalance. We have primarily confirmed that GA signaling in D was different from T and M plants (the data will be published elsewhere). The phytohormone regulation in D was very complex, but hormonal imbalance should be one aspect.
Carbohydrate transport and metabolism were significantly suppressed in D. Combined with photosynthesis, energy production and conversion, this was the biggest group (27 unigene, 18.88% of DEGs, S3 Table) of DEGs between T and D. Carbohydrate not only provides the energy for various biological processes, but also provides materials for cell wall construction. Carbohydrate transport and metabolism were significantly suppressed in D. For example, the expressions of putative alpha-L-fucosidase 1 and mitochondrial outer membrane protein porin 1-like genes were almost inhibited completely ( Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table). Alpha-Lfucosidase 1 is an element of carbohydrate metabolism. The transcripts of alpha-L-fucosidase 1-like gene were abundance. Its expression was sharply suppressed in D ( Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table), thus we considered it as an important gene influencing wheat floral organ differentiation. The cell wall biogenesis / degradation were more active and important for differentiating spikes [62,63]. The male sterility induced by a gametocide SQ-1 might be associated with imbalance of energy metabolism [78]. The wheat line 3558 has two type tillers of inhibited and normal. Development was blocked or delayed at the elongation stage in the inhibited tillers of 3558. Weakened energy metabolism might be one reason that the inhibited tillers could not joint and develop into spikes [79]. The plants of D were very short (<30 cm) and completely sterile. The suppressed carbohydrate transport and metabolism, photosynthesis, energy production and conversion, should have played a basic role for morphological abnormalities of D.
Significant DEGs need investigation in detail. Some genes significantly differently expressed in spikes of D and T had been picked out through mRNA sequencing and QRT-PCR analysis. The expressions of an unknown gene in 1Ds prolamin gene locus (retrotransposon protein), mitochondrial outer membrane protein porin 1, replication factor A protein 1, uncharacterized protein LOC284861, transcription-repair-coupling factor, lysine-rich arabinogalactan protein 19, ABA-inducible protein PHV A1, putative alpha-L-fucosidase 1, uncharacterized protein DKFZp434B061, splicing factor 3A subunit 1-like, setaria italica titin-like, hypothetical protein F775_20975, G-protein coupled receptor 114-like and others without any similar known gene or sequences in GenBank were almost completely inhibited in D. The expressions of a amidase 1, bowman-birk type wound-induced proteinase inhibitor WIP1, homeobox-leucine zipper protein HOX24, serine/arginine repetitive matrix protein 2, metallothionein-like protein 1, LEA protein 12, isoflavone 2'-hydroxylase, probable WRKY transcription factor 70-like gene, and others without any similar known gene or sequences in GenBank were significantly up-regulated in D (Table 4 Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table). Although some of them were absolutely novel genes in wheat without any putative functions, the entirely opposite expression patterns between T and D implied the important roles they played in wheat floral organ development, and they were the major causal factors for spike abnormality of D (Fig 7A-7L Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table). For example, the unknown gene T2-14438, similar to a hypothetical protein F775_20975 of Aegilops tauschii, whose expression was significantly suppressed in spikes of D ( Table 4. The actin gene was used as internal control. The number 1-16 indicate the sampling dates, and the details see S1 Table), that implied an important role in spike development. The functions of these genes should be investigated in detail.

The molecular bases for abnormal spikes of D
Gene function loss can result in plant dwarf, such as an extreme dwarf phenotype of the GAsensitive mutant of sunflower, dwarf2 [80]. Pistillody in alloplasmic wheats is caused by alterations to the expression pattern of class-B MADS-box genes, DL, AP1, AP2, FL etc [29,52,53]. Rice hybrid sterility is the most common form of postzygotic reproductive isolation in plants. The hybrid sterility genes is in very different functional categories, aspartic protease, small ubiquitin-like modifier E3 ligase-like protein, F-box protein, transcription factor, the catalytic subunit of a Na + /K + ATPase, polypeptide with a single MADF DNA-binding domain near its C terminus end [81]. These cases indicated that one gene mutation is sufficient to result in deleterious function and causing abnormal phenotype, even if a single amino acid exchange in the coding peptide. One phenotype may caused by different factors. Interestingly, the key genes involved in dwarfism, pistillody and sterility reported above were also important genes in wheat floral organ differentiation identified in this study, such as MADS-box, DL, AP1, AP2, FL [51,66], E3 ligase-like protein, transcription factor [81], energy metabolism [79].
Although the genetic changes of mutant dms has not been identified yet, significant DEGs between T and D were carefully studied. The molecular base for abnormal spike differentiation of dms was proposed here (Fig 8). TaFL gene involves in specifying floral organ. A YABBY gene TaDL (DL) involves in specifying the pistil identity. The expressions of class A genes TaAP1-2 and TaAP2 were inhibited in D. The DNA replication was affected, most transcription factors were repressed, but some are activated. The biological processes of translation, photosynthesis, glycolysis, cell division are significantly repressed in D, signal transduction mechanism is activated. Furthermore, some unknown genes are suppressed or activated. These data constitute a frame work for transcriptome profile in D. Certainly it needs to be confirmed step by step in future.

Conclusions
The wheat spike transcriptomes of D and T at early stage of floral organ differentiation were carefully analyzed. We identified 419 genes related to wheat spike differentiation, including typical homeotic genes, DNA replication, transcription and translation factors, carbohydrate transport and metabolism related genes. The identified homeotic genes were MADS-box (MFO1, FUL2, FUL3, PAP2, etc), AP2, FL and DL-like genes that played important roles in wheat floral organ differentiation. The expression levels of 143 DEGs in differentiating spikes were significantly different, the expression of AP1-like gene (TaAP1-2) was inhibited in D, and AP2-like gene (TaAP2) was significantly suppressed. Among the DEGs, 16.08% genes encoding transcription factors; 18.88% genes were of carbohydrate transport and metabolism, photosynthesis, energy production and conversion; more than 21% were novel and unknown function genes. Most biological processes such as transcription, translation, cell division, photosynthesis, carbohydrate transport and metabolism, energy production and conversion were repressed in D. Particularly most metabolism pathways of photosynthesis and carbohydrate metabolism were suppressed. Inhibition of homeotic genes, suppression of most biological processes and energy deficiency were considered as the major causal factors for abnormal spike development of D.
Supporting Information