A Transcriptomic Evaluation of the Mechanism of Programmed Cell Death of the Replaceable Bud in Chinese Chestnut

1. Introduction

Chinese chestnut (Castanea mollissima BL.) is an important nut-producing tree grown in temperate regions worldwide [1]. However, chestnut trees suffer reduced productivity over time due to branch growth characteristics. In chestnut, the apical bud of each existing fruiting branch will develop into a new fruiting branch in the following year, resulting in the majority of fruit development taking place along the peripheral crown. A notable exception, the spontaneous mutant chestnut cultivar (cv.) ‘Tima Zhenzhu’, was first described in 1979. In ‘Tima Zhenzhu’ trees, the apical bud (hereafter “replaceable bud”) senesces and dies, resulting in a compact crown [2]. Research indicates that the death of the replaceable bud involves programmed cell death (PCD), as evidenced by the presence of certain markers of PCD, including chromatin condensation, nuclear degradation, DNA laddering, tonoplast invagination, vacuolar rupture, and autophagy, among others [1].

PCD, an orderly process of cellular suicide mediated by intracellular death programs [3], is crucial for normal tissue development and stress response [4]. PCD occurs during both reproductive and vegetative tissue development, and is essential for trichome formation, sex determination, xylogenesis, leaf senescence, and aerenchyma formation, among other processes [5,6,7]. Each of these processes requires the destruction of targeted cells, for example in the replaceable bud, root cap, nucellar tissue, and tapetum. Due to a lack of apoptotic regulator and executor homologs of animal PCD, the process of PCD is distinctive in plants [8]. PCD regulation involves a network of molecules with prepare the cell for PCD, as well as molecules which act as initiators, effectors, and degraders [7]

Cellular preparation for PCD is coordinated primarily by transcriptional regulation of hormone signaling, most commonly by the ethylene signaling pathway [9,10]. Ethylene promotes several PCD processes, including leaf senescence [11], flower deterioration following pollination [12], and root decline under hypoxia [13]. Other phytohormones also play key roles in PCD-related processes: gibberellic acid (GA) promotes aleuronal and tapetal PCD; GA is antagonized by abscisic acid (ABA) [14,15]; auxin is important for xylogenesis, root cap morphogenesis, and leaf senescence [16]; brassinosteroid (BR), auxin, and cytokinin promote tracheary element (TE) PCD [17]; and salicylic acid (SA) and jasmonic acid (JA) act to regulate both aerenchyma formation and leaf senescence [18,19].

Transcription factors (TFs) link phytohormone signaling to PCD regulation. For example, several TFs appear to be central to both nucellar and tapetal PCD, including MYB, MADS-box, bHLH, and NAC [16,17,20]. Downstream of hormone signaling, PCD can be triggered in specific cell types by diverse cellular events such as changes in intracellular Ca2+ concentration, buildup of reactive nitrogen species (RNS) and reactive oxygen species (ROS), activation of protein kinases, acidification of the cytoplasm, and modification of the cytoskeleton [16,21,22]. In addition, a plethora of hydrolytic enzymes, including nucleases and proteases, have been suggested as putative PCD executors. For example, in zinnia, the nuclease ZEN1 is involved in DNA fragmentation in TEs [23]. In Arabidopsis, xylem cysteine proteinases (AtXCP1 and AtXCP2) and metacaspase 9 (AtMC9) function as autolysis executors in TEs [16]. In rice, the cysteine protease (CEP) gene promotes nucellar degeneration [24], and two aspartic proteases (OsAP37 and OsAP25) stimulate tapetal PCD [25]. In barley endosperm, a vacuolar protease (HvVPE1) is responsible for degradation of nuclear DNA [26].

In this work, we performed transcriptomic profiling of the replaceable bud of Chinese chestnut cv. ‘Tima Zhenzhu’ before (S20), during (S25), and after PCD (S30). We sought to (i) investigate the expression patterns, functions, and associated metabolic pathways of any differentially expressed genes (DEGs); and (ii) elucidate the regulatory and signaling pathways related to replaceable bud PCD in chestnut cv. ‘Tima Zhenzhu’. The results of this work will not only serve to bolster our understanding of the molecular mechanism of replaceable bud PCD, but may also provide a relevant resource for the future genetic improvement of chestnut.

2. Materials and Methods

2.1. Plant materials

Chestnut cv. ‘Tima Zhenzhu’ trees were kept at the Changli Institute of Pomology, Academy of Agricultural and Forestry Sciences, Hebei Province, China (118°51′ E, 39°53′ N). Replaceable buds were collected at 20 (S20), 25 (S25), and 30 (S30) days after flowers, corresponding to the developmental stage before, during, and after PCD (Wang et al. 2012). At S20, the buds are green and display a normal phenotype; at S25, the buds are yellow-green; at S30, the buds are brown-yellow and display a partial abscission layer (Figure 1 and Figure 2). Three biological replicates were collected at each phase, and all analyses were conducted using three biological replicates.

2.2. Methods

2.2.1.Transcriptomic Sequencing and Library Construction

