Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

RNA-Seq analysis reveals potential regulators of programmed cell death and leaf remodelling in lace plant (Aponogeton madagascariensis)

BMC Plant Biology volume 21, Article number: 375 (2021) Cite this article

Abstract

Background

The lace plant (Aponogeton madagascariensis) is an aquatic monocot that develops leaves with uniquely formed perforations through the use of a developmentally regulated process called programmed cell death (PCD). The process of perforation formation in lace plant leaves is subdivided into several developmental stages: pre-perforation, window, perforation formation, perforation expansion and mature. The first three emerging “imperforate leaves” do not form perforations, while all subsequent leaves form perforations via developmentally regulated PCD. PCD is active in cells called “PCD cells” that do not retain the antioxidant anthocyanin in spaces called areoles framed by the leaf veins of window stage leaves. Cells near the veins called “NPCD cells” retain a red pigmentation from anthocyanin and do not undergo PCD. While the cellular changes that occur during PCD are well studied, the gene expression patterns underlying these changes and driving PCD during leaf morphogenesis are mostly unknown. We sought to characterize differentially expressed genes (DEGs) that mediate lace plant leaf remodelling and PCD. This was achieved performing gene expression analysis using transcriptomics and comparing DEGs among different stages of leaf development, and between NPCD and PCD cells isolated by laser capture microdissection.

Results

Transcriptomes were sequenced from imperforate, pre-perforation, window, and mature leaf stages, as well as PCD and NPCD cells isolated from window stage leaves. Differential expression analysis of the data revealed distinct gene expression profiles: pre-perforation and window stage leaves were characterized by higher expression of genes involved in anthocyanin biosynthesis, plant proteases, expansins, and autophagy-related genes. Mature and imperforate leaves upregulated genes associated with chlorophyll development, photosynthesis, and negative regulators of PCD. PCD cells were found to have a higher expression of genes involved with ethylene biosynthesis, brassinosteroid biosynthesis, and hydrolase activity whereas NPCD cells possessed higher expression of auxin transport, auxin signalling, aspartyl proteases, cysteine protease, Bag5, and anthocyanin biosynthesis enzymes.

Conclusions

RNA sequencing was used to generate a de novo transcriptome for A. madagascariensis leaves and revealed numerous DEGs potentially involved in PCD and leaf remodelling. The data generated from this investigation will be useful for future experiments on lace plant leaf development and PCD in planta.

Background

Programmed cell death

Programmed cell death (PCD) is an essential developmental process that ensures the removal of cells, for tissue remodelling or in response to environmentally induced stress [1,2,3,4]. Plant developmental PCD is controlled by internal and external signals resulting in the organization of developing tissues [5,6,7]. Examples include embryonic suspensor deletion [8], aerenchyma formation [9], and xylem differentiation [10, 11]. Due to the experimental challenges associated with physically separating PCD destined cells spatially and temporally in plant model systems, there is presently little known about the genetic mechanisms that control developmental PCD. Plant systems with adjacent regions of differing cells fates arising from perforation formation during leaf morphogenesis can provide unique insight into the PCD process [12].

The lace plant model system

The formation of leaf perforations during development is unique and has been studied in plant families such as Araceae and Aponogetonaceae [12]. Monstera obliqua, M. deliciosa, and Aponogeton madagascariensis utilize PCD to dismantle and clear designated cells during early shoot development, creating perforations in the leaf blade.

The lace plant A. madagascariensis is an aquatic monocot that has recently emerged as a model system for studying PCD, due to the accessibility and predictability of PCD in developing leaves. Lace plant leaves are also thin and translucent, making them ideal for live-cell microscopy. Finally, the sterile propagation of whole lace plants in axenic environments create opportunities for pharmacological studies [3, 12].

The lace plant generates leaves with a perforated morphology, wherein specific cells bounded within the vasculature are removed via developmentally regulated PCD (Fig. 1A). Leaves of the lace plant emerge in a heteroblastic series from an underground corm, and while the first 3–4 leaves (termed imperforate leaves, Fig. 1B, C) to emerge do not form perforations at all during development, all successive emerging leaves become perforated (termed adult leaves). Early adult leaves in the “pre-perforation stage” (Fig. 1B, D) emerge from the corm in a furled form, with an abundance of anthocyanin and complete vein pattern [12, 13]. As pre-perforation leaves unfurl they transition to the “window stage” (Fig. 1B, E). Window stage is reached when cells in the central portion of areoles located between longitudinal and transverse veins undergo PCD, losing their anthocyanin and chlorophyll pigmentation (“PCD cells”, Fig. 1E). Cells proximal to the veins retain both anthocyanin and chlorophyll pigmentation and do not undergo PCD; and are therefore called “non-PCD cells” (NPCD cells, Fig. 1E). The perforations increase in size until halting 4–5 cell layers from the veins. Once perforation formation is nearly completed, the mesophyll cells at the NPCD cell edge of perforations transdifferentiate into epidermal cells. The leaves reach the mature stage once perforation expansion halts and NPCD cells achieve homeostasis (Fig. 1B, F).

Fig. 1
figure 1

The lace plant programmed cell death (PCD) model system. A-B Lace plant grown in axenic Magenta box culture with the pre-perforation stage (P), window stage (W), mature stage (M), and imperforate leaves (I). C Imperforate leaves emerge from the corm lacking anthocyanin and forming areoles with no perforations. D Successive pre-perforation stage leaves emerge from the corm with anthocyanin pigmentation (indicated by asterisks). E PCD can be seen actively occurring in the window stage of development. Between longitudinal and transverse veins, in spaces known as areoles, a gradient of cell death can be observed. Non-PCD cells (NPCD; bounded by white dashed lines) cells persist beyond maturity. PCD cells (bounded by black dashed lines) have lost their anthocyanins, are translucent and on the verge of death. PCD and perforation formation is complete in mature stage leaves (F) and anthocyanin pigmentation is visibly reduced, and homeostasis for NPCD cells is reached. Scale bars: A = 1 cm; B = 2 cm; C = 200 μm

Full size image

Lace plant PCD mechanism

The visible gradient of PCD within the areoles of window stage leaves represents a powerful model system for studying cellular changes across this gradient during developmental PCD [12]. The cellular changes that occur across the gradient of PCD in the lace plant have been well characterized, and some of the early events of lace plant PCD include the reorganization of the actin cytoskeleton, chloroplast ring formation around the nucleus, mitochondrial aggregation, and increased number of transvacuolar strands [12]. The manner in which organelles are taken up into the vacuole in membrane-bound vesicles suggest the involvement of autophagy [14]. Further examples of cellular changes include loss of mitochondrial membrane potential, DNA fragmentation, and activation of caspase-like proteases before vacuolar collapse, plasma membrane shrinkage and cell wall degradation [15, 16].

In spite of the well characterized progression of PCD, little is known about the molecular mechanisms that control lace plant PCD regulation and execution [17, 18], in part due to a lack of molecular information for the Aponogetonaceae family. The advancement of comparative RNA sequencing (RNA-Seq) analysis between PCD and NPCD-like cells in other plant models has helped characterize differentially expressed genes (DEGs) that resemble PCD regulators. RNA-Seq analysis of separated embryonal mass and suspensor cells of Picea abies has shown that a spruce homolog of bax inhibitor-1 transcript is upregulated in early PCD suspensor cells and plays a role in regulating vacuolar cell death in suspensor cells [19]. Moreover, RNA-Seq analysis of early and late cavern forming leaf aerenchyma cells of Typha angustifolia revealed expansins, calmodulin-like proteins and pectinases transcripts that were directly related to lysigenous aerenchyma induction [20].

To date there is only one transcriptome study, Rantong et al. (2016) [18], which investigated lace plant leaf stages using complementary DNA-amplified fragment length polymorphism (cDNA-AFLP). This study identified 79 annotated DEGs which are involved in processes such as photosynthesis, stress responses, pathogen defence, and PCD. Importantly, their results suggested that ubiquitin-proteosome machinery may be involved in lace plant PCD. However, as cDNA-AFLP captures only a fraction of the transcriptome, how expression patterns across the entire A. madagascariensis genome change during PCD remains unknown.

In this study, we used high-throughput RNA-Seq to compare global transcriptome expression profiles of different stages of lace plant development as well as PCD and NPCD cells. We separated PCD and NPCD cells within window stage leaves by laser capture microdissection, allowing us to identify DEGs that may be involved in PCD and survival. We also identified, through differential expression analyses of RNA-Seq data, genes highly expressed in PCD and NPCD cells that are potential PCD inductors, executors and/or regulators. To identify key regulators of lace plant leaf remodelling, we additionally tested for genes that are highly expressed among perforating and non-perforating leaves.

Our objectives are to identify and compare DEGs among different leaf developmental stages and between PCD and NPCD cells. We hypothesize that imperforate and mature leaves have significantly higher levels of expression of genes involved in photosynthesis and negative regulation of PCD while pre-perforation and window leaves will have significantly higher levels of expression of genes responsible for anthocyanin biosynthesis, caspase-like activity, cell wall organization, and pro-PCD regulation.

Results and discussion

RNA-Seq data overview

To identify potential regulators of lace plant developmental PCD and leaf remodelling, we generated a novel lace plant transcriptome and identified DEGs in comparisons of pre-perforation, window, mature and imperforate leaf stages, and NPCD and PCD cell types using RNA-Seq. Eighteen paired-end RNA-Seq libraries were generated from three biological replicates of each imperforate, pre-perforation, window, mature leaf stages, NPCD cells and PCD cells. The Illumina (San Diego, CA, USA) NovaSeq6000 sequencing platform was used for paired-end sequencing at Génome Québec (Montréal, QC, Canada), with a 100 bp read length. A total of 1,320,261,351 reads were generated, and data for individual biological libraries were deposited to the NCBI SRA database with the following SRA accession IDs: SRR10524134-SR10524151 and BIOPROJECTID: PRJNA591467. After read filtering, 1,288,318,561 reads remained and 1,102,201,639 reads aligned concordantly (Additional file 1).

We assembled 908,119 transcripts with an N50 length of 1,041 bp, and 49.9% GC content from eighteen RNA-Seq libraries. These transcripts translated to 437,825 protein coding genes (Table 1). Gene Ontology (GO) annotated DEGs across leaf stages and cell types accounted for 4,339 of the 106,222 (4.08%) A. madagascariensis GO annotated genes across all leaf stages and cell samples with a 1% adjusted P-value cut-off (Fig. 2A, B). Of the 10,416 DEGs, 2808, 313, 1541, and 1267 genes were upregulated exclusively in pre-perforation, window, mature, and imperforate leaves respectively (Fig. 2A). Between cell types, 482 and 166 genes were exclusively upregulated (at least twofold) in PCD and NPCD cells, respectively. Remaining DEGs were expressed mutually between different combinations of leaf stages (Fig. 2).

Table 1 Number of differentially expressed and GO annotated genes across leaf and cell-type samples
Full size table
Fig. 2
figure 2

