Quantitative Changes in the Transcription of Phytohormone-Related Genes: Some Transcription Factors Are Major Causes of the Wheat Mutant dmc Not Tillering

Tiller number is an important agronomic trait for grain yield of wheat (Triticum aestivum L.). A dwarf-monoculm wheat mutant (dmc) was obtained from cultivar Guomai 301 (wild type, WT). Here, we explored the molecular basis for the restrained tiller development of the mutant dmc. Two bulked samples of the mutant dmc (T1, T2 and T3) and WT (T4, T5 and T6) with three biological replicates were comparatively analyzed at the transcriptional level by bulked RNA sequencing (RNA-Seq). In total, 68.8 Gb data and 463 million reads were generated, 80% of which were mapped to the wheat reference genome of Chinese Spring. A total of 4904 differentially expressed genes (DEGs) were identified between the mutant dmc and WT. DEGs and their related major biological functions were characterized based on GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) categories. These results were confirmed by quantitatively analyzing the expression profiles of twelve selected DEGs via real-time qRT-PCR. The down-regulated gene expressions related to phytohormone syntheses of auxin, zeatin, cytokinin and some transcription factor (TF) families of TALE, and WOX might be the major causes of the mutant dmc, not tillering. Our work provides a foundation for subsequent tiller development research in the future.


Introduction
Wheat (Triticum aestivum L.) is one of the most important food crops in the world. Since tiller number is an important agronomic trait for grain yield [1,2], it has always been one of the key traits to select in breeding programs. Generally, low and high tillering wheat mutants do not have very high grain yield. Moreover, tiller number as well as grain number and weight affect yield.
In comparison to the model species Arabidopsis, only a few genes affecting tiller initiation and outgrowth have been cloned and described in crops. For example, the rice gene MONOCULM 1 (MOC1) encodes a GRAS (GIBBERELLIC-ACID INSENSITIVE (GAI), REPRESSOR of GAI (RGA) and SCARECROW (SCR)) domain-containing protein that affects the initiation and outgrowth of axillary meristems [27,28]. Mutant lines that do not have a functional moc1 gene exhibit a severe reduction in tiller number. The rice gene Tillering and Dwarf 1 (TAD1), which encodes a multi-subunit E3 ligase, regulates rice tillering by degenerating MOC1 [29]. The rice TEOSINTE BRANCHED 1 (OsTB1) gene encodes a TCP domain protein and negatively regulates lateral branching [30]. In barley, the INTERMEDIUM-C (INT-C) gene is an ortholog of maize and rice TB1 genes which has an effect on seedling tiller number [31]. The barley uzu gene encodes a putative brassinosteroid (BR) receptor HvBRI1 and regulates tiller number [32]. In bread wheat, overexpression of tae-miR156 significantly affects tillering, probably by regulating a group of SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) genes; furthermore, miR156-TaSPLs and strigolactone signaling pathways might have a potential association with tillering [33]. Overexpression of the maize tb1 gene in wheat results in reduced tillering [34]. The molecular mechanism of wheat tillering remains largely unknown.
At present, the RNA-Sequencing technique has been widely used to identify differentially expressed genes (DEGs) among various biological samples, so as to explore the possible mechanisms leading to various morphologies [35,36]. Transcriptome analyses regarding tiller developments in switchgrass (Panicum virgatum L.) and sorghum (Sorghum bicolor L.) have been reported [37][38][39]. In this study, a wheat mutant that does not develop tillers was obtained from the ethyl methanesulfonate (EMS)-treated wheat cultivar Guomai 301 (WT), which was designated as dmc (dwarf-monoculm). We characterized the mutant dmc and WT at the transcriptome level by RNA-sequencing technology.

