Transcriptome Analysis of a Premature Leaf Senescence Mutant of Common Wheat (Triticum aestivum L.)

Leaf senescence is an important agronomic trait that affects both crop yield and quality. In this study, we characterized a premature leaf senescence mutant of wheat (Triticum aestivum L.) obtained by ethylmethane sulfonate (EMS) mutagenesis, named m68. Genetic analysis showed that the leaf senescence phenotype of m68 is controlled by a single recessive nuclear gene. We compared the transcriptome of wheat leaves between the wild type (WT) and the m68 mutant at four time points. Differentially expressed gene (DEG) analysis revealed many genes that were closely related to senescence genes. Gene Ontology (GO) enrichment analysis suggested that transcription factors and protein transport genes might function in the beginning of leaf senescence, while genes that were associated with chlorophyll and carbon metabolism might function in the later stage. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that the genes that are involved in plant hormone signal transduction were significantly enriched. Through expression pattern clustering of DEGs, we identified 1012 genes that were induced during senescence, and we found that the WRKY family and zinc finger transcription factors might be more important than other transcription factors in the early stage of leaf senescence. These results will not only support further gene cloning and functional analysis of m68, but also facilitate the study of leaf senescence in wheat.


Introduction
Senescence plays a key role in the adaptability of plants. Effective senescence can enhance the adaptation of plants to the environment [1]. In crops, leaf senescence is an important agronomic trait that not only affects crop yield, but also affects crop quality [2,3]. Leaf senescence is the last stage of leaf development that can induce leaf cell death and nutrient transfer from senescing leaves to seeds or storage organs. This process is not a negative or unregulated degeneration step; rather, it is a positive step that involves fine gene regulation [4]. During senescence, the metabolism of leaf cells changes. Specifically, assimilation decreases the while catabolism is enhanced, e.g., chloroplast degradation occurs, the photosynthetic capacity decreases, and macromolecular material degrades [5]. Leaf senescence is affected by both internal and external factors. The former mainly include leaf age, foliar reactive oxygen species (ROS) levels, hormones, and sugar content. The latter mainly include lighting conditions, soil moisture, nutrient content, temperature, and pathogens [6,7]. Because the regulation of leaf senescence is so important and complicated, many studies on senescence-associated genes (SAGs) have been carried out in different plant species [5,8].
senescence and affects the protein, zinc, and iron contents of wild wheat grains [27]. TaNAMs were down-regulated by RNAi technology, which led to the stay-green phenotype and the reduction in nutritive material transport from flag leaves to seeds [38]. The difficulty that is associated with wheat gene studies is caused by the complexity of the huge wheat genome, with a high ratio of repeat sequences [39]. Recently, in the wake of the development of sequencing and assembly technologies, the wheat reference genome has been greatly improved to provide a favorable condition for the study of wheat genes.
Using RNA-seq to study leaf senescence in wheat has rarely been reported. In the present study, RNA-seq analysis was performed with the leaf senescence mutant m68 and the wild-type (WT) variety YZ4110 at different time points. Differentially expressed genes (DEGs) were identified between m68 and the WT, and the Gene Ontology (GO) enrichment, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of DEGs were analyzed at different time points. Analysis of transcriptome data will be helpful for understanding leaf senescence in wheat.