Overview of differentially expressed genes in lace plant leaves and cell types. Venn diagram showing the number of mutual and exclusive differentially expressed genes (DEGs) between lace plant leaf stages (A) and between NPCD and PCD cells (B); fold change > 2.0, P < 0.01. Heatmaps of gene expression levels of DEGs plotted as log2(FPKM + 1) for indicated samples and biological replicates for leaf stage comparison (C) and comparison of NPCD and PCD cells (D). Low and high gene expression levels are shown in purple and yellow, respectively. Below heatmaps, biological replicates 1–3 are indicated by P, pre-perforation; W, window; M, mature; and I, imperforate. Supporting data for gene expression found in Additional file 3 and 5

Full size image

Transcriptomic profiles of the lace plant developmental leaf stages

BLAST-based comparisons of assembled transcripts yielded 106,222 GO annotated genes that were homologous to sequences in public databases. Of the annotated genes, 30,932 (29.12%) were most similar to Arabidopsis thaliana, 6,509 (6.12%) to Oryza sativa and 692 (0.65%) to Zea mays.

A. madagascariensis leaf DEGs were divided into four main clusters generated through the tree cutting method of expression patterns across four leaf developmental leaf stages (Fig. 3, Additional file 3). RNA-Seq data were divided into the pre-perforation cluster (4544 DEGs), the window cluster (718), the mature cluster (1572), and the imperforate cluster (1387) (Fold change > 2.0, P < 0.01, false discovery rate (FDR) = 1%; Additional file 3). GO enrichment analyses identified biological functions enriched in the four main clusters based on expression patterns across the four stages of leaf development (FDR = 1%, Additional file 3).

Fig. 3
figure 3

Transcriptomic analysis of lace plant leaf developmental stages. Top: Four main clusters grouped by highest expression in respective pre-perforation (P), window (W), mature (M) and imperforate (I) leaf stage biological replicates (n = 3). For each cluster, individual DEG expression values (shown as the transformed log2(FPKM + 1) values) are plotted as grey lines and the mean expression profile is highlighted in blue. The total number of DEGs per cluster is shown below each plot (P < 0.01, fold change > 2.0). Bottom: Heatmaps of composite gene expression for indicated proteins, with green and red corresponding to high and low gene expression, respectively. Supporting data are found in Additional file 3

Full size image

The pre-perforation cluster, representing genes that were most highly expressed in the pre-perforation stage leaves and exhibited reduced expression in subsequent stages, included genes encoding for proteins involved in flavonoid biosynthesis, anthocyanin biosynthesis, ethylene-activated signalling pathway, endopeptidase activity, autophagosome formation, and regulation of PCD. The window stage leaf cluster represented genes that were most highly expressed in the pre-perforation and window stages and then reduced in subsequent stages. This cluster included genes encoding for proteins involved in ATP-binding, ATPase activity, ion binding, response to auxin, response to oxygen-containing compounds, peroxidase activity, and hydrolase activity. Both the mature and imperforate clusters represented genes that increased in expression in later leaf development stages, and both contained genes encoding for proteins involved in photosystem I and II, chlorophyll-binding, light-harvesting complex organization, catalytic activity, ion binding, and cell wall biosynthesis.

Taken together, clustering data support the hypothesis that growth and organizational processes are enriched in developing pre-perforation and window stage leaves where metabolic processes must occur to fuel development. Many energy-related metabolic processes occur in the mature and imperforate leaves where development is near completion and flavonoid synthesis is reduced (Fig. 3, Additional file 3). General patterns of gene expression for select biological functions across the leaf clusters showed that genes involved in photosynthesis and negative regulation of PCD are expressed at higher levels in mature and imperforate leaves where PCD is not as active and homeostasis is reached. This suggests that imperforate and mature leaves are the major site of photosynthesis, whereas pre-perforation and window leaves specialize in growth, responding to hormones, and executing PCD. All clusters demonstrated high expression of genes involved in cell wall modifying enzymes such as pectinesterases. However, pre-perforation and window leaves possess a greater number of highly expressed expansins, pectinesterases and subtilisin-like proteases than mature and imperforate leaves. Pre-perforation and window stage leaves represent the culmination of many developmental processes such as regulation of PCD, cell wall organization, lignin, and stomata development. Genes encoding proteins involved in hormone synthesis and transport were found to be differentially expressed between leaf stages. Pre-perforation and window stage leaves had high expression levels of genes encoding for auxin biosynthesis and transport, abscisic acid (ABA) biosynthesis, brassinosteroid (BR) biosynthesis, cytokinin (CK) biosynthesis, gibberellin (GA) biosynthesis, ethylene biosynthesis, ethylene receptor activity, ethylene signalling pathway, jasmonate biosynthesis and salicylic acid (SA) response. Mature and imperforate leaves expressed similar levels of ABA transport genes in comparison to pre-perforation and window stage leaves. In total, early developing leaves revealed an expression pattern of leaf development similar to other monocots like Agave deserti and Z. mays [21] where expression of most transcription factors (TFs) and hormone biosynthesis tend to be at their highest. Likewise, mature leaves express genes related to photosynthesis [22].

Insights from comparative transcriptomics of NPCD and PCD cells

Laser capture microdissection was used to separate NPCD and PCD cells from A. madagascariensis window stage leaves to reveal potential regulators of PCD by RNA-Seq analysis. NPCD and PCD transcriptomes were divided into two main clusters based on cluster analysis expression patterns (Fig. 4) using the tree cutting method. This cluster analysis identified approximately 431 genes that were differentially expressed in either NPCD (326 DEGs) or PCD (105 DEGs) samples with a minimum of a twofold change (P < 0.01). In comparison to upregulated genes found in leaf stages; this represents a small fraction of the transcriptome which may be a result of limited biological replicates used in this experiment. We tested for GO enrichment (FDR = 1%; Additional file 5) in each of the NPCD and PCD clusters. The NPCD cluster contained genes encoding for proteins involved in respiratory burst activity, leaf senescence (including negative regulation of senescence), protein autoubiquitination, and the ethylene-activated signalling pathway. The PCD cluster contained genes encoding for proteins involved in ethylene biosynthesis, cell wall modifiers, protease inhibitors, and ROS generation, and PCD regulation (Fig. 4).

Fig. 4
figure 4

Transcriptomic analysis of NPCD vs PCD cells of the lace plant. Top: Two main clusters grouped by highest expression in NPCD and PCD cell biological replicates (n = 3). For each cluster, individual DEG expression values (shown as the transformed log2(FPKM + 1) values) are plotted as grey lines and the mean expression profile is highlighted in blue. The total number of DEGs per cluster is shown below each plot (P < 0.01, fold change > 2.0). Bottom: Heatmaps of composite gene expression for indicated proteins, with green and red corresponding to high and low gene expression, respectively. Supporting data are found in Additional file5

Full size image

Comparing the relative GO counts between NPCD and PCD clusters (Fig. 5) we found that NPCD cells upregulated more genes encoding for proteins involved in flavonoid biosynthesis, cysteine protease activity, metalloprotease activity, positive regulation of autophagy, PCD regulation, and cell wall organization. Conversely, PCD cells upregulated more genes involved in ethylene biosynthesis, photosystem II and I, BR biosynthesis, and cutin biosynthesis. NPCD and PCD cells expressed similar numbers of genes in GO categories for homeobox, myeloblastosis (MYB), zinc finger TF families and serine endopeptidase activity, suggesting that these families act as regulatory elements across both cell types during differentiation.

Fig. 5
figure 5

Heat map of GO category orthologs differentially expressed in leaf stages, NPCD and PCD cells. Heatmaps show composite select GO term counts normalized by the cluster with highest count for specific biological processes. The number of GO category orthologs were compared among lace plant leaf stages (left) and between NPCD and PCD cells (right). Colour gradient ranges from white (zero genes upregulated) to dark blue (highest number of genes upregulated). Supporting data are found in Additional file 3 and 5

Full size image

Additionally, five DEGs were selected for fold-change expression validation by quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) across leaf stages and cell types (Additional file 4). The transcripts for Bag5, expansin A-29, aquaporin 4–4, anthocyanin regulatory protein C1, nuclear transcription factor YC-1, and α-tubulin were chosen as validation targets. The fold-change expression of the selected genes normalized by α-tubulin followed a similar pattern to that seen in the cluster analysis among different stages of leaf development and between cell types (Additional file 4) and were deemed valid.

Transcription factors

Most of the DEGs that encode for TFs were found in early developing leaves (Fig. 5, Additional file 3). TFs identified in the pre-perforation leaf cluster included YABBY, MYB, bHLH (basic Helix-Loop-Helix), Zn finger, GATA, C2C2 CO (constans), homeobox, DOF (DNA-binding with one finger), TCP (teosinte, branched1, cycloidea and PCF), and trihelix families (Fig. 5, Additional file 3). These identified families presumably underpin the expression of a broad array of genes important for early leaf development, in particular, the establishment of axial polarity, stomata development, and vascular tissue, as seen in O. sativa and A. thaliana [23,24,25,26,27,28,29,30,31]. Mature and imperforate leaves expressed a greater number of genes than pre-perforation and window leaves that encode for MADS-box and NAC (no apical meristem) family TFs which may underpin the promotion of photosynthetic development, lignin, wax and secondary cell wall development which is enhanced in these later developing leaves.

DEGs encoding for TFs were also detected in NPCD and PCD cell DEG expression profiles (Fig. 5, Additional file 5). Trihelix and WRKY families were expressed at a high level in NPCD cells. These identified TFs were most likely responsible for transcribing genes involved in stress response and are known to be upregulated in mesophyll cells of developing O. sativa leaves [21]. Two constan family TFs were highly expressed in PCD cells. Constan TF families may be involved in mediating PCD cell collapse or suppressing anti-PCD genes. Both NPCD and PCD cells upregulated an equal number of MYB, Zn finger, and homeobox family TFs. There may be MYB TFs in both types of cells but these two groups possess opposite responsibilities in terms of promoting and suppressing flavonoid biosynthesis. The differences in transcriptional regulation of gene expression between perforating and non-perforating lace plant leaves may control the key programming events that result in differential growth.

Plant hormones

We detected higher expression patterns for auxin, ABA, CK, GA, ethylene, and jasmonate hormone biosynthesis genes in the pre-perforation cluster (Fig. 5, Additional file 3). The window stage cluster contained the highest levels of expression for BR biosynthesis and ethylene receptor activity. The mature and pre-perforation stage clusters contained the highest levels of expression for ABA hormone transport (Fig. 5, Additional file 3). These results support previous findings of several plant-specific hormones being involved in PCD signalling, including SA, jasmonic acid, ABA, GA, and ethylene [32,33,34]. The hormone biosynthesis patterns observed in early lace plant development are similar to the monocot leaves of A. deserti, Agave tequilana and Z. mays [21, 22].

To date, several pharmacological whole plant experiments have revealed how lace plant perforation formation is dependent on auxin, and ethylene biosynthesis [35,36,37]. Further work is required to disentangle the roles of each plant hormone in mediating lace plant leaf development from perforation formation, outside of ethylene which has been studied extensively.