Morphology of the dmc Mutant
The mutant dmc was obtained from wheat cultivar "Guomai 301" treated with EMS. Mutant dmc basically does not tiller, and only has a main stem (Table 1, Figure 1A). Some individuals occasionally have a small tiller number. The plant height of the mutant dmc was significantly lower than that of the WT, i.e., 48.00 cm, which was 74.53% of the WT height (Table 1, Figure 1A). Both the spike length and the seed length of the mutant dmc were shorter than those of the WT (Table 1, Figure 1B-D). The diameters of the internodes of the mutant dmc were reduced ( Figure 1E). The internode number of the mutant dmc was four; however, that of the WT was five (Table 1, Figure 1F). Additionally, we continuously observed and compared the tillers of the mutant dmc and WT from the early tillering stage. The mutant dmc had no or one tiller at the early tillering stage ( Figure 2C), and it had no more tillers during the middle tillering period. Its tillering was significantly inhibited ( Figure 2D). Generally, the mutant dmc had only one main stem at the jointing stage ( Figure 2E). The tiller primordia ( Figure S1) were used as samples for transcriptomic analysis in this study. Tiller primordium samples were dissected at the three-leaf stage to four-leaf stage. Two super bulk samples of the mutant dmc (T1, T2, and T3) and WT (T4, T5, and T6) with three biological replicates were prepared. Each bulk sample included more than ten independent individuals. tiller primordia ( Figure S1) were used as samples for transcriptomic analysis in this study. Tiller primordium samples were dissected at the three-leaf stage to four-leaf stage . Two super bulk  samples of the mutant dmc (T1, T2, and T3) and WT (T4, T5, and T6) with three biological replicates were prepared. Each bulk sample included more than ten independent individuals.    tiller primordia ( Figure S1) were used as samples for transcriptomic analysis in this study. Tiller primordium samples were dissected at the three-leaf stage to four-leaf stage. Two super bulk samples of the mutant dmc (T1, T2, and T3) and WT (T4, T5, and T6) with three biological replicates were prepared. Each bulk sample included more than ten independent individuals.

Genetic Diversity between the WT and Mutant dmc
To validate whether the dmc was derived from random open pollination in the field, Polymerase Chain Reaction (PCR) amplification of the WT and mutant dmc were carried out with 431 primer pairs of wheat Simple Sequence Repeats (SSR) markers evenly distributed on 21 wheat chromosomes. There was no polymorphic SSR locus between WT and mutant dmc ( Figure S2), which demonstrated that their genetic backgrounds were highly consistent. This demonstrated that dmc was a real mutant derived from WT.

RNA Sequencing Data
Six libraries were analyzed using RNA sequencing (mutant dmc: T1, T2, T3; WT: T4, T5, T6). We obtained a total of 68.8 Gb clean bases and about 463 million reads (single-end reads) ( Table S1). The GC contents of the six libraries were 55.84-56.85%, and the average Q30 percentage was 89.25% (Table S1). Most transcripts were 100 to 300 bp in length ( Figure S3), and the biological replicates were highly consistent ( Figure S4). The reads were compared with the T. aestivum reference genome. There were 374 million (80.74%) reads mapped to the reference genome (Table S2). The average percentage of unique mapped reads was more than 71% (Table S2). It was clear that a high-quality transcriptome data set was obtained.

Annotation and Functional Classification of the Unigenes
In total, 113,619 unique genes (unigenes) were obtained from the six libraries, and 109,685 unigenes were annotated by BLAST in several databases (Table S3). Furthermore, 10,080 new genes among the 113,619 unigenes were obtained, and 7075 new genes were annotated (Table S4). Functional classification in GO showed that the DEGs were classified into cellular component, molecular function, biological process, and many subcategories. Within the cellular component, molecular function and biological process categories, the most represented DEGs were classified as "cell", "organelle" and "cell part", "catalytic activity", "binding"; "metabolic process", "cellular process" and "single-organism process" ( Figure S5).
According to homology unigene distribution of various species in the Nr database ( Figure S6), the order from highest to lowest abundance was Aegilops tauschii, H. vulgare, Triticum urartu, T. aestivum, Brachypodium distachyon, O. sativa, Z. mays, Setaria italica and S. bicolor.

DEGs between the WT and Mutant dmc
To investigate the gene expression profile variation, a total of 4904 differentially expressed genes (DEGs) between the WT and mutant dmc were identified. Among them, 1506 were expressed at a low level, and 3398 were highly expressed in dmc compared to the WT ( Figure 3A). The expression levels of DEGs in dmc and WT are shown as volcano plots ( Figure 3B). The expression patterns of DEGs were hierarchically clustered ( Figure 3C). The result showed there were significant DEGs, which may be the major genes related to wheat tillering.

