Sequencing data processing and TF genes identificationThe global gene expression profiles of 20 tobacco leaf sample collected from different developmental stages were downloaded from the SRA database and were processed as the methods previously described by Li, et al.15. Following data processing, 272,920 unigenes were identified for subsequent analysis. In this study, TF prediction server (http://planttfdb.gao-lab.org/prediction.php) was used to predict the ability of identified unigenes to encode transcription factor based on their sequence characteristics. As a result, a total of 3617 transcription factor genes were identified to be expressed in tobacco leaf and classified into 56 TF families using the classification criterion in PlantTFDB (Fig. 1, Table S2). The number of TFs belonging to bHLH, B3, MYB_related, C3H, bZIP, ERF, WRKY, C2H2, FAR1, and NAC families was very large, and they were widely expressed in tobacco leaves at different senescence stages.Fig. 1The distribution of tobacco transcription factors in different families.Among, some members of TF families expressed at specific stages during tobacco leaf senescence. As shown in Fig. 2, 239 senescence-associated TF genes (21 early senescence-associated TF genes, 35 middle senescence-associated TF genes, and 183 late senescence-associated TF genes) belong to 35 TF families (Fig. 2 and Table S2). Members of the C3H, ERF, HD-ZIP, and WRKY families were found to be associated with tobacco leaf senescence at all three senescing stages, suggesting that TF genes in these families play an extensive regulatory role in the process of tobacco leaf senescence.Fig. 2The distribution of early, middle, and late senescence-associated transcription factors in different families.Transcription factors are involved in the growth and development of plants by regulating target genes. The expression levels of TF genes have an important impact on the degree of regulation. Therefore, the expression trends of transcription factor genes during leaf senescence are crucial to the study of the regulatory role of transcription factors during leaf senescence. As shown in Fig. 3, most identified early, middle, and late senescence response TF genes, respectively, exhibit relatively higher expression levels at their corresponding stages than at the other two senescence stages, indicating the potential regulatory roles of these senescence response TF genes in the early, middle, and late stages of tobacco leaf senescence. These TF genes may be associated with the variations of leaf morphology and physiological metabolism in the tobacco leaf senescence process.Fig. 3Trend plots of the normalized expression levels of early, middle, and late senescence-associated TF genes in different leaf samples. The expression value in the broken line graphs were log2-transformed FPKM values.Weighted gene co-expression networks construction, identification, and expression analysis of senescence-related TF genesTo investigate the important regulatory genes of tobacco leaf senescence, we screened for genes that were differentially expressed during the senescence of tobacco leaves. As UL, ML, and LL samples were all collected at 15 DAT (the starting time point), the gene expressions at 15 DAT were defined as the reference. All subsequent time points were compared to the starting point. As a result, a total of 64,964 differentially expressed genes (DEGs) were identified through pair comparison between different stages of UL, ML, and LL, under the threshold of fold change ≥ 2.0 and FDR ≤ 0.001 (Table S3). Of all the identified senescence-associated DEGs, 1,359 were TF genes, accounting for 30% of the identified TF genes expressed in tobacco leaves, indicating tight regulation of the transcriptional activity (Table S4). Although the systemic analysis provides sufficient information about genes differentially expressed during tobacco senescence, the large number of identified DEGs resulted in the lack of specificity in distinguishing key senescence-related genes. To further excavate the co-expressed modules and key genes associated with leaf senescence, we constructed a co-expression network based on the identified DEGs. As a result, DEGs were clustered into 21 modules (module1-module 21), with modules ranging in size from 61 to 8876 unigenes (Table S5). Each module represented genes with highly correlated expression profiles. The correlation coefficients between modules and developmental stages were shown in Fig. 4. According to the senescence phenotype of the leaf and the expression feature of senescence marker genes (CP1, SAG12 homolog in tobacco) in all samples, the upper and middle leaves collected at 15 DAT, 25 DAT, and 35 DAT were considered to be in the early senescence stage. The upper and middle leaves collected at 45 DAT, 55 DAT, 65 DAT, and lower leaves collected at 15 DAT and 25 DAT were considered to be in the middle senescence stage. The upper and middle leaves collected at 75 DAT, 85 DAT, and lower leaves collected at 35 DAT and 45 DAT were considered to be in the late senescence stage. We found that most samples at the particular senescence stage were associated with the same module. The most relevant modules of samples at different senescence stages were different from each other. Therefore, the co-expression analyses also highlight notable differences among different leaf senescence stages. By module-trait correlation analysis, relative to other modules, 6 (modules 2, 3, 5, 6, 8, 9), 8 (modules 1, 4, 11, 12, 13, 14, 15, 16), and 6 (modules 7, 10, 17, 18, 19, 20) modules were found to be the most relevant modules with leaf samples at early, middle, and late senescence stage, respectively. 21, 35, and 183 TF genes were clustered into modules associated with early, middle, and late senescence stages, respectively. Hence, these TF genes were regarded as early, middle, and late senescence-associated TF genes.Fig. 4Network analysis of tobacco leaf senescence. (A) Gene clustering tree (dendrogram) obtained by hierarchical clustering of adjacency-based dissimilarity. The colored row below the dendrogram indicates different modules identified by WGCNA. (B) Correlations between module eigengenes and leaf sample. Biological traits (rows) were regrouped in 3 categories (early senescence trait, middle senescence trait and late senescence trait). In the heatmap, positive correlation and negative correlation were represented by red color and blue color.Gene ontology and pathway analysis of TF genes associated with different senescence stagesTo characterize the functional changes of tobacco leaf senescence-associated TFs, we conducted GO annotation of these TFs (Fig. 5). Based on the functional annotation of senescence-associated TF genes, we assigned 21 early senescence-associated TF genes, 35 middle senescence-associated TF genes, and 183 late senescence-associated TF genes with 39, 46, and 196 GO terms, respectively. These genes were grouped into three main GO categories: molecular function, cellular component, and biological process (Fig. 5A, Table S2). Among them, the most abundant GO terms assigned by these senescence-associated TF genes in the cellular component, molecular function, and biological process category were nucleus (GO:0005634), sequence-specific DNA binding transcription factor activity (GO:0003700), and regulation of transcription, DNA-dependent (GO:0006355), indicating that TFs in tobacco leaf at different senescence stages were active. In the biological process category, many senescence-associated TFs were correlated with leaf development, suggesting that these TFs are key regulators of leaf development (Fig. 5B). According to the GO enrichment analysis results, six GO terms, including transcription regulator activity (GO:0140110), organic cyclic compound binding (GO:0097159), heterocyclic compound binding (GO:1901363), binding (GO:0005488), biosynthetic process (GO:0009058), and nitrogen compound metabolic process (GO:0051171) (corrected p-value < 0.05) (Fig. 5C, Table S6), were significantly enriched by early, middle, and late senescence-associated TF genes.Fig. 5GO analysis of senescence-associated TFs. (A) GO classifications of early, middle, and late senescence-associated TF genes. (B) The enriched GO terms (corrected P-value ≤ 0.05) annotated by senescence-associated TF genes. Each dot indicates a GO term. The dot sizes indicate the number of genes. Colors indicate the significance of the GO term. (C) The results of enrichment analysis of GO terms annotated by senescence-associated TFs in biological process category are shown as REVIGO scatterplots in which similar GO terms are grouped in arbitrary two-dimensional space based on semantic similarity. Each dot indicates a GO term. The dot sizes indicate the number of genes. Colors indicate the significance of the GO term.To investigate which metabolic pathways the identified TFs were associated with, we obtained the pathway annotations of TF genes using the KEGG database23,24,25. As a result, we found that early, middle, and late senescence-associated TF genes were associated with 10, 2, and 13 metabolic pathways, respectively (Fig. 6, Table S7). According to the KEGG enrichment analysis results, no pathways were enriched by the early senescence-associated TF genes. Middle senescence-associated TF genes were enriched in the plant pathogen interaction pathway. Pathways of circadian rhythm plant, plant hormone signal transduction, MAPK signaling pathway plant, and plant pathogen interaction were enriched by late senescence-associated TF genes (Fig. 6, Table S7). These data suggest that complex regulatory mechanisms underlie tobacco leaf senescence via TF genes.Fig. 6Barcharts represent the KEGG pathway annotated by early, middle, and late senescence response TF genes. X axis is the normalized corrected p-value. Y axis is the KEGG pathway. The KEGG pathways enriched by early, middle, and late senescence response TF genes were marked by the red dot.Transcription factor regulatory network in senescent tobacco leavesTranscription factors are indispensable regulatory factors in the life activities of higher plants. Some transcription factor genes may regulate the expression of many key genes associated with important metabolic pathways. Therefore, identifying the target genes of senescence-associated transcription factors is important for unraveling regulatory networks in senescing tobacco leaves. PlantRegMap is a database that integrates a series of regulatory data, including regulatory information of transcription factors and other gene regulatory elements obtained from relative references and experimental data. Based on the information of TF-TF-target interactions in tobacco extracted from the PlantRegMap database, 797, 756, and 4552 unigenes were found to be regulated by early, middle, and late senescence-associated TF genes, respectively (Table S8). To explore the dynamics of TF regulatory networks across tobacco leaf senescence, we constructed early, middle, and late senescence-associated TF gene regulatory networks using Cytoscape software22. (The early, middle, and late senescence-associated TF gene regulatory networks were shown in Figure S1). Hub genes, which are strongly related to numerous genes, have been shown to play important regulatory roles in gene expression networks. Based on the degree of the node, 2, 1, and 6 hub genes in the early, middle, and late senescence-associated TF regulatory networks, respectively, were identified. Genes with high values of degree (> 300), including Unigene34594 (HSF family member), Unigene29760 (ERF family member), CL11135.Contig1 (WRKY family member), CL11180.Contig2 (TCP family member), CL20581.Contig2 (C3H family member), Unigene47914 (HSF family member), CL7305.Contig4 (HSF family member), Unigene47579 (CAMTA family member), and CL17451.Contig2 (ERF family member), were observed as hub genes in the network.To further excavate the TF regulating effects on the tobacco leaf senescence, we focused on analyzing the functions and biological pathways associated with targets of senescence-associated TFs at different developmental stages. As a result, 797, 756, and 4552 targets of early, middle, and late senescence-associated TFs were assigned to 38, 41, and 51 GO terms, respectively (Table S8). A high percentage of targets of senescence-associated TFs from GO terms of cell part (GO:004446), binding (GO:0005488), and response to stimulus (GO:0050896) (Fig. 7A). GO terms of response to abiotic stimulus (GO:0009628), response to stress (GO:0006950), and response to stimulus (GO:0050896) were significantly enriched (corrected p-value < 0.05, see Fig. 7B, Table S6) by the targets of early and late senescence-associated TFs, suggestive of protective mechanisms against potential stresses during tobacco leaf senescence. No significant enrichment of GO term was observed among the targets of middle senescence-associated TFs.Fig. 7GO analysis of senescence-associated TFs. (A) GO classifications of targets of early, middle, and late senescence-associated TFs. (B) The enriched GO terms (corrected P-value ≤ 0.05) annotated by targets of senescence-associated TF genes. Each dot indicates a GO term. The dot sizes indicate the number of genes. Colors indicate the significance of the GO term.Previous studies have revealed many biological pathways associated with leaf senescence, including stress resistance, plant hormone signal transduction, and protein metabolism pathways. In the present study, the targets of early, middle, and late senescence-associated TFs were assigned to 94, 85, and 130 metabolic pathways according to KEGG annotation (Table S7). Based on the KEGG enrichment analysis results, targets of early senescence-associated TFs were enriched in protein processing in the endoplasmic reticulum, zeatin biosynthesis, and spliceosome pathways. Middle senescence-associated TFs enriched pathways of flavone and flavonol biosynthesis, zeatin biosynthesis, and plant pathogen interaction. Late senescence-associated TFs were enriched in three pathways, including protein processing in the endoplasmic reticulum, ubiquitin-mediated proteolysis, and flavonoid biosynthesis (Fig. 8, Table S7). These results revealed the targets of senescence-associated TFs responsible for essential biological functions during leaf senescence. These senescence-associated TFs may control tobacco leaf senescence by regulating the expression of genes involved in these metabolic pathways.Fig. 8Heatmap of KEGG pathway which enriched by the targets of early, middle, and late senescence-associated TF genes. The color from red to blue represent the − Log10 corrected p-values.qRT-PCR analysis of hub TF genes selected from networkIn the current study, hub genes selected from senescence response TF-TF-target networks were subjected to qRT-PCR analysis using samples collected from mature leaves (ML) at different time points (15, 25, 35, 45, 55, 65, 75, and 85 (DAT)). As a result, the ex-pression patterns of these hub TF genes determined by qRT-PCR were basically consistent with the transcriptome sequencing data, confirming the reliability of the transcriptome data used in this study. The hub TF genes detected from different TF-TF-target networks exhibited different time-course expression characteristics during leaf senescence. The results will help infer their potential roles in regulating the senescence of tobacco leaves. The qRT-PCR results are shown in Fig. 9, and the primer sequences are available in Table S1.Fig. 9Expression analysis of hub genes associated with senescence response TF-TF-targets networks at indicated time points. The error bars above line present ± SD with three biological replicates and asterisks denote a statistically significant difference compared with expression level at 15 DAT according to one-way ANOVA analysis (***P < 0.01).
