Wild-type tobacco plants were grown from seed in a growth room at 25 °C with 16-h light periods, as previously described [17]. A. tumefaciens and Escherichia coli strains were grown in a Luria–Bertani medium at 30 and 37 °C, respectively. Antibiotics were used at the following concentrations: ampicillin (100 μg/mL), kanamycin (50 μg/mL), and rifampicin (100 μg/mL). VIGS of selected PCD-related genes were done in 4-week-old tobacco (N. benthamiana) seedlings as described by Velásquez et al. [13]. Primers used in constructing pTRV2 vectors of different gene are shown in Table S1. A. tumefaciens harboring the gateway compatible pTRV2 [16] and pTRV2::PDS [18] were used as negative and positive controls, respectively. Silencing of phytoene desaturase (or PDS) causes the plants to photo-bleach, hence, it can be used as a reporter gene to detect silencing efficiency.

Leaf discs were obtained from 7-week-old plants of WT or VIGS lines using a 10-mm diameter cork borer, then, submerged in OA (Sigma–Aldrich) at 20 mM for 0, 12, 24 and 48 h for the initial experiment to detect the time of PCD occurrence in the WT, and for 0 and 24 h for the VIGS lines experiment. Cell death was determined in each experiment and quantified by using Evans blue as previously described [19]. Measurements were expressed as relative values, with “one” corresponding to a maximum of a sample tested. DNA laddering was detected across time for the WT experiment, where total plant DNAs were extracted from leaf discs across time using a modified method [20]. DNA samples were separated on a 1.5% agarose gel in 0.5 × Tris–borate–EDTA, stained with ethidium bromide and visualized using the Gel Doc XR from Bio-Rad.