DEGs between the WT and Mutant dmc
To investigate the gene expression profile variation, a total of 4904 differentially expressed genes (DEGs) between the WT and mutant dmc were identified. Among them, 1506 were expressed at a low level, and 3398 were highly expressed in dmc compared to the WT ( Figure 3A). The expression levels of DEGs in dmc and WT are shown as volcano plots ( Figure 3B). The expression patterns of DEGs were hierarchically clustered ( Figure 3C). The result showed there were significant DEGs, which may be the major genes related to wheat tillering.
In order to further explore the key genes, significant DEGs (Log2FC ≥3 or ≤−2) were screened between the mutant dmc and WT (Table S5, Table S6). These DEGs were classified into several groups, such as transcription factor, signal transduction mechanism, and carbohydrate metabolism, etc. Chloroplastic D-3-phosphoglycerate dehydrogenase 2 was the most significant highly expressed gene in the mutant dmc; its Log2FC value was 11.43. Histone H2B.1 was the gene that was most significantly expressed at a low level in the mutant dmc; its Log2FC value was 12.28.

Functional Classification of the DEGs in GO
To further explore wheat tiller-related biological pathways or processes, 3991 DEGs were classified into 54 subcategories in the GO database ( Figure 4). According to the percentage of DEGs in all genes, the significant subcategories were biological phase (two DEGs), biological adhesion In order to further explore the key genes, significant DEGs (Log 2 FC ≥3 or ≤−2) were screened between the mutant dmc and WT (Table S5, Table S6). These DEGs were classified into several groups, such as transcription factor, signal transduction mechanism, and carbohydrate metabolism, etc. Chloroplastic D-3-phosphoglycerate dehydrogenase 2 was the most significant highly expressed gene in the mutant dmc; its Log 2 FC value was 11.43. Histone H2B.1 was the gene that was most significantly expressed at a low level in the mutant dmc; its Log 2 FC value was 12.28.

Functional Classification of the DEGs in GO
To further explore wheat tiller-related biological pathways or processes, 3991 DEGs were classified into 54 subcategories in the GO database ( Figure 4). According to the percentage of DEGs in all genes, the significant subcategories were biological phase (two DEGs), biological adhesion (three DEGs) (belonging to biological process), the membrane part (599 DEGs), the extracellular region (268 DEGs) and membrane-enclosed lumen (17 DEGs) (belonging to the cellular component), structural molecule activity (31 DEGs), molecular transducer activity (24 DEGs), receptor (13 DEGs), and protein binding transcription factor (two DEGs) (belonging to molecular function).

The DEGs Involved in Phytohormone Metabolisms in dmc
The DEGs involved in plant hormone metabolisms were enriched in the mutant dmc, and a total of 83 DEGs were obtained ( Figure 7, Table S9). Among them, 52 DEGs were expressed at a low level, and 31 DEGs were highly expressed in the mutant dmc compared to the WT. The DEGs of cytokinin hydroxylase, cytokinin dehydrogenase and cytokinin phosphoribohydrolase were all expressed at a low level ( Figure 7C, Table S9); three DEGs of abscisic acid 8 -hydroxylase and abscisic stress-response protein were expressed at a high or low level ( Figure 7B, Table S9). Two DEGs of gibberellin 20 oxidase and gibberellin 2-beta-dioxygenase were expressed at a low level ( Figure 7D, Table S9); 19 DEGs involved in auxin metabolism were expressed at a low level ( Figure 7A, Table S9), 17 DEGs involved in ethylene metabolism were expressed at a low level in the mutant dmc compared to the WT (Table S9). Abscisic stress-responding protein 3 (AA0320460) and auxin-induced protein X15 (AA1998680) were the DEGs that were most significantly expressed at high and low levels, respectively, in the mutant dmc compared to the WT (Table S9).

The DEGs Involved in Phytohormone Metabolisms in dmc
The DEGs involved in plant hormone metabolisms were enriched in the mutant dmc, and a total of 83 DEGs were obtained (Figure 7, Table S9). Among them, 52 DEGs were expressed at a low level, and 31 DEGs were highly expressed in the mutant dmc compared to the WT. The DEGs of cytokinin hydroxylase, cytokinin dehydrogenase and cytokinin phosphoribohydrolase were all expressed at a low level ( Figure 7C, Table S9); three DEGs of abscisic acid 8′-hydroxylase and abscisic stress-response protein were expressed at a high or low level ( Figure 7B, Table S9). Two DEGs of gibberellin 20 oxidase and gibberellin 2-beta-dioxygenase were expressed at a low level ( Figure 7D, Table S9); 19 DEGs involved in auxin metabolism were expressed at a low level ( Figure 7A, Table  S9), 17 DEGs involved in ethylene metabolism were expressed at a low level in the mutant dmc compared to the WT (Table S9). Abscisic stress-responding protein 3 (AA0320460) and auxin-induced protein X15 (AA1998680) were the DEGs that were most significantly expressed at high and low levels, respectively, in the mutant dmc compared to the WT (Table S9).