Expression patterns in the NPCD cell cluster included more highly expressed genes related to auxin transport, auxin signalling pathway, GA biosynthesis, ethylene activated signalling pathway, jasmonate biosynthesis and SA response in comparison to the PCD cell cluster (Fig. 5, Additional file 5). PCD cells upregulated 1 gene encoding for ethylene biosynthesis relative to NPCD cells. Ethylene biosynthesis and ethylene receptors are involved in promoting PCD in cells destined to die in the areoles of lace plant leaves [35, 37]. Additionally, PCD cells upregulate more genes associated with BR biosynthesis in comparison to NPCD cells. BRs are believed to mediate the timing of ROS production, and in turn, PCD execution in tapetal cells of Solanum lycopersicum [38]. BRs may also play a similar role in PCD cell triggering, as suggested by the higher expression of genes involved with BR synthesis. NPCD cells upregulated ethylene response factor (ERF) RAP2-3 in comparison to PCD cells. ERF-RAP2-3 has been identified as playing an important role in ethylene mediated hypoxia stress in A. thaliana seedlings [39] and may play a role in protecting NPCD cells from PCD cells which accumulate superoxide during PCD execution. Together, our results support the hypothesis that plant hormones are involved in PCD, and targeted approaches are needed to disentangle their specific roles and functions.

Anthocyanin biosynthesis enzymes

Forty-six upregulated genes were categorized as enzymes with flavonoid biosynthesis, isomerase, hydrolase, oxidoreductase or lyase activities in the pre-perforation leaf cluster while < 10 genes for these enzymes were found in each of the window, mature, and imperforate leaf clusters (Fig. 5, Additional file 3). Eight of the upregulated genes in the pre-perforation cluster encoded for the flavonoid biosynthesis pathway. Pre-perforation and window leaves upregulated genes that promote the biosynthesis of early stage and late-stage flavonoids as well as downstream transferase enzymes for promoting the synthesis of anthocyanins, flavonols and anthocyanidins. NPCD cells also upregulated 5 genes that encoded for flavonoid biosynthesis in comparison to PCD cells (Fig. 5, Additional file 5). Of the genes that are involved in flavonoid biosynthesis only flavonone 3-dioxygenase 2, which produces dihydroflavonol [40], was upregulated in mature leaves and imperforate leaves as well as PCD cells (Additional files 3 and 5). Upregulated genes that encode for enzymes that promote anthocyanin biosynthesis are summarized in Fig. 6.

Fig. 6
figure 6

DEGs involved in flavonoid biosynthesis. General late-flavonoid biosynthetic pathway (grey boxes) and genes (cyan boxes) expression in pre-perforation (P), window (W), mature (M) and imperforate (I) leaf stages (4-box expression comparison strings) or NPCD and PCD cell samples (2-box expression comparison strings). Detailed gene expression data provided in Additional file 3 and 5. Expression levels displayed in log2(FPKM + 1) across samples. Green indicates high expression values; red indicates low expression values. Pathway modified from Argout et al., (2008) [41]

Full size image

The most notable feature of lace plant leaf development is the visible gradient of anthocyanin pigmentation during PCD in early developing leaves. Previous studies have investigated the role of exogenous ROS and antioxidants on lace plant leaf development, and found that they are key regulators of the establishment of perforation formation in the lace plant model system [13]. Pre-perforation and window stage leaves have the highest levels of anthocyanin compared to imperforate and mature leaves [13, 42]. Our results are consistent with these findings, as pre-perforation leaves, window leaves, and NCPD cells expressed genes encoding enzymes that promote the biosynthesis of glucoside constituted anthocyanins (Fig. 6). The upregulation of anthocyanidin-3-O-glucosyltransferases (A3OGT), anthocyanidin-5,3-O-glucosyltransferases (A5,3OGT), and anthocyanin-5-aromatic acyltransferase (A5AAT) suggests that anthocyanidin-3-O-glucosides, anthocyanidin-5,3-O-glucosides, and anthocyanidin-3-glucoside-6-hydroxycinnamoyl glucosides are being actively synthesized in early-stage leaves, leading to the accumulation of anthocyanins. The presence of pink-red coloured anthocyanin vacuoles in window stage leaves is consistent with the hypothesis that the vacuole has a pH of ~ 5.5 pH which favours the formation of flavylium cations [43].

Recent attention has been drawn to the formation of anthocyanin vacuolar inclusions (AVIs) in window stage leaves, which then dissipate as leaves enter maturity (Gunawardena lab, 2020 unpublished data). Kallam et al. (2017) [43] recently reported that the composition of anthocyanin substituents or decorations determined the solubility of AVIs in tobacco lines. The aromatic acylation of anthocyanidin 3-O-glucoside by anthocyanidin 3-O-glucoside-6″-O-coumaroyltransferase (A3T) promotes the formation of AVIs, while the rhamnosyl decorations decrease the formation of AVIs by competition. The absence of upregulated aromatic rhamnosyltransferases or malonyltransferase DEGs (which inhibit AVI formation) in early leaves suggests a role for AVIs in mediating proper leaf development [43]. NPCD cells upregulated genes encoding for anthocyanidin 3-O-glucosyltransferases, consistent with the promotion of ROS scavengers and protection from PCD. The identification of anthocyanin compounds across developmental leaf stages and NPCD cells is currently under investigation using global mass spectrometry.

Potential regulators of programmed cell death

Cell wall modification enzymes and aquaporins

Lace plant perforation formation relies on thin cuticle erosion, removal of polysaccharides and degradation of cellulose [44] to mediate perforation expansion. The deposition of suberin is imperative to prevent infection and loss of nutrients in NPCD cells while neighbouring PCD cells collapse. We found that the window stage cluster had high expression of the greatest number of suberin biosynthesis-related genes (Fig. 5, Additional file 3), which is consistent with the hypothesis that suberin deposition is most active during the window stage of development [12].

Pre-perforation and window stage clusters contain high expression of 67 and 7 genes, respectively, that were categorized under either cell wall organization or cell wall biosynthesis. Mature and imperforate clusters showed high expression of 9 and 5 genes respectively. The NPCD cluster showed high expression of 4 genes, while the PCD cluster showed high expression of 1 gene (Fig. 5, Additional file 5). All clusters contained higher expression for orthologs of hydrolases, pectinesterases, glucosidases, and xyloglucan glycosyltransferase activities. Mature and imperforate leaves upregulated no genes that encode for expansin enzymes.

The facilitation of cell expansion and cell wall loosening is important for not only the growth of developing lace plant leaves but also the execution of PCD cells and reorganization of NPCD cells. Pre-perforation and window stage leaves upregulated several genes that encode for xyloglucan endotransglucosyltransferases, expansins (expansin-A29 detected by qRT-PCR, Additional file 4), and pectinesterases (Fig. 3) which are responsible for loosening and reorganizing the cell wall during growth [45,46,47]. Mature and imperforate stage leaves followed a similar pattern but did not upregulate genes for expansins (Fig. 3). This likely indicates that this activity is no longer needed once leaf maturity is reached.

Aquaporins play an important role in cell expansion by controlling water uptake [48]. Pre-perforation and window stage leaves upregulated several genes encoding for tonoplast intrinsic protein (TIP), nodule intrinsic protein (NIP), and plasma membrane intrinsic protein (PIP) aquaporins. NPCD cell upregulated a TIP4-4 aquaporin gene (Fig. 4). TIP4-4 was detected by qRT-PCR (Additional file 4). PCD cells upregulated a TIP2-3 gene and mature and imperforate leaves upregulated both. The expression of aquaporin genes in all samples suggests that cell expansion by vacuole enlargement is needed throughout leaf development, and as well as in differentiating NPCD and PCD cells to control the progression of expansion or cell burst [48].

Heat shock proteins

Heat shock proteins (Hsps) are synthesized in response to stress to maintain homeostasis by refolding proteins before they become irreversibly denatured [49, 50]. The pre-perforation leaf cluster contained 25 genes and the window leaf cluster contained 2 genes categorized as protein folding. Mature and imperforate leaf clusters contained 1 and 2 genes under unfolded protein folding, respectively (Fig. 3, Additional file 3). The NPCD cell cluster contained 1 gene and the PCD cell cluster contained 0 genes categorized under unfolded protein folding (Additional file 5). ATP-dependent molecular chaperones such as Hsp81-1, Hsp81-3, and Hsp70-15 were upregulated in pre-perforation and window stage leaves.

In O. sativa protoplast models, mtHsp70 overexpression prevents ROS burst and PCD in protoplasts under oxidative stress [51]. In contrast to suppressing stress induced PCD, Hsp70 of Capsicum annuum promotes the hypersensitive response in infected Nicotiana benthamiana leaves [52] by nuclear localization of an effector protein. Hsp70s and their respective Bcl-2-associated athanogene (Bag) proteins can modulate animal PCD and many cellular processes [1], warranting their investigation in lace plant PCD. Experimental treatment with the Hsp70 inhibitor chlorophenylethynylsulfonamide (PES-Cl) caused a significant decline in the number of perforations, caspase-like activity and anthocyanin levels in window stage leaves [42], suggesting Hsp70 plays a role in mediating lace plant PCD. Hsp70 proteins are developmentally regulated and significantly higher in pre-perforation and window stage leaves, which is consistent with the expression pattern of an Hsp70-15 gene (Fig. 3). While our results support the hypothesis that Hsp70 activity affects lace plant leaf development at a threshold level, where it localizes is still unknown.

Genes for lace plant homologs of O. sativa Hsp81-1 and Hsp81-3 were transcriptionally upregulated in pre-perforation stage leaves (Fig. 3, Additional file 3). Hsp81s promote salt stress tolerance in O. sativa and over-expression experiments in A. thaliana show that it promotes heat tolerance [53, 54]. Proteomic models predict that AtHsp81 can form a complex with AtHsp70 [55] most likely for protein quality control, and our results suggest that both lace plant homologs of Hsp70 [42] and Hsp81 are being synthesized for maintaining protein homeostasis during early leaf development. Additional pharmacological whole plant experiments are required to improve characterization of Hsps and Bag protein function in lace plant development.

Bag proteins

The Bag protein family has gained recent attention in the field of plant developmental biology for their role in mediating PCD, senescence, and autophagy in several plant systems [1, 55]. Bag proteins are major nucleotide exchange factors for Hsp70 [56] and help accelerate the Hsp70 protein fold cycle and modulate PCD pathways in animals and plants [57, 58]. The Hsp70 co-chaperone Bag7 was more highly expressed in pre-perforation and window leaves. Bag3 and mitochondrial (mt) Bag5 were upregulated in mature and imperforate leaves (mtBag5 detected by qRT-PCR, Additional file 4). Expression values for mtBag5 were significantly higher in NPCD cells in comparison to PCD cells (Fig. 4, Additional file 5).

In the lace plant, we found that pre-perforation and window stage leaves upregulated genes encoding homologs of A. thaliana for Bag7 (Fig. 3, Additional file 3). Bag7 is ER-localized in A. thaliana and is involved in the unfolded protein response, and can localize to the nucleus to interact with multiple proteins to modulate PCD pathways [59, 60].