Phenotypic Characterization and Genetic Analysis
At the seedling stage, there was no obvious difference between the WT and m68 mutant ( Figure 1A,B). Initial senescence only occurred in the bottom leaves at the shooting stage in m68 ( Figure 1C). Then, the leaves in the middle part of m68 plants showed senescence symptoms around the heading date. Approximately one week later, flag leaves started to senescence. During the filling stage, the leaf senescence process of the m68 mutant accelerated dramatically ( Figure 1D). The final spikelet number, grain number, floret number per spike, 1000-grain weight, plant height, and effective tiller number of m68 were significantly different from those of the WT (Table 1). To determine the genetic basis of these phenotypes, we conducted genetic analysis. Among a total of 318 F 2 plants, 73 showed the premature leaf senescence phenotype of m68, whereas the others showed the WT phenotype. The segregation ratio of the m68 phenotype to the WT phenotype was 1:3.36. Based on the chi-square (χ 2 ) test, the F 2 population segregation ratio was in accordance with the expected ratio of 1:3 (χ 2 = 0.71 < χ 2 0.05,1 = 3.84). The results indicated that premature leaf senescence of m68 is controlled by a single recessive nuclear gene. cassette (ABC) transporter that provides durable resistance to multiple fungal pathogens and induces flag leaf tip necrosis-like leaf senescence [37]. GPC-B1 is an NAC family transcription factor that accelerates leaf senescence and affects the protein, zinc, and iron contents of wild wheat grains [27]. TaNAMs were down-regulated by RNAi technology, which led to the stay-green phenotype and the reduction in nutritive material transport from flag leaves to seeds [38]. The difficulty that is associated with wheat gene studies is caused by the complexity of the huge wheat genome, with a high ratio of repeat sequences [39]. Recently, in the wake of the development of sequencing and assembly technologies, the wheat reference genome has been greatly improved to provide a favorable condition for the study of wheat genes.
Using RNA-seq to study leaf senescence in wheat has rarely been reported. In the present study, RNA-seq analysis was performed with the leaf senescence mutant m68 and the wild-type (WT) variety YZ4110 at different time points. Differentially expressed genes (DEGs) were identified between m68 and the WT, and the Gene Ontology (GO) enrichment, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of DEGs were analyzed at different time points. Analysis of transcriptome data will be helpful for understanding leaf senescence in wheat.

Phenotypic Characterization and Genetic Analysis
At the seedling stage, there was no obvious difference between the WT and m68 mutant ( Figure  1A,B). Initial senescence only occurred in the bottom leaves at the shooting stage in m68 ( Figure 1C). Then, the leaves in the middle part of m68 plants showed senescence symptoms around the heading date. Approximately one week later, flag leaves started to senescence. During the filling stage, the leaf senescence process of the m68 mutant accelerated dramatically ( Figure 1D). The final spikelet number, grain number, floret number per spike, 1000-grain weight, plant height, and effective tiller number of m68 were significantly different from those of the WT (Table 1). To determine the genetic basis of these phenotypes, we conducted genetic analysis. Among a total of 318 F2 plants, 73 showed the premature leaf senescence phenotype of m68, whereas the others showed the WT phenotype. The segregation ratio of the m68 phenotype to the WT phenotype was 1:3.36. Based on the chi-square (χ 2 ) test, the F2 population segregation ratio was in accordance with the expected ratio of 1:3 (χ 2 = 0.71 < χ 2 0.05,1 = 3.84). The results indicated that premature leaf senescence of m68 is controlled by a single recessive nuclear gene.   As the development of flag leaves is crucial to agronomic traits, we observed the flag leaf phenotype of the m68 mutant in detail beginning on the heading date, i.e., 19 April 2016. We separated the phenotype into four stages. In the first stage (stage 1, S1), on 19 April 2016, the flag leaves of m68 and the WT showed almost no differences (Figure 2A). In the second stage (stage 2, S2), on 25 April, the m68 flag leaf tip exhibited senescence ( Figure 2B). Then, the flag leaf of m68 underwent accelerated senescence, and senescence was obvious on 3 May 2016 (stage 3, S3) ( Figure 2C). In the final stage (stage 4, S4), the flag leaf of m68 had almost completely turned yellow, while the flag leaf of the WT had not senesced ( Figure 2D). Our measurements of maximum quantum efficiency of photosystem II (F v /F m ) and the malondialdehyde (MDA) contents indicated that in m68 leaves, F v /F m gradually decreased ( Figure 2E), but the MDA content increased significantly as the leaf senescence progressed ( Figure 2F). The expression of TaPAO, which is associated with chlorophyll degradation, increased rapidly during leaf senescence in m68 ( Figure 2G). As the development of flag leaves is crucial to agronomic traits, we observed the flag leaf phenotype of the m68 mutant in detail beginning on the heading date, i.e., 19 April 2016. We separated the phenotype into four stages. In the first stage (stage 1, S1), on 19 April 2016, the flag leaves of m68 and the WT showed almost no differences (Figure 2A). In the second stage (stage 2, S2), on 25 April, the m68 flag leaf tip exhibited senescence ( Figure 2B). Then, the flag leaf of m68 underwent accelerated senescence, and senescence was obvious on 3 May 2016 (stage 3, S3) ( Figure  2C). In the final stage (stage 4, S4), the flag leaf of m68 had almost completely turned yellow, while the flag leaf of the WT had not senesced ( Figure 2D). Our measurements of maximum quantum efficiency of photosystem II (Fv/Fm) and the malondialdehyde (MDA) contents indicated that in m68 leaves, Fv/Fm gradually decreased ( Figure 2E), but the MDA content increased significantly as the leaf senescence progressed ( Figure 2F). The expression of TaPAO, which is associated with chlorophyll degradation, increased rapidly during leaf senescence in m68 ( Figure 2G).