The DEGs Related to Carbohydrate Metabolism in the Mutant dmc
A total of 303 DEGs related to carbohydrate metabolism were obtained and they were divided into 19 groups ( Figure 8A, Table S10). Among these, the DEGs related to the metabolism of glucose, fructose, starch, fucose, ribulose, mannose, xlanase, glucoside, galactoside, galactinol and galacturonokinase were all highly expressed in dmc ( Figure 8A, Table S10). The DEGs related to metabolisms of glucan, xyloglucan, sucrose, sugar, amylase, glucosidase, glycosyltransferase, and polygalacturonase were simultaneously expressed at high or low levels in dmc. Furthermore, the percentages of DEGs expressed at low levels were lower than those of the highly-expressed DEGs in all carbohydrate metabolisms ( Figure 8A, Table S10).

The DEGs Related to Carbohydrate Metabolism in the Mutant dmc
A total of 303 DEGs related to carbohydrate metabolism were obtained and they were divided into 19 groups ( Figure 8A, Table S10). Among these, the DEGs related to the metabolism of glucose, fructose, starch, fucose, ribulose, mannose, xlanase, glucoside, galactoside, galactinol and galacturonokinase were all highly expressed in dmc ( Figure 8A, Table S10). The DEGs related to metabolisms of glucan, xyloglucan, sucrose, sugar, amylase, glucosidase, glycosyltransferase, and polygalacturonase were simultaneously expressed at high or low levels in dmc. Furthermore, the percentages of DEGs expressed at low levels were lower than those of the highly-expressed DEGs in all carbohydrate metabolisms ( Figure 8A, Table S10).

Transcription Factor Type DEGs in the Mutant dmc
A total of 35 significant DEGs of transcription factors (Log 2 FC ≥3 or ≤−2) were obtained in dmc ( Figure 8B, Table S5, Table S6), and they belonged to zinc finger protein, heat stress transcription factor, NAC domain-containing protein, WRKY transcription factor, ethylene-responsive transcription factor and MADS-box protein families, etc. Among these DEGs regulated by transcription factors ( Figure 8B, Table S11), an SPX domain-containing membrane protein Os02g45520 (AA1610920) was the most significantly highly expressed, and probable WRKY transcription factor 12 (AA0534130) was the DEG that was most significantly expressed at a low level in the mutant dmc compared to the WT. The DEGs belonging to families B3, DBB, Dof, GRF, LBD, LFY, MADS, SRS, and WOX were only expressed at low levels; the DEGs belonging to families GATA, GRA, and NF-YC were only expressed at high levels in the mutant dmc compared to the WT ( Figure 8B, Table S11). The most enriched transcription factor family was ERF, followed by MYB, bHLH and HSF ( Figure 8B, Table S11).

Transcription Factor Type DEGs in the Mutant dmc
A total of 35 significant DEGs of transcription factors (Log2FC ≥3 or ≤−2) were obtained in dmc ( Figure 8B, Table S5, Table S6), and they belonged to zinc finger protein, heat stress transcription factor, NAC domain-containing protein, WRKY transcription factor, ethylene-responsive transcription factor and MADS-box protein families, etc. Among these DEGs regulated by transcription factors ( Figure 8B, Table S11), an SPX domain-containing membrane protein Os02g45520 (AA1610920) was the most significantly highly expressed, and probable WRKY transcription factor 12 (AA0534130) was the DEG that was most significantly expressed at a low level in the mutant dmc compared to the WT. The DEGs belonging to families B3, DBB, Dof, GRF, LBD, LFY, MADS, SRS, and WOX were only expressed at low levels; the DEGs belonging to families GATA, GRA, and NF-YC were only expressed at high levels in the mutant dmc compared to the WT ( Figure 8B, Table S11). The most enriched transcription factor family was ERF, followed by MYB, bHLH and HSF ( Figure 8B, Table S11).