RNAseq was carried out with 1 μg of RNA per sample. A NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, MA, USA) was utilized to ascertain both the purity and concentration of the RNA. An RNA Nano 6000 Assay Kit and a 2100 Bioanalyzer (Agilent Technologies, CA, USA) were utilized to ascertain RNA integrity. Transcriptomic analysis was conducted by the Biomarker Technologies Company (Beijing, China). A NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, MA, USA) was utilized to generate sequencing libraries. The AMPure XP system (Beckman Coulter, CA, USA) was utilized to perform purification of library fragments. The Illumina NovaSeq 6000 platform (Illumina, CA, USA) was utilized to perform library sequencing.

2.2.2. Identification and Functional Annotation of DEGs

Raw sequence data was cleaned up by filtering out adapter sequences and reads with low quality. HISAT2 software [27] was utilized to map clean reads to the C. mollissima reference genome [28]. StringTie software [29] was utilized to quantify gene expression as fragments per kilo bases of exons for per million mapped reads (FPKM) values. Genes were considered as “expressed” if their FPKM was > 0. DESeq2 software [30] was utilized to discover DEGs between the S20, S25, and S30 timepoints, with a threshold false discovery rate (FDR) set at < 0.01 and fold-change ≥ 2. Finally, the Database for Annotation, Visualization, and Integrated Discovery (DAVID) was utilized for both analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and annotation of gene ontology (GO), as implemented in BMK Cloud provided by Biomarker Technologies Company (Beijing, China)

2.2.3. Protein-Protein Interaction Network Construction

Protein-protein interactions were queried using the STRING database (https://string-db.org/cgi). The proteins encoded by DEGs putatively involved in replaceable bud PCD were networked using their respective tobacco homologs. The Cytoscape software [31] was utilized to visualize the resulting protein-protein interaction network.

2.2.4. Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) Validation

An RNAqueous Total RNA Isolation Kit (Solebo Technology, China) was utilized to extract total RNA from the replaceable bud before (S20), during (S25), and after (S30) PCD. The HiScript II Q RT SuperMix for qPCR (Yisheng Biotechnology, China) was utilized to synthesize cDNA from the extracted RNA, with Actin serving as the internal reference gene. All qRT-PCR primers can be found in Table S7. A LightCycler 480II Real-Time PCR Detection System (Roche, Switzerland) was utilized to perform qRT-PCR. The reaction mixture contained 2 µL of cDNA template, 10 μL of ChamQ SYBR qPCR Master Mix, 6.8 μL of deionized water, 0.4 μL of ROX Reference Dye2, 0.4 µL of the forward primers, and 0.4 µL of the reverse primers. The amplification program was conducted as follows: 95 °C for 5 min; 40 cycles of 95 °C for 10 s and 60 °C for 30 s. To ensure data reproducibility, three independent experiments were conducted. The 2^-△△Ct method [17] was utilized to quantify the levels of relative gene expression.

3. Results

3.1. Transcriptomic Analysis

Samples of chestnut cv. ‘Tima Zhenzhu’ replaceable bud tissue were used to prepare RNA libraries before (S20), during (S25), and after (S30) PCD. After quality control, clean reads were mapped to the C. mollissima reference genome, producing 86.94% (37,370,780), 90.25% (40,861,735), and 89.77% (40,896,455) mapped reads in the S20, S25, and S30 libraries, respectively (Table 1). Of these, uniquely mapped reads accounted for 82.87% (35,620,696 reads) of the S20 library, 85.43% (38,678,320 reads) of the S25 library, and 83.21% (37,872,455) of the S30 library. The rate of multiple mapped reads across the three libraries was less than 6.56%, and the percentage of Q30 bases was 93.24%, 94.29%, and 93.31% in the S20, S25, and S30 libraries, respectively (Table 1), indicating that our data was of sufficiently high quality. To assess the reliability of the constructed libraries, heat mapping and principal component analysis (PCA) were used to examine relationships among different samples and replicates (Figure 3). Overall, the read counts between both samples and replicates showed high correlations, with correlation coefficients (r) ranging between 0.98-0.99, 0.775-0.945, and 0.703-0.872 in S20, S25, and S30, respectively (Figure 3a). Significant differences were observed between samples (S20, S25, and S30), with a variance of 66.2% between samples and a variance of 19.6% within the same sample (Figure 3b).

3.2. Global Transcriptional Changes before, during, and after PCD in Replaceable Buds of Chestnut cv. ‘Tima Zhenzhu’

In order to identify DEGs, a stringent threshold of FDR < 0.01 and |log2FC| > 1 was used. In total, 5779 DEGs were identified in S20 vs S25, 9,867 in S20 vs S30, and 2,674 in S25 vs S30 (Figure. 4). The large number of DEGs (10,870 in total; Table S1) suggested that replaceable bud PCD may be a complex process regulated by an extensive population of genes. Moreover, we found a significantly greater number of DEGs in S20 vs S30, compared to the other two timepoint comparisons, indicating that gene expression in replaceable bud PCD is highly variable through time. Clear differences were found in global gene expression patterns between the three timepoints (Figure 4). In general, the amount of both down- and up-regulated genes was large in all comparisons, although the amount of down-regulated genes tended to be slightly higher (Figure 4).

3.3. Functional Enrichment Analysis of the Common-DEGs

