- © 2019 American Society of Plant Biologists. All rights reserved.
Abstract
Premature senescence in annual crops reduces yield, while delayed senescence, termed stay-green, imposes positive and negative impacts on yield and nutrition quality. Despite its importance, scant information is available on the genetic architecture of senescence in maize (Zea mays) and other cereals. We combined a systematic characterization of natural diversity for senescence in maize and coexpression networks derived from transcriptome analysis of normally senescing and stay-green lines. Sixty-four candidate genes were identified by genome-wide association study (GWAS), and 14 of these genes are supported by additional evidence for involvement in senescence-related processes including proteolysis, sugar transport and signaling, and sink activity. Eight of the GWAS candidates, independently supported by a coexpression network underlying stay-green, include a trehalose-6-phosphate synthase, a NAC transcription factor, and two xylan biosynthetic enzymes. Source–sink communication and the activity of cell walls as a secondary sink emerge as key determinants of stay-green. Mutant analysis supports the role of a candidate encoding Cys protease in stay-green in Arabidopsis (Arabidopsis thaliana), and analysis of natural alleles suggests a similar role in maize. This study provides a foundation for enhanced understanding and manipulation of senescence for increasing carbon yield, nutritional quality, and stress tolerance of maize and other cereals.
INTRODUCTION
Agricultural productivity is essentially the amount of carbohydrates generated by photosynthetic assimilation of CO2 by plant leaves and stored into heterotrophic organs harvested for human/industrial use. Improving the photosynthetic assimilation of crop plants therefore is a viable approach to increase agricultural productivity. However, the photosynthetic ability is inherently linked with the plant developmental program, and for monocarpic (annual) plants that are the major source of food, starts to decrease during the reproductive phase. For instance, during the grain-filling period in maize (Zea mays), one of the most important staple foods on earth, photosynthesis rates decline by as much as 50% within a 6-week period after anthesis (Ying et al., 2000). The decline of photosynthesis is triggered upon the onset of senescence, a highly regulated, well-coordinated, and biologically active process that marks the end of the life cycle of the leaf and ultimately the whole plant (Quirino et al., 2000; Buchanan-Wollaston et al., 2003; Lim et al., 2007; Wingler et al., 2009; Thomas, 2013).
Senescence is important for recycling resources, particularly nitrogen and carbon, from old organs to either the newly developing organs or storage sinks, thereby contributing to the fitness of a plant (Lim et al., 2007). However, most of the dry matter accumulating in grains of modern maize hybrids is fixed during the grain fill, with very little attributed to remobilization from organs developed during the pre-flowering period (Below et al., 1981; Cliquet et al., 1990; Ciampitti and Vyn, 2013). Thus, given the potential benefits of delaying senescence without adversely affecting grain yield, developing late senescing or stay-green cultivars with delayed physiological maturity and prolonged grain-filling period is an important crop breeding strategy (Duvick et al., 2010; Gregersen et al., 2013). In maize, selection for stay-green trait accounts for up to 63% of the total increase in dry matter accumulation in newer era hybrids released after the 2000s compared with older hybrids released in the 1930s (Lee and Tollenaar, 2007). Stay-green has been associated with increased drought tolerance in multiple crop species including wheat (Triticum aestivum), sorghum (Sorghum bicolor), barley (Hordeum vulgare), and maize (Tuinstra et al., 1997; Vijayalakshmi et al., 2010; Emebiri, 2013; Trachsel et al., 2016). Maize and sorghum stay-green lines possess increased resistance to diseases, especially those caused by stalk-rotting pathogens (King and Frederickson, 1976; Bertolini et al., 1983). Stay-green also improves produce quality such as higher Suc and protein content in maize and higher stalk carbohydrates in sorghum (McBee et al., 1983; Gentinetta et al., 1986). For certain production systems, the stay-green trait can also have serious negative consequences such as delayed dry down resulting in higher moisture in the stover (leaves, stalks, and husk), which interferes with timely mechanical harvesting with combine harvesters (Yang et al., 2010).
Senescence is a complex and dynamic process that, while developmentally programmed, can be initiated and/or accelerated by a number of internal and external factors and overlaps with other processes including autophagy, aging, and death. Genetic dissection of this quantitative trait in Arabidopsis (Arabidopsis thaliana) has implicated a wide range of processes and genes including those involved in the synthesis of hormones such as cytokinin and ethylene, sugar sensing and signaling, source–sink translocation, and autophagy and aging (Schippers, 2015; Wingler, 2018). Senescence is also promoted by several environmental factors including lack of light (darkness), excessive light intensity/dosage, drought, salinity, nutrient (e.g., nitrogen, sulfur, and iron) starvation, and biotic stresses (Quirino et al., 2000; Lim et al., 2007; Thomas, 2013; Schippers, 2015). The partitioning of sugars to various sinks and the associated signals emerging from the interplay of the photosynthetic output of the plant (source strength) and the need for the photoassimilate (sink strength) play a key role in the regulation of senescence (Lee and Tollenaar, 2007; Wingler et al., 2009; Thomas, 2013; Thomas and Ougham, 2014). A weaker source induces early senescence due to increased mobilization of dry matter from leaves and other organs to the developing sinks, thereby starving the source organs (Rajcan and Tollenaar, 1999; Thomas, 2013). Weaker sink demand also accelerates senescence (Lee and Tollenaar, 2007), likely due to the accumulation of sugars in leaves beyond a certain critical threshold followed by sugar-mediated signaling of leaf senescence (Wingler and Roitsch, 2008; Thomas, 2013). While the role of source–sink interaction in the regulation of senescence is well established, the underlying mechanisms and genetic elements are largely unknown.
Valuable insights into the genetic architecture of senescence have been provided by the forward genetic screens in multiple plant species; however, these mostly represent the major effect of genes producing visible mutant phenotypes. Biparental mapping populations have also been used to identify quantitative trait loci in Arabidopsis (Wingler et al., 2010; Chardon et al., 2014), maize (Zheng et al., 2009; Wang et al., 2012; Belícuas et al., 2014; Trachsel et al., 2016), rice (Oryza sativa; Jiang et al., 2004; Abdelkhalik et al., 2005; Yoo et al., 2007), sorghum (Tuinstra et al., 1997; Xu et al., 2000; Harris et al., 2007), and wheat (Vijayalakshmi et al., 2010; Bogard et al., 2011; Pinto et al., 2016). These studies have been useful in understanding the genetic architecture of senescence but often lack the resolution to identify the underlying genes. Transcriptomic analyses in Arabidopsis have identified several nuclear senescence-associated genes (SAGs) and transcription factors, and chloroplast genes (Breeze et al., 2011; Guo and Gan, 2012; Woo et al., 2016). We have performed transcriptomic analysis of leaf and internode tissues in the maize B73 inbred line undergoing premature senescence induced by lack of seed sink and identified genes involved in source–sink-regulated senescence (Sekhon et al., 2012). Differentially expressed (DE) genes during senescence have also been identified in sorghum (Johnson et al., 2015; Wu et al., 2016), wheat (Gregersen and Holm, 2007), and rice (Ma et al., 2005). While these studies provide general insights into global senescence-related transcriptome, identification of the underlying causal genes is difficult due to a large number of DE candidates.
Given the diversity of internal and external regulators of senescence, a systems genetics (SysGen) approach (Civelek and Lusis, 2014) is needed for a comprehensive understanding of the genes and genetic networks underlying this complex developmental process. We characterized a large proportion of natural variation captured in a diversity panel and constructed gene coexpression networks based on comparative transcriptomic analysis of stay-green and natural senescence. Investigation and integration of these system-level data revealed both known and novel genes/processes and provided several new and valuable insights into the regulation of senescence in maize and other grasses.
RESULTS
Characterization of the Natural Variation for Leaf Senescence in Maize
Maize diversity characterized in this study was captured in a panel that (1) includes maize lines from two major heterotic groups (Stiff Stalk Synthetic and Non-Stiff Stalk Synthetic), (2) includes only those inbred lines that flowered within a reasonable window in the midwestern United States, and (3) was genotyped by RNA sequencing (RNA-seq) of whole seedling mRNA (Hirsch et al., 2014). We evaluated natural variation for onset of leaf senescence in a subset of this panel in three environments by growing 328 inbred lines during summer 2011 in Verona, Wisconsin (W11); 364 inbred lines during summer 2013 in Arlington, Wisconsin (W13); and 223 inbred lines during summer 2017 in Pendleton, South Carolina (S17). The maximum quantum efficiency of PSII photochemistry (Fv/Fm), a reliable indicator of monocarpic senescence (Wingler et al., 2004, 2010), was recorded in dark-adapted, field-grown plants at 36 d after anthesis (DAA). We found substantial phenotypic variation ranging from complete senescence (Fv/Fm = 0) to lack of senescence (Fv/Fm = 0.82), with the latter signifying stay-green genotypes (Supplemental Figure 1; Supplemental Table 1). The genotypic variance was significant for all three environments and for the combined W11_13 analysis (Table 1). Consistent with the significant impact of environment on senescence, the genotype ☓ environment component was significant for the combined W11_13 analysis. The broad sense heritability of Fv/Fm observed in each environment was 37.3% (W11), 72.6% (W13), 41.2% (W11_13), and 43.7% (S17).
Identification of Loci Associated with Senescence by Exploiting Natural Allelic Diversity
To identify the genes and genetic elements underlying senescence, we tested association of 438,161 RNA-seq single-nucleotide polymorphism (SNP) markers reported for the maize diversity panel (Hirsch et al., 2014) with senescence, recorded as Fv/Fm, for individual environments (W11, W13, and S17) and for combined WI environments (W11_13) (Supplemental Figure 2). In total, 62 unique SNPs were identified to be associated with senescence in one or more environments, of which five were highly significant (P ≤ 5.45E−9), 11 were significant (P ≤ 2.73E−7), and 46 met a suggestive threshold (P ≤ 5.45E−6) based on a proposed criteria (see “Methods”; Table 1; Lander and Kruglyak, 1995; Li et al., 2012).
Based on the evaluation of a large collection of more than 2800 inbred lines, linkage disequilibrium in maize drops within ∼10 kb with an r2 of 0.2, but the extent of linkage disequilibrium is dependent on the germplasm included in the analysis and the region of the genome (Romay et al., 2013). Therefore, only the SNP with the lowest P-value within a 200-kb window (100 kb ± SNP position) is reported. Two SNPs are represented twice due to the presence of two candidate genes in the vicinities of each of these SNPs.
Total phenotypic variance (R2) explained by all SNPs ranged between 5.1 and 13.5%, with a median value of 6.4%. To investigate the role these SNPs in senescence, we defined the candidate regions represented by each SNP to span 100 kb on each flank of the SNP. The 200-kb window size is based on similar studies in maize (Li et al., 2013b, 2016; Diepenbrock et al., 2017) and aimed to facilitate the identification of the correct candidate gene associated with the SNP. A search for all annotated genes within these 200-kb windows surrounding the identified SNPs yielded 411 candidate genes, with a range of 1 to 17 candidate genes per window and a median of 6.5 (Supplemental Data Set 1). These genes were further scrutinized using a SysGen approach to identify the most likely candidates and their role in senescence.
Physiological and Metabolic Analysis of Senescence and Stay-Green
We chose two inbred lines for detailed analysis of senescence and stay-green: PHG35 and B73. PHG35, a Non-Stiff Stalk inbred line from Oh07-Midland/Iodent background, is a stay-green line (Popelka, 2012), while B73, a Stiff Stalk inbred line from Iowa Stiff Stalk Synthetic background, is a non-stay-green line. Photosynthetic activity, as indicated by Fv/Fm of dark-adapted ear leaf, started to decline in both inbred lines from 30 DAA, albeit the rate of decline was different (Figure 1A). B73 showed a sharper decline after 33 DAA, with complete loss of photosynthetic activity at 39 DAA, while PHG35 had very high activity up to 39 DAA but declined quickly thereafter. Both chlorophyll a (Chl a) and Chl b followed a similar trend, with a sharper decline in B73 after 33 DAA and a gradual loss in PHG35 (Supplemental Figure 3). Overall, the 33 to 39 DAA window represents a transition zone separating the non-stay-green and stay-green phenotype displayed by the two inbred lines (Figure 1A).
Physiological and Metabolic Differences between a Non-Stay-Green (B73) and a Stay-Green (PHG35) Maize Inbred Line during Senescence.
(A) Onset and progression of senescence deduced from the maximum quantum efficiency of PSII (Fv/Fm). Gray box shows the transition zone separating non-stay-green and stay-green phenotype.
(B) to (D) Accumulation of Glc (B), Fru (C), and Suc (D) in leaf.
(E) to (G) Accumulation of Glc (E), Fru (F), and Suc (G) in internodes. Numbers on x axis denote DAA. Orange and dark green lines in all graphs represent B73 and PHG35, respectively. The error bars at each data point indicate se obtained from three biological replicates. One and two asterisks denote the significance of differences at P ≤ 0.05 and 0.01, respectively.
Since sugar status of the source and the sink tissues has a major role in senescence, we examined the accumulation of major nonstructural carbohydrates (NSCs) in leaf and internode tissues of B73 and PHG35 during post-flowering development. Up to 33 DAA, accumulation of hexoses (Glc and Fru) was higher in B73 leaves, but this trend reversed at 36 DAA with PHG35 accumulating more hexoses (Figures 1B and 1C). Accumulation of Suc showed an opposite trend, with higher accumulation in PHG35 up to 33 DAA, indicating active transport to sink (Figure 1D). The content of hexoses and Suc showed a sharp decline at 42 DAA, reflecting advanced senescence and death of the leaf tissue. The hexose-to-Suc ratio was higher in B73 during early stages but showed a reversed trend in PHG35 before the onset of senescence of this inbred line (Supplemental Figure 3). In the internodes, hexose and Suc concentration remained mostly unchanged in B73 except for a slight increase in Glc at 42 DAA (Figures 1E to 1G). PHG35 had lower amounts of hexose and Suc at 33 DAA, indicating remobilization of sugars to grain, but showed increased accumulation of all three sugars at later stages, eventually accumulating significantly higher Glc and Suc compared with B73. To summarize, the pattern of hexose accumulation in source tissue (i.e., leaf) is consistent with the onset of senescence in both inbreds. Glc and Fru hyperaccumulated in leaf of B73 before the observed decrease in Fv/Fm at 36 DAA, while such hyperaccumulation was later in PHG35, consistent with a drop in Fv/Fm at 42 DAA. Finally, PHG35 had higher soluble sugars in internodes at the end of the season compared with B73, indicating stronger alternative sink activity in the stay-green inbred.
Characterization of Transcriptome Underlying Natural Senescence
To understand the dynamics of the transcriptome during the onset and progression of senescence, we performed RNA-seq on the leaf at the ear-bearing node at key developmental stages selected based on the physiological and metabolic data (Sekhon et al., 2012). We focused on the normally senescing B73 inbred line and divided the leaf lifespan into the nonsenescent phase represented by 9 DAA and the senescent phase from 15 to 42 DAA (Supplemental Data Set 2). By comparing the expression at 9 DAA to each of the stages representing the senescence phase, we identified 14,518 (33.1%) genes showing differential expression at one or more stages during senescence (Supplemental Data Set 3). Distribution of up- and downregulated genes at different stages revealed that, except at 15 DAA, more genes were upregulated than down-regulated during senescence (Figure 2A). These DE genes were grouped into 22 clusters based on the kinetics of expression that included six clusters with 4005 (27.6%) upregulated genes, six clusters with 3043 (21.0%) downregulated genes, four clusters with 1346 (9.3%) genes upregulated during phenotypically visible senescence phase, and two clusters with 782 (5.37%) genes downregulated during phenotypically visible senescence phase (Supplemental Figure 4A). The other four clusters included remaining genes with other expression profiles.
Transcriptome Dynamics Associated with Natural Senescence and Stay-Green Trait Deduced from Comparison of a Non-Stay-Green (B73) and a Stay-Green (PHG35) Inbred Line of Maize.
(A) Differentially expressed genes during the progression of natural senescence in B73.
(B) Differentially expressed genes underlying stay-green trait identified from comparison of B73 with PHG35.
(C) Genes unique to progression of natural senescence (B73_U), unique to stay-green (SGR_U), and common to both phenomena.
(D) Significantly overrepresented GOBPs for each of the categories described in (C).
Gene ontology (GO) enrichment analysis of the six upregulated and six downregulated clusters identified several were significantly enriched biological processes (GOBPs) and molecular functions (GOMFs; Supplemental Figure 4B). Notable upregulated GOBPs included autophagy, lipid metabolisms, carbohydrate metabolism, carbohydrate transport, and amino acid transport. Key downregulated GOBPs included chlorophyll biosynthesis chloroplast organization, multiple GOBPs related to photosynthesis and functioning of photosystem I (PSI) and PSII, and, surprisingly, abscisic acid biosynthesis. Notable upregulated GOMFs included amino acid transporters, kinases, glycosyltransferase activity, and transcription factor activity, while notable downregulated GOMFs included electron transport, chlorophyll binding, and carbon-nitrogen ligases. These results highlight the magnitude of functional transitions in the leaf during senescence and identify upregulated catabolic processes/activities and downregulated biogenesis processes.
Characterization of Transcriptome Underlying Stay-Green
To understand the regulation of stay-green, we compared the transcriptome of B73 (non-stay-green) and PHG35 (stay-green) at specific developmental stages relevant to the onset and progression of senescence. In total, 4713 (10.7%) of the genes were DE between these two genotypes at four selected stages (Figure 2B; Supplemental Data Set 4). The majority (70.7%) of the DE genes were observed at the 33 DAA stage, which precedes the appearance of visual symptoms of senescence in B73, but not in PHG35 (Figure 1A). Clustering of the DE genes based on expression kinetics produced one cluster (Supplemental Figure 5A), further supporting a synchronized transcriptional control of stay-green. Predominant enriched GOBPs included cell wall biogenesis, xylan biosynthesis, carbon fixation, regulation of autophagy, and lipid and amino acid metabolisms (Supplemental Figure 5B), while enriched GOMFs included abscisic acid binding, glycosyltransferase activity, carbohydrate and polysaccharide binding, and electron transport (Supplemental Figure 5C). Overall, this analysis suggests an important role for carbon fixation, transport, and incorporation of sugars into cell wall in the determination of stay-green trait in maize.
To understand the molecular activities specific to either progression of natural senescence or stay-green, and those shared between these two phenomena, we further examined the DE genes. Of 14,688 unique DE genes described in Supplemental Data Sets 3 and 4, 10,062 were unique to natural senescence, 170 were unique to stay-green, and 4456 were common (Figure 2C). Remarkably, response to brassinosteroids was one of the most significantly enriched GOBPs in stay-green set, and other notable terms included oligosaccharide synthesis, negative regulation of kinase activity, Trp metabolism, and histone Lys methylation (Figure 2D). The GOBPs enriched during natural senescence included photosynthesis, redox homeostasis, nucleosome assembly, carbon utilization, and aging. Finally, GOBPs enriched in the common group included defense response, lipid transport, glutathione metabolism, and xylan biosynthesis.
Coexpression network analysis of the 4713 DE genes underlying stay-green produced a predominant subnetwork with 663 unique genes (Figure 3; Supplemental Data Set 5). Highlighting a major role for transcriptional regulation in senescence, 45 genes (6.8%) in this subnetwork are putative transcription factors, with a majority (8 genes each) representing NAC and MADS families. Putative Arabidopsis orthologs of 109 of the genes (16.4%) in this subnetwork are DE during senescence and designated as SAGs (Li et al., 2017).
A Major Coexpression Subnetwork Underlying Stay-Green.
Coexpression analysis of transcriptome deduced from comparison of B73 and PHG35 inbred lines produced a major subnetwork (n = 663 genes). Diamonds represent transcription factors, inverted arrows denote putative maize orthologs of SAGs in Arabidopsis, and the circles represent the rest of the annotated genes. Blue shapes are the genes that were also independently detected by GWAS.
Identification of Candidate Genes and Networks
The data generated from this study and previously published data on the mechanisms underlying senescence allowed us to implement a SysGen framework for discovery and prioritization of candidate genes. We used multiple criteria: (1) location of a SNP identified in the GWAS analysis, (2) annotated genes within a 100-kb window on each side of the SNP (200-kb window total), (3) relative expression of the 200-kb window genes in B73 and PHG35 inbred lines, (4) presence of a candidate gene in the stay-green coexpression network, and (5) deduced role of a candidate gene in senescence. For criterion 5, annotation of the 200-kb window genes, if available, was used as a term to search the literature and, based on role of a gene or the gene family in literature and/or the predicted orthology to SAGs in Arabidopsis (Li et al., 2014, 2017), the gene was considered to have a possible role in senescence. In the absence of any evidence based on criteria 3 to 5, the gene closest to the identified GWAS SNP was reported.
Using this strategy, the 62 unique GWAS associations identified in the study were first linked to 64 candidate genes (Table 1; Supplemental Data Set 1). Two SNPs, one each on chromosome 2 and chromosome 6, each had two candidate genes, one based on putative role in senescence based on the literature archive and the other indicated by the presence in the coexpression network. For SNP number 11 (Table 1), a Cys protease encoded by maize insect resistance3 (mir3) was considered a candidate due to reported role of these proteases in senescence in other species (Velasco-Arroyo et al., 2016; Pružinská et al., 2017). Likewise, for SNP number 43a, ZmSWEET1b was deemed a candidate based on reported role a SUGARS WILL EVENTUALLY BE EXPORTED TRANSPORTERS (SWEET) gene in senescence (Zhou et al., 2014).
Of the 64 candidate genes, 48 (75%) are DE during progression of senescence in B73, or in the comparison of B73 with PHG35 and therefore underlie stay-green, or during both processes (Supplemental Figure 6). Thirty (46.9%) of the candidate genes have been annotated as SAGs in Arabidopsis (Li et al., 2014, 2017), indicating both a substantial overlap and unique mechanisms used by cereals during senescence. Remarkably, 14 (21.9%) of the candidate genes are supported by more than one form of evidence and thus provide a very high-quality set of targets for further investigations (Figure 4; Table 1). Notable candidates in this set include a trehalose-6-phosphate (T6P) synthase (trps13), a NAC transcription factor (nactf9), an abscisic acid–responsive Group 3 late embryogenesis abundant protein (mlg3), a Cys protease (mir3), a β-glucosidase (ZmBGLU42), two cell wall–related genes (ZmIRX15-L and ZmGUX1), two sugar transporters (ZmSWEET1b and ZmMST4), a cell wall invertase (incw4), and a mitochondrial protein (dek10). Of these 14 candidates, 8 genes represent processes or functions relevant to senescence but lack conclusive experimental support, and the rest of the genes are novel. Eight genes of the 14 high-quality set are supported both by GWAS and coexpression network analysis and thus represent a core set of candidates (Figure 4).
Systems Genetic Analysis of Maize Leaf Senescence by Combining GWAS, Coexpression Networks, Pairwise Epistatic Interactions, and Prior Information.
Outer arcs denote 10 maize chromosomes along with 200-kb window around 62 GWAS SNPs shown as yellow boxes. Numbers in the boxes correspond to the serial number of the SNPs and candidate genes listed in Table 1. Selected candidate genes lying within the windows are shown underneath the arcs. Genes in blue font indicate those detected both by GWAS and coexpression network analysis (indicated by bold italicized font in Table 1) and in black font represent those detected only by GWAS but have a perceived role in senescence. Fine gray lines represent coexpression network connections of candidate genes with other genes in the genome. Thick light blue lines denote epistatic interaction between pairs of GWAS SNPs. Chr, chromosome.
The coexpression network and epistasis analysis revealed several putative interactions among the candidate genes and regions (Figure 4). One of the candidate SNPs (no. 36 in Table 1) was epistatic to three GWAS SNPs, while five others (nos. 3, 5, 18, 29, and 59) had two epistatic interactions, which suggests genetic interaction among the candidate regions or genes identified by GWAS. For instance, SNP number 5 harbors a candidate FLOWERING LOCUS K gene known to control flowering via suppression of FLOWERING LOCUS C transcription in Arabidopsis (Ripoll et al., 2009) and since FLOWERING LOCUS C also regulates senescence (Wingler, 2011), FLOWERING LOCUS K is expected to have a similar role. Interestingly, SNP number 5 is epistatic to SNP number 18 that lies near a candidate gene encoding for molybdenum cofactor sulfurase implicated in the regulation of senescence by controlling synthesis of xanthine dehydrogenase (Brychkova et al., 2008). Thus, while both candidates likely have a direct role in senescence independent of each other, our analysis suggests that these genes may be acting together. Likewise, several interactions were identified by coexpression networks underlying the stay-green trait (see “Discussion”). To summarize, a SysGen framework combining multiple pieces of evidence provides several highly valuable targets for future investigations to generate a novel mechanistic understanding of senescence and stay-green in maize and other cereals.
Functional Validation of Candidate Genes
To examine the role of mir3 (Zm00001d002065), associated with a highly significant SNP (no. 11) on chromosome 2 (Figures 5A and 5B), we compared natural alleles of this gene in B73 and PHG35 that are polymorphic at 20 nucleotide positions throughout the coding region (Figure 5C). Five of these SNPs have a moderate-to-high predicted effect and one SNP has high predicted impact as this SNP creates a putative splice acceptor site (Figure 5C; Supplemental Table 2). For both inbred lines, expression of mir3 increased at 30 DAA, and B73 had higher transcript accumulation compared with PHG35 for most later developmental stages (Figure 5D). Cys protease activity of Mir3, which has a predicted molecular mass of 37.2 kD upon maturity and removal of signal peptide (Pechan et al., 1999), increased dramatically in B73 at 36 DAA as indicated by presence of an intense band of expected size in the presence of gelatin loss of the band upon addition of E-64 Cys protease inhibitor (Figure 5E). Importantly, the protease activity in PHG35 was delayed by 6 d as evident from significantly lower activity at 36 DAA but dramatic increase afterward such that the activity in PHG35 at 42 DAA was same as that in B73 at 36 DAA (Figures 5E and 5F). Overall, these data suggest that delayed Cys protease activity contributes to prolonged photosynthetic activity and stay-green in PHG35. Besides Mir3, several other Cys proteases are active in both inbred lines at multiple developmental stages as indicated by the presence of larger bands and disappearance of those bands upon addition of the inhibitor. A similar variety of Cys proteases have also been reported in senescing wheat leaves wherein proteases of 36 to 46 kD were associated with leaf senescence (Martínez et al., 2007). These data indicate that the onset of senescence is associated with natural allelic variation of mir3 and Cys protease activity in the leaves. Further studies will confirm the identity of Mir3 as the causal Cys protease in maize leaf senescence and stay-green.
Functional Analysis of Natural Alleles of a Cys Protease Encoded by mir3.
(A) Manhattan plot of GWAS for natural senescence. The horizontal line indicates the suggested threshold (P ≤ 5.45E−6), and black and gray arrows indicate a Cys protease and a β-glucosidase gene discussed in the text, respectively.
(B) Regional association plot showing SNP used to identify mir3. Chr, chromosome.
(C) Sequence polymorphisms in alleles of mir3 in B73 and PHG35 inbred lines. One polymorphism with high predicted impact due to the generation of a putative splice acceptor site is marked by a black oval.
(D) Expression of mir3 in B73 and PHG35. Error bars are the se based on three biological replicates.
(E) Cys protease activity in leaves of B73 and PHG35 during the onset and progression of senescence. The protease activity was assayed in gels containing gelatin as a substrate in the absence (top) and presence (bottom) of E-64 protease inhibitor. The expected Mir3 band is marked by an arrow.
(F) Relative protease activity quantified based on signal intensity. Error bars are se obtained from biological replicates. One and two asterisks denote the significance of differences at P ≤ 0.05 and 0.01, respectively.
To further investigate the role of mir3, we annotated C1A papain-like Cys proteases (PLCPs) in four plant species based on the presence of the characteristic papain domain (Beers et al., 2004; van der Hoorn et al., 2004). Genome-wide search identified 63, 48, 36, and 30 PLCPs in maize, sorghum, Arabidopsis, and rice, respectively. Phylogenetic analysis of these protein sequences revealed that the maize Mir3 formed a close subgroup with Arabidopsis RD21A (AT1G47128), indicating a close evolutionary relationship between these two proteins (Supplemental Figure 7). Sequence analysis of the10 PLCPs in this subgroup revealed that RD21A, which encodes for a Cys protease (Koizumi et al., 1993), shares 66% amino acid sequence identity with maize Mir3 for complete protein sequence and 84% sequence identity for the papain domain (Supplemental Figure 8).
Based on this predicted orthology, we performed functional analysis of RD21A in Arabidopsis. A T-DNA insertion resulted in a complete loss of RD21A transcripts in Arabidopsis leaves (Figures 6A and 6B; Supplemental Figure 9). The mutant plants showed delayed onset of senescence at 24 DAA as evident from yellowing of the leaves in the wild-type Columbia (Col), and in the wild-type segregants derived from plants heterozygous for the insertional allele (Figure 6C). The Col and wild-type plants showed significantly advanced senescence in second and fourth rosette leaves at 28 DAA, in the seventh rosette leaf at 31 DAA, and complete senescence in all three leaves at 37 DAA (Figures 6D to 6F). However, the mutant plants maintained significantly higher photosynthetic activity at 37 DAA albeit at a lower level compared with earlier developmental stages, indicating delayed onset and a slower rate of senescence. The genotypic, transcriptional, and phenotypic characterization of RD21A was further confirmed with Cys protease zymography. The RD21A encoded protease, which has a predicted molecular mass of 33 kD upon maturity and removal of the signal peptide (Pružinská et al., 2017), was undetected in Col, the wild type, and the mutant plants at 10 DAA (Figures 6G and 6H). However, while RD21A abundance increased dramatically in the Col and wild-type plants, the protein was barely detectable in the mutants. Overall, these data further support the role of the mir3 ortholog in stay-green.
T-DNA Insertional Mutation in the Putative mir3-Ortholog Leads to Loss of Protease Activity and Delayed Onset of Senescence in Arabidopsis.
(A) Gene structure of Arabidopsis Cys protease AT1G47128 (RD21A), showing exons (boxes), untranslated regions (open boxes), translated regions (dark shaded boxes), introns (lines connecting boxes), location of T-DNA insertion (vertical arrowhead), primers used for genotyping (horizontal arrows), and conserved peptidase C1A (papain) domain (hatched boxes).
(B) RT-qPCR data showing expression of RD21A in the mutant, wild-type segregants from the initial, heterozygous mutant in Col background (Wt.), and Col plants at two developmental stages in the combined tissue of leaf 4 and leaf 7. Error bars are se obtained from three biological replicates.
(C) Phenotypes of representative Col, wild-type segregant (Wt.), and mutant plants. Numbers in the white box in represent DAA.
(D) to (F) Fv/Fm of three individual rosette leaves with the number representing the location of the leaf in the rosette.
(G) Comparison of Cys protease activity in same tissue and genotypes used for RT-qPCR. The protease activity was assayed in gels containing gelatin as a substrate in the absence (top) and presence (bottom) of E-64 protease inhibitor with the expected RD21A band marked by an arrow.
(H) Relative protease activity quantified based on signal intensity. Error bars are se obtained from three biological replicates. One and two asterisks denote the significance of differences between and mutant and wild-type segregant (Wt.) at P ≤ 0.05 and 0.01, respectively.
A highly significant SNP (no. 53) on chromosome 8 lies in gene ZmBLGU42 (Zm00001d012159; Figure 7A) which putatively encodes for a β-glucosidase orthologous to rice Os1bglu4 (Rouyi et al., 2014). We found that, compared with PHG35, transcript levels of ZmBLGU42 were higher in B73 at 33, 39, and 42 DAA but identical in both inbred lines at 36 DAA (Figure 7B). While the expression data were not conclusive, the assay of β-glucosidase activity in leaves of these inbred lines in zymograms by staining with 4-methylumbelliferyl-β-d-glucuronide revealed differential activity displayed by the B73 and PHG35 alleles. The onset of activity was earlier in B73 at 36 DAA, while PHG35 showed the activity at 42 DAA (Figures 7C and 7D). These data suggest that differential posttranscriptional activity of ZmBLGU42 alleles underlies differences in senescence phenotype.
Functional Analysis of Natural Alleles of a GWAS Candidate Encoding for a β-Glucosidase.
(A) Regional association plot showing SNP used to identify the β-glucosidase gene ZmBLGU42 (Zm00001d012159).
(B) Expression of ZmBLGU42 in B73 and PHG35 based on RNA-seq data. Error bars are se obtained from biological replicates.
(C) Comparison of β-glucosidase activity in leaves of B73 and PHG35 at post-flowering stages covering the onset and progression of senescence. The β-glucosidase activity was assayed by staining with 4-methylumbelliferyl-β-d-glucuronide.
(D) Relative β-glucosidase activity quantified based on signal intensity. Error bars are se obtained from three biological replicates. One and two asterisks denote the significance of differences at P ≤ 0.05 and 0.01, respectively.
DISCUSSION
Through systematic genomic and transcriptomic analyses, implementation of SysGen framework, and experimental validations, we provide a broad model of senescence control and progression in a global staple crop. Sixty-two unique associations identified by the characterization of natural variation in a diversity panel of U.S. dent maize were linked to 64 candidate genes. These genes implicate diverse processes in senescence including proteolysis (a Cys protease), sugar signaling (a T6P synthase), sugar transport (two sugar transporters), unloading of sugars to sink (a cell wall invertase), and, remarkably, synthesis of cell wall polysaccharides. Most of the genes reported in this study are novel candidates that have not been experimentally validated despite their logical connection with senescence based on published literature. Half of the candidate genes have been designated as SAGs in Arabidopsis (Li et al., 2014, 2017), thus highlighting both mechanistic similarities and differences in monocot and eudicot senescence.
A small number of previous inheritance studies on maize senescence used biparental populations and identified valuable quantitative trait loci underlying the stay-green trait (Zheng et al., 2009; Wang et al., 2012; Belícuas et al., 2014). Our study expands this knowledge by evaluating the role of a larger number of alleles captured in a diversity panel encompassing a substantial portion of the genetic variation in U.S. midwestern dent maize.
The previous phenotyping has relied on leaf greenness produced by chlorophyll through green leaf area (Zheng et al., 2009; Wang et al., 2012), visual scoring of leaf greenness (Belícuas et al., 2014), or chlorophyll content (Trachsel et al., 2016). However, greenness can indicate either functional stay-green plants that maintain high levels of photosynthetic activity or visual stay-green plant that maintain high chlorophyll, but not the photosynthetic activity (Thomas and Howarth, 2000). For estimation of functional stay-green, we recorded maximum quantum efficiency of PSII photochemistry (Fv/Fm) in dark-adapted, field-grown plants. Aging and high chlorophyll content coupled with lower CO2 assimilation could cause oxidative damage to the leaf cells, and plants prevent such damage by systematic dismantling of photosynthetic apparatus and by nonphotochemical quenching (NPQ; Lu and Zhang, 1998; Wingler et al., 2004). The PSII centers are dismantled before PSI centers, resulting in a greater decrease in PSII activity in senescing leaves, while NPQ levels increase during this period (Ghosh et al., 2001; Wingler et al., 2004). Thus, NPQ, PSII activity of non-dark-adapted leaves (Fv′/Fm′), and PSII activity of dark-adapted leaves (Fv/Fm) are all reasonable parameters for measuring photosynthetic activity. Our approach to record Fv/Fm of dark-adapted plants, a widely used parameter to access photosynthetic activity and senescence (Maxwell and Johnson, 2000; Wingler et al., 2004, 2010; Bresson et al., 2018), provides a better estimation of photosynthetic activity and therefore senescence.
Three spatially and temporally distinct environments allowed us to capture a significant portion of the genetic architecture of senescence that may otherwise go unnoticed. Various environmental factors including heat, drought, nutrient status, and shading influence senescence (Woo et al., 2013; Schippers, 2015), and therefore lead to significant phenotypic variation across environments. Biotic stresses, mainly the prevalence of diseases that may also be accompanied by hypersensitive response and death of leaf tissue, may also affect senescence in a given environment. Variation in heritability across environments reflects the role of environment in regulation on senescence. Similar variation in heritability of Fv/Fm was observed in a rice diversity panel with the estimates in two distinct environments being 12% and 34% (Qu et al., 2017). Dramatic fluctuations in heritability across environment have been reported for various quantitative traits in maize (Asaro et al., 2016).
High levels of biological activity during the onset and progression of senescence are evident from a massive change in transcriptional activity in one-third of the annotated maize genes. By contrast, a much smaller fraction of genes (9.2%) underlie stay-green, indicating fewer and focused processes. The majority (92.4%) of these genes were DE at 33 DAA, indicating a swift and massive transcriptome changes that differentiate a stay-green and a normally senescing phenotype. Overrepresentation of transcripts related to cell wall biogenesis and, in particular, hemicellulose synthesis in the stay-green transcriptome indicates a role of plant cell walls as an alternative sink for storage of excess carbohydrates. Remarkably, identification of genes that are unique to the B73-PHG35 comparison, and hence likely determine presence or absence of stay-green in these two inbred lines, produced an even smaller number (170) of genes.
Response to brassinosteroids is a prominent GOBP enriched in this set, bringing a known but mechanistically less understood regulator of senescence to light. Impairment of brassinosteroids biosynthesis, as observed in det2 mutant, or signaling, as seen in bri1 and bes1 mutants, results in stay-green phenotype in Arabidopsis (Chory et al., 1994; Clouse, 1996; Yin et al., 2002). Interestingly, a putative maize ortholog of bri1 was also identified in our GWAS study (Table 1), providing additional support to the role of brassinosteroids in stay-green.
The role of epigenetic mechanisms, specifically histone Lys methylation, in stay-green is another important finding of this analysis. In Arabidopsis, histone demethylase activity of a JmjC domain-containing protein JMJ16 is responsible for suppression of age-dependent senescence (Liu et al., 2019). Coexpression network analysis of stay-green also provided some very intriguing insights. Substantial representation of transcription factors (6.8%) in the subnetwork demonstrates a major role of transcriptional regulation, particularly the NAC and MADS families, in specifying stay-green. Besides the genes that appeared in the SysGen analysis, this analysis provided other interesting candidates. For instance, a sorghum NAC ortholog of a coexpression subnetwork member nactf36 is a key determinant of pithy or juicy stems (Zhang et al., 2018) and hence likely plays a vital role in determining the sink strength and therefore senescence. These findings provide some very interesting and promising lines of investigation to be pursued in the future.
Identification of mir3, a putative Cys protease (Pechan et al., 1999) in our screen, highlights the role of proteolysis in senescence. Proteases are required for hydrolysis of ribulose-1,5-bisphosphate carboxylase/oxygenase and other chloroplastic proteins for recycling nitrogen and other nutrients (Guo et al., 2004; Breeze et al., 2011). A large number of proteases existing in plants respond to specific senescence-triggering factors (Martínez et al., 2007), and based on the activity of two natural alleles (B73 and PHG35) in zymograms, mir3 appears to have a key role in stay-green. T-DNA insertional knockout of mir3 ortholog RD21A in Arabidopsis also resulted in a modest, albeit significant, delay in leaf senescence in our study and another recent study (Pružinská et al., 2017). Suppression of a Cys protease HvPap-1 in barley also delayed abiotic stress-induced senescence, thus underscoring the importance of these proteases in exacerbating biotic stresses faced by crop plants (Velasco-Arroyo et al., 2016). Despite the well-established role of proteolysis in senescence, the extent to which proteases play a role in stay-green is poorly understood and, if anything, suppression of proteolytic activities is generally believed to accelerate senescence (Roberts et al., 2012; Thomas and Ougham, 2014). Could the observation that certain natural alleles of mir3 condition stay-green without any obvious adverse effects on plants fitness be due to the precise localization of this protein in the chloroplast? Future experiments involving Mir3-specific antibodies for identification of Mir3 targets and interacting partners, and transgenics in maize and other plant species will confirm the role of this protease in senescence and stay-green.
Key insights into the role of ZmBLGU42, a β-glucosidase orthologous to Arabidopsis BGLU42 (Zamioudis et al., 2014) and rice Os1bglu4 (Rouyi et al., 2014), in senescence comes from the early onset of glucosidase activity in the B73 compared with PHG35. The β-glucosidase activity also increases during Arabidopsis leaf senescence (Patro et al., 2014), although the exact role of this process remains to be determined. β-Glucosidase hydrolyzes cell wall oligosaccharides to release Glc and, during sugar starvation of senescing leaf cells, maintains energy homeostasis in these cells (Patro et al., 2014). Since cell walls remain intact until the late stages of senescence, hexoses generated by β-glucosidase may be important for energy-dependent senescence processes (Patro et al., 2014). Since hexoses act as signal molecules through the hexokinase1 (HXK1) pathway to control senescence (Nuccio et al., 2015; Wingler, 2018), β-glucosidase may have a role in sugar signaling. Indeed, leaf hexoses were increased in B73 compared with PHG35, while Suc levels remained stable, resulting in a higher hexose:Suc ratio during senescence, and this trend was delayed in PHG35 (Figure 1). Hexose levels also rise dramatically in senescing Arabidopsis leaves, while Suc levels remain stable (Quirino et al., 2001; Wingler et al., 2006). Plants have separate sensors for Suc and hexoses and any changes in Suc:hexose ratio initiate distinct signal transduction, leading to specific changes in transcription (Smeekens, 2000). Since the senescing leaves will produce less photoassimilate, the rise in hexose levels in senescing leaves instead of a decrease is inexplicable (Wingler and Roitsch, 2008). Cell wall hydrolysis by ZmBLGU42 encoded β-glucosidase may be at least one of the mechanism responsible for generating these hexoses.
A loss of carbon demand from the primary sink upon grain fill combined with the absence of a secondary sink(s) will also lead to accumulation of hexose in leaf cells as seen B73 compared with PHG35 at ∼33 DAA. Hexoses act as positive regulators of the hexose sensor HXK1, which promotes senescence either by suppressing cytokinin signaling or through the conversion of Glc to Glc 6-phosphate, which then inhibits Suc nonfermenting-1-related kinase1 (SnRK1), a suppressor of senescence (Zhang et al., 2009; Smeekens et al., 2010; Wingler, 2018). Our study supports a direct role of sugar-mediated signaling in stay-green by identification of trps13 that synthesizes T6P, a Glc disaccharide, from Glc 6-phosphate and UDP-Glc. T6P is a negative regulator of SnRK1 and required for the onset of senescence associated with high sugar accumulation in leaves (Wingler et al., 2012). T6P is crucial for maintaining optimum Suc levels in plants by controlling starch breakdown at the source in response to the sink demand, and by regulating Suc consumption and growth at the sink through SnRK1 and other unknown mechanisms (Figueroa and Lunn, 2016). Attenuation of T6P by expressing a rice T6P phosphatase in maize resulted in prolonged and higher photosynthetic activity (Oszvald et al., 2018). T6P is an important target for biotechnological manipulations to enhance crop productivity and stress tolerance, and our data linking this molecule to stay-green further highlight the need for future investigations to fully understand such roles.
Given the role of sugars in senescence signaling, transporting these sugars from leaf cells to an alternative sink, evident from higher sugars in PHG35 internodes, should weaken such signaling. Identification of two sugar transporters (ZmSWEET1b and ZmMST), a cell wall invertase (incw4), and a mitochondrial pentatricopeptide repeat protein (dek10) further supports a key role for sugar partitioning in senescence. Disruption of Suc transport from leaves often leads to senescence-like symptoms in maize as seen in sugar transport mutants tdy1, tdy2, and sut1 (Baker and Braun, 2007; Slewinski et al., 2009, 2012). Arabidopsis hexose transporter STP13, a putative ortholog of ZmMST4, is associated with senescence and nitrogen use efficiency (Nørholm et al., 2006; Schofield et al., 2009), and putative rice ortholog, OsMST4, transports hexoses generated by cell wall invertases (Wang et al., 2007).
Extracellular invertases hydrolyze and remove Suc from the vascular system into a sink tissue and therefore regulate senescence by determining sink strength, source–sink translocation of sugars, and intracellular hexose:Suc ratio (Balibrea Lara et al., 2004; Rolland et al., 2006; Ruan et al., 2010). Enhanced expression of a cell wall invertase encoded by incw2 (also know as miniature1) resulted in higher starch accumulation and grain yield in maize (Li et al., 2013a). While incw2 is expressed in the basal endosperm transfer layer connecting the maternal and the endosperm tissue, incw4 is expressed in both vegetative and reproductive organs (Kim et al., 2000), thus making this invertase a good candidate for unloading sugars to a vegetative sink. A coordinated induction of an extracellular invertase and a hexose transporter by cytokinin treatment in cell culture has been suggested to enhance the uptake of sugars by a sink (Ehness and Roitsch, 1997). dek10, a ubiquitously expressed protein in maize plants, controls sugar uptake by the grain sink by orchestrating the development of the basal endosperm transfer layer (Qi et al., 2017). Our study suggests a key role of the cell wall invertases, hexose transporters, cytokinin, and Dek10 nexus in regulating the sink strength of internode parenchyma or other vegetative tissues and thus enhancing the overall sink strength of plants, and further studies will explore such a role. Interestingly, dek10 also affects nucleosome assembly during senescence (Gal et al., 2015), thus warranting investigations into the role of epigenetic mechanisms in stay-green.
Support for a model that transport and storage of sugars in an alternative sink contribute to stay-green comes from the discovery of two hemicellulose synthesis genes (ZmIRX15-L and ZmGUX1) and a NAC transcription factor (nactf9) in our study. Putative Arabidopsis orthologs of ZmIRX15-L and ZmGUX1 are involved in the synthesis of the major hemicellulose xylan in the primary and secondary cell walls (Jensen et al., 2011; Bromley et al., 2013). Remarkably, Arabidopsis IRX15 (but not IRX15-L) has been shown to be coexpressed with a NAC transcription factor SECONDARY WALL-ASSOCIATED NAC DOMAIN2 that regulates genes involved in cellulose and hemicellulose synthesis (Hussey et al., 2011). The NAC transcription factors can also enhance or suppress senescence by regulating the synthesis of abscisic acid, an important regulator of senescence, and also contribute to transition from carbon-capture to the nitrogen-recycling phase (Liang et al., 2014; Thomas and Ougham, 2014). Such a role is also supported by GWAS candidate mlg3, a late embryogenesis abundant protein, which is regulated by abscisic acid (White and Rivin, 1995). Therefore, nactf9 appears to be one of the master regulators of stay-green, acting in conjunction with ZmIRX15-L, ZmGUX1, mlg3, and other candidate genes identified in our study, and future work will explore such a role.
While both GWAS and transcriptome analyses produce false positives, our SysGen approach, which relies on combining GWAS, coexpression networks, and prior body of knowledge on senescence, provided a set of high-confidence candidate genes underlying stay-green and senescence (Figure 4). The SysGen framework identified source–sink communication as one of the key mechanisms underlying stay-green trait as 9 of the 14 high-confidence candidate genes play an important putative role in sugar transport and/or signaling (Figure 8). Cell wall appears to have an important role in senescence both as a source of sugars that can act in signaling at the source (leaf) and as a sink for excess sugars accumulating due to inactivity of the primary sink (in internode and other similar organs). In such a model, excess hexoses accumulating in the leaf cell due to loss of primary sink after completion of gain fill and the breakdown of cell wall polysaccharides by β-glucosidase (ZmBLU42), signal senescence through hexokinase (HXK1), and T6P (trps13).
A Model for Regulation of Senescence by Source–Sink Communication.
Nine candidate genes (shown in blue font) identified from system genetic analysis are involved in sugar transport, signaling, and storage and thus likely regulate senescence. Suc generated via photosynthesis during the reproductive phase of the maize plant is destined for the primary sink (grain). Upon loss of sugar demand by the primary sink after completion of grain fill (indicated by a red cross mark), some of the candidate genes likely play role in transporting sugars to an alternative sink exemplified by stalk parenchyma or cob pith cells, while others are involved in the deposition of sugars into the secondary cell wall. Other candidate genes (shown within the light green oval) are likely important for regulation of senescence via sugar-mediated signaling. CK, cytokinin.
In stay-green genotypes, senescence signaling is attenuated by transport of Suc out of leaf cells by sugar transporters (ZmMST4, ZmSWEET1b) to alternative sinks. These sinks are activated by dek10, hydrolyzed by cell wall invertase (incw4), and incorporated as cell wall polysaccharides by the activity of a NAC (nactf9) and xylan synthesis genes ZmIRX15-L and ZmGUX1. Effective partitioning of sugars to alternative sinks achieved by the combined activity of these genes therefore is a major determinant of stay-green. The most significant finding of this study is however the novel candidate genes not placed in this model due to the lack of any experimental evidence. Further examination of the role of these genes in determining the sink strength and, in particular, the strength of alternative sink to allow diversion of photosynthate after grain fill will be valuable toward engineering stay-green genotypes with higher net energy content.
With ever-worsening climate for agriculture and decreasing availability of freshwater, the stay-green trait has proven to be effective in improving drought resistance. Another less explored but equally important advantage of stay-green is increased overall carbon yield in the form of sugars accumulating in multiple sinks. Extra carbon captured in stay-green could be diverted to stover for developing dual-purpose (grain/biomass) crops for food and biofuel and/or as animal feedstock. However, scant data are available on the molecular mechanisms that would enable such a diversion, and this study provides promising candidates. In certain agricultural production systems, stay-green has also been found to be an undesirable trait due to, for instance, high moisture in stover and grain (delayed dry down) at harvesting in stay-green maize hybrids grown in the midwestern United States. Identification of the genetic determinants of senescence holds the key to interventions such as genetically uncoupling stay-green from higher grain/stover moisture. Deeper insights into sources–sink aspects of senescence are also needed for understanding the determinants of annuality and perenniality of cereals. The detailed genetic architecture of senescence holds a great promise for both applied and fundamental aspects of senescence research. With an appreciable amount of conserved gene and genome functionality among cereals, the information presented here will potentially translate to other important grasses including sorghum and rice.
METHODS
Plant Materials and Phenotypic Data Collection on the Diversity Panel
A maize (Zea mays) diversity panel (Hirsch et al., 2014) was grown in three environments for phenotypic data collection on senescence (Supplemental Data Set 6). In summers 2011 and 2013, the panel was grown at the University of Wisconsin West Madison Agricultural Research Station in Verona and at the University of Wisconsin Arlington Research Station in Arlington, respectively. In summer 2017, the panel was grown at Clemson University Simpson Research Farm in Pendleton, South Carolina. At each location, the experiment consisted of two replications planted in a randomized complete block design. In each one-row plot, five to seven competitive plants were open pollinated and used for phenotypic data.
To record senescence at a uniform developmental stage, senescence data were recorded at 36 d after 50% of plants in a plot showed silk (elongated style) emergence. Judging by green color, chlorophyll content, and sugar accumulation in the plant, maize genotypes generally start to show the onset of senescence ∼30 DAA (Swank et al., 1982; Crafts-Brandner et al., 1984); thus, we recorded data at 36 DAA to ensure that a large portion of the variation is captured.
The maximum quantum efficiency of PSII photochemistry (Fv/Fm) was used as an indicator of senescence (Wingler et al., 2004, 2010). To ensure that the photosynthetic reaction centers are fully oxidized and available for photochemistry, data were recorded after allowing dark adaption for 2 h past sunset. Fv/Fm was recorded with a Pocket PEA fluorimeter (Hansatech Instruments) by exposing the leaves to a light intensity of 3500 µmol m−2 s−1 for 1 s. Data were recorded two nodes above the ear-bearing node on the leaf blade at the midpoint between the midrib and the leaf edge ∼12 inches from the collar (joint of leaf sheath and blade).
GWAS Analysis
Best linear unbiased predictions of the phenotypic values were calculated for each genotype in the individual environments (2011, 2013, and 2017) using the mixed linear model capability of JMP Pro v13.2.0 (SAS Institute). Best linear unbiased predictions were also calculated for the 2011 and 2013 environments combined since those were grown in similar geographical locations. The environments are named W11, W13, W11_13, and S17. GWAS was performed for all these environments using mixed linear model (Yu et al., 2006) implemented in Genome Association and Prediction Integrated Tool (Lipka et al., 2012).
A minor allele frequency cutoff of 0.01 was used. In total, 438,161 polymorphic markers were used for the GWAS analyses; however, since dense marker panels likely exhibit high levels of linkage disequilibrium, the effective number of markers for significance testing was determined by Genetic Type I Error Calculator (Li et al., 2012). This analysis showed that the effective number of markers is 183,449.
Three significance thresholds were calculated according to a Bonferroni correction based on the effective number of markers. Based on Lander and Kruglyak (1995), we classified SNPs according to three P-value thresholds. The three thresholds, referred to as suggestive, significant, and highly significant, are associated with a false positive result being expected to occur 1.0, 0.05, and 0.01 times, respectively, in a genome scan. The reason for the three thresholds is to attempt to control for both Type I (false positives) and Type II (false negatives) errors. The suggestive threshold allows an increased overall chance of false positives but greatly reduces the overall chance of false negatives, thereby reducing the chance of completely missing an important region. We follow the recommendation (Lander and Kruglyak, 1995) that markers with a P-value less than the suggestive threshold are only speculative and require further study. The highly significant threshold allows an increased overall chance of false negatives but greatly reduces the overall chance of false positives. The significant threshold attempts to balance both false positives and false negatives. The threshold P-values for this study were 5.45E−6 for suggestive, 2.73E−7 for significant, and 5.45E−9 for highly significant.
The SNP set used for GWAS is based on maize RefGen_v2, and the coordinates of the significant SNPs reported here were converted to RefGen_v4 using CrossMap v0.2.7 (Zhao et al., 2014). All SNP variants in the diversity panel were evaluated for their potential effect on gene function using snpEff (Cingolani et al., 2012).
Tissue and Phenotypic Data Collection for Transcriptomic Experiments
For transcriptomic analyses, inbred lines B73 and PHG35 lines were grown in a randomized complete block design with three replications (i.e., three blocks) at Clemson University Simpson Research Farm in Pendleton in summer 2017. In each block, a six-row plot of each inbred line was grown with row length and row-to-row distance being 4.57 and 0.762 m, respectively.
In each six-row plot, the primary and all subsequent ears of all plants were covered with shoot bags before silk emergence. Upon silk emergence in at least 80% of the plants in a plot, shoot bags were removed on the same day to allow open pollination. Supplemental hand pollinations were performed the same day to ensure maximum seed set, and this day was designated as 0 DAA. Only those plants that showed silk emergence at 0 DAA were included in the experiment.
Two tissues, leaf and internode, were collected at 3-d intervals from 0 to 45 DAA. At each stage and for each inbred line, three biological replicates of each tissue were collected. A biological replicate consisted of pooled tissue from two different, competitive (nonborder and well-spaced from neighbors) plants harvested from the same block and dissected to separate the leaf and internode tissue. The leaf samples consisted of a 30-cm section of the leaf blade at the primary ear-bearing node excluding ligule midrib and sheath. The internode samples consisted of a 10-cm section of the internode above the primary ear-bearing node excluding the nodes on either end. Samples were chopped into small pieces, flash frozen in liquid nitrogen, and stored at −80°C. All samplings were performed between 8:00 am and 12:00 pm.
Senescence data, recorded as Fv/Fm, were collected on 7 to 15 randomly chosen plants from the three blocks with the same approach as described above in the Methods section for the phenotypic data collection on the diversity panel.
RNA-Seq and Data Analysis
RNA extraction and sample preparation were performed as described previously (Sekhon et al., 2013). mRNA isolation, cDNA synthesis, and library preparation were performed by Novogene following the standard Illumina protocol. Sequencing on the Illumina HiSeq platform generated 150-bp paired-end reads. Reads were trimmed of low-quality bases and adapter sequences with the trimmomatic software package (Bolger et al., 2014). Reads were aligned to the reference maize genome (Zea_mays.AGPv4.40.gtf) (Jiao et al., 2017) with the Tophat2 v2.1.1 package (Trapnell et al., 2012) and Bowtie2 v.2.3.4.1 (Supplemental Table 3; Langmead and Salzberg, 2012).
One biological replicate (replicate 1) was dropped for both B73 and PHG35 at 42 DAA stage for all subsequent analyses due to a low number of reads resulting from the turnover of RNA at this late developmental stage. Fragments per kilobase million (FPKM) were determined with Cufflinks v.2.2.1 (Trapnell et al., 2012). The overlap of sequence reads with genes (exons) was determined with the HTSeq-count v.0.10.0 software package (Anders et al., 2015). Counts for each replicate and genotype were normalized to remove sample-specific effects with EdgeR (Robinson et al., 2010).
Relative pairwise changes in gene level expression were computed using the EdgeR package. Genes with a P ≤ 0.05 were corrected by the false discovery rate (FDR) ≤ 0.05 (Benjamini and Hochberg, 1995). To identify DE genes during natural senescence, log2 (fold change [FC]) was first calculated by subtracting the log2(FPKM) in B73 at 9 DAA from all other time points in B73 and, at each stage, genes with FDR-corrected P ≤ 0.05 and a log2(FC) of ≥ 1 were considered to be DE. The unique DE genes for each stage were combined to develop a list of all genes involved in natural senescence for all downstream analyses. To identify DE genes underlying stay-green, log2(FC) in B73 relative to PHG35 was calculated by subtracting log2(FPKM) in PHG35 at a stage from log2(FPKM) in B73 at that stage, and DE genes at each stage were those with P ≤ 0.05 and a log2(FC) of ≥ 1. All unique DE genes underlying stay-green were used for subsequent analyses.
For GO enrichment, candidate genes above the assigned thresholds of FC were clustered into functional pathways with the GOSeq software tool (Young et al., 2010) and filtered by P ≤ 0.05. JMP Pro v.13.0 was used to evaluate candidate genes via k-means clustering (Wang et al., 2011). Circular figures were prepared with the Circos plotting tool (Krzywinski et al., 2009). Genes within 200-kb region around the GWAS SNPs were identified with the BEDTOOLS package (https://bedtools.readthedocs.io/en/latest/).
Coexpression networks
Coexpression networks were constructed with DE genes underlying stay-green using our random matrix theory based approach (Luo et al., 2007). We filtered out genes that have FPKM value of zero in more than half of 40 samples as well as genes with maximum FPKM value > 1.0 across all. After filtering, we have 4474 DE genes for the coexpression network construction. The random matrix theory based method output a coexpression network with 856 genes represented in 48 subnetworks with size greater than two genes.
Estimation of Epistatic Interactions
All possible two-way interaction terms (additive by additive epistasis terms) were included in a model along with all the marker main effect terms (additive terms) chosen in the initial GWAS screening. The interaction terms to include in the final model were selected using the least absolute shrinkage and selection operator methodology (Tibshirani, 1996) implemented in JMP Pro v.13.2.0 (SAS Institute). The value of the Wald chi-square statistic for each included interaction term was used to determine the width of the connecting line in the network graph. Only those interactions with the Wald chi-square ≥ 1.0 are presented.
Annotation of Cys Proteases
A Pfam domain pepetidase_C1 (papain C1A, Pfam 00,112) was used as a query in EnsemblPlants (http://plants.ensembl.org/index.html) Biomart tool (Kinsella et al., 2011) to identify the PLCPs in maize, sorghum (Sorghum bicolor), rice (Oryza sativa), and Arabidopsis (Arabidopsis thaliana). The retrieved sequences were annotated by National Center for Biotechnology Information Conserved Domain Database (Marchler-Bauer and Bryant, 2004) to confirm the presence of papain_C1A. The protein alignments for a final set of 177 curated PLCPs were performed by ClustalW implemented in MEGA7 (Kumar et al., 2016). The phylogenetic tree was inferred using the neighbor-joining method (Saitou and Nei, 1987), and evolutionary distances were calculated using Poisson correction method (Zuckerkandl and Pauling, 1965). The Bootstrap consensus tree inferred from 1000 replicates was taken to represent the evolutionary history (Felsenstein, 1985).
Cys Protease Zymography
Protease zymography was based on previously published protocols (Wagstaff et al., 2002; Martínez et al., 2007) and a soluble protein extraction buffer (Eason and De Vré, 1995). Briefly, each zymogram consisted of an SDS-PAGE gel with an 11% (v/v) separating gel containing 0.05% (v/v) gelatin (Type I from porcine skin, Sigma-Aldrich) and a 4% (v/v) stacking gel. Soluble protein (5 µg) was loaded in each lane in an equal volume representing 3 volumes of protein extract:1 volume of 4× protein loading buffer (100 mM Tris-HCl, pH 6.8, 20% [v/v] glycerol, 2% [w/v] SDS, and 0.04% [w/v] bromophenol blue) and electrophoresed at 100 V. Following electrophoresis, gels were washed for 1 h at room temperature in wash buffer (75 mM sodium-acetate, pH 4.5, 2% [v/v] Triton X-100, and 8 mM DTT) on an orbital shaker at 50 rpm to remove SDS. Gels were then developed for 16 h in incubation buffer (75 mM sodium-acetate, pH 4.5, and 8 mM DTT) at 37°C, washed briefly (3 × 5 min) in reverse osmosis water, stained for 1 h in LabSafe Gel Blue (G Biosciences), and imaged using a ChemiDoc MP Imaging System (Bio-Rad). To confirm Cys protease activity, gels were developed for 16 h in incubation buffer containing 20 µM E-64 Cys protease inhibitor and then processed in the same manner as done after the elecetrophoresis.
To quantify relative levels of Cys protease activity from the gel images, the peak area for the Cys protease was determined using ImageJ software (https://imagej.nih.gov/ij/docs/guide/index.html). The percentage of activity per sample was determined by calculating each time point and expressed that as a percentage of the total measured activity (peak area per time point/total measured peak area), allowing a comparison between the two genotypes. For both maize and Arabidopsis, three biological replicates were used for zymography, and each biological replicate consisted of pooled leaf tissues from two plants.
Arabidopsis Analysis
A mutant allele (stock: CS351853) of RD21A (AT1G47128) carrying a T-DNA insertion in the fourth exon of the gene in Col background was obtained from the Arabidopsis Information Resource (www.arabidopsis.org) and confirmed by sequencing and RT-PCR (Figures 6A and 6B; Supplemental Figure 9). Arabidopsis plants were subjected to 4 d of cold treatment and later transferred to a growth room at 22°C with 16-h day and 8-h night cycle. DNA from young leaves was isolated as described previously (Kim et al., 2016). Genotyping was performed using gene-specific forward (GSF, 5′-GGACCCTATGTTTCATTCAGTTTCA-3′), gene-specific reverse (GSR, 5′-TCAAAATTTTAGGCGGGATGACT-3′), and left border (LB, 5′-TAGCATCTGAATTTCATAACCAATCTCGATACAC-3′) primers with three combinations (GSF+GSR, LB+GSR, and GSF+GSR+LB) as described previously (O’Malley et al., 2015). PCR was performed in 20-μL reaction volume with reaction conditions as initial denaturation at 95°C for 3 min, 32 cycles of 30 s at 95°C, 30 s for 58°C and 1 min at 72°C, and the final extension at 72°C for 5 min.
For time-sensitive analyses, key plant developmental stages and specific leaves were identified using the published description of Arabidopsis ecotype Col development (Boyes et al., 2001). We used DAA, the number of days after the appearance of first flower (also known as stage 6), for comparative analysis. Total RNA was isolated from the second rosette leaf after 10 and 20 DAA in two biological replicates using TRIzol (Thermo Fisher Scientific) according to the manufacturer’s instructions. The first-strand cDNA synthesis was performed using SensiFAST cDNA synthesis kit (Bioline) using 1 µg of RNA. RT-qPCR was performed using SensiFAST kit (Bioline) according to manufacturer’s instructions using GSF and GSR primers. Normalization was based the Arabidopsis housekeeping gene ACTIN (AtACT-RT_F-GTTGGGATGAACCAGAAGGA; AtACT-RT_R-GAACCACCGATCCAGACACT). FCs were calculated using 2^(−ΔΔCT) method.
The onset of senescence was based on Fv/Fm data on one representative leaf from second, fourth, and seventh rosette using a light intensity of 3500 µmol m−2 s−1 for 1 s. DAA was used as a benchmark to ensure data collection at the same developmental stage in all plants. Data were collected on six plants belonging to each of the genotypic classes.
Annotation of Maize SAGs
A list of Arabidopsis SAGs was obtained from the leaf senescence database (Li et al., 2014), and the maize orthologs were identified as reported previously (Jiao et al., 2017) and accessed from MaizeGDB (Andorf et al., 2016).
Metabolic Analyses
NSCs were extracted from 50 mg of freeze-dried powdered leaf and internode tissue by adding 800 μL of 0.7 M perchloric acid (Tobias et al., 1992). Samples were incubated on ice with shaking for 20 to 30 min and centrifuged for 15 min at 12,000 rpm at 4°C. The supernatant was transferred to a new tube, neutralized to pH 7.0 by adding 240 μL of solution containing 2 M KOH and 0.4 M MES, and centrifuged for 10 min at 12,000 rpm at 4°C. The clear supernatant was transferred to a new tube and used for NSC analysis. Glc, Fru, and Suc were measured in the sequential reaction according to a published method (Zhao et al., 2010), with minor modifications. A Glc reagent containing 1 mM ATP, 1 mM NAD, 1 mM MgCl2, 2 U/mL hexokinase, and 4 U/mL Glc-6 phosphate dehydrogenase was added for Glc measurement, followed by addition of 2 U/reaction phosphoglucose isomerase for Fru measurement and, finally, the addition of 25 U/reaction of invertase for Suc measurements. The incubation time after adding enzymes was 30 min at 30°C for each reaction, and absorbance was taken at 340 nm. Purified Glc, Fru, and Suc (Sigma-Aldrich), each with at 0.5 mg/mL concentration, were used as standards with sequential dilutions. The amounts of NSC were estimated as Glc equivalents after comparison with the standards. Chl a and Chl b were measured and quantified following a previously published protocol (Lichtenthaler, 1987), with minor modifications as reported by Sekhon et al. (2012). For all metabolites, analyses were performed on three biological replicates.
Statistical Analysis
The metabolites, protease activity, and physiological data were analyzed using two-tailed Student’s t test using JMP Pro 14 (SAS Institute), and statistical significance was based on P ≤ 0.05. The values are presented as means ± se of the biological replicates.
Accession Numbers
Sequence data from this article can be found at the National Center for Biotechnology Information Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA505891.
Supplemental Data
Supplemental Figure 1. Distribution and correlation of Fv/Fm data across environments.
Supplemental Figure 2. Manhattan and quantile-quantile plots for GWAS analysis.
Supplemental Figure 3. The senescence-related decline in the content of chlorophyll a and chlorophyll b in leaves of B73 and PHG35 inbred lines.
Supplemental Figure 4. Dynamics of the transcriptome during the onset and progression of natural senescence in B73 inbred line of maize.
Supplemental Figure 5. Transcriptome dynamics associated with stay-green trait deduced from comparison of a stay-green (B73) and a stay-green (PHG35) inbred line of maize.
Supplemental Figure 6. Expression pattern of 48 differentially expressed GWAS candidate genes.
Supplemental Figure 7. Phylogeny of the papain-like Cys proteases (PLCP) in maize, sorghum, Arabidopsis, and rice.
Supplemental Figure 8. Multiple sequence alignment of predicted protein sequences belonging to a subgroup of 10 closely related PLCPs highlighted in Supplemental Figure 7.
Supplemental Figure 9. T-DNA insertion in Arabidopsis Cys protease AT1G47128 (RD21A) resulted in the loss of transcription.
Supplemental Table 1. Descriptive statistics and heritability of senescence measured as maximum quantum efficiency of photosystem II photochemistry (Fv/Fm) in the maize diversity panel.
Supplemental Table 2. SNP variants in maize Cys protease gene mir3 present in B73 and PHG35 and predicted consequences of these variants on gene activity.
Supplemental Table 3. Mapping summary of RNA-seq data.
Supplemental Data Set 1. List of candidate genes with 200-kb window around each selected GWAS SNP.
Supplemental Data Set 2. FPKM of biological replicates of B73 and PHG35 at various developmental stages included in the study.
Supplemental Data Set 3. Differentially expressed genes during at least one of the developmental stages covering natural senescence.
Supplemental Data Set 4. Differentially expressed underlying stay-green identified based on differential expression B73 compared with PHG35.
Supplemental Data Set 5. Genes represented in a major co-expression subnetwork underlying stay-green trait.
Supplemental Data Set 6. Inbred lines evaluated in the study.
Dive Curated Terms
The following phenotypic, genotypic, and functional terms are of significance to the work described in this paper:
Acknowledgments
We thank Jakyung Yi, Shannon Alford, Taylor DeHart, and many talented undergraduate student helpers for assistance in the collection of data and samples and Chris McMahan for helpful discussions. Special thanks to Yalan Huo for help in preparing graphics. We also thank anonymous reviewers for valuable suggestions to improve this article. This work was supported by Clemson University and by the U.S. Department of Energy Great Lakes Bioenergy Research Center, Department of Energy Office of Science (BER-DE-FC02-07ER64494).
AUTHOR CONTRIBUTIONS
R.S.S. and S.M.K. designed the research; R.S.S., R.K., B.S.F., A.J.A., and M.W.B. performed the research; R.S.S., C.S., F.L., T.M.B., W.C.B., and N.d.L. analyzed the data; and R.S.S., C.S., R.K., and S.M.K. wrote the article.
Footnotes
- Received December 10, 2018.
- Revised May 22, 2019.
- Accepted June 19, 2019.
- Published June 25, 2019.
References