Expression Profiles of Twelve Genes in Wheat Tiller Primordia
To demonstrate the reliability of the sequencing results, we selected twelve of the significant DEGs to perform real-time qRT-PCR. Linear regression analysis showed that the results of the transcriptome analysis were reliable (Figure 9, Figure S7).   A histone H2B.1 gene (CS42-U-AA2080400, Figure 9A), a PGR5-like protein 1A gene (CS42-1BS-AA0175910, Figure 9B) and a WRKY transcription factor 12 gene (CS42-2DL-AA0534130, Figure 9C) were significantly suppressed in tiller primordia of dmc, and these three genes were also significantly suppressed in the leaves of dmc compared to the WT. Mutation significantly suppressed the expression of an ATP-dependent zinc metalloprotease FTSH 5 gene (CS42-5BL-AA1349930, Figure 9D), a BOI-related E3 ubiquitin-protein ligase 3 gene (CS42-U-AA2143280, Figure 9E), an acid phosphatase 1 gene (CS42-2DS-AA0603970, Figure 9F) and a bidirectional sugar transporter SWEET3a gene (CS42-1BS-AA0159600, Figure 9G) in tiller primordia of dmc; these four genes were specifically expressed in tiller primordia compared to leaves except the bidirectional sugar transporter SWEET3a gene (CS42-1BS-AA0159600, Figure 9G). The expression of an arginine decarboxylase gene (CS42-3DL-AA0834780, Figure 9H); a putative F-box/FBD/LRR-repeat protein At5g22670 gene (CS42-1BS-AA0166120, Figure 9I); a SPX domain-containing membrane protein Os02g45520 gene (CS42-6BL-AA1610920, Figure 9J); and a GDSL esterase/lipase At1g28600 gene (CS42-2AL-AA0292740, Figure 9K) were significantly activated in tiller primordia of dmc compared to WT, and their expression profiles in leaves of dmc were also similar. The expression of an abscisic stress-response protein 1 gene (Wheat-newGene-4506, Figure 9L) was irregular. The twelve genes were either only just expressed or completely suppressed in tiller primordia of dmc. The results suggested that these genes played important roles in the tiller development of wheat.

Lack of Vitality in the Mutant dmc
The morphology of the mutant dmc was dwarfish, with yellowish cotyledons and almost no tillers, which demonstrated its lack of vitality. The phenotype of the mutant dmc was novel. In wheat, the leaves of the tiller inhibition mutants tin3 and ftin were much darker than the WT plants [1,4]; however, the leaves of the mutant dmc were yellowish. Both spikes and seed were much larger in the mutant tin3 [1] compared with the mutant dmc. Normally, there is a highly negative correlation between tiller number and plant height in rice and wheat [6]. In wheat, NAUH167, a high-tillering mutant was dwarfish compared to the WT [6], the plant height of the tiller inhibition mutant ftin was slightly higher than that of the WT [4], the phenotype of the mutant dmc was not the same as NAUH167 or ftin, and the mutant dmc was not only dwarfish but also had almost no tillers. The mutant dmc had fewer tiller primordia. Most of the primordia could not grow or stopped developing, which resulted in no tillers. Though most DEGs related to photosynthesis were highly expressed in dmc, this occurred in contrast to its lower vitality and biomass.