RNA-Seq Analysis and DEGs
In total, 24 sample libraries were constructed and sequenced. We obtained 40.79 to 41.51 million reads from each sample. The GC content of the raw reads ranged from 53.66 to 56.33% in different libraries, and the Q30 percentage exceeded 94.5%. These libraries contained 85.61-89.34% of mapped reads and 77.81-82.53% of unique reads mapped to the T. aestivum _CS42_TGAC_v1 ( Table 2). GAPDH was selected as the internal standard. Error bars represent the means ± SD (n = 3). ** indicates p < 0.01; * indicates p < 0.05.

RNA-Seq Analysis and DEGs
In total, 24 sample libraries were constructed and sequenced. We obtained 40.79 to 41.51 million reads from each sample. The GC content of the raw reads ranged from 53.66 to 56.33% in different libraries, and the Q30 percentage exceeded 94.5%. These libraries contained 85.61-89.34% of mapped reads and 77.81-82.53% of unique reads mapped to the T. aestivum _CS42_TGAC_v1 ( Table 2). The differentially expressed genes (DEGs) between m68 and the WT at the four development stages were identified to investigate the genes that might be responsible for leaf senescence in the m68 mutant (Supplementary Table S1). At S1 (WT_S1 vs. M_S1), although the flag leaf of the m68 mutant had not started senescing, 2396 genes were significantly differentiated. These included 2030 up-regulated and 366 down-regulated genes in m68 ( Figure 3A). When the m68 flag leaf tip exhibited senescence at S2, 8383 DEGs were identified (WT_S2 vs. M_S2), including 5241 up-regulated and 3142 down-regulated genes in m68 ( Figure 3B). During S3, the leaf senescence phenotype of m68 was obvious, and more DEGs between WT_S3 and M_S3 were found, including 6903 up-regulated and 6452 down-regulated genes in m68 ( Figure 3C). At S4, the m68 mutant flag leaf had turned completely yellow. The number of DEGs between WT_S4 and M_S4 increased to 16130, which included 7369 up-regulated and 8761 down-regulated genes in m68 ( Figure 3D). Moreover, we also detected the specific DEGs at different stages ( Figure 3E,F); the Venn diagrams showed that 94, 698, 1652, and 3911 genes were down-regulated at S1, S2, S3, and S4, respectively, while 367, 1358, 1616, and 2842 genes specific were up-regulated at S1, S2, S3, and S4, respectively. The numbers of consistently down-and up-regulated genes in all of the stages were 140 and 805, respectively.

