Comparative Transcriptome-Based Mining of Senescence-Related MADS, NAC, and WRKY Transcription Factors in the Rapid-Senescence Line DLS-91 of Brassica rapa

Leaf senescence is a developmental process induced by various molecular and environmental stimuli that may affect crop yield. The dark-induced leaf senescence-91 (DLS-91) plants displayed rapid leaf senescence, dramatically decreased chlorophyll contents, low photochemical efficiencies, and upregulation of the senescence-associated marker gene BrSAG12-1. To understand DLS molecular mechanism, we examined transcriptomic changes in DLS-91 and control line DLS-42 following 0, 1, and 4 days of dark treatment (DDT) stages. We identified 501, 446, and 456 DEGs, of which 16.7%, 17.2%, and 14.4% encoded TFs, in samples from the three stages. qRT-PCR validation of 16 genes, namely, 7 MADS, 6 NAC, and 3 WRKY, suggested that BrAGL8-1, BrAGL15-1, and BrWRKY70-1 contribute to the rapid leaf senescence of DLS-91 before (0 DDT) and after (1 and 4 DDT) dark treatment, whereas BrNAC046-2, BrNAC029-2/BrNAP, and BrNAC092-1/ORE1 TFs may regulate this process at a later stage (4 DDT). In-silico analysis of cis-acting regulatory elements of BrAGL8-1, BrAGL42-1, BrNAC029-2, BrNAC092-1, and BrWRKY70-3 of B. rapa provides insight into the regulation of these genes. Our study has uncovered several AGL-MADS, WRKY, and NAC TFs potentially worthy of further study to understand the underlying mechanism of rapid DLS in DLS-91.