Auxin and Cytokinin Metabolisms Were Suppressed in dmc
Phytohormones play essential roles in plant growth and development throughout the vegetative to reproductive stages [40][41][42][43]. Tillering is a key developmental event in wheat, which needs a balanced phytohormone metabolism to maintain normal tillering [44][45][46].
In Arabidopsis, auxin can inhibit bud outgrowth in the highly branched shoot mutant axr1, and the effect of auxin on the wild type is more obvious [47]. In rice, a knock-down mutant of OsIAA6 produces more tillers due to the regulation of the auxin transporter OsPIN1 [48]. However, a further study found that auxin does not enter the lateral buds [49,50], and this phenomenon can be supported by cytokinin and strigolactone [50,51]. In Arabidopsis, cytokinin acts to overcome auxin-mediated apical dominance, allowing buds to evade apical dominance [52]. In this study, we obtained 29 DEGs related to auxin metabolism, and 19 DEGs were expressed at a low level ( Figure 6, Table S9): the level of auxin in the mutant dmc might be low; however, this might not be consistent with the premise that auxin inhibits bud outgrowth. The KEGG pathway of zeatin biosynthesis (ko00908; Figure S8) was suppressed, and the suppressed differentially expressed zeatin biosynthesis genes were adenylate isopentenyltransferase 1 (chloroplastic), three cytokinin hydroxylases, and three cytokinin dehydrogenases ( Figure S8, Table S9). Furthermore, the plant hormone signal transduction genes (ko04075; Figure S9) in KEGG showed that the expression of zeatin biosynthesis can affect cytokinin biosynthesis: the expression of downstream signal transduction components and the two-component response regulator (ARR) were down-regulated ( Figure S9). This result is in agreement with a previous study on rice, which indicated that a high level of cytokinin can increase the tiller number [53].
Interestingly, overexpressing gibberellin 2-oxidases (GA2ox) can produce more tillers in rice [54]. Gibberellin localization in vascular tissue is required to control auxin-dependent bud outgrowth in hybrid aspen (Pterocarya stenoptera) [55]. Gibberellin is also required for cytokinin-mediated axillary bud outgrowth in Jatropha curcas [56]. In this study, two DEGs encoding gibberellin 20 oxidase 2, and DEGs encoding gibberellin 2-beta-dioxygenase 1 and gibberellin 3-beta-dioxygenase 2-1, which were expressed at low levels, were obtained ( Figure 6, Table S9), and the possible high level of gibberellin could be consistent with the mutant dmc not tillering. The hormone regulation of wheat tillers is very complex, and further research is needed in the future.
The plant height of the mutant dmc was decreased significantly. It was about 74.53% of the WT height ( Figure 1A). Many plant hormones can promote plant growth, like auxin [42], gibberellin [57,58], and cytokinin [59]. The suppressed expression of genes involved in auxin, zeatin and cytokinin metabolisms ( Figure S9) may be the key causes accounting for dwarfism of the mutant dmc.

Carbohydrate Metabolism and Phenylpropanoid Biosynthesis Were Active in dmc
Many studies have demonstrated that sugar can regulate bud growth [60][61][62]. In wheat, the phenotype of the tiller inhibition mutant (tin) was associated with precocious internode elongation; further study showed that the mutant tin transferred sucrose away from buds to elongating internodes [61]. In pea (Pisum sativum), total sucrose levels were significantly increased in node 2 buds after decapitation; the axillary buds were rapidly released in intact plants supplemented with external sucrose [60]. Sugar also affects bud outgrowth in rose (Rosa sp.) [62,63]. In this study, the enriched pathways of starch and sucrose metabolism (ko00500) in KEGG were significantly activated ( Figure 5, Figure 8A). Two DEGs related to sucrose synthase were expressed at a low level, and nine DEGs related to sucrose 1-fructosyltransferase, which were expressed at a high level, were obtained ( Figure 7, Table S10). The potentially higher activities of the corresponding enzymes may decrease the content of small molecular carbohydrates, which may inhibit wheat tillering.
Interestingly, phenylpropanoid biosynthesis (ko00940; Figure S10) was the most enhanced pathway in KEGG, which regulates the biosyntheses of different lignins such as syringyl, guaiacyl, p-hydroxy-phenyl, and 5-hydroxy-guaiacyl lignins. A total of 93 annotated DEGs belonged to the phenylpropanoid biosynthesis ( Figure S10), and among them 56 were various peroxidases (Table S10) which likely suppressed the biosyntheses of lignins, and the concentration of lignins in the mutant might be low. This result might be consistent with the dwarf phenotype of the mutant dmc. 1-aminocyclopropane-1-carboxylic acid oxidase (ACCO) catalyzes the final step of ethylene biosynthesis. In this study, seven DEGs of ACCO and seven DEGs of ACCO homologs were obtained, and they might affect the biosynthesis of ethylene in dmc.