GO Enrichment Analysis of DEGs
The enriched GO terms of down-and up-regulated DEGs at different time points were analyzed ( Figure 4 and Supplementary Table S2). A total of 141 down-regulated DEGs were significantly enriched in GO terms. Among these terms, none were found at S1 ( Figure 4A). However, 52 were found at S2, such as the chloroplast thylakoid membrane, photosynthesis, light harvesting, photosystem II, and chlorophyll binding. Ninety-one terms, such as ribulose-bisphosphate carboxylase activity, fructose 1,6-bisphosphate 1-phosphatase activity, ribosome, carbohydrate metabolic process, and chloroplast, were detected at S3. At S4, there existed 108 GO terms, among which, 29 were simultaneously enriched at S2 and S3 ( Figure 4A). In terms of the up-regulated DEGs, there were 143 significantly enriched GO terms ( Figure 4B). At S1, there were 48 including thirteen stage-specific GO terms such as lyase activity phosphate ion transport, transcription factor activity, sequence-specific DNA binding and inorganic phosphate transmembrane transporter activity. Sixty-one terms, such as integral component of membrane, response to stress, sucrose alpha-glucosidase activity, chlorophyll catabolic process, and oxidation-reduction process, were detected at S2. Sixty and 80 were found at S3 and S4, respectively. Among the 143 significantly enriched GO terms, 16 were enriched at all four development stages ( Figure 4B), including terms in calcium ion binding, carbohydrate metabolic process, oxidation-reduction process, and integral component of the membrane.    The Y-axis indicates the number of genes belonging to the GO terms below. The X-axis indicates GO terms. The DEGs at S1, S2, S3, and S4 are shown in red, blue, purple and yellow, respectively.

KEGG Pathway Analysis of DEGs
To further investigate the leaf senescence pathway in the m68 mutant, KEGG analysis of DEGs was performed at different stages. Among the four stages, DEGs were significantly enriched in 18 KEGG pathways (Figure 5). At S1, DEGs were enriched in 5 KEGG pathways, and four of them, i.e., the plant-pathogen interaction (ko04626), glutathione metabolism (ko00480), the calcium signaling pathway (ko04020), and the cyclic adenosine 3′,5′-monophosphate (cAMP) signaling pathway (ko04024), were only enriched in this stage. DEGs at S2 and S3 shared 6 enriched KEGG pathways, The X-axis indicates GO terms. The DEGs at S1, S2, S3, and S4 are shown in red, blue, purple and yellow, respectively.