The DEGs shared by at least two comparisons (“common-DEGs”; 6,137 in total) were further evaluated to assess their core biological functions related to replaceable bud PCD (Figure 4d; Table S2). The common-DEGs were annotated with 35 GO terms (P ≤ 0.001) belonging to 14 molecular functions (MFs), 15 cellular components (CCs), and 19 biological processes (BPs) (Figure 5; Table S3). Annotated BP GO terms included “microtubule cytoskeleton organization” (GO:0000226),“oxidation-reduction process” (GO:0055114), “lignin catabolic process” (GO:0046274), “protein phosphorylation” (GO:0006468), “cell wall macromolecule catabolic process” (GO:0016998 “metabolic process” (GO:0008152), “microtubule-based movement” (GO:0007018), “carbohydrate metabolic process” (GO:0005975), “DNA replication initiation” (GO:0006270), ), “secondary metabolic process” (GO:0019748), and “regulation of transcription, DNA-templated” (GO:0006355). Annotated CC GO terms included “integral component of membrane” (GO:0016021), “anchored component of plasma membrane” (GO:0046658), “microtubule” (GO:0005874), “intracellular membrane-bounded organelle” (GO:0043231), “kinesin complex” (GO:0005871), “apoplast” (GO:0048046), and “plasma membrane” (GO:0005886). Annotated MF GO terms included “transferase activity, transferring acyl groups other than amino-acyl groups” (GO:0016747), “microtubule motor activity” (GO:0003777), “polysaccharide binding” (GO:0030247), “protein serine/threonine kinase activity” (GO:0004674), “hydrolase activity, hydrolyzing O-glycosyl compounds” (GO:0004553), “hydroquinone: oxygen oxidoreductase activity” (GO:0052716), “2-alkenal reductase [NAD(P)] activity” (GO:0032440), “microtubule binding” (GO:0008017), “serine-type carboxypeptidase activity” (GO:0004185), and “DNA binding” (GO:0003677).

Several GO terms associated with either autophagy or PCD were identified (Additional File 3: Table S3), including “negative regulation of programmed cell death” (GO:0043069), “singlet oxygen-mediated programmed cell death” (GO:0010343), “positive regulation of programmed cell death” (GO:0043068), “poly(A)-specific ribonuclease activity” (GO:0004535), “positive regulation of autophagy” (GO:0010508), “regulation of programmed cell death” (GO:0043067), “autophagy” (GO:0006914), “regulation of autophagy” (GO:0010506), “single-stranded DNA endodeoxyribonuclease activity” (GO:0000014), “autophagy of mitochondrion” (GO:0000422), “autophagy of nucleus” (GO:0044804), “single-stranded DNA 5’-3’ exodeoxyribonuclease activity” (GO:0045145), “endoribonuclease activity, producing 5’-phosphomonoesters” (GO:0016891), “5’-3’ exodeoxyribonuclease activity” (GO:0035312), “ribonuclease activity” (GO:0004540), “RNA-DNA hybrid ribonuclease activity” (GO:0004523), and “3’-5’-exoribonuclease activity” (GO:0000175).