Some Transcription Factors Are Suppressed in dmc
Transcription factors (TFs) play essential roles in plant leaf [64], flower [65], and branch [66] growth and development.
In maize (Zea mays), the Knotted1-like homeobox (KNOX) TFs up-regulates GA2ox1 [67]. In this study, four DEGs of knotted-1-like 1, three knotted-1-like 12 and three KNOX3 (TALE family) were obtained, and they were all expressed at a low level, and these might down-regulate the expression of GA2ox1 and inhibit tillering in the mutant dmc ( Figure 8, Table S11). This observation is in agreement with the observation that overexpressing GA2ox produced more tillers in rice [54]. In rice, overexpression of growth-regulating factors OsGRF3 and OsGRF10 reduces tillers. The expression of Oskn2, one rice KNOX gene, was down-regulated by overexpression of OsGRF3 [68]. In this study, one DEG homologous to TFs GRF2, GRF10, and GRF12; two DEGs homologous to GRF1 and GRF5; three DEGs homologous to GRF9 were obtained (GRF family), however, they were all expressed at a low level in dmc (Figure 8, Table S11). The relationship between GRFs expressed at a low level and the phenotype of the mutant dmc needs further study.
DWARF TILLER1, a WUSCHEL-related homeobox (WOX) TFs, is a positive regulator of tiller growth in rice [66]. In rice, the completely sterile and reduced tillering 1 mutant (srt1) was caused by a mutation in WUSCHEL (OsWUS, one member of the WOX gene family). The homeobox domain of SRT1 is essential for tillering [69]. In this study, two DEGs homologous to TFs WOX 4 (WOX family) were expressed at a low level in dmc (Figure 8, Table S11), which might affect the tillering of dmc. Furthermore, The DEGs of TF families B3, DBB, Dof, GRF, LBD, LFY, MADS and SRS were only expressed at a low level (Figure 8, Table S11); they might affect the tillering of dmc.

Signal Transduction and Photosynthesis in dmc
Signal transduction is related to many pathways and metabolic processes. In Arabidopsis, short internodes (SHI) is a suppressor of GA responses [70]. In this study, six DEGs associated with short internodes, expressed at a low level, were obtained, and they might promote GA responses in the dmc mutant that does not produce tillers. The larger percentage of DEGs related to E3 ubiquitin-protein ligase, NRT1/PTR FAMILY, F-box protein, pentatricopeptide repeat-containing protein and E3 ubiquitin-protein ligase (Table S12) might also affect the tillering and plant growth of the mutant dmc.
Photosynthesis (ko00195; Figure S11) was the third most enhanced pathway in the mutant dmc (Figure 4), and photosystem I, photosystem II, cytochrome b6/f complex, photosynthethic electron transport and F-type ATPase related to photosynthesis were all highly expressed ( Figure S11).
The enhanced expression of the photosynthetic components did not correspond with the phenotype of the mutant dmc, which requires further study.

Plant Materials and Growth Conditions
The wheat cultivar 'Guomai 301' was bred in our laboratory. Wheat cultivar "Guomai 301" has medium tillers, but a high percentage of earbearing tillers. The tillering ability of the WT conformed to the requirements of high-yield wheat breeding, which ensures a relatively dense population (600-700 spikes per hectare) and large spikes with more grains (36.1-37.4 grains per spike). In October 2012, the seeds of the WT were treated with EMS and planted at the Shangqiu Experimental Farm, Henan Province, China (34 • [71]. The samples for transcriptome sequencing analysis were prepared in 2017.

Trait Measurements, Morphology Observations and SSR Analysis
The tiller number, spike number, plant height, and internode number of the main stem, among other traits (Table 1), of dmc and the WT were observed and measured. Each sample was prepared by random selection of more than 20 individuals of dmc and the WT. The statistical tests were performed using Student's t-test, and the variation was expressed as the standard deviation (SD).
The tillers of the mutant dmc and WT were observed from the early tillering stage with an inverted microscope (Olympus 3111286), and the images were captured by a camera (Nikon Coolpix 4500).
The outer leaves and sheaths were removed to show the tiller primordia and very small tillers at the seedling stage with an anatomical needle. The large tillers at the adult plant stage were observed by the naked eye.

RNA Extraction, Library Preparation and Sequencing
Tiller primordial samples ( Figure S1) were dissected at the three-leaf stage to four-leaf stage. Two super bulk samples of the mutant dmc (T1, T2, and T3) and WT (T4, T5, and T6) with three biological replicates were prepared. Each bulk sample included more than ten independent individuals.
RNA was extracted using Trizol reagent (TransGen Biotech, Beijing, China) according to the manufacturer's protocol. RNA concentration was measured using a NanoDrop 2000 (NanoDrop Technologies, Wilmington, DE, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit on the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).
The RNA bulk samples were prepared with 1 µg RNA per individual. Sequencing libraries were constructed using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) following the manufacturer's recommendations. The preferential cDNA fragments were 240 bp in length. The correlation among biological replicates was analyzed using Pearson's Correlation Coefficient (r) [76]. The six cDNA libraries were sequenced with Illumina HiSeq Xten from Biomarker Biotechnology Corporation (Beijing, China).