KEGG Pathway Analysis of DEGs
To further investigate the leaf senescence pathway in the m68 mutant, KEGG analysis of DEGs was performed at different stages. Among the four stages, DEGs were significantly enriched in 18 KEGG pathways (Figure 5). At S1, DEGs were enriched in 5 KEGG pathways, and four of them, i.e., the plant-pathogen interaction (ko04626), glutathione metabolism (ko00480), the calcium signaling pathway (ko04020), and the cyclic adenosine 3 ,5 -monophosphate (cAMP) signaling pathway (ko04024), were only enriched in this stage. DEGs at S2 and S3 shared 6 enriched KEGG pathways, including starch and sucrose metabolism (ko00500), plant hormone signal transduction (ko04075), carbon fixation in photosynthetic organisms (ko00710), amino sugar and nucleotide sugar metabolism (ko00520), glyoxylate and dicarboxylate metabolism (ko00630), and photosynthesis (ko00195). At the late stage of leaf senescence, i.e., S4, DEGs were enriched in seven KEGG pathways, involving carbon metabolism (ko01200), lysosomes (ko04142), glycolysis (ko00010), plant hormone signal transduction (ko04075), starch and sucrose metabolism (ko00500), porphyrin and chlorophyll metabolism (ko00860) and alanine, aspartate, and glutamate metabolism (ko00250). including starch and sucrose metabolism (ko00500), plant hormone signal transduction (ko04075), carbon fixation in photosynthetic organisms (ko00710), amino sugar and nucleotide sugar metabolism (ko00520), glyoxylate and dicarboxylate metabolism (ko00630), and photosynthesis (ko00195). At the late stage of leaf senescence, i.e., S4, DEGs were enriched in seven KEGG pathways, involving carbon metabolism (ko01200), lysosomes (ko04142), glycolysis (ko00010), plant hormone signal transduction (ko04075), starch and sucrose metabolism (ko00500), porphyrin and chlorophyll metabolism (ko00860) and alanine, aspartate, and glutamate metabolism (ko00250). At S2, the m68 flag leaf tips exhibited senescence, and many DEGs were significantly enriched in plant hormone signal transduction; therefore, genes that were involved in the plant hormone signal transduction pathway were selected and analyzed further. We identified 24 DEGs orthologs of the wheat genes related to ABA, IAA, ET, JA, and SA responses ( , which were also upregulated, belonged to the SA signaling pathway. We checked the expression of six SA response gene orthologs using qRT-PCR ( Figure 6). The results showed that the expression of these genes was obviously increased, which is consistent with the transcriptome analysis. At S2, the m68 flag leaf tips exhibited senescence, and many DEGs were significantly enriched in plant hormone signal transduction; therefore, genes that were involved in the plant hormone signal transduction pathway were selected and analyzed further. We identified 24 DEGs orthologs of the wheat genes related to ABA, IAA, ET, JA, and SA responses ( Table 3) , which were also up-regulated, belonged to the SA signaling pathway. We checked the expression of six SA response gene orthologs using qRT-PCR ( Figure 6). The results showed that the expression of these genes was obviously increased, which is consistent with the transcriptome analysis.

Expression Pattern Clustering Analysis of DEGs
To identify genes closely that were related to senescence, we performed expression pattern cluster analysis of DEGs at different time points. All of the DEGs were used for cluster analysis. Hierarchical clustering based on DEGs ( Figure 7A) and a 10% tree cutoff was performed using the R package cutree. In total, 4708 genes were classified into 393 clusters, which showed different

Expression Pattern Clustering Analysis of DEGs
To identify genes closely that were related to senescence, we performed expression pattern cluster analysis of DEGs at different time points. All of the DEGs were used for cluster analysis. Hierarchical clustering based on DEGs ( Figure 7A) and a 10% tree cutoff was performed using the R package cutree. In total, 4708 genes were classified into 393 clusters, which showed different expression patterns between the WT and m68. We further analyzed the genes of 27 subclusters, which contained 1012 genes ( Figure 7B and Supplementary Table S3). The expression levels of these cluster genes were almost unchanged in the WT at all four time points, but in the m68 mutant, they were obviously up-regulated during leaf senescence. There were 75 genes encoding kinase family proteins, 23 of which were receptor-like kinase genes. Additionally, 56 genes were annotated as transcription factors or DNA binding factors, 24 of which were orthologs of WRKY transcription factors, and 26 were orthologs of zinc finger transcription factors. In addition, 42 genes encoded different types of transporters, including inorganic phosphate transporters, high-affinity potassium transporters, boron transporter proteins, potassium transporters, ABC transporters, and amino acid transporters. expression patterns between the WT and m68. We further analyzed the genes of 27 subclusters, which contained 1012 genes ( Figure 7B and Supplementary Table S3). The expression levels of these cluster genes were almost unchanged in the WT at all four time points, but in the m68 mutant, they were obviously up-regulated during leaf senescence. There were 75 genes encoding kinase family proteins, 23 of which were receptor-like kinase genes. Additionally, 56 genes were annotated as transcription factors or DNA binding factors, 24 of which were orthologs of WRKY transcription factors, and 26 were orthologs of zinc finger transcription factors. In addition, 42 genes encoded different types of transporters, including inorganic phosphate transporters, high-affinity potassium transporters, boron transporter proteins, potassium transporters, ABC transporters, and amino acid transporters.

Discussion
In this study, we identified a wheat mutant, m68, which had many effects on plant development. Agronomic traits, such as grain weight, kernel number per spike, and floret number decreased significantly in m68. Because the leaves of m68 started to senescence much earlier than those of the WT, we found some indicators that were detected in previous senescence studies. These senescence indicators include MDA, which is toxic to cells and can lead to cell death; Fv/Fm, which affects leaf photosynthesis and the accumulation of photosynthetic products; and, TaPAO, which participates in the degradation of chlorophyll [8,14]. We examined those indicators and found that the MDA content and the expression level of TaPAO were significantly up-regulated, and Fv/Fm decreased remarkably in m68 (Figure 2E-G). Therefore, m68 is a premature senescence mutant. Using genetic analysis, we further determined that the premature senescence phenotype was controlled by a single recessive nuclear gene.
We randomly selected 16 genes among the DEGs to test their expression at different stages using qRT-PCR. Our qRT-PCR analysis results showed that those genes were indeed differentially expressed, which is in accordance with our transcriptome analysis results (Supplementary Figure S1). Therefore, our genome-wide transcriptome profiling analysis is reliable. The detection of specific DEGs at different stages indicated that the numbers of both the specific down-regulated and upregulated genes were increased from S1 to S4. Through the GO enrichment analysis of down-and up-regulated DEGs at different time points, we detected a total of 141 and 143 enriched GO terms in down-and up-regulated DEGs at different time points, respectively. We found that the downregulated DEGs did not have significantly enriched GO term at S1, while for the up-regulated DEGs, significantly enriched GO terms appeared at S1. For up-regulated DEGs, thirteen of 143 GO terms

Discussion
In this study, we identified a wheat mutant, m68, which had many effects on plant development. Agronomic traits, such as grain weight, kernel number per spike, and floret number decreased significantly in m68. Because the leaves of m68 started to senescence much earlier than those of the WT, we found some indicators that were detected in previous senescence studies. These senescence indicators include MDA, which is toxic to cells and can lead to cell death; F v /F m , which affects leaf photosynthesis and the accumulation of photosynthetic products; and, TaPAO, which participates in the degradation of chlorophyll [8,14]. We examined those indicators and found that the MDA content and the expression level of TaPAO were significantly up-regulated, and F v /F m decreased remarkably in m68 ( Figure 2E-G). Therefore, m68 is a premature senescence mutant. Using genetic analysis, we further determined that the premature senescence phenotype was controlled by a single recessive nuclear gene.
We randomly selected 16 genes among the DEGs to test their expression at different stages using qRT-PCR. Our qRT-PCR analysis results showed that those genes were indeed differentially expressed, which is in accordance with our transcriptome analysis results (Supplementary Figure S1). Therefore, our genome-wide transcriptome profiling analysis is reliable. The detection of specific DEGs at different stages indicated that the numbers of both the specific down-regulated and up-regulated genes were increased from S1 to S4. Through the GO enrichment analysis of down-and up-regulated DEGs at different time points, we detected a total of 141 and 143 enriched GO terms in down-and up-regulated DEGs at different time points, respectively. We found that the down-regulated DEGs did not have significantly enriched GO term at S1, while for the up-regulated DEGs, significantly enriched GO terms appeared at S1. For up-regulated DEGs, thirteen of 143 GO terms were specifically enriched at S1, which may play important roles in induced leaf senescence or onset of leaf senescence-related signal transduction, such as transcription factor activity and sequence-specific DNA binding. AtWRKY53, AtWRKY70, and AtWRKY40 transcription factors participated in SA signaling transduction pathways [30,31,56], the expression levels of their orthologous genes were up-regulated in m68. AtWRKY67 might participate in stress resistance-induced leaf senescence [57], and the expression levels of its orthologous genes were also up-regulated in m68. Therefore, the functions of these WRKY transcription factors might be very conserved in regulated leaf senescence and may include functions as upstream regulators of senescence. The mutant gene in m68 might function upstream of the senescence regulatory network.
We also found that down-regulated DEGs were enriched in GO terms such as chlorophyll binding, chloroplast, chloroplast organization, photosystem I, and photosystem II at S2, S3, and S4, which were defined as the degenerative phase in previous studies [58]. In accordance with the results from previous studies, many chlorophyll metabolism genes were differentially expressed in m68, which might play important roles in wheat leaf senescence. For example, the expression of orthologs of chloroplast precursor-relative and TPR (tetratricopeptide repeat) genes, such as MSTRG.35978.1, MSTRG.34045.1, and MSTRG.3773.1, were reduced in m68. Chlorophyll synthesis-related genes, such as orthologs of At5G54190, At1G44446 and At3G51820 genes, were reduced in m68. GATA transcription factors can influence leaf senescence by regulating the chlorophyll content [17][18][19][20]. In our study, the GATA zinc finger family genes MSTRG.14332.1 and MSTRG.41348.1 were down-regulated during leaf senescence, which agrees with the function of GATA transcription factors. PPR proteins can influence the chloroplast content by regulating chloroplast gene expression, such as OspTAC2 in rice [59] and PPR53 can also affect chloroplast gene expression and leaf development in maize [16]. In our study, PPR-like genes MSTRG.37348.1 and MSTRG.23240.1 were down-regulated in m68. The above mentioned genes might play crucial roles in leaf senescence in wheat at the S2, S3 and S4 stages.
KEGG pathway analysis showed that plant hormone signaling severely affected the leaf senescence process. A previous study showed that ET, SA, JA, and ABA promoted leaf senescence, while CK, GA and IAA were observed to delay the progression of leaf senescence [21]. In our study, we found that the orthologs of the SA receptor gene NPR3 and its downstream genes, CaMs and PR-related genes, such as MSTRG.13794.1, MSTRG.27143.1, MSTRG.29737.1, MSTRG.32097.1, and MSTRG.26416.1, were up-regulated in the m68 mutant ( Table 3). The JA response gene JAZ2 was also up-regulated ( Table 3). The expression levels of ABA signal transduction pathway genes, such as ABF1, ABF4, and HAB1 significantly increased in m68 ( Table 3). The IAA signaling transduction-related genes IAA7, IAA13, and IAA11 were down-regulated in the m68 mutant (Table 3). However, ET receptor gene EIN2 expression was down-regulated in the m68 mutant, which is supposed to promote leaf senescence (Table 3). These results seem to imply that leaf senescence in wheat is regulated by a variety of hormones and might have a unique regulatory mechanism.
The clusters of the expression pattern of DEGs provided us with a good opportunity to identify SAGs in wheat. We mainly paid attention to the clusters whose expression pattern had no obvious change during WT leaf development but showed a rapid increase during the leaf development of the mutant. There were 27 subclusters that contained 1012 genes showing this expression pattern ( Figure 7B and Supplementary Table S3). All of the genes in this cluster were induced significantly during senescence and might be essential to wheat leaf senescence. Among these genes, many types of receptor-like protein kinases, transporter family proteins and transcription factors were enriched. The number of WRKY family and zinc finger transcription factors was much higher than that of the other transcription factors in these subclusters. It is suggested that they may play much more important roles than the other transcription factors in the regulation of leaf senescence in wheat. In addition, we noticed that flag leaf senescence occurred after the heading date, which was much later than that in the lower leaves of m68. Therefore, we deduced that oriented signal transduction and substance transport occur in plants around the heading date and are important in leaf senescence. We identified receptor genes and transporter genes that might participate in signal transduction and substance transport during leaf senescence of wheat.

Plant Materials and Genetic Analysis
The m68 mutant was generated from an EMS mutant library of the common wheat cultivar Yanzhan 4110 (YZ4110), which was constructed and conserved in our laboratory. We planted m68 and the YZ4110 in 2013-2016 in Beijing, and investigated agronomic traits, such as plant height and spike length each year.
We obtained inheritance-stable homozygote m68 mutant seeds and crossed them with the YZ4110 to generate a heterozygote F 1 population. F 1 seeds were produced from self-crossed F 1 plants and were sown on 7 October 2015, in Beijing. In each row, 15 seeds were evenly planted along a length of 2 m. We observed the flag leaves of 318 F 2 plants on 3 May 2016. The number of plants exhibiting premature senescence was counted. A chi-square (χ 2 ) test was used to analyze the segregation ratios.

Plant Materials Used in RNA-Seq
The WT and m68 mutant were sown on 7 October 2015, i.e., 30 rows each, in Beijing (China), and were grown under field conditions. In each row, 15 seeds were evenly planted along a length of 2 m. The flag leaves of the WT and m68 mutant were collected from the heading stage (19 April 2016) and were continuously sampled every 6 days until 9 May 2016. In each of the four stages, leaves from 10 plants were pooled to produce three biological replicates. All of the samples were immediately frozen in liquid nitrogen and were stored at −80 • C in a refrigerator.

F v /F m and MDA Content Measurement
The F v /F m and the MDA content of flag leaves of the WT and m68 were measured at different sampling times. F v /F m was measured using a handheld fluorometer (FluorPen FP100, Photon Systems Instruments, Drasov, Czech Republic), according to the manufacturer's instructions. At each time point, 20 plants were measured; an average of three measurements were collected per plant. The MDA content was measured by utilizing a 5% thiobarbituric acid reaction, according to previous study methods [60].

RNA Extraction and Library Construction
Total RNA was isolated using Trizol reagent (Promega, Madison, WI, USA) following the manufacturer's instructions. RNA purity was tested using a NanoDrop 2000 (Thermo, Waltham, MA, USA), RNA concentration was measured using an Agilent 2100 Bioanalyzer, and RNA integrity was evaluated with an Agilent RNA6000 NanoKit according to the manufacturer's protocol.
Each sample required approximately 4 to 8 µg of total RNA to construct the RNA-Seq libraries. The method consisted mainly of the following steps: (1) mRNA was first purified from total RNA using polyA magnetic bead enrichment. (2) Chemical fragmentation was conducted. (3) Random hexamer priming was used to convert mRNA into first single-strand cDNA and to synthesize second-strand cDNA. (4) Sequencing adaptors were added to cDNA fragments. Suitable fragments were selected based on the separation by agarose gel electrophoresis. (5) PCR amplification construction of RNA libraries was completed. In total, 24 RNA libraries were constructed, and the quality and quantity of each RNA library were assessed using Nanodrop ND-1000 spectroscopy and an Agilent 2100 Bioanalyzer according to the protocol guide.

Illumina Sequencing and Read Mapping
RNA-Seq was performed using an Illumina HiSeq 4000 sequencing platform by OnMath Technologies, Chengdu, Sichuan, China. The raw reads were cleaned by removing adapter and low-quality sequences. Mapping of 24 sample reads to wheat genome (T. aestivum. TGACv1) was conducted with Tophat (version 2.1.0, default parameters). The new transcripts were assembled with StringTie (version 1.3.3b, default parameters). New gtf assembling and merging for each sample were used known TGAC gtf (version: TGACv1).
All of the sequencing data in this study are stored at the Sequence Read Archive under BioProject ID PRJA431543.
GO enrichment analysis of DEGs was performed using the Goseq R package. All of the significantly DEGs between compare groups (both known genes and new assembled genes) were blasted in Pfam database (version 31.0, March 2017) and extracted all gene's GO id for enrichment. GOseq was applied for enrichment and topGO was used for plot directed acyclic graph base on significantly enrichment gene (p value < 0.05 and q value < 0.05) [63]. We used KOBAS software 2.0 to determine significant enrichment of DEGs in KEGG pathways (p-value < 0.05 and corrected p-value < 0.05) [64]. Hierarchical clustering based on DEGs and a 10% tree cutoff was performed using the R package cutree.

Quantitative Real-Time PCR
To validate the RNA-seq data, we randomly selected 16 genes in DEGs to test their expression at different stages using qRT-PCR, which were the same RNA samples used for RNA-seq. According to the sequences, primer version 5.0 was used to design gene-specific primers. The primer sequences are listed in Supplementary Table S4. qRT-PCR was performed using a Roche LightCycler 480 Real-Time System (Roche, Switzerland). The wheat gene GAPDH (NCBI accession: AF251217.1) was used as an internal control for normalization of expression. Each experiment was performed with three biological replicates. qRT-PCR was performed using a 15-µL reaction volume that consisted of 7.5 µL of SYBR Mix (Toyobo, Japan), 1 µL of cDNA, 2 µL of gene-specific primers (2 µmol·L −1 ), and 4.5 µL of ddH 2 O. The PCR program parameters were 95 • C for 1 min, followed by 40 cycles at 95 • C for 20 s, 60 • C for 20 s, and 72 • C for 40 s. The relative gene expression levels were evaluated according to a previous study [65].