Leaf discs of 7-week-old wild-type tobacco (N. benthamiana) plants were submerged in OA at 20 mM for 2, 6, 12 and 24 h. Then, discs were frozen in liquid nitrogen to extract total RNAs from similar-sized leaf samples using Trizol (Invitrogen, Life Tech, Grand Island, NY, USA) and treated with RNase-free DNase (Promega Corporation, Madison, WI, USA) in the presence of 1 U/μL of RNasin^® Plus RNase Inhibitor (Promega) for 2 h at 37 °C. Untreated leaf discs were used as a negative control. RNAs were quantified, and 30 μg (300 ng/μL) were used for RNA-Seq. Purified RNA samples were shipped to Genomics Unit, MSU, USA in three replicates for deep sequencing and generation of datasets (at least 100 million reads per sample). Filtered reads were aligned with up to two mismatches to the tobacco (N. benthamiana) draft genome sequences (v0.4.2 Scaffolds + NrContigs) of the Sol Genomics Network: (ftp.solgenomics.net/genomes/Nicotiana_benthamiana/annotation/Niben.genome.v0.4.4.transcripts.annotated.fasta). RSEM v1.1.6, an RNA-Seq quantification tool, was used to estimate the relative abundances and expected read counts for the transcripts. RSEM uses the Bowtie aligner (Bowtie v0.12.1) to map the reads against the transcripts. Transcript quantification of the reference-aligned reads was performed with RSEM, which allowed for the assessment of transcript abundances based on the mapping of RNA-Seq reads to the assembled transcriptome. Expected read counts were used as input to differential expression analysis by EdgeR (version 3.0.0, R version 2.1.5). The median was used as the common dispersion factor for differential expression (DE) analyses. The remaining unmapped sequences were re-aligned against the contigs collectively de novo assembled using the Trinity RNA-Seq Assembly package (r2013-02-25) from total unmapped sequences of all treatments and replicates [21]. Blastx was performed against the NCBI non-redundant protein database with an E-value cut off of 1e⁻⁵. Fold change of DE transcripts was calculated against the published tobacco actin gene that exhibited relatively constant expression levels. Accordingly, FC of ≥ 5 was selected to detect DE transcripts that might significantly have a role during PCD. To identify clusters with functional enrichment, we determined a significant Pearson correlation through permutation analysis [22]. The resulting clusters were refined by visual inspection and analyzed for GO term enrichment using Blast2GO (http://www.blast2go.org/). We also clustered the RPKM data to provide a representation of absolute abundance of the transcripts.

This procedure was done to validate RNA-Seq data analysis and to validate the occurrence of VIGS in tobacco. In the first experiment, eight differentially expressed (DE) transcripts representing selected clusters were validated. In the second experiment, VIGS lines of the studied PCD-related genes were validated with the WT and VIGS of empty pTRV2 used as controls. Leaf discs of 7-week-old N. benthamiana plants of the different genotypes were submerged in OA at 20 mM and samples were taken after 0, 2, 6, 12 and 24 h for the first validation, while 0 and 2 h for the second validation. The latter indicates the time point where the selected genes were highly expressed. Then, discs were frozen in liquid nitrogen to extract RNAs. Total RNA was extracted from three similar-sized (10 mm²) leaf discs per plant (approximately 50 mg tissue) collected from upper leaves using Trizol (Invitrogen) and treated with RNase-free DNase (Promega Inc.). First-strand cDNA was synthesized using 2.5 μg of total RNA, 0.5 μg oligo (dT) primer and Superscript II reverse transcriptase (Invitrogen) to a final volume of 20 μL. sqRT-PCR was performed in 20-μL reactions using 1 μL cDNA, 1 × PCR buffer (with 1.5 mM MgCl₂), 200 μM dNTPs, 200 nM of each gene-specific primers (Table S1) and 0.2 U of Taq DNA polymerase (Promega Inc.). Forty PCR cycles for each gene product include denaturation at 94 °C for 15 s, annealing at 55 °C for 30  s, and extension at 72 °C for 45 s. Amplicons were analyzed on a 1.2% agarose gel stained with ethidium bromide and visualized using the Gel Doc XR from Bio-Rad Laboratories (Hercules, CA, USA). The actin gene (372 bp) was used in either experiment as a housekeeping control. Primers were designed using Netprimer software (http://www.premierbiosoft.com/netprimer/index.html) with the following criteria: length ∼ 20 bases, GC content ∼ 50%, minimal secondary structures, comparable annealing temperatures (55 °C) of the primer pairs, and PCR products of 275–375 bp.

For detecting the TRV virus in the newly emerged leaves [23], a first-strand from TRV RNA in leaf RNA extracts of VIGS lines was synthesized using TRV1 primer (5′-CAGTCTATACACAGAAACAGA-3′). Then, PCR amplification was performed using primers of TRV1 and TRV2 (5′-GACGTGTGTACTCAAGGGTT-3′). VIGS of empty TRV2 was used as a positive control. Initial denaturation (2 min at 94 °C) was first done, followed by 40 cycles (denaturation for 15 s at 94 °C, annealing for 2 min at 58 °C, extension for 1 min at 72 °C; and a final extension step for 10 min at 72 °C).

An experiment was conducted to detect the performances of five T-DNA insertion KO SALK lines at 0, 100, 150 or 200 mM NaCl as compared to the WT plants (Col). Mutant seeds (Col) were provided by the SALK Institute, Genomic Analysis Laboratory (SIGnAL) (http://signal.salk.edu/tdnaprimers.2.html). The lines were screened for homozygous segregants by PCR using different PCD-related gene-specific primers and left border (LB) T-DNA primers. To score the germination rate at day 6, 30 seeds/replicate/treatment/line were surface sterilized, sown in Petri dishes containing MS medium [24] with different concentrations of NaCl (0, 100, 150 or 200 mM). Experiments were conducted in three replicates and analyzed using a randomized complete block design. The plates were kept in the dark at 4 °C for 2 days and then shifted to 21 ± 2 °C (day/night) under a 16-h light/8-h dark cycle. Emergence of the radicle from seed coat was considered as completion of germination. To score other characteristics under salt stress, two-week-old plants of different mutants, grown at 0 mM NaCl under a 16-h light/8-h dark cycle at 21 ± 2 °C (day/night), were transferred to the soil (4 plants/pot/replicate) with different concentrations of NaCl (0, 100, 150 or 200 mM) to grow for 2 more weeks. Light intensity was ∼ 175 μmoles/m²/s¹ for the 16-h photoperiod. Then, the individual plant root length, the number of leaves per plant, and the rosette area were scored for 4-week-old plants of different SALK lines as well as the control WT plants. The experiment was conducted using a randomized complete block design with three replicates. Statistical analyses of different experiments were performed following the procedure outlined by Gomez and Gomez [25], and multiple comparisons were performed based on the Duncan's New Multiple Range test [26].

Qualitative and quantitative time-course analysis of cell death upon treatment with OA in tobacco WT leaf discs were measured using Evans blue staining. Gradual cell death was clearly observed across time of OA treatment (Fig. S1). Our results indicated that treatment beyond 24 h is not necessary. Longer treatment time (48 h) resulted in recovering poor quality RNAs for further gene expression analysis. Kim et al. [19] indicated substantial increases in tobacco cell death 48 through 96 h after OA treatment with concentration (e.g. 20 mM) similar to this in the present study. However, they indicated that a 96-h treatment resulted in extensive tissue collapse. This reflects the difficulty of isolating good quality RNA for gene expression analysis. They also indicated that lower concentration of OA (10 mM) yielded similar but delayed results.

Based on the kinetics of the cell death in response to OA, tobacco gDNAs were isolated at several time points. DNA cleavage resulting in a characteristic DNA laddering pattern associated with apoptotic cell death started to appear 12 h after OA treatment (Fig. S2). OA was reported by Kim et al. [19] as a non-specific multifunctional phytotoxin in plant tissue that enhances fungal pathogenicity and induces PCD by inhibiting the oxidative burst of the plant.

This analysis was conducted in order to detect new PCD-related genes induced by OA treatment in tobacco (N. benthamiana). Involvement of these genes in conferring salt stress tolerance in Arabidopsis was done. These genes will be candidates for future improvement of plant salt tolerance via genetic transgenesis. Sequencing of cDNAs recovered from RNA samples of tobacco collected during the course of PCD yielded between 51 and 55 million reads corresponding to over 5 billion nucleotides per sample (Table 1). The raw sequencing reads will eventually be deposited in the Small Read Archive (SRA) at GenBank. Between ∼ 60 and 65% of tobacco reads were mapped to the draft genome sequences in exonic regions (v0.4.2 Scaffolds + NrContigs) of the Sol Genomics Network (http://solgenomics.net/), while the remaining ∼ 35–40% matched no sequences (Table 1). The unmapped reads of all treatments and replicates were collectively de novo assembled and the resulting contigs were used as a new reference genome. Then, the unmapped reads were aligned to this new reference genome and the resulted transcripts were added to those mapped on the draft genome for subsequent analysis.

RPKM-derived read counts were compared using a likelihood ratio test [21]. Statistical analysis was reliable when applied to genes with an RPKM value ≥ 2. To determine DE transcripts, a two-fold (or greater) change in expression and false discovery rate (FDR) of 10⁻³ or less was required. The resulting number of DE transcripts from the alignment with tobacco (N. benthamiana) reference genome as well as the de novo assembled transcripts was ∼ 2650 (data provided upon request). DE transcripts were, then, utilized in generating the expression clusters across time of OA treatment. Expression profiles of the DE transcripts were determined by a cluster analysis based on the k-means method using Pearson's correlation distance so that the similarity in relative change for each transcript or among transcripts across time of OA treatment was determined. These data were, then, subjected to hierarchical clustering using the Pearson correlation as the distance metric (Fig. 1). Two-dimensional hierarchical clustering classified DE profiles into 43 expression clusters according to the similarity of their expression profiles in which 16 of them (Table S2) fitted within the required expression profiles of upregulation after 2 h. The selected clusters might include genes with major effects on triggering or retarding PCD. A number of 23 PCD-related genes were selected for further analysis (Fig. 2) in which eight of them successfully validated the RNA-Seq data via sqRT-PCR (Fig. 3) across different time points of OA treatment. Data across the different selected clusters indicated the regulation of a number of 23 PCD-related genes (Table 2 and Fig. 2). These genes encode Bax inhibitor-1 proteins, Mlo, metacaspase 1, programmed cell death protein 2-like, apoptosis-inducing factors, plasminogen activator inhibitor, Stm1, xylosyltransferase, and exostosin.

Effectiveness of VIGS in tobacco was proven via the use of pTRV2 plasmid carrying a fragment of the PDS gene, e.g., TRV2::PDS. Silencing of the PDS gene has resulted in photo-bleaching of the new leaves 10 days after infiltration (Fig. 4). TRV virus was successfully traced in the 23 VIGS lines via sqRT-PCR with primers specific to the 16-kDa gene (463 bp) of TRV (Fig. 5). sqRT-PCR was also used to confirm occurrence of knockdown in different genes via VIGS (Fig. 6). The results of PCD assay, in terms of Evans blue staining after 24 h, in VIGS of different PCD-related genes are shown in Table 3. Two controls were used in the study, e.g., WT and VIGS with empty pTRV2. The results indicated that VIGS of five PCD-related genes (lines VIGS-G8, VIGS-G13, VIGS-G15, VIGS-G18 and VIGS-G21) responded differently from WT and VIGS with empty pTRV2. One of them, namely Bax Inhibitor-1 (BI-1), has an antiapoptotic effect where the mean relative cell death value due to the exposure OA in its knockdown VIGS lines was significantly higher than the two controls. The other four genes have an apoptotic effect as the mean relative cell death values under treatment in their knockdown VIGS lines were significantly lower than for the two controls. These four genes encode Mlo-like protein, apoptosis-inducing factor (AIF) homolog a, Stm1 and RCD1-SRO-TAF4 (RST). Analogs of these five genes in Arabidopsis were further analyzed for their role during salt stress.

We speculated to detect the gene analogues in Arabidopsis and use their knockouts (Table S3) available at the SALK Institute, Genomic Analysis Laboratory (SIGnAL, http://signal.salk.edu/tdnaprimers.2.html) in order to detect the influence of these genes under salt stress. The SALK lines were cultivated in order to recover the mutants in homozygous state on an MS medium with kanamycin (100 mg/L) used as a selectable marker for T-DNA insertions. It is worth noting that the overall performance of mutant plants under normal growth conditions was similar to the WT plants. Seeds of five lines of T-DNA insertion homozygous KO mutants as well as those of the WT were first examined at different concentrations of NaCl (0, 100, 150 and 200 mM) in order to detect the day where shoot emergence starts. This sign is considered as seed germination. We observed that seed germination was initiated at day 3 under different salt concentrations, while day 4 for seeds put to germinate under control condition. We observed no further changes in germination rate across genotypes and treatments after day 5 of seed germination. Therefore, we decided to score germination percentages of different lines under different salt concentrations at day 6 (Fig. 7a). The results indicated almost no significant differences among genotypes up to 150 mM NaCl. However, significant differences in terms of germination percentages across genotypes were found under 200 mM NaCl. The germination percentage of KO mutant of Mlo gene (G18) was significantly higher than the other mutants as well as the WT under 200 mM NaCl. On the other hand, the germination percentage of BI-1 mutant (G13) was significantly lower than the other genotypes under the highest concentration.

Other sets of plants of different genotypes were grown on the MS medium for two weeks under normal condition and transferred to the soil to grow for two more weeks under different salt treatments. Performance of 4-week-old plants for two mutants with contrasting genotypes (G8 referring to the apoptotic gene with RST domain and G13 referring to the antiapoptotic gene BI-1) based on PCD assay in tobacco, grown under 0 and 200 mM NaCl as compared to the WT control (Col) is shown in Fig. S3. The characteristics statistically analyzed for the 4-week-old plants of different genotypes are shown in Fig. 7b–d. The results for the mean root length (mm) and rosette area (cm²) indicated that BI-1 mutant had the highest performance under 200 mM NaCl as compared to the WT control, while Mlo mutant has the lowest. There are no significant differences among mutant and WT plants under different salt concentrations for the number of leaves/plant. We can conclude that the results of the salt stress experiment in Arabidopsis are in harmony to those of PCD assay for two gene analogs in tobacco.

The gene encoding the inhibitor of apoptosis-promoting Bax1 (BI-1) protein interacts with Bcl-2 protein to block the action of Bax in inducing PCD. In plant, BI-1 gene confers resistance against stress-induced metabolic perturbations or PCD [27,28] by blocking Bax gene product in orchestration with Bcl-2 gene. The gene encoding the functional barley MLO protein is a negative control element of cell death and of defence responses to the fungus Blumeria graminis [29]. Interestingly, it was reported that the gene encoding Bax inhibitor-1 induces breakdown of Mlo-mediated penetration resistance to the fungus in barley. The recessive mutant mlo2 alleles resulted in spontaneous leaf cell death [30] in pepper plants and enhanced resistance against the virulent fungus Xanthomonas campestris pv. vesicatoria (Xcv) [31]. In other words, the gene encoding BI-1 might block PCD either indirectly or directly. The first approach involves the interaction with Bcl-2 gene to suppress Bax encoding gene and the second involves cross talking with Mlo gene to block PCD. Accordingly, the reduced levels of PCD of knockdown line of Mlo encoding gene in the present study (Table 3) might be due to the accumulative effects of BI-1 gene (natural) and VIGS (artificial).

In conclusion, we can either retard apoptotic genes or overexpress antiapoptotic genes in plant breeding programs as a useful approach to improve the salt or drought stress tolerance of agronomically important field crops. This conclusion does not apply to all the studied genes, as their function is much more complicated to be adopted without experimentation.