Transcriptome Analyses
The clean reads were compared with the T. aestivum reference genome [77] using TopHat2 [78]. The unigenes were annotated by BLAST [79] in Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (protein family), KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KO (KEGG Ortholog database), and GO (Gene Ontology) databases. To explore new transcripts and genes in the mutant dmc, the reads were assembled using Cufflinks software [80] based on the reference genome sequence, and the transcriptional intervals that had not been annotated were searched for and compared with the original genome annotation information. The KEGG Orthology of the new genes was analyzed using KOBAS2.0 [81], and the amino acid sequences of the new genes were predicted; finally, the new genes were annotated using HMMER [82], referring to the numerous databases. The standard gene expression levels were expressed as fragments per kilobase of transcript per million fragments mapped (FPKM) [83]. Differentially-expressed genes (DEGs) were identified by DEseq [84] with a false discovery rate (FDR) <0.01 and a fold change value (FC) ≥2.

qRT-PCR
The tiller primordia and leaf samples of the mutant dmc and guomai 301 were prepared at 10 time points for real-time qRT-PCR; the intervals of the 10 time points were seven days from 23 December 2016 to 24 February 2017. Twelve DEGs were selected to test and verify the sequencing data by qRT-PCR. The primers were designed using Primer Premier 5.0 (Table S13). The wheat actin gene was used as an internal control gene. The qRT-PCR reactions were performed in 20 µL volumes containing 10 µL TransStart ® Top Green Qpcr SuperMix (2×) (TransGen Biotech, Beijing, China), 2 µL primer mix (10 µM), 1 µL cDNA (50 ng) and 7 µL ddH 2 O. The PCR parameters were: 94 • C for 30 s, then 42 cycles of 94 • C for 5 s, 60 • C for 30 s. All qRT-PCR reactions were replicated three times. The gene expression levels were calculated according to the 2 −∆∆Ct method [85].

Ethical Approval
This article does not contain any studies with human participants or animals performed by any of the authors.
Wheat (Triticum aestivum L.) plants were used in this study. The wheat cultivar 'Guomai 301' was bred in our laboratory.
The wheat mutant dmc was selected from EMS treated Guomai 301 in our laboratory.

Conclusions
Transcriptomes of the tiller primordia from the wheat non-tiller mutant dmc and WT Guomai 301 were integratively analyzed. We identified a set of genes related to wheat tiller differentiation. Sixty-nine percent of significant DEGs (FC ≥ 2) were highly expressed and 31% were expressed at a low level in dmc. Phenylpropanoid biosynthesis (ko00940; Figure S10) had the most highly expressed DEGs (9.50%). Carbohydrate-related metabolisms consisted of the largest highly-expressed DEG group, including photosynthesis (6.74%), starch and sucrose metabolism (6.54%), photosynthesis-antenna proteins (6.23%), carbon fixation in photosynthetic organisms (5.52%), carbon metabolism (7.76%), and galactose metabolism (4.09%). The majority of genes that were expressed at a low level belonged to the classes of DNA replication, transcription, and translation, as well as zeatin synthesis. Functional model analysis indicated that variations in carbohydrate metabolism, phytohormones and transcription factors were the major causes of non-tillering in dmc. Synthesis and lower expression of genes related to protein synthesis, auxin, zeatin and cytokinine; syntheses enhanced the expression of genes involved in the biosynthesis of abscisic acid, gibberellins, and ethylene syntheses, consistent with the phenotype of dmc. The expression profiles of TF homologs, knotted-1-like 1, homeobox protein knotted-1-like 12 and KNOX3, and WOX 4 were also consistent with the phenotype of dmc. Other issues need further study, such as the relationship between phenylpropanoid biosynthesis, photosynthesis, variations of different phytohormone concentrations and the phenotype of dmc.