Introduction
Leaf senescence is a complex, programmed series of events involving various physiological and molecular processes, including chlorophyll degradation and a decline in photochemical efficiency [1,2]. Senescence in plants, which occurs through nutrient mobilization and recycling [3][4][5], is induced by biotic factors, such as aging and pathogen attack, and by abiotic factors, including extreme temperatures, nutrient deficiency, insufficient light, darkness, and drought [6].
Chinese cabbage (Brassica rapa ssp. pekinensis) is a leafy vegetable that generally grows best from late summer to early winter but that can also be cultivated year-round depending on the genotype. Chinese cabbage is readily affected by various external environmental factors that lead to multiple nutrient losses, cell death-induced leaf yellowing, weight loss due to reduced relative humidity, and decomposition [7,8]. In Chinese cabbage, the leaf is the economically most important part of the plant; consequently, leaf senescence is responsible for major crop yield losses. Nevertheless, studies on leaf senescence in Chinese cabbage are limited or lacking [9]. To study senescence, dark conditions are often used because they induce rapid, simultaneous senescence in detached leaves, similar to that observed during normal plant development [10,11]. In addition, many functional senescence associated genes (SAGs) are involved in some of the same mechanisms during both age-triggered and dark-induced senescence [12].
NGS technologies, such as transcriptome sequencing and microarray-based analysis of global gene expression, have been successfully applied for understanding the molecular mechanisms of leaf senescence in different crop species. For example, Lu et al. [25] identified 881 transcription factors (TFs), mainly belonging to bHLH, MYB, ERF, MYB-related, NAC, and WRKY families, in the seasonal leaf-senescence transcriptome of black cottonwood. In sunflower, a comparative transcriptome study uncovered genotype-specific genes and TFs related to leaf senescence [26]. In monocot plants such as maize, members of CO-like, NAC, ERF, GRAS, WRKY, and ZF-HD TF families have been found to be involved in early leaf senescence [27].
Multiple studies have described the functions of leaf senescence-associated TFs, which mainly include NACs and WRKYs along with a few MADS family members. For example, ORE1/ANAC092 plays crucial roles in aging leaf and flower tissues [28,29], and the NAP2 gene positively regulates leaf senescence by directly binding to SAGs in tomato plants [30]. In addition, Oda-Yamamizo and colleagues have shown that ANAC046 promotes leaf senescence in Arabidopsis by activating Chl catabolic genes and SAGs [31]. Other studies have been conducted on the role of NAC TFs in the regulation of leaf senescence [4,32].
At present, only a few MADS-box TFs, which are primarily known as flowering developmental regulators [44], have been functionally associated with leaf senescence. A few studies have identified MADS-box genes with functions related to leaf senescence. For instance, overexpression of the FYFL MADS-box gene delays leaf and sepal senescence in tomato [45]. A global gene expression analysis during different leaf senescent stages revealed the differential expressions of AtAGL8 (AT5G60910) and AtAGL20 (AT2G45660) [46]. Likewise, Fang et al. [47] have reported that overexpression of AGL15 delays SAG12 promoter activity and slows the senescence of sepal, petal, and silique tissues in Arabidopsis. Furthermore, TFs characterized as having a role in various trans-regulation, signaling, metabolic, and hormonal responses of senescence can be found in the LSD leaf senescence database (https://bigd.big.ac.cn/lsd/index.php, accessed on 17 March 2021) [48].
Leaf senescence is a complex process arising from cell death in leaf tissues. In our initial study reported here, we adopted different experimental approaches to identify factors causing leaf yellowing in Chinese cabbage.
We selected two genotypes exhibiting different senescence during dark-treatment senescence: DLS-42 (used as a control) and the rapid-senescence line DLS-91. To identify the cause of early senescence in DLS-91, we mainly focused on physiological factors, such as chlorophyll content and photochemical efficiency, as well as possible genetic factors identified by RNA-seq analysis.

DLS-91 Exhibited More Rapid Leaf Senescence than the Control Line during Dark Treatment
Because yellowing is the most distinct visual phenotype associated with leaf senescence, we used leaf color to compare different stages of the dark-induced leaf senescence process. To identify phenotypic changes occurring in B. rapa under dark treatment, we selected two inbred lines, the control line DLS-42 and the rapid-senescence line DLS-91, according to criteria specified in Materials and Methods. Similar to the protocol used in an earlier study in Arabidopsis [49], fourth leaves were detached from the two genotypes and incubated in darkness for up to 9 days ( Figure 1). Leaf yellowing was first observed on the third day of dark treatment (phase I) in DLS-91 and on the sixth day in DLS-42 ( Figure 1). By day 9, leaves of DLS-91 were entirely yellow; in contrast, the leaf color of DLS-42 on the final day was similar to that of DLS-91 on 5-DDT ( Figure 1). Therefore, to further detect leaf senescence symptoms, we selected 0-, 3-, and 7-DDT samples based on leaf color and analyzed chlorophyll content and photochemical efficiency (Fv/Fm) ( Figure 1). After 3 DDT, the reduction in the chlorophyll content of DLS-91 (ca. 40%) was higher than that of DLS-42 (ca. 25%); by day 7, chlorophyll was totally depleted from DLS-91 (ca. 97%) compared with the control (ca. 71%) ( Figure 2B). Similarly, phytochemical efficiency began to drop at 3 DDT and was completely reduced in DLS-91 leaves compared with those of DLS-42 at 7 DDT ( Figure 2C). According to these phenotypical and physiological parameters, leaf senescence occurs earlier in DLS-91 than in DLS-42 (Figures 1 and 2A,B). In addition, analysis of the expressions of six senescence-related genes [32,50], namely, Ethylene Insensitive 2 (BrEIN2), Senescence-associated gene 12 (BrSAG12-1 and BrSAG12-2), and ORESARA 1 genes BrORE1-1 and BrORE1-2, revealed that BrSAG12-1 was highly differentially expressed at 2 DDT between DLS-91 and DLS-42 ( Figures 2C and S3     , and (C) expression of the senescence-related marker SAG12-1 (with EF1-α as a reference gene). Data is the mean of three independent experiments. Error bars represent standard deviation of three biological replicates. Different letters indicate significant differences. One-way ANOVA test is performed using SPSS software. DAT, days of dark treatment.

Assessment of DLS-42 and DLS-91 RNA Sequence Data and Identification of Stage-and Genotype-Specific Genes
In this study, the first to address dark-induced leaf senescence at the transcriptomic level in B. rapa, we identified the complete set of differentially expressed genes (DEGs) during dark treatment in DLS-91. Previous transcriptomic analyses of leaf senescence in Figure 2. Phenotypic changes during dark-induced leaf senescence. Variation in (A) chlorophyll content, (B) photochemical efficiency (F v /F m ), and (C) expression of the senescence-related marker SAG12-1 (with EF1-α as a reference gene). Data is the mean of three independent experiments. Error bars represent standard deviation of three biological replicates. Different letters indicate significant differences. One-way ANOVA test is performed using SPSS software. DAT, days of dark treatment.

Assessment of DLS-42 and DLS-91 RNA Sequence Data and Identification of Stage-and Genotype-Specific Genes
In this study, the first to address dark-induced leaf senescence at the transcriptomic level in B. rapa, we identified the complete set of differentially expressed genes (DEGs) during dark treatment in DLS-91. Previous transcriptomic analyses of leaf senescence in Arabidopsis used 0-3-DDT leaf tissues and were able to identify important genes involved in leaf senescence [49,51]. In contrast, we carried out RNA sequencing, with three biological replicates, of 0-, 1-, and 4-DDT leaves of DLS-42 and DLS-91 ( Figure 1), which allowed genes involved in early and late stages of DLS to be identified in the rapid-senescence line. On average, 18 to 20 million raw reads were generated from both genotypes (DLS-42 and DLS-91). Q20 and Q30 percentages across all reads were more than 97% and 94%, respectively, with a sequencing error rate of 0.06%, and the average mapping rate for all samples of DLS-42 and DLS-91 lines was between 87% and 89%. Other important aspects of the sequencing data are statistically summarized in Table S1.

Functional Enrichment Analysis of DEGs
The genes and TFs found to be differentially expressed at each treatment stage were subjected to gene ontology (GO) annotation using the agriGO database (http:// systemsbiology.cau.edu.cn/agriGOv2/, accessed on 17 March 2021). The top 10 GO terms from each of three categories (biological process, molecular function, and cellular component) were used to generate the graphs shown in Figure 5. GO terms were assigned to 362, 419, and 413 DEGs from 0-, 1-, and 4-DDT samples, respectively, and are listed in Table S6. In the biological process category, most DEGs were related to "metabolic process", followed by "cellular process" and "organic substance metabolic process". Similarly, most genes related to molecular function were involved in "catalytic activity", followed by "binding" and "organic and heterocyclic compound binding". Under the cellular component category, many DEGs were localized to "cell", "cell part", and "intracellular" components ( Figure 5A-C). TFs at the three leaf senescence stages were associated with the 84, 77, and 72 GO terms listed in Table S7. Under the biological process category, most TFs detected at 0 DDT were assigned to "metabolic process", followed by "cellular and cellular metabolic process". At 1 and 4 DDT, in contrast, TFs were enriched in 10 biological process terms, such as "regulation of biosynthetic process and regulation of cellular biosynthetic process". In regard to molecular function, most TFs were enriched in "binding" and "organic and heterocyclic compound binding" at all sampled stages. The most heavily enriched cellular component terms were "cell" "cell part" and "organelle" ( To identify pathways associated with genes upregulated and downregulated at each leaf senescence stage, we performed a KEGG pathway enrichment analysis on the entire set of DEGs (Table S8). As a result, genes upregulated at 0, 1, and 4 DDT were assigned to 67, 48, and 59 pathways, respectively. Most of these genes were involved in metabolic To identify pathways associated with genes upregulated and downregulated at each leaf senescence stage, we performed a KEGG pathway enrichment analysis on the entire set of DEGs (Table S8). As a result, genes upregulated at 0, 1, and 4 DDT were assigned to 67, 48, and 59 pathways, respectively. Most of these genes were involved in metabolic and secondary metabolite biosynthetic pathways, followed by glycolysis/gluconeogenesis and DNA replication-related pathways. In addition, a few genes differentially expressed at 0 and 1 DDT were related to terpenoid and chlorophyll biosynthesis. Similarly, several DEGs at 1 and 4 DDT were involved in fatty acid biosynthesis and degradation, while others up-or downregulated at 0 and 4 DDT were associated with amino acid degradation (Table S8). In addition, genes downregulated at 0, 1, or 4 DDT were found to participate in 61, 65, and 65 pathways, respectively, whereas genes expressed continuously at all three sampled stages were involved in pathways related to metabolism, secondary metabolite biosynthesis, and purine, carbon, and propanoate metabolisms (Table S8).

Transcript Abundance of Genes with Various Functions during Senescence
In general, leaf senescence comprises a series of molecular reactions involving various genes and TFs at different stages [52]. To excavate the entire set of senescence-associated genes (SAGs) in DLS-91, we assembled a comprehensive set of genes and TFs with functions related to leaf senescence in different plant species [1,4,48,51,[53][54][55][56][57][58][59], including leaf coloration, chlorophyll biosynthesis and degradation, photosystems I and II, phytohormone signal components, chromatin mediated regulation, transcription and post-transcription, translation and post-translation mediated regulators, leaf senescence and plant immunity, and genes that are upregulated and downregulated during DLS. In total, 3003 orthologous genes with a wide range of expression levels among different time frames were identified from the complete transcriptome data (Table 1 and Table S9).  [57,59] Regulators that act at the translational level 1 1.54 −0.90 −0.65 [59] Regulators that act at the post-translational level 25  Out of 55 chlorophyll-related genes identified from the transcriptome data, 44 were related to chlorophyll biosynthesis, and 11 were chlorophyll-degrading genes. Five of the chlorophyll biosynthetic genes, PORB, PORC, GSA1, HEMB1, and HEMF1, were highly expressed (log 2 FC > 1) at 0 DDT in DLS-91 compared with DLS-42, with PORC and PORB also highly expressed at 1 DDT and 4 DDT, respectively ( Figure 6A). Among chlorophylldegrading genes, CLH1 was highly expressed at 0 DDT in DLS-91 compared with the remaining two stages and with DLS-42. In contrast, NYE1 was more strongly expressed during later stages (1 and 4 DDT) in DLS-91 ( Figure 6A). In addition, expression levels of chlorophyll-degrading genes and the rate of chlorophyll content reduction in leaves were higher in DLS-91 than in DLS-42. Overall, these results indicate that chlorophyll biosynthetic and degradation functions were activated earlier in DLS-91 than in DLS-42 ( Figure 2A; Table S9).
phyll-degrading genes, CLH1 was highly expressed at 0 DDT in DLS-91 compared with the remaining two stages and with DLS-42. In contrast, NYE1 was more strongly expressed during later stages (1 and 4 DDT) in DLS-91 ( Figure 6A). In addition, expression levels of chlorophyll-degrading genes and the rate of chlorophyll content reduction in leaves were higher in DLS-91 than in DLS-42. Overall, these results indicate that chlorophyll biosynthetic and degradation functions were activated earlier in DLS-91 than in DLS-42 ( Figure 2A; Table  S9).
A total of 55 genes related to photosystems I and II were also identified from the transcriptome data. Among them, 11 genes-PSAH2, LHCA1, LHCA4, and PSAG of PSI and two LHCB6, two LHB1B1, two CAB1 and one LHCB3 gene(s) of PSII-were highly expressed (log2 FC > 1) at the late stage (4 DDT) in DLS-91 ( Figure 6B). In addition, the PSAH2 gene of PSI was highly expressed (log2 FC > 1) at all treatment stages (0-4 DDT), and BraAnng000560 of PSII exhibited stage-specific expression (at 1 DDT) in DLS-91 compared with DLS-42 ( Figure 6B). We also identified 29 genes related to phytohormones, auxin, abscisic acid, brassinosteroids, cytokinins, ethylene, abscisic acid, jasmonic acid, and salicylic acid. In regard to the three sampled stages, the largest number of genes (JAZ7-2, 2 SID2, SID2-2, 3 PAD4, PAD4-2, and PAD4-3) were more highly expressed in 1-DDT samples in DLS-91 compared with DLS-42 ( Figure 6C), with JAZ7-1 and BRI1 more strongly expressed at 0 DDT as well (Table S9). This result indicates that hormones such as jasmonic acid, abscisic acid, and brassinosteroids are produced early in DLS-91-at 0 DDT, whereas salicylic acid is produced during dark treatment (1 DDT) in this genotype. Consistent with these findings, Yin et al. [60] have suggested that SA delays leaf senescence by regulating concentrations of methyl jasmonates. In Arabidopsis, ABI5 has been determined to delay DLS [61]. In our study, the expression FC of the ortholog of ABI5/GIA1 (BraA05g009140) was higher at all sampled treatment stages in DLS-91. A total of 55 genes related to photosystems I and II were also identified from the transcriptome data. Among them, 11 genes-PSAH2, LHCA1, LHCA4, and PSAG of PSI and two LHCB6, two LHB1B1, two CAB1 and one LHCB3 gene(s) of PSII-were highly expressed (log 2 FC > 1) at the late stage (4 DDT) in DLS-91 ( Figure 6B). In addition, the PSAH2 gene of PSI was highly expressed (log 2 FC > 1) at all treatment stages (0-4 DDT), and BraAnng000560 of PSII exhibited stage-specific expression (at 1 DDT) in DLS-91 compared with DLS-42 ( Figure 6B).
We also identified 29 genes related to phytohormones, auxin, abscisic acid, brassinosteroids, cytokinins, ethylene, abscisic acid, jasmonic acid, and salicylic acid. In regard to the three sampled stages, the largest number of genes (JAZ7-2, 2 SID2, SID2-2, 3 PAD4, PAD4-2, and PAD4-3) were more highly expressed in 1-DDT samples in DLS-91 compared with DLS-42 ( Figure 6C), with JAZ7-1 and BRI1 more strongly expressed at 0 DDT as well (Table S9). This result indicates that hormones such as jasmonic acid, abscisic acid, and brassinosteroids are produced early in DLS-91-at 0 DDT, whereas salicylic acid is produced during dark treatment (1 DDT) in this genotype. Consistent with these findings, Yin et al. [60] have suggested that SA delays leaf senescence by regulating concentrations of methyl jasmonates. In Arabidopsis, ABI5 has been determined to delay DLS [61]. In our study, the expression FC of the ortholog of ABI5/GIA1 (BraA05g009140) was higher at all sampled treatment stages in DLS-91.
Genes shown to regulate senescence at the RNA/protein level (transcription, translation, and post-translation) in a previous investigation [56] were highly expressed in DLS-91 during at least two sampled stages in our study ( Figure 6D-F). In addition, a large number of genes related to transcriptional regulation during leaf senescence were highly expressed at 4 DDT in DLS-91 samples, which indicates that these genes might be involved, even in this genotype, in regulating senescence by controlling the transcription of their targets [56].
In addition to the above findings, the TFs shown in Figure 6G have been reported to play a role in leaf senescence and plant immunity [54]. Most of these orthologous genes showed enhanced expression after dark treatment, thus indicating that these genes were activated at early stages of DLS in DLS-91.
Totally 194 genes are downloaded from the LSD leaf senescence database included genes related to senescence and leaf color, most of which had their maximum expressions at 0 DDT (Table S9). Using age-triggered and dark-induced leaf senescence genes previously identified from the Arabidopsis transcriptome [1], we identified 2528 orthologs from our transcriptome data. Similar to previous expression data, the genes in our study exhibited a broad range of expression levels at different stages and between DLS-42 and DLS-91. In addition, these genes were up-or downregulated in a genotype-dependent manner (Table S9).

qRT-PCR Validation of MADS, NAC, and WRKY TFs
To further validate gene transcript abundances, we selected 16 genes (seven AGL, six NAC, and three WRKY genes) (Table S11) on the basis of their expression levels (FPKM values) and roles in leaf color change and leaf senescence ( Table 2). The transcript abundance of BrAGL8-1 was at least seven times higher at 0, 1, and 4 DDT in DLS-91 compared with DLS-42, and a similar expression pattern was observed in the qRT-PCR analysis ( Figure 7A). Two paralogs of BrAGL15 (BrAGL15-1 and BrAGL15-2) had contrasting expression patterns based on the RNA-seq data. Likewise, the qRT-PCR analysis indicated that the expression of BrAGL15-1 was genotype specific (i.e., limited to DLS-91) at all stages, whereas BrAGL15-2 was slightly more highly expressed in DLS-91 at 0 DDT than at the other two sampled stages ( Figure 7A and Table 2). Similarly, BrAGL20-2 and BrAGL20-3 paralogs of BrAGL20 exhibited genotype (DLS-42)-specific expression at all stages according to the transcriptome data, whereas BrAGL20-1 was slightly more highly expressed at 0 DDT in DLS-42 and at 1 and 4 DDT in DLS-91. The overall results of the qRT-PCR analysis of BrAGL20 genes were well correlated with the transcriptome expression data ( Figure 7A and Table 2). Finally, BrAGL42-2 exhibited higher transcript abundance in DLS-42 than in DLS-91 at all stages, and qRT-PCR also revealed its DLS-42-specific expression ( Figure 7A and Table 2). The expression patterns of AGL genes uncovered in this study suggest that these genes have completely genotype (DLS-42/DLS-91)-specific functions during DLS.

qRT-PCR Validation of MADS, NAC, and WRKY TFs
To further validate gene transcript abundances, we selected 16 genes (seven AGL, six NAC, and three WRKY genes) (Table S11) on the basis of their expression levels (FPKM values) and roles in leaf color change and leaf senescence ( Table 2). The transcript abundance of BrAGL8-1 was at least seven times higher at 0, 1, and 4 DDT in DLS-91 compared with DLS-42, and a similar expression pattern was observed in the qRT-PCR analysis ( Figure 7A). Two paralogs of BrAGL15 (BrAGL15-1 and BrAGL15-2) had contrasting expression patterns based on the RNA-seq data. Likewise, the qRT-PCR analysis indicated that the expression of BrAGL15-1 was genotype specific (i.e., limited to DLS-91) at all stages, whereas BrAGL15-2 was slightly more highly expressed in DLS-91 at 0 DDT than at the other two sampled stages ( Figure 7A and Table 2). Similarly, BrAGL20-2 and BrAGL20-3 paralogs of BrAGL20 exhibited genotype (DLS-42)-specific expression at all stages according to the transcriptome data, whereas BrAGL20-1 was slightly more highly expressed at 0 DDT in DLS-42 and at 1 and 4 DDT in DLS-91. The overall results of the qRT-PCR analysis of BrAGL20 genes were well correlated with the transcriptome expression data ( Figure 7A and Table 2). Finally, BrAGL42-2 exhibited higher transcript abundance in DLS-42 than in DLS-91 at all stages, and qRT-PCR also revealed its DLS-42-specific expression ( Figure 7A and Table 2). The expression patterns of AGL genes uncovered in this study suggest that these genes have completely genotype (DLS-42/DLS-91)-specific functions during DLS.   Among the analyzed NAC genes, BrNAC029-2 was approximately five times more highly expressed in 0-and 1-DDT samples of DLS-42 than in DLS-91 but was more strongly expressed in 4-DDT samples of DLS-91 than those of DLS-41. The results of qRT-PCR validation of this gene were correlated with the transcriptome data ( Figure 7B and Table 2). The BrNAC046-2 gene was highly expressed at 1 and 4 DDT in DLS-91, and qRT-PCRbased relative expression levels of this gene were consistent with its transcript abundance ( Figure 7B and Table 2). Although transcripts of BrNAC082-3 accumulated at all stages in both genotypes, this gene was more highly expressed in DLS-42 than in DLS-91. In the qRT-PCR analysis, the expression of this gene gradually increased from 0 to 4 DDT in both genotypes, but the strongest expression was detected in DLS-42-similar to the transcriptome data ( Figure 7B and Table 2). Two duplicates of BrNAC90 (BrNAC090-1 and BrNAC090-2) displayed high stage-specific expression in DLS-91: BrNAC090-1 at 4 DDT and BrNAC90-2 at 1 DDT. In addition, BrNAC090-1 was 10 times more highly expressed in DLS-42 than in DLS-91 at 0 DDT, an observation consistent with the qRT-PCR results ( Figure 7B and Table 2). Moreover, BrNAC092-1 was more highly expressed at 0 and 1 DDT in DLS-42 and at 4 DDT in DLS-91, which was similar to the qRT-PCR data ( Figure 7B and Table 2). The varied expressions of BrNAC genes at 4 DDT in DLS-91 suggest that these genes are involved in DLS regulation in a stage-specific manner.
In regard to WRKY genes, BrWRKY6-1 FPKM-based read counts and qRT-PCR expression levels were high at 1 and 4 DDT in DLS-91. In contrast, BrWRKY45 displayed high transcript abundance and qRT-PCR expression at 0 and 1 DDT in DSL-42 and at 4 DDT in DLS-91. Finally, BrWRKY70-1 was more strongly expressed in DLS-91 than in DLS-42 at all stages according to both transcriptome and qRT-PCR data ( Figure 7C and Table 2). We therefore speculate that WRKYs regulate leaf senescence in both genotypeand stage-specific manners.

Genomic Analysis of MADS, NAC, and WRKY TFs
We analyzed evolutionary relationships, intron-exon structures, and conserved motifs of AGL, NAC, and WRKY genes as specified in Materials and Methods (Figure 8). In the phylogenetic tree of AGL genes, two sister relationships were uncovered: one between BrAGL20 genes and one between BrAGL15 genes ( Figure 8A). In terms of AGL intron-exon organization, the large intron region and intron phases 0 and 2 were found to be conserved. BrAGL8-1 and BrAGL42-2, which were distinct in the phylogenetic tree, have the highest number of introns ( Figure 8B). Five conserved motifs were identified: a SRF-TF motif, two K-box motifs, and two unknown or low-complexity motifs ( Figure 8C). Two sister relationships were apparent in the NAC gene phylogenetic tree: one involving BrNAC090 genes, and the other between BrNAC046-2 and BrNACL092-1 ( Figure 8D). All NAC genes except for BrNAC082-3 were found to have the same number of introns and similar intron phases (0 and 1); BrNAC082-3, which did not group with any other NAC genes in the phylogenetic tree, possesses three introns ( Figure 8E). The six NAC proteins contain four NAM motifs and one unknown or low complexity motif ( Figure 8F). One pair of WRKY genes, BrWRKY6-1 and BrWRKY45, were sister to one another in the WRKY phylogenetic tree ( Figure 8G). The three analyzed WRKY genes have different intron-exon distributions, different numbers of introns, and different intron phase patterns ( Figure 8I). The three predicted WRKY proteins contain two WRKY motifs and three unknown or low-complexity motifs ( Figure 8H).
( Figure 8F). One pair of WRKY genes, BrWRKY6-1 and BrWRKY45, were sister to one another in the WRKY phylogenetic tree ( Figure 8G). The three analyzed WRKY genes have different intron-exon distributions, different numbers of introns, and different intron phase patterns ( Figure 8I). The three predicted WRKY proteins contain two WRKY motifs and three unknown or low-complexity motifs ( Figure 8H).

Prediction of MADS, NAC, and WRKY Binding Sites in DEGs and Identification of Other Cis-Elements in Promoter Regions of Selected TFs
To identify the gene targets of MADS, NAC, and WRKY TFs, we analyzed the 2-kb upstream promoter region in all DEGs ( Figure 9A). As a result, we identified 361 genes containing a putative MADS-box TF-binding site (CArG-box). Some of these genes, such as
We also searched for other important cis-elements with various functions in the promoter regions of 16 selected MADS, NAC, and WRKY DEGs ( Figure 9B). We identified various combinations of cis-regulatory elements (CREs) and/or gene-specific CREs. The functionally most important of these CREs were related to light, abscisic acid, methyl jasmonates, gibberellin, auxin, salicylic acid, circadian clock, and other CREs (Table S13). Light-responsive elements (LREs) were also commonly conserved (Table S13). The highest number of LREs in the seven AGL genes was found in BrAGL8-1 and BrAGL42-2, most of which were Box 4 (ATTAAT) or G-box ([C/T] ACGT[T/G]) elements. Three AGL genes (BrAGL8-1, BrAGL15-1, and BrAGL20-1) contain a circadian motif (CAAAGATATC), with BrAGL15-1 additionally possessing some gene-specific motifs, such as the ATCT motif (AATCTAATCC) and the LAMP element (CCTTATCCA) (Table S13). In contrast, all NAC TF genes except for BrNAC090-1 and BrNAC092-1 have an equal number of light-and hormone-responsive elements. Interestingly, we observed some motifs (GTGGC motif and chs-CMA2a) that are conserved between BrAGL8-1 and BrNAC029-2. A few other LREs and circadian clock-related CREs are also conserved between NAC and AGL family TFs (Table S13). Finally, we found that the three WRKY TFs have different combinations of functional CREs: BrWRKY6-1 possesses only LREs, whereas BrWRKY45 has more hormonerelated elements than LREs, and BrWRKY70-1 has an equal number. Circadian clock-related CREs are also conserved in BrWRKY70-1. Information on other related conserved CREs is provided in Tables S12 and S13. Detailed information on depicted light-related, hormone-related, and other conserved motifs is provided in Tables S12 and S13.

Discussion
Various anabolic and catabolic processes take place during dark-induced leaf senescence. In this study, we compared two genotypes, the control DLS-42 and the rapid-senescence line DLS-91, and identified the physiological and molecular changes causing DLS-91 to more rapidly senesce during dark treatment. We determined the genetic factors-DEGs and TFs-through comparative transcriptome profiling and identified physiological characters by analyzing chlorophyll and photochemical contents. In regard to all these factors, DLS-91 exhibits a broader spectrum of variation compared with DLS-42. Our analysis of cis-elements also revealed that AGL genes contain both NAC-and WRKY-binding sites. Although many previous studies have focused on the role of NAC and WRKY TF genes, our results more strongly suggest the involvement of AGL genes in the rapid senescence of DLS-91. We thus speculate that AGL genes, along with NACs or WRKYs, regulate dark-induced rapid senescence in B. rapa DLS-91. The physiological and molecular parameters analyzed here are likely to be critical functional components involved in the rapid leaf senescence of DLS-91 [ Figure 9. Prediction of cis-elements in the 2-kb upstream promoter region of differentially expressed transcription factor genes. (A) Proportion of differentially expressed genes from the entire transcriptome data set containing MADS, NAC, and WRKY conserved motifs. The table details the binding sequence of each cis-element. (B) Schematic representation of conserved motifs identified among AGAMOUS-LIKE, NAC, and WRKY differentially expressed transcription factors. Detailed information on depicted light-related, hormone-related, and other conserved motifs is provided in Tables S12 and S13.

Discussion
Various anabolic and catabolic processes take place during dark-induced leaf senescence. In this study, we compared two genotypes, the control DLS-42 and the rapidsenescence line DLS-91, and identified the physiological and molecular changes causing DLS-91 to more rapidly senesce during dark treatment. We determined the genetic factors-DEGs and TFs-through comparative transcriptome profiling and identified physiological characters by analyzing chlorophyll and photochemical contents. In regard to all these factors, DLS-91 exhibits a broader spectrum of variation compared with DLS-42. Our analysis of cis-elements also revealed that AGL genes contain both NAC-and WRKY-binding sites. Although many previous studies have focused on the role of NAC and WRKY TF genes, our results more strongly suggest the involvement of AGL genes in the rapid senescence of DLS-91. We thus speculate that AGL genes, along with NACs or WRKYs, regulate dark-induced rapid senescence in B. rapa DLS-91. The physiological and molecular parameters analyzed here are likely to be critical functional components involved in the rapid leaf senescence of DLS-91 [46,47,49,54,56,65,85].

Physiological and Molecular Responses of DLS-91 during DLS
Leaf yellowing, a distinct initial symptom of leaf senescence, is caused by chlorophyll depletion [51]. This reduction of chlorophyll is a characteristic reaction in the breakdown of photoreceptors during light depletion [1,86]. In the present study, we compared the chlorophyll content and photochemical efficiency (Fv/Fm) of leaf samples of DLS-42 and DLS-91 at three different phases of dark treatment (Figure 2A,B). The rapid loss of chlorophyll and reduction in photochemical efficiency observed in DLS-91 may have been due to the breakdown of photosynthetic components leading to diminished photosynthesis during dark treatment [1,86]. Similarly, leaves of Arabidopsis and barley subjected to dark treatment have been reported to exhibit a decline in chlorophyll and the Fv/Fm ratio according to the progression of dark exposure [11,85]. Genes related to chlorophyll biosynthesis (PORC, PORB, GSA1, HEMF1, and HEMB1) and degradation (CLH1) were also highly expressed in DLS-91 before dark treatment (i.e., 0 DDT) and then displayed reduced expression, which indicates that the high rate of chlorophyll reduction resulted in the rapid senescence of DLS-91 (Figures 1, 2A and 5A). We then used the leaf-senescence markers EIN2, ORE1, and SAG12 to analyze DLS-91-specific expression patterns during dark treatment [1,32,50]. According to the expression analysis of orthologous genes, the expression of BrSAG12-1 at 2 DDT was higher in DLS-91 compared with DLS-42. These data were correlated with physiological parameters, thus indicating that BrSAG12-1 is a genotype-specific early leaf senescence regulator in B. rapa ( Figure 2C and Figure S3). Because SAG12 is controlled by developmental signals and not induced by various stress conditions in Arabidopsis, this gene is believed to be a reliable marker for natural leaf senescence [87,88]. As shown in Figure 2C, however, the BrSAG12 gene is strongly induced during dark treatment in B. rapa. Previous studies have also reported that high accumulation of SAG12 in early-senescent leaves is due to dark treatment [1,89]. Given the phenotypic, physiological, and molecular evidence, we therefore hypothesize that leaf senescence can vary according to genotype in B. rapa.

Roles of MADS, NAC, and WRKY TFs in Various Processes Related to Senescence
In general, TFs act as central regulators of various processes, including leaf senescence [4,89,90], but no transcriptomic studies focused on TFs related to leaf senescence in B. rapa have been reported. In the present study, we identified 848 genes, including 15.6% encoding TFs, that were differentially expressed in DLS-91 at different stages of DLS (Figure 3; Tables S4 and S5). Many TFs, including NAC (21.8%), WRKY (18%), and MADS (9.8%) TFs (Table S5), were differentially expressed, which indicates their significant roles in various physiological and biochemical processes during leaf senescence [89][90][91] and confirms the necessity of studying TFs regulating leaf senescence in B. rapa.
Our transcriptome analysis uncovered 10 AGL-MADS-box genes (BrAGL8, BrAGL15, BrAGL20, and BrAGL42) that are orthologs to Arabidopsis genes AtAGL8 [46], AtAGL15 [47], AtAGL20 [64], and AtAGL42 [65] that function in senescence according to mutant and transgene expression analyses. Among the three copies of BrAGL8, BrAGL8-1 had the highest transcript abundance (FPKM values) at all three sampled stages in DLS-91 compared with DLS-42. We further validated the expression of this gene by qRT-PCR and obtained results similar to those of the transcriptome analysis ( Figure 6 and Table S10). In Arabidopsis, the flowering time regulator AGL-8 is also expressed in rosette and cauline leaves [52,92,93]. Interestingly, the promoter regions of BrAGL8-1 and BrNAP/BrNAC029-2 both possess a common LRE (the GTGGC motif), and we also discovered putative NAC (AGATTCGT) and MADS-box (CCTTTTTTGG) binding sequences in the promoter regions of both genes (Tables S12 and S13). Previous investigations on SEPALLATA3 have revealed that this gene along with AGL-9 TF binds to the ANAC092 promoter, which suggests that NAC genes can control AGL-MADS TF regulation of various functions [29,94]. In the case of BrAGL15, BrAGL15-1 was only expressed at all stages in DLS-91 according to qRT-PCR and transcriptome data. In contrast, the expression of BrAGL15-2 did not differ between DLS-42 and DLS-91. The variation in the expressions of these two gene duplicates may be due to structural differences, as indicated by their phylogenetic positions as well as the arrangement of their motifs, introns, and exons (Figures 6 and 7A). Previous studies of Arabidopsis have shown that AGL15 increases the longevity of floral organs by delaying their senescence and that overexpression of this gene does not delay leaf senescence [47,95]. Consequently, BrAGL15-1 and BrAGL15-2 may have genotypeand non-genotype-specific functions during dark-induced senescence. The three paralogs of BrAGL-20/SOC1 (BrAGL20-1, BrAGL20-2, and BrAGL20-3) exhibited a wide range of expression patterns (Figures 6 and 7A). Among them, BrAGL20-2 was highly expressed in DLS-42 at all stages. Sequence analysis uncovered differences in the structure, phylogenetic placement, and first-and second-intron sizes of BrAGL20-2 compared with the other two paralogs. Interestingly, the BrAGL20-2 gene promoter region contains a putative WRKY binding sequence (Table S12). Arabidopsis SOC1 mutants display inhibited Chl degradation and retain maximum photochemical efficiency, which suggests that SOC1 negatively regulates the general senescence process [65]. BrAGL42 gene duplicates were highly specifically expressed in DLS-42 according to our qRT-PCR and transcriptome analyses. Interestingly, the BrAGL42-1 gene includes both NAC (ACGTATCT) and WRKY(AGTCAA) binding sequences in its promoter region. Arabidopsis AGL42/FYF mutants exhibit delayed floral senescence/abscission due to repression of ethylene responses [65]. Taken together, these results suggest that the B. rapa AGL orthologs BrAGL-20-2, BrAGL-42-1, and BrAGL-42-2 play crucial, genotype-dependent roles during DLS.
The NAC TF genes AtNAP/ANAC029 [96], ANAC046 [31], VNI1/ANAC082 [66], and ORE1/ANAC092 [80] have various functions related to leaf senescence in Arabidopsis. In the present study, all of the orthologs of these genes, except for BrNAC082-3, were more highly expressed in DLS-91 than in DLS-42 at late stages of dark treatment (4 DDT). In addition, BrNAC029/NAP and BrNAC092/ORE1 were highly expressed before (0 DDT) and at the early stage (1 DDT) of dark treatment in DLS-91, whereas BrNAC082-3 was more highly expressed in DLS-42 than in DLS-91 at all dark-treatment stages. According to our analysis of BrNAC082 sequences, BrNAC082-3 is phylogenetically distinct and has more introns than other NAC genes (Figures 6 and 7B). Similarly, BrNAC090-1 and BrNAC090-2 were highly expressed at 0 DDT in DLS-42 and at 1 DDT in DLS-91. The variation in the expressions of these two genes may be due to structural differences (Figures 6 and 7B). An earlier study of Arabidopsis mutants of ANAC017, ANAC082, and ANAC090 demonstrated that these "NAC troika" genes negatively regulate leaf senescence by controlling salicylic acid and reactive oxygen species responses [66]. As mentioned above, conserved binding sites of MADS-box and WRKY TFs are present in the promoters of these NAC genes. Furthermore, most of these genes are highly expressed in DLS-91 leaves at the late stage (4 DDT) of senescence (Table S10). These results indicate that NAC genes have a regulatory role in DLS at the late stage of senescence in DLS-91.
Among the three validated BrWRKY6 orthologs, the expression of BrWRKY6-1 was slightly elevated at 1 and 4 DDT in leaves of DLS-91 than DLS-42. A study involving an Arabidopsis mutant has shown that WRKY6 positively regulates senescence through interaction with the senescence-induced receptor-like kinase (SIRK) gene [76]. BrWRKY6 may function during dark-induced senescence in B. rapa, but further functional analysis is needed to support our finding of stage/genotype-specific expression. In contrast, the expression of BrWRKY45 varied only slightly between the two genotypes. As a critical component of the GA-mediated signaling pathway, WRKY45 functions to positively regulate age-triggered leaf senescence [35]. In regard to BrWRKY70, BrWRKY70-1 was highly expressed in DLS-91 at all sampled stages. Similar to BrWRKY70-1, BrWRKY70-2 was strongly expressed at all stages, whereas BrWRKY70-3, whose promoter contains two NACbinding sequences (ATGAATCT and AGGTACAT; Table S12), was only highly expressed at 1 DDT in DLS-91. Among the three paralogous WRKY54 genes, BrWRKY54-2 was highly expressed in DLS-91. Previous studies in Arabidopsis have shown that WRKY54 and WRKY70 work together to negatively regulate leaf senescence [78]. Finally, other WRKY genes were more strongly expressed in DLS-91 during dark treatment (1 and 4 DDT) (Table  S10). Although these results suggest that these genes also have major roles during DLS in DLS-91, additional functional studies are needed to support the expression data ( Figure 6 and Table S10).

Estimation of Chlorophyll Content and Selection of Plant Materials
A core collection of 177 B. rapa lines was cultivated in a growth chamber at 23 • C under a 16 h light photoperiod. The B. rapa seedlings were grown in soil for 4 weeks, and their fourth leaves were then detached for the dark treatment. The detached leaves were incubated for 1 to 9 days in Petri dishes wrapped in aluminum foil and containing a floral foam block soaked in distilled water ( Figure S1). For chlorophyll extraction, the detached leaves were subsequently incubated in 95% ethanol for 17 h in darkness. The extracts were then centrifuged at 12,000× g for 5 min at 4 • C. Chlorophyll was quantified by measuring absorbance at 645 and 663 nm, and chlorophyll content was calculated as the ratio of ([20.23 × A 645 ] + [8.023 × A 657 ]) g −1 fresh weight [1,51]. On the basis of their leaf chlorophyll content at 3 and 7 DDT (as explained in results) the control line DLS-042 and the rapid-senescence line DLS-091 were selected for further analysis ( Figure S2). At least four biological replicates were used for the chlorophyll estimation.

Photochemical Efficiency
To estimate DLS-42 and DLS-91 photochemical efficiencies, healthy, fully grown fourth leaves were detached from 4-week-old plants and subjected to dark treatment. A Handy PEA meter (Hansatech, Norfolk, UK) was used to measure Fv/Fm values of the dark-treated leaves according to the manufacturer's instructions. At least five biological replicates were used for the estimation of Fv/Fm values.

RNA Sequencing and Transcript Annotation
DLS-42 and DLS-91 plants were grown for 1 month, with fourth leaves then detached from three biological replicates and subjected to dark treatment. After 0, 1, and 4 DDT, total RNA was extracted from leaves of 18 samples (2 lines × 3 stages × 3 biological replicates) using an RNeasy Mini kit (Qiagen, Seoul, Korea) and sequenced on a Illumina Hi-seq2000 platform by Seeders, Inc. (Daejeon, Korea). The sequence data have been submitted to the NCBI database and are available under the accession number PRJNA723680. The resulting RNA sequences were analyzed according to Rameneni et al. [97] After raw reads were trimmed using the Dynamic-Trim and Length-Sort tools in the Solexa QA package [98], clean reads were assembled with Velvet v1.2.08 and Oases v0.2.08 [99]. The assembled transcripts were functionally annotated using BRAD [100] and KEGG databases. Bowtie2 v2.1.0 software was used for mapping of transcripts [101]. To calculate gene expression values, transcript normalization was performed with the R package DESeq [102].

Identification of DEGs and TFs and GO Annotation
DEGs were selected according to their log 2 FC values. In this study, genes with log 2 FC > 1 or log 2 FC < −1 were considered to be upregulated or downregulated, respectively [97]. For identification of TFs, we downloaded all 4127 B. rapa TFs in the Plant Transcription Factor Database (http://planttfdb.cbi.pku.edu.cn/, accessed on 17 March 2021) [103]. All of the assembled transcripts were compared against the downloaded TF data using BlastX software.
GO alignment [104] of the complete set of transcripts was carried out using the agriGO database. Graphs were generated based on default parameters, and the transcripts were classified into three functional categories: biological process, cellular component, and molecular function. To predict the functional characteristics of upregulated and downregulated genes, p-values computed from Fisher's exact test were used to identify significant associations between genes and GO terms, and a KEGG enrichment analysis was performed using the KOBAS online tool (http://kobas.cbi.pku.edu.cn/kobas3, accessed on 17 March 2021) [105] to obtain pathway annotations for each transcript.

Validation of Transcript Abundance by qRT-PCR
To validate gene expression levels estimated from the transcriptome data, we selected 16 functionally characterized genes exhibiting variation in expression. Total RNA was isolated using an RNeasy mini kit (Qiagen) from 0-, 1-, and 4-DDT leaves of DLS-42 and DLS-91. A DNAse kit (Qiagen) was used to digest genomic DNA. cDNA was synthesized using a ReverTra Ace-α kit (Toyobo, Japan) and used as a template for qRT-PCR analyses on a CFX96 Real-Time system (Bio-Rad, Hercules, CA, USA). The qRT-PCR analysis was conducted on 13 SAGs with three biological replicates and three technical replicates using the following cycling conditions: 95 • C for 3 min, followed by 39 cycles of 95 • C for 15 s and 58 • C for 20 s [97]. Relative expression levels were determined by normalizing the data according to the comparative 2 −∆∆Ct method [106] using actin and EF-1α genes as a reference.

Promoter Analysis
The 2-kb region upstream of the transcription start site of selected DEGs was extracted, and all conserved cis-regulatory motifs were predicted using the plantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/, accessed on 17 March 2021) [107].

Genomic Analysis
A multiple alignment of the 16 protein sequences was generated in Clustal X [108] and used to generate an unrooted tree by the distance-based neighbor-joining method in MEGA X [109] with 1000 bootstrap iterations. Conserved motifs were identified using the MEME suite v5.3.3 (https://meme-suite.org/meme/, accessed on 17 March 2021) [110]. Analysis of gene intron-exon structural variation based on genomic and coding sequences was conducted using the GSDS v2.0 online tool (http://gsds.gao-lab.org/Gsds_about.php, accessed on 17 March 2021) [111].

Conclusions
Here, we showed for the first time, to our knowledge, the possible DLS underlying mechanism (s) in B. rapa. According to phenotypic, chlorophyll content, and photochemical efficiency data, DLS-91 undergoes more rapid leaf senescence than DLS-42 during dark treatment. Analysis of transcriptome-and qRT-PCR-based expression data demonstrated that the DLS-91-specific expression of most AGL-type MADS TFs varies at different stages of leaf senescence. In addition, NAC and WRKY binding sequences are conserved in the promoter regions of AGL genes. The discovery of an expression pattern specific to the rapid-DLS line DLS-91 provides essential information for identifying genes involved in dark-induced leaf senescence in B. rapa. Further functional analyses of these candidate genes are required to reveal the molecular mechanism underlying DLS.

Conflicts of Interest:
The authors declare no conflict of interest.