Mature and imperforate lace plant leaves upregulated Bag3 and Bag5 homologs relative to pre-perforation and window leaves (Fig. 3), whereas NPCD cells upregulated a Bag5 homolog compared to PCD cells (Fig. 4). The role of Bag3 involvement in human cargo mediated autophagy has made it a potential cancer therapy target, but its role in plant PCD is unknown [61]. In A. thaliana leaf systems, Bag5 has been found to bind heat shock cognate 70 (Hsc70) and localize to the mitochondria to promote ROS generation and clearing of chlorophyll while leaves are under senesence [58]. The fact that Bag5 genes were upregulated in NPCD cells provides an opportunity to investigate the possible role of this gene in separating NPCD cells from PCD cells or preventing mitochondrial burst [62]. Plant Bag5 plays a role in regulating Ca2+ in the mitochondria and possibly the outcomes of mitochondrial stress response in NPCD cells. Bag family proteins seem to play an important regulator role in plant PCD but their precise biochemical roles are still unknown in plants [61].

Autophagy-related proteins

Dauphinee et al. (2019) found that autophagy predominately contributes to cell survival and that there is no clear evidence for the direct involvement of autophagy and the induction of PCD during the perforation formation of lace plant leaves [14]. Only pre-perforation and window stage leaves were found to upregulate genes in the GO category of autophagosome assembly such as homologs of AuTophaGy-related protein 16 (Atg16) and Atg18a (Fig. 3 and 5, Additional file 3). The pre-perforation cluster contained a SNF1-related protein kinase catalytic subunit alpha KIN10 and a lysosomal amino acid transporter 1, both of which are involved in autophagy regulation. Mature leaves, imperforate leaves, and NPCD cells (in comparison to PCD cells, Fig. 3 and 4) upregulated WRKY33, a TF involved in the positive regulation of autophagy [63]. WRKY33 was upregulated in NPCD cells relative to PCD cells, suggesting a requirement for stress resistance in NPCD cells during perforation development. WRKY33 is important for plant resistance to necrotrophic pathogens and interacts with the Atg18a in the nucleus [64] in response to biotic stress.

We detected higher expression of genes for Atg18a and Atg16 in pre-perforation and window stage leaves, suggesting changes in the regulation of autophagosome formation and mitophagy occur during plant developmental PCD. There are eight members in the AtATG18 gene family (AtATG18a–AtATG18h) [63]; each member has a different expression pattern, and only AtATG18a shows an increased transcript level under starvation conditions and during senescence in A. thaliana. AtATG18a expression is also upregulated and is required for autophagy during oxidative, salt, and osmotic stresses. RNA interference (RNAi) of AtATG18a makes plants autophagy-defective and more sensitive to various stress conditions [63, 65,66,67]. Atg16 oligomerizes to form an Atg12-Atg5·Atg16 complex that is essential for autophagy [68].

The A. thaliana protein kinase SNRK KIN10 has been shown to induce several autophagy genes such as Atg8, a protein found to be developmentally regulated in lace plant leaves [69, 70]. During early leaf development, the level of expression of genes involved in photosynthesis is low in comparison to mature and imperforate leaves. Prolonged stress often causes decreased photosynthesis and increased ROS accumulation, which can trigger PCD pathways [71]. Autophagy is used in early leaf development to provide energy while photosynthesis and chlorophyll genes are not upregulated until maturity and optimal photosynthesis activity is reached. It is likely WRKY33 is transcribed to increase autophagosomes formation during early leaf development and to promote survival in stressed NPCD cells.

Regulators of programmed cell death

Genes falling under the GO category “PCD regulation” were differentially expressed across lace plant leaf stages and between NPCD and PCD cells. We found that the pre-perforation leaf cluster contained homologs for L-type lectin domain kinase IX.I and mechanosensitive ion channel protein 10 (MSL10). MSL10 promotes PCD in response to pathogen invasion, and mechanical stress in A. thaliana [72,73,74]. Additionally, the pre-perforation cluster contained genes encoding for a BOI-related E3-ubiquitin ligase 2. BOI E3-ubiquitin ligases are capable of inhibiting PCD by limiting α-picolinic acid generation and by ubiquitination of apoptotic inhibitors [75, 76]. The upregulation of MYB33 in pre-perforation and window stage leaves (Fig. 3, Additional files 3) is of particular interest for further lace plant investigations, due to its role in the promotion of PCD in anthers and seeds of other monocots like Hordeum vulgare and O. sativa [77].

The mature and imperforate leaf stage clusters contained genes encoding a BOI related E3 ubiquitin ligase 2 and a respiratory burst oxidase homolog-F (Rboh-F) respectively. BOI related E3 ubiquitin ligase 2 genes negatively regulate PCD by suppressing ROS generation [75, 76] which is consistent with mature and imperforate leaves where perforation formation and superoxide accumulation is less active [13, 42]. Rboh-F is mostly responsible for ROS generation by ABA signaling in A. thaliana systems [78] and is implicated in immunity.

We also observed differential expression patterns for PCD regulation genes across NPCD and PCD cells (Fig. 5, Additional file 5). NPCD upregulated genes encoding for aspartyl protease AED3 and ERF-RAP2-3 which have been previously described as pro-PCD contributors. PCD cells upregulated a gene encoding for primary amine oxidase 1 (PAO1) relative to NPCD cells. PAOs have been shown to play a role in generating ROS in differentiating tissue during organ development and during PCD [79, 80]. Pro-PCD genes such as aspartyl protease AED3 upregulated in NPCD was an unexpected result and may potentially be explained by these genes having a function to promote a stress response or senescence in NPCD cells during PCD activation.

Plant proteases and programmed cell death

The pre-perforation, window, mature and imperforate clusters were all found to contain upregulated genes encoding for enzymes with endopeptidase activity. The pre-perforation leaf cluster shows 59 genes encoding for enzymes with aspartic, serine, and cysteine endopeptidase activity (Fig. 5, Additional file 3). NPCD cells upregulated 27 genes encoding for enzymes with aspartic endopeptidase activity (2 for PCD cells), 1 gene with serine activity for both cell types, 4 with cysteine activity (versus for 0 in PCD cells), and 1 with metalloprotease (0 for PCD; Fig. 4 and 5, Additional file 5).

Previous research has pinpointed the roles of caspase-like enzymes in plant development as PCD initiators or executors [81,82,83]. Subtilisin-like proteases have potential PCD regulation roles with an autocatalysis activity-containing domain that re-enters the cell once the prodomain is removed to execute PCD, and all lace plant leaf stage clusters upregulated several subtilisin-like proteases. The similar expression of subtilisin-like proteases across leaf stages is consistent with a role in leaf remodelling and homeostasis. All leaf stages and NPCD cells upregulated a Bowman-birk serine protease inhibitor. Serine protease inhibition activity may indicate a role for coordinated inhibition of serine endopeptidase activity for proper PCD execution.

Multiple cysteine proteases have documented roles in developmental plant PCD. For example, cysteine protease 1 mediates tapetal PCD in A. thaliana [82] and cysteine protease R2D1A is found to enhance PCD in innate immunity of A. thaliana [83, 84]. Cysteine protease activity can be tightly controlled by the activity of serpin and Kunitz protease inhibitors (reversible inhibition [83]), which are upregulated in mature and imperforate leaves in the form of Kunitz protease inhibitor 2 and cysteine protease inhibitor A. Protease cascades may trigger lace plant leaf PCD. Inhibitors such as the Kunitz and cysteine inhibitors in mature leaves could play a role in stopping the effector phase of cysteine protease activity and preventing PCD from becoming active again [85]. We found that both genes for cysteine protease 1 and RD21A previously mentioned are transcriptionally upregulated in NPCD cells. The reason(s) for the accumulation of mRNA for these proteases in cells destined to survive is still uncertain, this could suggest they play an important role in mediating PCD cell collapse.

Aspartyl proteases also regulate plant developmental PCD; specifically, the deletion of reproductive tissues [86, 87]. In lace plants, pre-perforation, window stage leaves, and NPCD cells all had high expression of genes encoding for aspartyl protease AED3, which is involved in PCD (Fig. 3 and 4; Additional file 3 and 5). Aspartyl protease AED3 may be upregulated for transdifferentiating NPCD cells into endodermis during PCD progression.

A metalloproteinase 2-MMP gene was upregulated in NPCD cells in comparison to PCD cells (Fig. 4). SI2-MMPs have been found to inhibit epidermal cell death in S. lycopersicum L and, in contrast, 2-MMP promotes early senescence in A. thaliana [88, 89]. The upregulation of 2-MMP in NPCD cells may indicate their role in differentiating NPCD from PCD cells, or inhibiting PCD execution.

The protease substrate phosphoenolpyruvate carboxykinase 1 (PEPCK1) can be cleaved by A. thaliana metacaspase 9 (MC9) and MC9 in turn promotes the clearance of root xylem tissue [90] and possibly gluconeogenesis [91]. A gene encoding an A. thaliana PEPCK1 homolog [85] was upregulated in our imperforate leaf cluster (Fig. 3, Additional file 3). The cleavage of PEPCK1 by MC9 may promote gluconeogenesis in imperforate leaves and supports the hypothesis that imperforate leaves serve to generate energy for the lace plant before undergoing rapid senescence.

Interestingly, our RNA-Seq analysis did not reveal differential expression across leaf stages or cell types of any lace plant genes encoding for vacuolar processing enzymes (VPEs), which are known to play a major role in developmental PCD of lace plant leaves [35, 92]. Using qRT-PCR methods, Rantong & Gunawardena (2018) showed that transcript levels of AmVPE-1 and -2 are significantly upregulated in early developing leaves and window leaves (normalized to AmActin). Previous work has highlighted the fact that VPE activity is important for vacuolar collapse in lace PCD cells. The absence of this process from results of our study and that of Rantong & Gunawardena (2018) suggests that autoprocessing and post-translational modification of the VPE protein might explain its functional activity, rather than just higher accumulation of mRNA [92].

Conclusions

The cellular dynamics and chronological events of lace plant leaf PCD are well documented. Here, we investigate the molecular basis of this process by characterizing the transcriptomic profiles of different stages of leaf development and PCD and NPCD cells isolated from window stage leaves. We profiled DEGs to summarize genes controlling the mechanism of developmental PCD and leaf remodelling (Fig. 7).

Fig. 7
figure 7

Summary of differentially expressed genes involved in lace plant leaf remodelling, based on RNA-Seq data. Summary of transcriptionally upregulated genes involved in lace plant leaf remodelling and differentiation of NPCD and PCD cells. GO terms bounded in white shaded boxes are most highly expressed GO counts in respective sample cluster (Additional files 3 and 5)

Full size image