KEGG enrichment analysis was also conducted on the common-DEGs in order to ascertain which biological pathways they may be involved in (Table S4). Among the 125 enriched KEGG pathways found, the most significant pathways (P ≤ 0.001) were “plant hormone signal transduction” (ko04075), “peroxisome” (ko04146), “linoleic acid metabolism” (ko00591), “alpha-linolenic acid metabolism” (ko00592), “nicotinate and nicotinamide metabolism” (ko00760), “phenylpropanoid biosynthesis” (ko00940), “valine, leucine, and isoleucine degradation” (ko00280), “glutathione metabolism” (ko00480 “arginine and proline metabolism” (ko00330), “ABC transporters” (ko02010), “cutin, suberine, and wax biosynthesis” (ko00073), “fatty acid degradation” (ko00071), “glycerophospholipid metabolism” (ko00564), and “ether lipid metabolism” (ko00565) (Table 2 and Figure 6).

3.4. DEGs Related to Phytohormone Signal Transduction Pathways

In total, the “plant hormone signal transduction” (ko04075) KEGG pathway was found to contain 93 DEGs (Table 3 and Table S5), the bulk of which were related to ethylene, BR, SA, GA, ABA, cytokinin, or auxin signaling pathways (Figure 6). Within the ethylene pathway, the following DEGs were up-regulated at S25 and S30, compared to S20: EVM0010215 (encoding mitogen-activated protein kinase like d5), EVM0016645 (encoding ethylene response sensor 1), EVM0019242 (encoding ethylene receptor 2-like), EVM0007260 (encoding EIN3-binding F-box protein), three DEGs coding for ethylene-insensitive protein 3 (EVM0013705, EVM0023854, and EVM0025637), and three DEGs for ethylene-responsive transcription factors (EVM0030131, EVM0003409, and EVM0022486). In the auxin pathway, the following DEGs were down-regulated at S25 and S30, compares to S20: 15 DEGs encoding auxin-responsive proteins (IAA) (e.g., EVM0003597, EVM0004089, and EVM0005109), four DEGs for auxin response factors (EVM0015180, EVM0028394, EVM0032825, and EVM0002460), two DEGs coding for auxin-induced proteins (EVM0027083 and EVM0014325), and two DEGs coding for auxin transporter-like proteins (EVM0011123 and EVM0030472). Additionally in the auxin pathway, EVM0019998, coding for auxin influx transport protein, was up-regulated at S25 and S30, compared to S20. Moreover, five up- and four down-regulated DEGs belonged to the auxin-responsive protein SAUR genes (e.g., EVM0003671, EVM0004337, and EVM0007905). In the cytokinin pathway, one DEG (EVM0008286) coding for histidine kinase and one DEG (EVM0003407) coding for two-component response regulators were up-regulated, while three DEGs (EVM0013339, EVM0019639, and EVM0020813) coding for histidine kinases, eight DEGs (e.g., EVM0000777, EVM0008810, EVM0016676, and EVM0017362) coding for two-component response regulators, and two DEGs (EVM0015700 and EVM0018230) coding for histidine-containing phosphotransfer proteins were down-regulated at S25 and S30, compared to S20. In the ABA pathway, 12 DEGs were up-regulated and seven DEGs were down-regulated at S25 and S30, compared to S20, among which four DEGs code for ABA-insensitive 5-like proteins (two up- and two down-regulated; EVM0006332, EVM0010332, EVM0019421, and EVM0031140), six DEGs for protein phosphatase 2C (one down- and five up-regulated; e.g., EVM0003340, EVM0007208, and EVM0009898), and nine DEGs for serine/threonine-protein kinases (four down- and five up-regulated; e.g., EVM0002534, EVM0003135, and EVM0003803). In the GA pathway, two DEGs encoding gibberellin receptors (EVM0018995 and EVM0023348) and two DEGs encoding DELLA proteins (EVM0002951 and EVM0009195) were down-regulated, whereas two DEGs encoding transcription factor PIF3 (EVM0021995 and EVM0023743) were up-regulated at S25 and S30, compared to S20. In the SA pathway, three DEGs encoding pathogenesis-related protein 1-like (EVM0020541, EVM0033371, and EVM0001078) and EVM0018658, encoding regulatory protein, were consistently up-regulated. In the BR pathway, several DEGs coding for brassinosteroid insensitive 1-associated receptor kinase 1 were up- (e.g., EVM0005680) or down-regulated (e.g., EVM0008845, EVM0022513, and EVM0027385 for cyclin-D3).

3.5. DEGs Related to Replaceable Bud PCD Regulation in Chestnut cv. ‘Tima Zhenzhu’

In order to identify DEGs central to the regulation of replaceable bud PCD, the functional annotations of all DEGs were compared to the NCBI non-redundant protein sequences(NR) database. A total of 441 DEGs were found to be related to PCD processes, including both the initiation and execution of PCD as well as phytohormone signaling (Table S6). In particular, the ethylene signaling pathway was found to be highly enriched, involving DEGs such as EVM0010215 (coding for mitogen-activated protein kinase like d5), EVM0016645 (coding for ethylene response sensor 1), EVM0019242 (coding for ethylene receptor 2-like), EVM0007260 (coding for EIN3-binding F-Box protein), three DEGs coding for ethylene-insensitive protein 3 (EVM0013705, EVM0023854, and EVM0025637), and three DEGs coding for ethylene-responsive transcription factors (EVM0030131, EVM0003409, and EVM0022486). Moreover, a total of 16 DEGs related to ethylene biosynthesis exhibited significantly higher expression levels at S25 and S30 than at S20, including three 1-aminocyc lopropane-1-carboxylate synthases (EVM0001055, EVM0011816, and EVM0029965), one s-adenosylmethionine synthase (EVM0032964), and 12 1-aminocyclopropane-1-carboxylate oxidases (e.g., EVM0000849, EVM0006088, EVM0006661, and EVM0010524).

Both PCD initiation and execution require the temporal and spatial coordination of many precisely-regulated processes, which work downstream of ethylene signaling (Li et al. 2019). We found that DEGs related to PCD initiation were primarily TFs, including 43 NAC TFs (e.g., EVM0010588, EVM0016926, and EVM0018368), 7 MADS-box TFs (e.g., EVM0019787, EVM0021142, and EVM0023741), 31 bHLH TFs (e.g., EVM0012487, EVM0013783, EVM0016064, and EVM0019464), and 49 MYB TFs (e.g., EVM0021418, EVM0024602, EVM0033473, and EVM0006679) (Table S6). Other functional components were also identified, including cytochrome c (EVM0024076), MAPK (EVM0012332, EVM0022934, EVM0010215, and EVM0030224), MAPKK (EVM0004329, EVM0000127, and EVM0025130), MAPKKK (e.g., EVM0015125, EVM0014496, and EVM0033128), and Ca2+ uniporters (EVM0018926 and EVM0024216). A total of 36 DEGs were associated with Ca2+ signaling (Table S6), including four Ca2+ channels (EVM0026184, EVM0007035, EVM0016547, and EVM0023315), nine Ca2+-dependent protein kinases (e.g., EVM0022389, EVM0025960, EVM0030576, and EVM0031549), nine Ca2+ B-like proteins (e.g., EVM0000079, EVM0000956, EVM0004877, and EVM0010534), and 14 calmodulins (e.g., EVM0030734, EVM0033478, EVM0002917, and EVM0004863). DEGs related to PCD execution were primarily enriched in autophagy and proteases. Protease genes included one metacaspase (EVM0011448), eight cysteine proteinases (e.g., EVM0007262, EVM0007929, EVM0020838, and EVM0025298), two vacuolar-processing enzymes (EVM0002134 and EVM0012232), 20 vacuolar protein sorting-associated proteins (e.g., EVM0012501, EVM0012829, EVM0013147, and EVM0015898), nine endonucleases (EVM0006680, EVM0015983, EVM0021509, and EVM0022275), 19 aspartic proteases (e.g., EVM0005045, EVM0005145, EVM0009402, and EVM0012500), five endoglucanases (EVM0016472, EVM0017671, EVM0022249, EVM0032182, and EVM0003400), two senescence-associated proteins (EVM0009758 and EVM0031003), 17 pectinesterases (e.g., EVM0005897, EVM0021830, and EVM0032186), 13 xyloglucan endotransglucosylase/hydrolases (e.g., EVM0023103, EVM0028198, and EVM0031923), three exosome complex components (EVM0015970, EVM0019411, and EVM0029579), one Werner syndrome ATP-dependent helicase (EVM0029274), 13 lipoxygenases (e.g., EVM0032274, EVM0000380, EVM0000461, and EVM0002948), and 16 acid phosphatases (e.g., EVM0006765, EVM0014183, and EVM0022260). Autophagy-related DEGs included 18 autophagy-related proteins (e.g., EVM0032306, EVM0002553, EVM0003991, and EVM0004113), seven protein-transport proteins (e.g., EVM0001633, EVM0007575, and EVM0011846), one syntaxin-related protein KNOLLE (EVM0031095), two vesicle transport proteins (EVM0014337 and EVM0016934), and one calreticulin (EVM0014858). Additionally, some DEGs involved in longevity regulation and stress response were identified, including seven chaperonins (e.g., EVM0008263, EVM0008832, and EVM0009026), six BAG family molecular chaperone regulators (e.g., EVM0001410, EVM0009410, and EVM0024418), two catalases (EVM0024729 and EVM0009488), three Cu-Zn superoxide dismutases (EVM0015967, EVM0017860, and EVM0032156), one Mn superoxide dismutase (EVM0018305), and 37 heat shock proteins (e.g., EVM0015018, EVM0015246, EVM0016538, and EVM0017835).

To further investigate the regulatory pathways related to replaceable bud PCD of chestnut cv. ‘Tima Zhenzhu’, a protein-protein interaction network was completed for the DEGs associated with plant PCD, according to their tobacco homologs (Figure 7). The resultant network included 114 DEGs, grouped into three distinct modules. The first network module was clustered into two groups: one consisting of heat shock proteins (e.g., Nthsp70, A0A1S4B2J0, Nthsp26a, A0A1S4CVK6, A0A1S3XCI4, A0A1S4CCK2, and A0A1S4CB32), and the other consisting of DEGs involved in either the ethylene signaling pathway (e.g., EIL4, EIL1, A0A1S4DNE7, acs2, A0A1S3YAY2, A0A1S3ZD07, and A0A1S3ZDB8) or the MAPK signaling cascade (e.g., A0A1S4AN76, NtMEK2, A0A1S4ADQ5, NPK2, MPK11, A0A1S4B6U4). Notably, between the two groups, A0A1S4B6Y7 (metacaspase-9-like) showed strong interactions with A0A1S3Z1L4 (Ca²⁺-transporting ATPase), A0A1S4DPX4 (Ca²⁺-transporting ATPase), and A0A1S4BNP6 (chaperonin 60 subunit), suggesting that these genes may interact at the nexus of stress response, and ethylene and MAPK signaling. The second network module consisted of autophagy-related proteins (A0A1S3Y2M7, ATG8e, A0A1S4A7C8, ATG4, A0A1S4D5L9, ATG5, A0A1S4C6J8, ATG9, ATG2, A0A1S4DLP2, ATG13a, and A0A1S4BPT6) and vacuolar protein sorting-associated proteins (A0A1S4DPC1, A0A1S4DKE1, A0A1S4CYV6, A0A1S4BA99, A0A1S3WXN7, and A0A1S4DLQ1). The third network module consisted of DNA repair proteins (e.g., A0A1S3ZD09, A0A1S3ZXG5, A0A1S4BBM1, and A0A1S3YFL3) and exonucleases (A0A1S4AJ03, A0A1S3ZDV9, and A0A1S4AUU7). Of these, A0A1S4D430, a polygalacturonase-like protein, showed strong interactions with eight pectinesterase proteins (PPME1, A0A1S4DRD3, A0A1S3XCR1, A0A1S4A2I1, A0A1S4C5C8, A0A1S3ZXA8, A0A1S4DB14, and A0A1S3YLM9); and A0A1S3Y4S2, a respiratory burst oxidase homolog protein, showed strong interactions with several Ca²⁺-dependent protein kinases (e.g., A0A1S4C4C3, CDPK14, A0A1S4DBM4, cdpk3, and A0A1S4ASR5).

3.6. qRT-PCR Validation of the Transcriptomic Data

Nine of the PCD-related DEGs were chosen for further qRT-PCR analysis to confirm their expression levels during replaceable bud PCD of chestnut cv. ‘Tima Zhenzhu’. Overall, the RNA-seq results were largely in agreement with the qRT-PCR results (Figure 8), suggesting that our transcriptomic analysis was reasonable and accurate.

4. Discussion

4.1. The role of the Ethylene Signaling in Replaceable Bud PCD of Chestnut cv. ‘Tima Zhenzhu’

Phytohormone signal transduction is likely the prevailing mechanism responsible for upstream regulation of PCD processes [7]. For example, ethylene is responsible for leaf senescence and root deterioration in maize [11,12]; auxin is responsible for developmentally-regulated PCD in lace plant [32]; GA is related to both aleuronal and tapetal PCD [14,15]; the combination of auxin, cytokinin, and BR are related to PCD of TEs [33]; and JA and SA are related to both aerenchyma formation and leaf senescence [18,19].

Here too, KEGG analysis indicated that the “plant hormone signal transduction” (ko04075) pathway was particularly enriched. Several key components of the ethylene signaling pathway were found to be up-regulated at S25, including mitogen-activated protein kinase like d5 (EVM0010215), EIN3 (EVM0013705, EVM0023854, and EVM0025637), ERS1 (EVM0016645), and EBF1 (EVM0007260) (Figure 6). In particular, EIN3 exerts control over cellular senescence and death through its involvement in a trifurcate feed-forward pathway [34]. Within the EIN2–EIN3–NAC regulatory cascade controlling leaf senescence-associated PCD, EIN3 acts to transcriptionally activate NAC TFs (ORE1 and AtNAP), which themselves act to positively regulate leaf senescence [35]. Additionally, EIN3 functions as a direct repressor of miR164 [34], which acts to post-transcriptionally negatively regulate ORE1. In order to efficiently regulate leaf senescence, EIN3 appears to simultaneously regulate the expression of both ORE1 and its particular negative regulator (miR164). By analogy, the up-regulation of EIN3 at S25 likely promoted replaceable bud PCD. Taken together, these findings suggest that ethylene signaling is activated during PCD preparation in replaceable buds of chestnut cv. ‘Tima Zhenzhu’.

In our previous study, we reported that the ethylene content was higher in replaceable buds at S25 than at S20 [36]. Our present study bolsters this observation, as we found that genes related to ethylene biosynthesis were consistently up-regulated at S25 (Table S6), including 1-aminocyclopropane-1-carboxylate oxidases (EVM0000849, EVM0006088, EVM0006661, EVM0010524, EVM0016161, EVM0017576, EVM0017856, EVM0019065, EVM0019151, EVM0019624, EVM0023867, and EVM0031148), s-adenosylmethionine synthase 2 (EVM0032964), and 1-aminocyclopropane-1-carboxylate synthases (EVM0001055, EVM0011816, and EVM0029965). We additionally identified DEGs involved in other hormone signaling pathways (Table 3), suggesting that several phytohormones may interact during replaceable bud PCD. For example, both GA and ABA antagonistically regulate aleuronal PCD [5,19]. Additionally, leaf senescence is promoted by a combination of ethylene, SA, and ABA, while this process is delayed by a combination of auxin, cytokinin, and GA [37,38]. As shown in Table 3 and Figure 6, DEGs involved in ethylene signaling, such as MAPK (EVM0010215), ERS1 (EVM0016645), EIN3 (EVM0013705 and EVM0025637), and ethylene-responsive TF 1B (EVM0030131 and EVM0003409), were up-regulated at S25 vs. S20. Contrarily, DEGs involved in the auxin pathway, such as AUX28 (EVM0014325), SAUR-like auxin-responsive proteins (EVM0004337, EVM0007905, EVM0019733, and EVM0021087), and auxin-responsive proteins (e.g., EVM0003597, EVM0004089, EVM0005109, EVM0007851, EVM0017536, and EVM0017744), or in the cytokinin pathway, such as AHP5 (EVM0015700) and two-component response regulator ARR family members (EVM0008810, EVM0017362, EVM0020552, and EVM0022760), were down-regulated at this same time point. In our previous study, we reported that the content of both auxin and cytokinin were significantly lower in replaceable buds at S25 than at S20 [36]. Taken together, it appears that replaceable bud PCD of chestnut cv. ‘Tima Zhenzhu’ is primarily controlled by high ethylene and low auxin and cytokinin signaling, suggesting that ethylene may antagonistically interact with auxin and cytokinin in this process.

The MAPK signaling cascade acts to transduce extracellular signals into biochemical and physiological responses [39]. Each MAPK signaling cascade contains at least three essential members: MAPK, MAPKK, and MAPKKK [40]. There is evidence that the MAPK signaling cascade may be essential for PCD in plants. For example, in mutant rice, several genes encoding MAPKs (e.g. MPK3, MPK5, and MPK13) are up-regulated during leaf PCD [41]. MAPK appears to work to initiate self-incompatibility (SI)- induced PCD when poppy flowers are exposed to incompatible pollen [42]. In Arabidopsis, MAPK cascades acts to regulate leaf senescence [43], mediate cell death, and induce ethylene production through the regulation of ACS2, ACS6, and ACS8 gene expression [40]. Here, we identified 16 DEGs (e.g., EVM0012332, EVM0022934, EVM0004329, EVM0000127, and EVM0025130) associated with MAPK signaling, with most of these exhibiting up-regulating during replaceable bud PCD (Table S6). Furthermore, a MAPK cascade involving A0A1S4ADQ5 (MAPKKK)-NtMEK2 (MAPKK)-A0A1S4AN76 (MAPK) was found in the protein-protein interaction network (Figure 7), and A0A1S4AN76 showed strong interactions with ethylene biosynthesis-related genes A0A1S4DNE7 (1-aminocyclopropane-1-carboxylate synthase), A0A1S3ZDB8 (S-adenosylmethionine synthase), A0A1S3YAY2 (1-aminocyclopropane-1-carboxylate oxidase-like), A0A1S3ZD07 (1-aminocyclopropane-1-carboxylate oxidase), and acs2 (1-aminocyclopropane-1-carboxylate synthase-like). Taken together, we suggest that the MAPK signaling cascade controls ethylene biosynthesis by transcriptional regulation of related genes during replaceable bud PCD of chestnut cv. ‘Tima Zhenzhu’.

4.2. PCD Initiation in Replaceable Buds of Chestnut cv. ‘Tima Zhenzhu’

Several TFs may act as bridges linking phytohormone signaling with PCD regulation [44,45]. As an example, ORESARA1 (ANAC092, a NAC TF) regulated leaf senescence both downstream of ethylene signaling and upstream of senescence-inducing genes, including NAC TFs such as BIFUNCTIONAL NUCLEASE 1 (BFN1) [46]. Homologs of MYB [47] and bHLH [48] TFs act to promote tapetal PCD. In Arabidopsis, the endospermous MADS-box TF AGAMOUS-LIKE 62 promotes PCD in nucellar tissue after fertilization through activation of the PCD-promoting MADS-box TF TRANSPARENT TESTA 16 [49]. In thus study, we discovered a large number of DEGs encoding TF homologs, including 48 MYB TFs (e.g., EVM0013996, EVM0000707, EVM0009774, EVM0011416, and EVM0021418), seven MADS-box TFs (e.g., EVM0019787, EVM0021142, and EVM0023741), 31 bHLH TFs (e.g., EVM0012487, EVM0013783, EVM0016064, and EVM0019464), and 43NAC TFs (e.g., EVM0010588, EVM0016926, EVM0018368, and EVM0018380) (Table S6). These results suggest that several TFs may act to ensure initiation of bud PCD of chestnut cv. ‘Tima Zhenzhu’.

Cytochrome c plays a part of the PCD signaling network, with the up-regulation of both cytochrome c oxidase and cytochrome c being an early event during the process of PCD [50,51]. The process of PCD in several plant systems is accompanied by the dispensation of mitochondrial cytochrome c, including tapetal PCD in sunflower and terminally differentiated suspensor PCD in runner bean cotyledons [52,53]. Here, we found that the expression levels of genes encoding cytochrome c (EVM0024076) and cytochrome c oxidases (EVM0029309, EVM0032493, EVM0012209, EVM0004365, EVM0022679, and EVM0012699) were notably up-regulated at S25 compared to S20 (Table S6). Therefore, we suggest that cytochrome c and cytochrome c oxidases may act to modulate replaceable bud PCD initiation in chestnut cv. ‘Tima Zhenzhu’.

Substantial evidence suggests that Ca²⁺ is crucial for PCD regulation in plants [54]. The loss of intracellular Ca²⁺ stores and/or alteration of intracellular Ca²⁺ levels is associated with several PCD processes, including the regulation of Ca²⁺-dependent hydrolytic enzymes and endonucleases [5]. In plants, the cytosolic increase of Ca²⁺ is an early event during PCD. Specifically, Ca²⁺-permeable channels, Ca²⁺ sensor calcineurin B-like proteins (CBLs), calmodulins (CAMs), and Ca2+-dependent protein kinases (CDPKs) are necessary for Ca²⁺ signal transduction and PCD [55], and stimulus-induced Ca²⁺ elevation is interpreted by downstream Ca²⁺ sensors [56]. In tomato, both CIPK6 and CBL10 are necessary for Pto kinase-triggered PCD upon perception of Pseudomonas syringae effectors [57]. In both yeast and transgenic plants, AtBAG6, a CaM-binding protein, is up-regulated by stress and promotes PCD [56]. CDPK3 positively regulates PCD in Arabidopsis [58]. Additionally, CDPK proteins are up-regulated during the process of artificial ageing in maize seeds [59]. In our previous study, we reported that the content of Ca²⁺ in the cytoplasm, nucleus, and vacuole increased significantly during replaceable bud PCD of chestnut cv. ‘Tima Zhenzhu’ [60]. Here, we identified 36 DEGs associated with Ca²⁺ signaling (Table S6), including four Ca²⁺ channels (EVM0026184, EVM0007035, EVM0016547, and EVM0023315), nine Ca²⁺-dependent protein kinases (e.g., EVM0022389, EVM0025960, EVM0030576, and EVM0031549), nine calcineurin B-like proteins (e.g., EVM0000079, EVM0000956, EVM0004877, and EVM0010534), and 14 calmodulins (e.g., EVM0030734, EVM0033478, EVM0002917, and EVM0004863). The bulk of these DEGs were significantly up-regulated at S25 and/or S30 compared to S20 (Table S6). Moreover, in the constructed protein-protein interaction network (Figure 7), A0A1S3Z7E2 (calcineurin B-like protein 9) showed strong interaction with A0A1S3YJV4 (calmodulin-binding receptor-like cytoplasmic kinase 3), A0A1S3YYX5 (calcineurin B-like protein 3 isoform X2), and MPK11 (mitogen-activated protein kinase 9). Taken together, these results suggest that the Ca²⁺-dependent signaling cascade may be related to the modulation of replaceable bud PCD in chestnut cv. ‘Tima Zhenzhu’.

4.3. PCD Execution in Replaceable Buds of Chestnut cv. ‘Tima Zhenzhu’

The reception of specific signal triggers initiates PCD execution and cellular corpse clearance, resulting in the activation and release of myriad hydrolytic enzymes, including various nucleases and proteases [44]. Our previous study found no evidence of DNA degradation at S20, with some DNA degradation occurring at S25, and heavy DNA degradation occurring at S30 [1]. Here, we found several DEGs encoding nucleases, including endonuclease V isoform X1 (EVM0021509) and endonuclease 1-like isoform X2 (EVM0022275), which were significantly up-regulated at S25 and S30 compared to S20 (Figure 8; Table S6). We also found a variety of DEGs encoding proteases, including cysteine proteinases RD21A-like (EVM0007929), metacaspase-9-like (EVM0011448), aspartic proteinases (e.g., EVM0005045, EVM0005145, EVM0009402, and EVM0012500), vacuolar-processing enzyme-like (EVM0002134 and EVM0012232), senescence-associated proteins (EVM0009758 and EVM0031003), endoglucanase 1-like, e1, 9-like, and 7 (EVM0016472, EVM0017671, EVM0022249, EVM0032182, and EVM0003400), pectinesterase 2 and QRT1 (EVM0005897, EVM0021830, and EVM0032186), xyloglucan endotransglucosylase/hydrolases 30, 23, and 28 (EVM0023103, EVM0028198, and EVM0031923), xyloglucan galactosyltransferase GT11 (EVM0022425), and exosome complex component RRP45A-like (EVM0015970). Previous studies on homologs of our identified cysteine proteinases RD21A-like, metacaspase-9-like, aspartic proteinases, vacuolar-processing enzyme-like, and senescence-associated proteins suggest that these enzymes are crucial for effective PCD by degrading many essential cellular targets [26,31,44,61,62]. Endoglucanases, pectinesterases, xyloglucan endotransglucosylase/hydrolases, and xyloglucan galactosyltransferases may act to recycle carbohydrates as cells undergo PCD. Finally, homologs of exosome complex component RRP45A-like likely participate in DNA degradation [17].

Autophagy is a metabolic mechanism whereby cytoplasmic substances and organelles are degraded by lysosomes or vacuoles [63,64], and is a requirement for efficient PCD in many plant systems. Autophagic cell death is indicated by several distinct morphological features in plants, such as an increase in vacuolar and cellular size, movement of organelles into the vacuole for destruction, and subsequent cell death resulting from the vacuolar lysis [65,66]. Autophagic plays both pro-survival and pro-death roles in plants PCD [66]. The core autophagy mechanism is organized by an evolutionarily-conserved population of ATG (AuTophaGy-related) genes [65]. In our previous study, we found that during replaceable bud PCD of chestnut cv. ‘Tima Zhenzhu’, small vacuoles appearing in the cytoplasm fused to form a large vacuole, resulting in the eventual degradation of other cellular components [1]. Here, in accordance with this observation, we found that “regulation of autophagy” pathway genes were up-regulated at S25 and/or S30, compared to S20, including autophagy-related proteins 8f, 16, 8i-like, 8C-like, 3 and 13a (EVM0027149, EVM0009961, EVM0003991, EVM0004768, EVM0015502, EVM0002553, EVM0027459, and EVM0017527), cysteine protease ATG4-like (EVM0001368), ubiquitin-like modifier-activating enzyme ATG7 (EVM0017393), ubiquitin-like protein ATG12 (EVM0018088), and CBL-interacting serine/threonine-protein kinase 1-like and 6-like (EVM0001436 and EVM0010887) (Table S9). Additionally, the highly selective autophagy of soluble proteins is mediated by the ESCRT (endosomal sorting complex required for transport) complex, which requires vacuolar protein sorting (VPS) [67]. We found six DEGs (five up- and one down-regulated at S25 and/or S30, compared to S20)encoding VPS-associated proteins (A0A1S4DPC1, A0A1S4DKE1, A0A1S4CYV6, A0A1S4BA99, A0A1S3WXN7, and A0A1S4DLQ1) which interacted with each other and with the network of autophagy-related proteins (A0A1S3Y2M7, ATG8e, A0A1S4A7C8, ATG4, A0A1S4D5L9, ATG5, A0A1S4C6J8, ATG9, ATG2, A0A1S4DLP2, ATG13a, and A0A1S4BPT6) (Figure 7; Table S8). It appears that autophagic processes are necessary for the timely progression of replaceable bud PCD of chestnut cv. ‘Tima Zhenzhu’, and that replaceable bud PCD may be dependent on VPS-associated proteins.

5. Conclusions

Through transcriptomic profiling, we identified the DEGs and signaling pathways responsible for regulating replaceable bud PCD in chestnut cv. ‘Tima Zhenzhu’. Based on our cumulative results, we offer a hypothetical model of replaceable bud PCD consisting of three overlapping processes (Figure 9). First, ethylene signaling is activated during preparation for PCD, in order to regulate the activity of downstream targets. Next, during PCD initiation, the up-regulation of several TFs (including MYB, MADS-box, bHLH, and NAC TFs) induces an increase in cytochrome c expression and in the cytosolic Ca²⁺ content, activating the Ca²⁺-dependent signaling cascade. Finally, during PCD execution, the process of autophagy and the activity of proteases (i.e., cysteine proteinases RD21A-like, metacaspase-9-like, vacuolar-processing enzyme-like, and senescence-associated proteins for hydrolysis) work synergistically to clear the cell of cellular components. When this process is complete, the replaceable bud senesces and dies.