Based on comparative transcriptomics our results support the hypothesis that NPCD and PCD cell differentiation is mediated by a differential balance of plant hormones and TF activities that both promote and limit the PCD pathway. GO enrichment analyses of DEGs suggest that autophagy, cell expansion, protease activity, ROS generation, and flavonoid biosynthesis work in concert to ensure promotion of perforation expansion during lace plant leaf development. The high level of expression of genes involved in these diverse biological functions differed significantly between early and late lace plant leaf developmental stages, indicating their involvement in regulating perforation initiation, execution, and leaf growth. The results of our investigation into the lace plant transcriptome and expression patterns reveal a variety of candidate genes with possible involvement in the initiation and progression of lace plant leaf cell death, generating new hypotheses and providing novel insights into plant developmental PCD. Future experiments on candidate DEGs will be required moving forward to characterize and confirm protein functionality in lace plant leaf perforation formation.

All things considered, the A. madagascariensis transcriptome data generated and analysed herein exemplifies the power of de novo Illumina-based RNA-Seq. Our transcriptomes serve as both a high-quality gene discovery resource and a framework for the detection of physiological changes through gene expression profiling. Combined with additional transcriptome annotation tools, experimental observations from model plant species will undoubtedly facilitate deeper insights into the biology of PCD and remodelling of lace plant leaves in the future.

Methods

Plant tissue culturing

Lace plant cultures were propagated aseptically as described in Gunawardena et al. (2006) [93]. Lace plant corms were cultured in Magenta GA-7 boxes, embedded in 100 mL of solid MS media made of 1.5% plant tissue culture agar (w/v, Phytotechnology Laboratories) in liquid MS [3% sucrose (w/v), 0.01% Myo-inositol (w/v), 0.215% MS basal salts (w/v, Phytotechnology Laboratories), 0.0025% thiamine-HCl (v/v), pH 5.7] and then submerged under 150 ml of liquid MS. Plant cultures were grown at 24 °C and exposed under light levels of 125 μmol m−2 s−1 on 12 h light/dark cycles with daylight deluxe fluorescent light bulbs (Philips). A. madagascariensis (Mirbel) H. Bruggen corms were originally purchased from The PlantGuy (AB, Canada).

Cultures grew for 30 days or until each plant produced 3 perforated mature leaves. One of each imperforate, pre-perforation, window stage leaves, and one of the most recently developed mature leaf were separated from one whole plant culture and washed thoroughly with distilled water before leaf tissue was excised from the midrib and flash frozen for downstream molecular work. Three independent experiments were carried out. For each experiment, one of each leaf stage was taken from one whole plant culture.

PCD and NPCD cell preparations

A Zeiss PALM Laser Capture Microdissection and Imaging System (North York, ON, Canada) was used to separate PCD and NPCD cells. The cells were collected from the areoles of window stage leaves, where each type is visibly distinguishable by colour. Cell type population samples were collected from 4 separate areoles (approximately 2000 cells per areole) per window stage leaf before being flash-frozen. Three independent experiments per cell type were carried out. For each experiment, one window stage leaf was taken from one whole plant culture. A sample diagram of the laser capture and catapult process of the leaf tissue sample is provided in Additional file 2.

RNA extraction and quality control

RNA was extracted from leaves of four developmental stages, and from the two different cell types (NPCD and PCD cells). RNA was extracted from 40 mg of flash-frozen, midrib free leaf lamina tissue from one of each imperforate, pre-perforation, window or mature stage leaves from 3 different whole-plant cultures as per instructions for the ReliaPrep RNA Kit (Promega). RNA samples were treated with DNAse I (Thermo Fisher). Eluted RNA quantity was estimated using a Nanodrop spectrophotometer (Thermo Fischer) and a RNA integrity number (RIN) was determined using a Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA). Only RNA samples with a RIN ≥ 6.5 were approved for cDNA conversion.

cDNA library preparation and Illumina sequencing

cDNA library preparation and sequencing were performed by Génome Québec (Montréal, QC, Canada). Eighteen paired-end RNA-Seq libraries of length 100 bp were generated on an Illumina NovaSeq6000 (CA, USA) using strand-specific Trueseq protocols. The raw read data obtained were deposited to NCBI and are accessible under the SRA. SRA accession IDs: SRR10524134-SR10524151 and BIOPROJECTID: PRJNA591467. Data were first inspected for quality by analysing FastQ files with FastQC [94]. Reads of low quality and containing adapter contaminations were trimmed with Trimmomatic v.0.35 [95] with a k-mer size of 25 and with parameters: ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:5:0 MINLEN:50. The quality of trimmed reads was assessed using FastQC v0.11.2 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

Transcriptome de novo assembly

High-quality adapter free reads were used to construct a de novo assembly with Trinity v2.3.1 [96] with default settings. Quality evaluation of assemblies was considered with major bioinformatics indicators such as transcript mean length, GC%, and an N50 value (Table 1). The Trinity pipeline clusters de novo assembled transcripts into genes and isoforms, and we worked only with ‘genes’ datasets, with the highest expressed isoform chosen as the representative for each gene. This Transcriptome Shotgun Assembly project has been deposited at DDBJ/ENA/GenBank under the accession NO. GJFM00000000. The version described in this paper is the first version, GJFM01000000.

Transcript quantification and identification of differentially expressed genes

The abundance of each gene was calculated by aligning reads from each sample to our de novo transcriptome with RNA-Seq by Expectation–Maximization (RSEM) [97] for each library. The trimmed mean of M-values (TMM) method [98] was used to calculate the normalization factors (one calculation for NPCD vs PCD cells and one calculation for comparisons among the leaf stages). Using Trinity [99], expression normalization was performed using TMM, following fragments per kilobase of exon model per million reads (FPKM) calculations. DEGs among the leaf stage and between the cell type libraries were identified using the Empirical Analysis of Digital Gene Expression (edgeR) statistical package [100] (http://bioconductor.org/packages/release/bioc/html/edgeR.html) performed with R (v3.3.2; R Core Team 2015). Genes that were more than twofold differentially expressed with an FDR of 1% were defined as differentially expressed [101].

Cluster analysis

Expression patterns of genes among leaf stage samples and between NPCD and PCD cells were separated using cluster analysis of DEGs. Hierarchical clustering of normalized gene expression was achieved using centralized and log2(FPKM + 1) transformation [99] and tree cutting at a depth of 40%, with heatmap visualization performed using R and the package Superheat [102].

Annotation and GO enrichment analysis

To identify the genes and functions associated with each transcript, assembled transcripts were annotated using Trinotate [103] and public genome and functional repositories. We annotated transcripts based on top matching sequence similarity to orthologs in public databases, including GO, the Kyoto Encyclopedia of Genes and Genomes (KEGG), the euKaryotic Ortholog Groups database (KOG), the Swiss Protein Resource (Swiss-Prot), and the Panther Database, using the BLASTX method with a cut-off E-value of 10−5 [104,105,106]. To eliminate transcripts derived from foreign organisms and lab contaminants, genes of non-plant origin were removed. Selected annotation data was included with Additional file 3 and 5.

GO-enrichment analysis was carried out using the Plant Transcription Factor Database v5.0 [107] program based on Fisher’s Exact Test with multiple testing correction of FDR = 1%. GO analysis was performed by comparing the GO terms in the test sample to the GO terms in the background reference of the entire lace plant de novo transcriptome generated from all samples.

Validation with qRT-PCR

Five selected DEGs (Bag5, expansin A-29, aquaporin 4–4, anthocyanin regulatory protein C1, nuclear transcription factor YC-1, and α-tubulin) were used to verify the expression results of RNA-Seq by using the ΔΔCT method. RNA from lace plant leaf stages and isolated cell types were collected and extracted as described above. Single-strand cDNA was synthesized using SuperScript®III First-Strand Synthesis System for qRT-PCR (Invitrogen, Burlington, ON, Canada) and oligo dT20 following the manufacturer’s instructions. qRT-PCR was conducted on a Rotor-Gene RG-3000 system (Corbett Research, Sydney, NSW, Australia) using 0.5 μl cDNA as a template and 0.4 mmol l−1 primers for all selected genes (Additional file 6) under the following conditions: 5 min at 94 °C, 35 cycles of 30 s at 94 °C, 30 s at 54 °C for all chosen genes and 1 min at 72 °C, followed by 72 °C. qPCR was conducted using a QuantiFast® SYBER® Green PCR Kit (Qiagen, Mississauga, ON, Canada). Melt curve analysis was conducted by Rotor-Gene 6 Software and experiments with at least 90% efficiency were used for analysis (Corbett Research). The experiment was performed in triplicate using three biological replicates of imperforate, pre-perforation, window, and mature stage leaves as well as NPCD cells and PCD cells. cDNA copy numbers for five chosen genes (Additional file 4) were determined from a standard curve of Ct values (R2 > 0.99) and normalized against the α-tubulin isoform [37].

Image analysis and processing

Images of leaf layouts were obtained using a Nikon L110 digital camera. Photoshop and Illustrator (Adobe Creative Cloud; Adobe Systems Inc.) were used to prepare images and remove backgrounds of detached leaves for publication. A Nikon AZ100 microscope acquired micrographs of leaf stages. Image adjustments were made evenly within and consistently across figures and included background removal, as well as fine tuning of brightness, contrast and colour balance using Photoshop.

Statistical analysis and data representation

One-way ANOVA followed by a Tukey test was used to identify significant differences among leaf stage means for qRT-PCR validation experiments, and a Student’s t-test was used to identify significant differences between means of cell types (Additional file 4). All data are illustrated as the mean ± standard error. Analyses were conducted using GraphPad Prism 5 software (GraphPad Software Inc.).

Availability of data and materials

Raw sequencing data files are available in the NCBI SRA database with project accession NO. PRJNA591467. This Transcriptome Shotgun Assembly project has been deposited at DDBJ/ENA/GenBank under the accession NO. GJFM00000000. The version described in this paper is the first version, GJFM01000000.

Abbreviations

DEG:

Differentially expressed gene

BP:

Base pair

GO:

Gene Ontology

nt:

Nucleotide

FPKM:

Fragments per kilobase of exon model per million mapped reads

RNA-Seq:

RNA sequencing

PCD:

Programmed cell death

NPCD:

Non-PCD

ROS:

Reactive oxygen species

qRT-PCR:

Quantitative reverse transcriptase polymerase chain reaction

References

  1. Kabbage M, Kessens R, Bartholomay LC, Williams B. The Life and Death of a Plant Cell Programmed cell death (PCD): a genetically regulated program of cellular suicide employed by both unicellular and multicellular organisms. Annu Rev Plant Biol. 2017;68:375–404. https://doi.org/10.1146/annurev-arplant-043015.

    Article  CAS  PubMed  Google Scholar 

  2. Kacprzyk J, Daly CT, McCabe PF. The botanical dance of death. Programmed cell death in plants. Adv Botanical Res. 2011;60:169–261. https://doi.org/10.1016/B978-0-12-385851-1.00004-4.

    Article  CAS  Google Scholar 

  3. Kacprzyk J, Dauphinee AN, Gallois P, Gunawardena AHLAN, McCabe PF. Methods to Study Plant Programmed Cell Death. Methods Mol Biol. 2016;1419:145–60. https://doi.org/10.1007/978-1-4939-3581-9_12 (New York: Humana Press).

    Article  PubMed  Google Scholar 

  4. Van Hautegem T, Waters AJ, Goodrich J, Nowack MK. Only in dying, life: Programmed cell death during plant development. Trends Plant Sci. 2015;20(2):102–13. https://doi.org/10.1016/j.tplants.2014.10.003.

    Article  CAS  PubMed  Google Scholar 

  5. Daneva A, Gao Z, Van DM, Nowack MK, Daneva A. Functions and Regulation of Programmed Cell Death in Plant Development. Annu Rev Cell Dev Biol. 2016;32:1–7. https://doi.org/10.1146/annurev-cellbio-111315-124915.

    Article  CAS  Google Scholar 

  6. Gunawardena AHLAN, McCabe PF. Plant Programmed Cell Death. Cham: Springer; 2015. https://doi.org/10.1007/978-3-319-21033-9.

    Book  Google Scholar 

  7. Kuriyama H, Fukuda H. Developmental programmed cell death in plants. Curr Opin Plant Biol. 2002;5(6):568–73. https://doi.org/10.1016/S1369-5266(02)00305-9.

    Article  CAS  PubMed  Google Scholar 

  8. Bozhkov PV, Filonova LH, Suarez MF, Helmersson A, Smertenko AP, Zhivotovsky B, et al. VEIDase is a principal caspase-like activity involved in plant programmed cell death and essential for embryonic pattern formation. Cell Death Differ. 2004;11(2):175–82. https://doi.org/10.1038/sj.cdd.4401330.

    Article  CAS  PubMed  Google Scholar 

  9. Gunawardena AHLAN, Pearce DM, Jackson MB, Hawes CR, Evans DE. Characterisation of programmed cell death during aerenchyma formation induced by ethylene or hypoxia in roots of maize (Zea mays L.). Planta. 2001;212(2):205–14. https://doi.org/10.1007/s004250000381.

    Article  CAS  PubMed  Google Scholar 

  10. Han JJ, Lin W, Oda Y, Cui KM, Fukuda H, He XQ. The proteasome is responsible for caspase-3-like activity during xylem development. Plant J. 2012;72(1):129–41. https://doi.org/10.1111/j.1365-313X.2012.05070.x.

    Article  CAS  PubMed  Google Scholar 

  11. Ohashi-Ito K, Oda Y, Fukuda H. Arabidopsis VASCULAR-RELATED NAC-DOMAIN6 Directly Regulates the Genes That Govern Programmed Cell Death and Secondary Wall Formation during Xylem Differentiation. Plant Cell. 2010;22(10):3461–73. https://doi.org/10.1105/tpc.110.075036.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Gunawardena AHLAN. Programmed cell death and tissue remodelling in plants. J Exp Bot. 2008;59:445–51. https://doi.org/10.1093/jxb/erm189.

    Article  CAS  PubMed  Google Scholar 

  13. Dauphinee AN, Fletcher JI, Denbigh GL, Lacroix CR, Gunawardena AHLAN. Remodelling of lace plant leaves: antioxidants and ROS are key regulators of programmed cell death. Planta. 2017;246(1):133–47. https://doi.org/10.1007/s00425-017-2683-y.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Dauphinee AN, Denbigh GL, Rollini A, Fraser M, Lacroix CR, Gunawardena AHLAN. The function of autophagy in lace plant programmed cell death. Front Plant Sci. 2019;10:1198. https://doi.org/10.3389/fpls.2019.01198.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Dauphinee AN, Warner TS, Gunawardena AHLAN. A comparison of induced and developmental cell death morphologies in lace plant (Aponogeton madagascariensis) leaves. BMC Plant Biol. 2014;14(1):1–13. https://doi.org/10.1186/s12870-014-0389-x.

    Article  Google Scholar 

  16. Lord CEN, Dauphinee AN, Watts RL, Gunawardena AHLAN. Unveiling interactions among Mitochondria, caspase-Like proteases, and the actin cytoskeleton during plant programmed cell death (PCD). PLoS One. 2013;8:3. https://doi.org/10.1371/journal.pone.0057110.

    Article  CAS  Google Scholar 

  17. Rantong G, Gunawardena AHLAN. Programmed cell death: genes involved in signaling, regulation, and execution in plants and animals. Botany. 2015;93(4):193–210. https://doi.org/10.1139/cjb-2014-0152.

    Article  CAS  Google Scholar 

  18. Rantong G, Van Der Kelen K, Van Breusegem F, Gunawardena AHLAN. Identification of Differentially Expressed Genes during Lace Plant Leaf Development. Int J Plant Sci. 2016;177(5):419–31. https://doi.org/10.1086/684748.

    Article  Google Scholar 

  19. Reza SH, Delhomme N, Street NR, Ramachandran P, Dalman K, Nilsson O, et al. Transcriptome analysis of embryonic domains in Norway spruce reveals potential regulators of suspensor cell death. PLoS ONE. 2018;13(3):e0192945. https://doi.org/10.1371/journal.pone.0192945.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Du XM, Ni XL, Ren XL, Xin GL, Jia GL, Liu HD, et al. De novo transcriptomic analysis to identify differentially expressed genes during the process of aerenchyma formation in Typha angustifolia leaves. Gene. 2018;662:66–75. https://doi.org/10.1016/j.gene.2018.03.099.

    Article  CAS  PubMed  Google Scholar 

  21. Li P, Ponnala L, Gandotra N, Wang L, Si Y, Tausta SL, et al. The developmental dynamics of the maize leaf transcriptome. Nat Genet. 2010;42(12):1060. https://doi.org/10.1038/ng.703.

    Article  CAS  PubMed  Google Scholar 

  22. Gross SM, Martin JA, Simpson J, Abraham-Juarez MJ, Wang Z, Visel A. De novo transcriptome assembly of drought tolerant CAM plants, Agave deserti and Agave tequilana. BMC Genomics. 2013;14(1):563. https://doi.org/10.1186/1471-2164-14-563.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Dai M, Hu Y, Zhao Y, Liu H, Zhou DX. A WUSCHEL-LIKE HOMEOBOX gene represses a YABBY gene expression required for rice leaf development. Plant Physiol. 2007;144(1):380–90. https://doi.org/10.1104/pp.107.095737.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Guo Z, Fujioka S, Blancaflor EB, Miao S, Gou X, Li J. TCP1 modulates brassinosteroid biosynthesis by regulating the expression of the key biosynthetic gene DWARF4 in arabidopsis thaliana. Plant Cell. 2010;22(4):1161–73. https://doi.org/10.1105/tpc.109.069203.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Husbands AY, Chitwood DH, Plavskin Y, Timmermans MCP. Signals and prepatterns: New insights into organ polarity in plants. Genes Dev. 2009;23(17):1986–97. https://doi.org/10.1101/gad.1819909.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Kaplan-Levy RN, Brewer PB, Quon T, Smyth DR. The trihelix family of transcription factors - light, stress and development. Trends Plant Sci. 2012;17(3):163–71. https://doi.org/10.1016/j.tplants.2011.12.002.

    Article  CAS  PubMed  Google Scholar 

  27. Kim HS, Kim SJ, Abbasi N, Bressan RA, Yun DJ, Yoo SD, et al. The DOF transcription factor Dof5.1 influences leaf axial patterning by promoting revoluta transcription in arabidopsis. Plant J. 2010;64(3):524–35. https://doi.org/10.1111/j.1365-313X.2010.04346.x.

    Article  CAS  PubMed  Google Scholar 

  28. Liu T, Ohashi-Ito K, Bergmann DC. Orthologs of Arabidopsis thaliana stomatal bHLH genes and regulation of stomatal development in grasses. Development. 2009;136(13):2265–76. https://doi.org/10.1242/dev.032938.

    Article  CAS  PubMed  Google Scholar 

  29. Reyes JC, Muro-Pastor MI, Florencio FJ. The GATA family of transcription factors in arabidopsis and rice. Plant Physiol. 2004;134(4):1718–32. https://doi.org/10.1104/pp.103.037788.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Sawa S, Watanabe K, Goto K, Kanaya E, Morita EH, Okada K. Filamentous flower, a meristem and organ identity gene of Arabidopsis, encodes a protein with a zinc finger and HMG-related domains. Genes Dev. 1999;13(9):1079–88. https://doi.org/10.1101/gad.13.9.1079.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Song YH, Lee I, Lee SY, Imaizumi T, Hong JC. CONSTANS and ASYMMETRIC LEAVES 1 complex is involved in the induction of FLOWERING LOCUS T in photoperiodic flowering in Arabidopsis. Plant J. 2012;69(2):332–42. https://doi.org/10.1111/j.1365-313X.2011.04793.x.

    Article  CAS  PubMed  Google Scholar 

  32. Gadjev I. Transcriptomic Footprints Disclose Specificity of Reactive Oxygen Species Signaling in Arabidopsis. PLANT Physiol. 2006;141(2):436–45. https://doi.org/10.1104/pp.106.078717.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Gechev TS, Van Breusegem F, Stone JM, Denev I, Laloi C. Reactive oxygen species as signals that modulate plant stress responses and programmed cell death. BioEssays. 2006;28(11):1091–101. https://doi.org/10.1002/bies.20493.

    Article  CAS  PubMed  Google Scholar 

  34. Petrov V, Hille J, Mueller-Roeber B, Gechev TS. ROS-mediated abiotic stress-induced programmed cell death in plants. Front Plant Sci. 2015;6:1. https://doi.org/10.3389/fpls.2015.00069.

    Article  Google Scholar 

  35. Dauphinee AN, Wright H, Rantong G, Gunawardena AHLAN. The involvement of ethylene in programmed cell death and climacteric-like behaviour during the remodelling of lace plant (Aponogeton madagascariensis) leaves. Botany. 2012;90:1237–44. https://doi.org/10.1139/b2012-093.

    Article  CAS  Google Scholar 

  36. Denbigh GL, Dauphinee AN, Fraser MS, Lacroix CR, Gunawardena AHLAN. The role of auxin in developmentally regulated programmed cell death in lace plant. Am J Bot. 2020;107(4):577–86. https://doi.org/10.1002/ajb2.1463.

    Article  CAS  PubMed  Google Scholar 

  37. Rantong G, Evans R, Gunawardena AHLAN. Lace plant ethylene receptors, AmERS1a and AmERS1c, regulate ethylene-induced programmed cell death during leaf morphogenesis. Plant Mol Biol. 2015;89(3):215–27. https://doi.org/10.1007/s11103-015-0356-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Yan MY, Xie DL, Cao JJ, Xia XJ, Shi K, Zhou YH, et al. Brassinosteroid-mediated reactive oxygen species are essential for tapetum degradation and pollen fertility in tomato. Plant J. 2020;102(5):931–47. https://doi.org/10.1111/tpj.14672.

    Article  CAS  PubMed  Google Scholar 

  39. Limami AM, Diab H, Lothier J. Nitrogen metabolism in plants under low oxygen stress. Planta. 2014;239(3):531–41. https://doi.org/10.1007/s00425-013-2015-9.

    Article  CAS  PubMed  Google Scholar 

  40. Bu C, Zhang Q, Zeng J, Cao X, Hao Z, Qiao D, et al. Identification of a novel anthocyanin synthesis pathway in the fungus Aspergillus sydowii H-1. BMC Genomics. 2020;21(1):29. https://doi.org/10.1186/s12864-019-6442-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Argout X, Fouet O, Wincker P, Gramacho K, Legavre T, Sabau X, et al. Towards the understanding of the cocoa transcriptome: Production and analysis of an exhaustive dataset of ESTs of Theobroma cacao L. generated from various tissues and under various conditions. BMC Genomics. 2008;9(1):512. https://doi.org/10.1186/1471-2164-9-512.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Rowarth NM, Dauphinee AN, Denbigh GL, Gunawardena AHLAN. Hsp70 plays a role in programmed cell death during the remodelling of leaves of the lace plant (Aponogeton madagascariensis). J Exp Bot. 2020;71(3):907–18. https://doi.org/10.1093/jxb/erz447.

    Article  CAS  PubMed  Google Scholar 

  43. Kallam K, Appelhagen I, Luo J, Albert N, Zhang H, Deroles S, et al. Aromatic Decoration Determines the Formation of Anthocyanic Vacuolar Inclusions. Curr Biol. 2017;27(7):945–57. https://doi.org/10.1016/j.cub.2017.02.027.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Gunawardena AHLAN, Greenwood JS, Dengler NG. Cell wall degradation and modification during programmed cell death in lace plant, Aponogeton madagascariensis (Aponogetonaceae). Am J Bot. 2007;94:1116–28. https://doi.org/10.3732/ajb.94.7.1116.

    Article  CAS  PubMed  Google Scholar 

  45. Franciosini A, Rymen B, Shibata M, Favero DS, Sugimoto K. Molecular networks orchestrating plant cell growth. Curr Opin Plant Biol. 2017;35:98–104. https://doi.org/10.1016/j.pbi.2016.11.010.

    Article  CAS  PubMed  Google Scholar 

  46. Landi S, Hausman JF, Guerriero G, Esposito S. Poaceae vs. Abiotic stress: Focus on drought and salt stress, recent insights and perspectives. Front Plant Sci. 2017;8:1214. https://doi.org/10.3389/fpls.2017.01214.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Tenhaken R. Cell wall remodeling under abiotic stress. Front Plant Sci. 2015;5:771. https://doi.org/10.3389/fpls.2014.00771.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Tyerman SD, Niemietz CM, Bramley H. Plant aquaporins: Multifunctional water and solute channels with expanding roles. Plant, Cell Environ. 2002;25(2):173–94. https://doi.org/10.1046/j.0016-8025.2001.00791.x.

    Article  CAS  Google Scholar 

  49. Kaimal JM, Kandasamy G, Gasser F, Andréasson C. Coordinated Hsp110 and Hsp104 Activities Power Protein Disaggregation in Saccharomyces cerevisiae. Mol Cell Biol. 2017;37(11):e00027-e117. https://doi.org/10.1128/MCB.00027-17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Shin Y, Klucken J, Patterson C, Hyman BT, McLean PJ. The Co-chaperone carboxyl terminus of Hsp70-interacting protein (CHIP) mediates α-synuclein degradation decisions between proteasomal and lysosomal pathways. J Biol Chem. 2005;280(25):23727–34. https://doi.org/10.1074/jbc.M503326200.

    Article  CAS  PubMed  Google Scholar 

  51. Qi Y, Wang H, Zou Y, Liu C, Liu Y, Wang Y, Zhang W. Over-expression of mitochondrial heat shock protein 70 suppresses programmed cell death in rice. FEBS Lett. 2011;585(1):231–9. https://doi.org/10.1016/j.febslet.2010.11.051.

    Article  CAS  PubMed  Google Scholar 

  52. Kim NH, Hwang BK. Pepper Heat Shock Protein 70a Interacts with the Type III Effector AvrBsT and Triggers Plant Cell Death and Immunity. Plant Physiol. 2015;167(2):307–22. https://doi.org/10.1104/pp.114.253898.

    Article  CAS  PubMed  Google Scholar 

  53. Dametto A, Buffon G, Blasi ÉA dos R, Sperotto RA. Ubiquitination pathway as a target to develop abiotic stress tolerance in rice. Plant Signal Behav. 2015;10(9):e1057369. https://doi.org/10.1080/15592324.2015.1057369.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Yabe N, Takahashi T, Komeda Y. Analysis of tissue-specific expression of arabidopsis thaliana HSP90-family gene HSP81. Plant Cell Physiol. 1994;35(8):1207–19. https://doi.org/10.1093/oxfordjournals.pcp.a078715.

    Article  CAS  PubMed  Google Scholar 

  55. Chen J, Gao T, Wan S, Zhang Y, Yang J, Yu Y, Wang W. Genome-wide identification, classification and expression analysis of the HSP gene superfamily in tea plant (Camellia sinensis). Int J Mol Sci. 2018;19(9):2633. https://doi.org/10.3390/ijms19092633.

    Article  CAS  PubMed Central  Google Scholar 

  56. Behl C. Breaking BAG: The Co-Chaperone BAG3 in Health and Disease. Trends Pharmacol Sci. 2016;37(8):672–88. https://doi.org/10.1016/j.tips.2016.04.007.

    Article  CAS  PubMed  Google Scholar 

  57. Kabbage M, Dickman MB. The BAG proteins: A ubiquitous family of chaperone regulators. Cell Mol Life Sci. 2008;65(9):1390–402. https://doi.org/10.1007/s00018-008-7535-2.

    Article  CAS  PubMed  Google Scholar 

  58. Li L, Xing Y, Chang D, Fang S, Cui B, Li Q, et al. CaM/BAG5/Hsc70 signaling complex dynamically regulates leaf senescence. Sci Rep. 2016;6(1):1–2. https://doi.org/10.1038/srep31889.

    Article  CAS  Google Scholar 

  59. Li Y, Williams B, Dickman M. Arabidopsis B-cell lymphoma2 (Bcl-2)-associated athanogene 7 (BAG7)-mediated heat tolerance requires translocation, sumoylation and binding to WRKY29. New Phytol. 2017;214(2):695–705. https://doi.org/10.1111/nph.14388.

    Article  CAS  PubMed  Google Scholar 

  60. Williams B, Kabbage M, Britt R, Dickman MB. AtBAG7, an Arabidopsis Bcl-2-associated athanogene, resides in the endoplasmic reticulum and is involved in the unfolded protein response. Proc Natl Acad Sci. 2010;107(13):6088–93. https://doi.org/10.1073/pnas.0912670107.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Thanthrige N, Jain S, Bhowmik SD, Ferguson BJ, Kabbage M, Mundree S, Williams B. Centrality of BAGs in Plant PCD, Stress Responses, and Host Defense. Trends Plant Sci. 2020. https://doi.org/10.1016/j.tplants.2020.04.012.

    Article  PubMed  Google Scholar 

  62. Lord CEN, Wertman JN, Lane S, Gunawardena AHLAN. Do mitochondria play a role in remodelling lace plant leaves during programmed cell death? BMC Plant Biol. 2011;11(1):102. https://doi.org/10.1186/1471-2229-11-102.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Xiong Y, Contento AL, Bassham DC. AtATG18a is required for the formation of autophagosomes during nutrient stress and senescence in Arabidopsis thaliana. Plant J. 2005;42(4):535–46. https://doi.org/10.1111/j.1365-313X.2005.02397.x.

    Article  CAS  PubMed  Google Scholar 

  64. Lai Z, Wang F, Zheng Z, Fan B, Chen Z. A critical role of autophagy in plant resistance to necrotrophic fungal pathogens. Plant J. 2011;66(6):953–68. https://doi.org/10.1111/j.1365-313X.2011.04553.x.

    Article  CAS  PubMed  Google Scholar 

  65. Hanaoka H, Noda T, Shirano Y, Kato T, Hayashi H, Shibata D, et al. Leaf senescence and starvation-induced chlorosis are accelerated by the disruption of an Arabidopsis autophagy gene. Plant Physiol. 2002;129(3):1181–93. https://doi.org/10.1104/pp.011024.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Liu Y, Xiong Y, Bassham DC. Autophagy is required for tolerance of drought and salt stress in plants. Autophagy. 2009;5(7):954–63. https://doi.org/10.4161/auto.5.7.9290.

    Article  CAS  PubMed  Google Scholar 

  67. Xiong Y, Contento AL, Bassham DC. Disruption ol autophagy results in constitutive oxidative stress in Arabidopsis. Autophagy. 2007;3(3):257–8. https://doi.org/10.4161/auto.3847.

    Article  CAS  PubMed  Google Scholar 

  68. Yang Z, Klionsky DJ. An Overview of the Molecular Mechanism of Autophagy. In: Autophagy in infection and immunity. Berlin: Springer; 2009. p. 1–32. https://doi.org/10.1007/978-3-642-00302-8_1.

  69. Baena-González E, Rolland F, Thevelein JM, Sheen J. A central integrator of transcription networks in plant stress and energy signalling. Nature. 2007;448(7156):938–42. https://doi.org/10.1038/nature06069.

    Article  CAS  PubMed  Google Scholar 

  70. Liu Y, Bassham DC. Autophagy: Pathways for Self-Eating in Plant Cells. Annu Rev Plant Biol. 2012;63:215–37. https://doi.org/10.1146/annurev-arplant-042811-105441.

    Article  CAS  PubMed  Google Scholar 

  71. Gechev TS, Dinakar C, Benina M, Toneva V, Bartels D. Molecular mechanisms of desiccation tolerance in resurrection plants. Cell Mol Life Sci. 2012;69(19):3175–86. https://doi.org/10.1007/s00018-012-1088-0.

    Article  CAS  PubMed  Google Scholar 

  72. Huang P, Ju HW, Min JH, Zhang X, Kim SH, Yang KY, et al. Overexpression of L-type lectin-like protein kinase 1 confers pathogen resistance and regulates salinity response in Arabidopsis thaliana. Plant Sci. 2013;203:98–106. https://doi.org/10.1016/j.plantsci.2012.12.019.

    Article  CAS  PubMed  Google Scholar 

  73. Maksaev G, Shoots JM, Ohri S, Haswell ES. Nonpolar residues in the presumptive pore-lining helix of mechanosensitive channel MSL10 influence channel behavior and establish a nonconducting function. Plant Direct. 2018;2(6):e00059. https://doi.org/10.1002/pld3.59.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Wang Y, Bouwmeester K. L-type lectin receptor kinases: New forces in plant immunity. PLoS Pathog. 2017;13(8):e1006433. https://doi.org/10.1371/journal.ppat.1006433.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Luo H, Laluk K, Lai Z, Veronese P, Song F, Mengiste T. The arabidopsis botrytis susceptible1 interactor defines a subclass of RING E3 ligases that regulate pathogen and stress responses. Plant Physiol. 2010;154(4):1766–82. https://doi.org/10.1104/pp.110.163915.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Yang YL, Li XM. The IAP family: Endogenous caspase inhibitors with multiple biological activities. Cell Res. 2000;10(3):169–77. https://doi.org/10.1038/sj.cr.7290046.

    Article  CAS  PubMed  Google Scholar 

  77. Millar AA, Gubler F. The Arabidopsis GAMYB-like Genes, MYB33 and MYB65, are microRNA-regulated genes that redundantly facilitate anther development. Plant Cell. 2005;17(3):705–21. https://doi.org/10.1105/tpc.104.027920.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Kwak JM, Mori IC, Pei ZM, Leonhardt N, Torres MA, Dangl JL, et al. NADPH oxidase AtrbohD and AtrbohF genes function in ROS-dependent ABA signaling in Arabidopsis. EMBO J. 2003;22(11):2623–33. https://doi.org/10.1093/emboj/cdg277.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Tavladoraki P, Cona A, Angelini R. Copper-containing amine oxidases and FAD-dependent polyamine oxidases are key players in plant tissue differentiation and organ development. Front Plant Sci. 2016;7:824. https://doi.org/10.3389/fpls.2016.00824.

    Article  PubMed  PubMed Central  Google Scholar 

  80. Galsurker O, Doron-Faigenboim A, Teper-Bamnolker P, Daus A, Fridman Y, Lers A, Eshel D. Cellular and molecular changes associated with onion skin formation suggest involvement of programmed cell death. Front Plant Sci. 2017;7:2031. https://doi.org/10.3389/fpls.2016.02031.

  81. Sueldo DJ, van der Hoorn RAL. Plant life needs cell death, but does plant cell death need Cys proteases? FEBS J. 2017;284(10):1577–85. https://doi.org/10.1111/febs.14034.

    Article  CAS  PubMed  Google Scholar 

  82. Zhang D, Liu D, Lv X, Wang Y, Xun Z, Liu Z, et al. The cysteine protease CEP1, a key executor involved in tapetal programmed cell death, regulates pollen development in Arabidopsis. Plant Cell. 2014;26(7):2939–61. https://doi.org/10.1105/tpc.114.127282.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Rustgi S, Boex-Fontvieille E, Reinbothe C, Von Wettstein D, Reinbothe S. Serpin1 and WSCP differentially regulate the activity of the cysteine protease RD21 during plant development in Arabidopsis thaliana. Proc Natl Acad Sci U S A. 2017;114(9):2212–7. https://doi.org/10.1073/pnas.1621496114.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Shindo T, Misas-Villamil JC, Hörger AC, Song J, van der Hoorn RAL. A role in immunity for arabidopsis cysteine protease RD21, the ortholog of the tomato immune protease C14. PLoS ONE. 2012;7(1):e29317. https://doi.org/10.1371/journal.pone.0029317.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Balakireva AV, Zamyatnin AA. Cutting out the gaps between proteases and programmed cell death. Front Plant Sci. 2019;10:704. https://doi.org/10.3389/fpls.2019.00704.

    Article  PubMed  PubMed Central  Google Scholar 

  86. Gao H, Li R, Guo Y. Arabidopsis aspartic proteases A36 and A39 play roles in plant reproduction. Plant Signal Behav. 2017;12(4):e1304343. https://doi.org/10.1080/15592324.2017.1304343.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Gao H, Zhang Y, Wang W, Zhao K, Liu C, Bai L, et al. Two membrane-anchored aspartic proteases contribute to pollen and ovule development. Plant Physiol. 2017;173(1):219–39. https://doi.org/10.1104/pp.16.01719.

    Article  CAS  PubMed  Google Scholar 

  88. Li D, Zhang H, Song Q, Wang L, Liu S, Hong Y, et al. Tomato Sl3-MMP, a member of the Matrix metalloproteinase family, is required for disease resistance against Botrytis cinerea and Pseudomonas syringae pv. tomato DC3000. BMC Plant Biol. 2015;15(1):1–8. https://doi.org/10.1186/s12870-015-0536-z.

    Article  CAS  Google Scholar 

  89. Zimmermann D, Gomez-Barrera JA, Pasule C, Brack-Frick UB, Sieferer E, Nicholson TM, et al. Cell death control by matrix metalloproteinases. Plant Physiol. 2016;171(2):1456–69. https://doi.org/10.1104/pp.16.00513.

    Article  PubMed  PubMed Central  Google Scholar 

  90. Bollhöner B, Zhang B, Stael S, Denancé N, Overmyer K, Goffner D, et al. Post mortem function of AtMC9 in xylem vessel elements. New Phytol. 2013;200(2):498–510. https://doi.org/10.1111/nph.12387.

    Article  CAS  PubMed  Google Scholar 

  91. Tsiatsiani L, Timmerman E, De Bock PJ, Vercammen D, Stael S, Van De Cotte B, et al. The Arabidopsis metacaspase9 degradome. Plant Cell. 2013;25(8):2831–47. https://doi.org/10.1105/tpc.113.115287.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  92. Rantong G, Gunawardena AHLAN. Vacuolar processing enzymes, AmVPE1 and AmVPE2, as potential executors of ethylene regulated programmed cell death in the lace plant (Aponogeton madagascariensis). Botany. 2018;96(4):235–47. https://doi.org/10.1139/cjb-2017-0184.

    Article  CAS  Google Scholar 

  93. Gunawardena AHLAN, Navachandrabala C, Kane M, Dengler NG, Teixeira da Silva JA. Lace plant: a novel system for studying developmental programmed cell death. Floriculture, ornamental and plant biotechnology: advances and tropical issues volume 1. 2006, Global Science Books, Ltd, Middlesex, 157–62.

  94. Andrews S. FastQC: A quality control tool for high throughput sequence data. Available online. http://www.bioinformatics.babraham.ac.uk/projects/fastqc

  95. Bolger AM, Lohse M, Usadel B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doi.org/10.1093/bioinformatics/btu170.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  96. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52. https://doi.org/10.1038/nbt.1883.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  97. Li B, Dewey CN. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323. https://doi.org/10.1186/1471-2105-12-323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):1–9. https://doi.org/10.1186/gb-2010-11-3-r25.

    Article  CAS  Google Scholar 

  99. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512. https://doi.org/10.1038/nprot.2013.084.

    Article  CAS  PubMed  Google Scholar 

  100. Robinson M, McCarthy DJ, Smyth GK. edgeR: differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40. https://doi.org/10.1093/bioinformatics/btp616.

    Article  CAS  PubMed  Google Scholar 

  101. Benjamini Y, Gavrilov Y. A simple forward selection procedure based on false discovery rate control. Ann Appl Stat. 2009;3(1):179–98. https://doi.org/10.1214/08-AOAS194.

    Article  Google Scholar 

  102. Barter RL, Yu B. Superheat: An R Package for Creating Beautiful and Extendable Heatmaps for Visualizing Complex Data. J Comput Graph Stat. 2018;27(4):910–22. https://doi.org/10.1080/10618600.2018.1473780.

    Article  PubMed  PubMed Central  Google Scholar 

  103. Bryant DM, Johnson K, DiTommaso T, Tickle T, Couger MB, Payzin-Dogru D, et al. A Tissue-Mapped Axolotl De Novo Transcriptome Enables Identification of Limb Regeneration Factors. Cell Rep. 2017;18(3):762–76. https://doi.org/10.1016/j.celrep.2016.12.063.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  104. The UniProt Consortium. UniProt: The universal protein knowledgebase. Nucleic Acids Res. 2018;46(5):2699. https://doi.org/10.1093/nar/gky092.

    Article  CAS  PubMed  Google Scholar 

  105. Kanehisa M. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000;28(1):27–30. https://doi.org/10.1093/nar/28.1.27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  106. Tatusov RL, Koonin EV, Lipman DJ. A genomic perspective on protein families. Science. 1997;278(5338):631–7. https://doi.org/10.1126/science.278.5338.631.

    Article  CAS  PubMed  Google Scholar 

  107. Jin J, Tian F, Yang DC, Meng YQ, Kong L, Luo J, Gao G. PlantTFDB 4.0: Toward a central hub for transcription factors and regulatory interactions in plants. Nucleic Acids Res. 2017;45:D1040. https://doi.org/10.1093/nar/gkw982.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We would like to thank Michaela Kember (Dalhousie University) for critically reviewing the manuscript and helping to provide lace plant leaves for RNA extraction.

Funding

This research was funded by a Natural Science and Engineering Research Council of Canada (NSERC) Discovery Grant (# 2017–04299) and Accelerator Supplements (# 2017–507825) awarded to AHLANG. NMR was supported by a NSERC Post-graduate scholarship, Nova Scotia Graduate Scholarship, and AHLANG’s Discovery Grant. There is no role of the funding body in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations

  1. Department of Biology, Dalhousie University, Halifax, NS, Canada

    Nathan M. Rowarth, Anthony L. Einfeldt & Arunika H. L. A. N. Gunawardena

  2. Department of Biochemistry & Molecular Biology, Dalhousie University, Halifax, NS, Canada

    Bruce A. Curtis & John M. Archibald

  3. Department of Biology, University of Prince Edward Island, Charlottetown, PEI, Canada

    Christian R. Lacroix

Authors
  1. Nathan M. Rowarth
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Bruce A. Curtis
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Anthony L. Einfeldt
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. John M. Archibald
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Christian R. Lacroix
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Arunika H. L. A. N. Gunawardena
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceived the study, participated in its design and coordination, secured funding and supervised all experimental work: AHLANG. Carried out all experiments, participated in its design, analysed data, drafted the first manuscript: NMR. Help with bioinformatics: BAC. Participated in study design and assisted with bioinformatics: ALE. Contributed to the manuscript revisions: NMR, ALE, BAC, JMA, CRL and AHLANG. All authors reviewed and approved the final manuscript.

Corresponding author

Correspondence to Arunika H. L. A. N. Gunawardena.

Ethics declarations

Ethics approval and consent to participate

Not Applicable.

Consent for publication

Not Applicable.

Competing interests

The authors do not declare any competing interests with the work.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Dataset S1.

Mapping rates of Trinity assembly.

Additional file 2: Figure S1.

Zeiss laser capture microdissection experiment. (A) Cell population separated by laser and (B) areole tissue remaining after cell population extracted.

Additional file 3: Dataset S2.

Leaf stage cluster annotation, median expression values, GO term counts and GO enrichment results.

Additional file 4: Figure S2.

qRT-PCR validation of leaf stages and cell types experiments. All copy numbers of probed genes were normalized by copy numbers of α-tubulin.

Additional file 5: Dataset S3

. NPCD and PCD cell cluster annotation, median expression values, GO term counts, and GO enrichment results.

Additional file 6: Dataset S6.

qRT-PCR primer information.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rowarth, N.M., Curtis, B.A., Einfeldt, A.L. et al. RNA-Seq analysis reveals potential regulators of programmed cell death and leaf remodelling in lace plant (Aponogeton madagascariensis). BMC Plant Biol 21, 375 (2021). https://doi.org/10.1186/s12870-021-03066-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-021-03066-7

Keywords

BMC Plant Biology

ISSN: 1471-2229

Contact us