Commercial brown onion cv. Orlando was grown in sandy soil in the northwestern Negev desert in the years 2013–2015. Onions were not sprayed with maleic hydrazide before leaf drop, as per common agricultural practice, and did not undergo field curing. Onions were harvested manually at 80–100% fallen leaves (top–down) and the leaves were removed with a sharp knife, leaving a ca. 10-cm long neck above the bulb, as previously described (Eshel et al., 2014), and omitting the postharvest fast-curing process. The onions were placed in a dark storage room at 2°C and 70% relative humidity until use. Experiments were conducted with undamaged bulbs of regular shape. Bulbs freshly harvest from the field contain a single completely dry skin (“skin”) and underneath it several scales characterized as thin yellowish scales. These yellowish scales were numbered as scales 1–4 (as illustrated in Figure 1A), have the ability to undergo cell death and form additional dry skins. As shown in our previous publication, scales 1–4 will be developed into skin during postharvest long cold storage or in several days during Fast Curing (Eshel et al., 2014). For morphological analysis, roots were removed from the bulbs which were separated into different successive scales from exterior to interior scales of bulb (Figure 1B).

The TUNEL reaction (TdT-mediated deoxy-uracil nick end labeling) is used to analyze DNA fragmentation by labeling the 3′-OH ends of the DNA-strand breaks. Histological analyses were performed on 10-μm-thick onion-scale sections cut by microtome, according to the method described by Teper-Bamnolker et al. (2010). For TUNEL staining, fixed tissues were rehydrated with Histoclear and decreasing concentrations of ethanol (100, 70, and 30%). Tissue permeabilization was performed with 20 μg/ml Proteinase K (Gibco BRL) in 10 mM Tris, pH 7.5, and 5 mM EDTA, pH 8.0, at 37°C for 30 min. After washing the tissue twice with phosphate-buffered saline (PBS), lysing enzyme (4 mg/ml) in 5 mM EDTA, pH 8.0, was added for 20 min with incubation at 37°C. TUNEL reaction was performed on slides using the In Situ Cell Death Detection kit with fluorescein (Roche Applied Science) according to the manufacturer’s instructions. To visualize nuclear morphology in onion cells, samples were stained with 4′-6-diamidino-2-phenylindole (DAPI; Sigma) at 1 μg/ml in PBS buffer for 10 min. Negative and positive controls were treated identically except for the omission of the enzyme solution (terminal deoxynucleotidyl transferase), or after incubation with DNAse I (Roche) for 30 min, respectively. DAPI and TUNEL-positive staining were observed with an IX81/FV500 confocal laser-scanning microscope (Olympus) equipped with a 488-nm argon ion laser and a 405-nm diode laser. DAPI was excited with the 405-nm diode laser, and the emission was filtered with BA 430–460 nm filters. TUNEL was excited with 488 nm of light, and the emission was filtered with a BA505IF filter. The transmitted light images were obtained using Nomarski differential interference contrast, and then images were subjected to Fluo View 500 software supplied with the confocal laser-scanning microscope. Percentage of TUNEL-positive cells is determine as compared to in the DAPI staining vs. cells number in the same tissue.

The first four scales of the onion bulb (from the outside) were analyzed by SEM. The scales were cut into small sections (circa 5 × 5 mm), fixed in 70% ethanol overnight at room temperature and then dehydrated in ethanol solutions (90% for 1 h, 95% for 1 h, and 100% for 2 h, twice). Sections were mounted by their longest axis in vertical orientation on a slotted cryo-SEM stub and held in position using a cryo-glue (OCT) compound (Tissue Tek^®, Sakura, Tokyo, Japan). Stubs were plunged into pre-frozen liquid nitrogen, and transferred under vacuum to a cryochamber (Alto 2100 Gatan, Abingdon, UK) with the stage temperature cooled to -180°C. Sections were fractured using a cold blade and the stage temperature was raised to 95°C to sublimate any contaminating ice and to enhance features on the fractured surface. The stage heater was turned off and once the temperature recovered to approximately 160°C, sections were coated with gold for 60 s. Prepared sections were then placed in the SEM (JSM LV6360 SEM, Jeol, Tokyo, Japan) and mounted on the stage with the temperature maintained at 160°C. The fractured surfaces were examined at ×500 magnification with an accelerating voltage of 12 kV before setting parameters for microanalysis.

The 1st thin outer scale and the 5th inner fleshy scale of the onion bulb were cut into small sections (circa 3 × 3 mm) and fixed in 2.5% (w/v) glutaraldeyde in 0.1 M cacodylate buffer (pH 7.4) for 2.5 h at room temperature, and then moved to 4°C for an additional 16 h. The tissues were rinsed four times, 10 min each, in cacodylate buffer and postfixed and stained with 2% (w/v) osmium tetroxide and 1.5% (w/v) potassium ferricyanide in 0.1 M cacodylate buffer for 2 h. Tissues were then washed four times in cacodylate buffer followed by dehydration in increasing concentrations of ethanol—30, 50, 70, 80, 90, and 95%—for 10 min each step followed by 100% anhydrous ethanol three times, 20 min each, and propylene oxide twice, 10 min each. Following dehydration, the tissues were infiltrated with increasing concentrations of Agar 100 resin in propylene oxide, consisting of 25, 50, 75, and 100% resin for 16 h each step. The tissues were then embedded in fresh resin and allowed to polymerize in an oven at 60°C for 48 h. Tissues embedded in blocks were sectioned with a diamond knife on an LKB 3 microtome and ultrathin sections (80 nm) were collected onto 300-mesh, carbon/formvar-coated copper grids. The sections on the grids were sequentially stained with uranyl acetate and lead citrate for 10 min each and viewed with a Tecnai 12 TEM 100kV (Phillips, Eindhoven, the Netherlands) equipped with a MegaView II CCD camera and Analysis^® version 3.0 software (SoftImaging System GmbH, Munster, Germany). Representative photographs are presented.

Tissue sections of onion scales from the first outer thin scale toward the inner fleshy scales of the bulb were hand-cut to determine cell viability. These small sections were immediately submerged in 50 μM FDA for 10 min at room temperature in the dark to maximize fluorescein formation. FDA fluorescence, indicating cell viability, was observed in a confocal laser-scanning microscope (Olympus). Excitation and emission wavelengths were 493 and 510 nm, respectively. Only cells that exhibited bright green fluorescence in the cytosol were considered viable.

Three different onion scales: first outer thin scale (1st), intermediate fleshy scale (5th) and inner fleshy scale (8th), were sampled in two biological replicates and used for RNA isolation, cDNA synthesis and sequencing. Scale samples were frozen in liquid nitrogen and stored at -80°C until RNA extraction. Total RNA of each sample was extracted using the CTAB protocol (Chang et al., 1993). Samples were treated with DNase (Epicentre, Madison, WI, USA) according to the supplier’s instructions. RNA purity and integrity were verified by RNA 6000 Nano Assay on an Agilent 2100 BioAnalyzer with a minimum RNA integrity number value of 7. Library preparation and sequencing were performed at the Genome Center, Life Sciences and Engineering, Technion, Israel. Fourteen single-end RNA-Seq libraries with a length of 100 nucleotides were prepared using Illumina Hiseq2000 and Trueseq protocols.

Raw reads were subjected to filtering and cleaning as follows: SortMeRNA tool was used to filter out rRNA (Kopylova et al., 2012) (doi: 10.1093/bioinformatics/bts611.) Then the FASTX Toolkit was used¹ (version 0.0.13.2) for: (i) trimming read-end nucleotides with quality scores < 30 using fastq_quality_trimmer; (ii) removing reads with less than 70% base pairs with quality score ≤ 30 using fastq_quality_filter. A total of ∼270 million cleaned reads, obtained after processing and cleaning, were assembled de novo using Trinity software (version trinityrnaseq_r20140717 2.1.1 (Grabherr et al., 2011) with the trimmomatic option to remove adaptors (Bolger et al., 2014) and 25 mer k-mer size. The assembled sequences that shared a number of k-mers (the set of isoforms of a gene) were referred to as “contigs.” The sets of all sequences that shared at least one k-mer were referred to as components. Filtering of the likely contig artifacts and low expressed contigs was carried out as follows: (i) abundance estimates were calculated for each contig using the RSEM software (Li and Dewey, 2011); (ii) only contigs representing more than 1% of the per-component (IsoPct) expression level were retained. The resulting de novo assembly generated a transcriptome consisting of 45,892 contigs with N50 of 1,694, median contig length of 965 bp, and average contig length of 1,208.62 bp.

Sequencing data were deposited in the NCBI Sequence Read Archive (SRA) database as bioproject PRJNA326316 (SRX1959523, SRX1959528, SRX1959536).

To assess the similarity of the onion transcriptome to those of other model and closely related species, analysis of sequence similarity was performed using the BLAST (Basic Local Alignment Search Tool) algorithm with an E-value cut-off of 10^-5 (Altschul et al., 1990). The BLASTX algorithm was used to search protein databases with a translated nucleotide query for comparison of the assembled contigs with sequences deposited in the databases of Arabidopsis², Oryza sativa³ and SwissProt proteins^4,5.

Transcript quantification (the number of reads per gene) from the RNA-Seq data was performed using the Bowtie aligner (Langmead et al., 2009) and the RNA-Seq by Expectation-Maximization (RSEM), which handles read-mapping uncertainty with a statistical model by estimating maximum-likelihood expression levels (Li and Dewey, 2011). Differential expression analysis was performed with the edgeR software suite (Robinson et al., 2010). Transcripts that were more than fourfold differentially expressed with false discovery-corrected statistical significance of at most 0.001 and log 2 of the fold change lower than -2 or greater than 2 were considered differentially expressed (Benjamini and Hochberg, 1995). The expression patterns of the transcripts in the different samples were studied using cluster analysis of the differentially expressed transcripts in at least one pairwise sample comparison. Then following the Trinity protocol (Haas et al., 2013), expression normalization was designed using TMM (trimmed mean of M-values), following FPKM (fragments per feature kilobase per million reads mapped) calculations. Hierarchical clustering of the normalized gene expression [using centralized and log 2 transformation; (Haas et al., 2013)] and heatmap visualization were performed using R Bioconductor (Gentleman et al., 2004). We used the “Venny” tool (Oliveros, 2007) for Venn diagram construction.

Gene ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations were performed using the GSEA server⁶. GO-enrichment analysis was carried out using the Blast2GO (Conesa et al., 2005) program based on Fisher’s Exact Test (Upton, 1992) with multiple testing correction of false discovery rate (FDR) (Benjamini and Hochberg, 1995). The threshold was set to a FDR with corrected P-value of less than 0.05. GO analysis was performed by comparing the GO terms in the test sample to the GO terms in a background reference. GO provides a structured and controlled terminology to define gene products according to three domains: molecular function (the biochemical activity of a gene product), biological process (operations or sets of molecular events to which the gene product contributes), and cellular component (cell parts in which a gene product is active).

Peripheral onion scales are chronologically older than those located toward the center of the bulb (Figure 1A). To characterize the morphological differences between scales located at different positions in the bulb, we analyzed the thickness and color of sequential scales. Scale thickness increased from the outer toward the inner bulb position (Figures 1B,C). The thinnest tissue was the skin, and scale thickness increased gradually with a fold change of sequential scales of around 1–1.6; maximum thickness was found in the 8th inner fleshy scale (Figures 1B,C). Scales also differed in color. The measured hue angle (H°) increased dramatically from the skin toward the inner scales of the bulb. The skin had a typical brown color with H° value of 58, whereas the following scales had a light yellow to white color with high H° values of around 115 (Figures 1A,D). The chroma (C) which, according to McGuire (1992) represents color saturation, decreased gradually from 24.7 in the skin to close to 1 in the three innermost scales of the bulb (Figure 1D).

To characterize the differences in tissue anatomy between the different outer scales compared to the skin, SEM analysis was performed for the four chronologically ordered outer scales. The analysis revealed gradual deterioration of cells progressing from the 3rd inner scale, which included only what seemed to be intact cells, toward the skin which contained only a thin layer of collapsed cell residues (Figure 2). Unexpectedly, cell deterioration in the 1st outer scale initiated from the inner epidermis and progressed to the outer epidermis (Figure 2). Cell viability in the different scales was monitored by FDA staining. All of the cells in the 3rd inner scale were viable, whereas viability declined successively toward the skin (Figure 2). Within the 1st outer scale, a progressive reduction in cell viability was observed from the inner epidermis toward the parenchyma cells, while only the outer epidermis remained viable (Figure 2).

To detect fragmented nuclear DNA in situ, a TUNEL assay was applied to six scales, from outer to inner, in the same bulb. TUNEL-positive cells could be detected extensively in the 1st outer scale (around 99% of the nuclei), but their density gradually decreased from the periphery toward the center to the 6th inner fleshy scale, where no DNA fragmentation was detected (Figure 3). The average percentage of TUNEL-positive cells in the 1st, 2nd, 3rd, 4th, 5th, and 6th scales were 99, 93, 64, 55, and 4% respectively.

To examine the ultrastructural changes in the cell more closely, onion parenchyma cells of the outer thin (1st) and inner fleshy (5th) scales were analyzed by TEM. Unusual structures were found in the parenchyma cells of the outer scale which were not detected in the inner scale cells. In some parenchyma cells, there was visible formation of vesicles in the cytoplasm (Figures 4A,B), knob-like bodies on the surface of the tonoplast, and rupture of the tonoplast (Figures 4C–E). Disappearance of cytoplasm and organelles was clearly detected in the outer scale cells (Figures 4C–E). In addition, condensed granular substances were observed in the vacuoles in the outer scale cells (Figures 4C,D), but were not detected in the inner scale cells. In the latter, the organelles were clearly observable in the cytoplasm, the tonoplast seemed intact and the vacuole content was clear (Figures 4F–J).

Six selected pools of mRNA samples, representing three scale positions: outer scale (1st), intermediate fleshy scale (5th) and inner scale (8th) of the onion bulb served for the construction of high-throughput parallel RNA-Seq libraries (Supplementary Figure S1). Each of the cDNA libraries yielded 16.4–20.9 million 100-bp one-end reads (Table 1). Quality trimming and filtration of all of the libraries resulted in ∼270 million cleaned reads that were assembled using Trinity, generating 45,891 contigs for the transcriptome catalog. The average contig length was 1208.61 bp; half of these (N50) were at least 1,694 bp long. The transcriptome catalog of onion was compared with the database of Oryza sativa, the most sequenced and annotated monocot species. BLASTX search against the rice database resulted in at least one significant hit for 25,703 contigs of the onion transcriptome (56%). Comparison against the TAIR (Arabidopsis thaliana) database resulted in significant hits for 25,398 contigs (54%). GO terms were assigned to 25,340 contigs and comparison to the KEGG database⁶ resulted in hits for 18,305 contigs. Databases of Arabidopsis, rice and SwissProt showed a general similarity of 54–56% with our data. InterProScan via the Blast2GO tool found 25,648 contigs with known protein motifs in the multiple databases (PROSITE, PRINTS, Pfam, ProDom, SMART, TIGRFAMs, PIR superfamily, SUPERFAMILY, Gene3D, PANTHER and HAMAP).

The outer scale (1st), intermediate fleshy scale (5th) and inner scale (8th) transcriptomes were analyzed for differential gene expression. About 85% of the cleaned reads could be mapped to the de novo transcript catalog of onion (Table 1). Three pairwise combinations of the scales were compared and 2,542 significantly differentially expressed genes (DEGs) were identified. The comparison of 5th vs. 1st scales revealed 722 DEGs, 183 upregulated and 539 downregulated (Figure 5A). Out of 2,382 DEGs in the comparison of 8th vs. 1st scales, 827 were upregulated and 1,555 were downregulated. Finally, we identified only 509 DEGs in the comparison of 8th vs. 5th scales, 257 upregulated and 252 downregulated (Figure 5A). These results showed the highest proportion of DEGs between the 8th and 1st scales, an intermediate proportion between the 5th and 1st scales, and the lowest differential gene expression for the 8th vs. 5th scales (Figure 5A). Thus, the number of DEGs tended to increase as the distance between scales increased. Moreover, in the comparisons of either 5th or 8th scale to the 1st scale, there were more downregulated than upregulated DEGs (Figure 5A).

To identify shared and unique DEGs among the scales, we generated a Venn diagram based on the three comparisons. Overlap of the 5th vs. 1st scale comparison with the 8th vs. 1st scale comparison revealed 556 shared DEGs (Supplementary Figure S2). Thus most of the DEGs from the comparison of 5th vs. 1st scale (556 out of 722 DEGs, 70%) were included in the 8th vs. 1st scale comparison. Out of 509 DEGs in the comparison of 8th vs. 5th scale, 345 DEGs (59%) were included in the comparison of the 8th vs. 1st scale. However, only 22 DEGs overlapped between the 5th-to-1st scale comparison and the 8th-to-5th scale comparison (Supplementary Figure S2). These results indicated a higher difference between more distant scales. In addition, in the comparison of 8th to 1st scale, 59% of the DEGs (1,407 out of 2,832) were specific to this comparison. Only 74 DEGs were shared by all three groups.

Hierarchical cluster analysis of gene expression revealed clearly differentiated patterns of gene expression among the three investigated scales. The heatmap analysis revealed five main clusters of coexpressed genes (Figure 5B), each labeled in a different color according to the pattern of the difference among the three scales: cyan and red clusters describe a low expression level in the 1st scale compared to the 5th and 8th scales, whereas the blue, green and pink clusters describe high expression in the 1st scale compared to the 5th and 8th scales. We focused on four of the five clusters: cyan, blue, green, and pink, which had a common GO-term profile. The cyan cluster (776 genes) profile was the only one that included genes with higher expression in the inner 8th scale (Figures 5B and 6A). GO-enrichment analysis of this cluster revealed gene-expression patterns reflecting metabolic activity, with most of the genes associated with cell growth, metabolism and developmental process (Figure 6A). The blue cluster (963 genes), the biggest cluster of coexpressed genes, was characterized by substantially higher expression levels in the 1st scale compared to the 5th scale, and a further moderate decrease in the inner 8th scale (Figures 5B and 6B). Genes in this cluster were related to defense response, respiratory burst, and response to stress, PCD and diversity of signal transduction processes (Figure 6B). The green cluster (335 genes) was also characterized by decreased expression levels with a slight decrease in the 5th scale compared to the 1st, and a more significant decrease in the 8th scale (Figures 5B and 6C). The most significant genes in this cluster were associated with peptide catabolic process, and sucrose and oligosaccharide transport (Figure 6C). These findings suggest changing carbohydrate contents in the different scales. The pink cluster (381 genes) was characterized by decreased expression levels from the 1st scale to the two inner scales (Figures 5B and 6D). Most of these genes were annotated as biological processes related to pigmentation, such as flavonoid biosynthetic processes and phenylpropanoid metabolic processes (Figure 6D). The final, red cluster (Figure 5B) was relatively small, consisting of 87 genes, and lacked any common biological role. The results of this cluster are therefore not shown.

We analyzed the expression of PCD-related genes in the transcriptome data in the three different scales. Most of these genes were induced in the 1st outer scale and had lower expression in the 5th and 8th inner scales as shown (Figure 7). The PCD-related genes were subdivided into several classes; oxidases, mitogen-activated protein kinase (MAPK) cascade, Ca²⁺/Calmodulin cascade, transcription factors, proteases, hormonal regulation, signal cascades, cell wall degradation and nucleases. The oxidases genes encoding to copper amine oxidases (CuAO) and amine oxidases (POAs) were highly expressed in the 1st scale compared to the two inner scales.

Among the MAPK cascades, genes encoding to MAP kinase 4 (MPK4), MAP kinase kinase 5 (MAPKK5), MAP kinase kinase 15 (MAPKKK15) and MAP kinase substrate 1 (MKS1) were highly expressed in the 1st scale (Figure 7; Table 2). This is probably followed by expression of genes members of the Ca²⁺/calmodulin cascade such as calcium-binding EF-hand family (CBP-EF) and calmodulin-domain protein kinase (CDPK) which also were highly expressed in the 1st scale (Figure 7; Table 2). This cascades then targets various effector proteins in the cytoplasm or nucleus, which include other kinases, enzymes, or transcription factors (TFs) (Khokhlatchev et al., 1998).

The targets of the MAPK signal cascade include various types of TFs. Among these, many of the genes are expressed in the 1st scale, including genes encoding WRKY DNA-binding protein 30 and 75 (WRKY30, WRKY75), MYB-domain family protein (MYB), NAC-domain family (NAC) zinc finger family such as the C3HC4-type RING, which presumably play a role in regulating specific protein levels via the ubiquitination pathway (Figure 7; Table 2). These are likely candidates for the control of later PCD processes. Proteases encoding genes were overrepresented in the 1st scale mainly due to expression of genes homologous to genes encoding papain family cysteine protease (RD19C), cysteine proteinase (CYP), Aspartic proteinase A1 (ASPA1) and metacaspases 1 and 9 (AMC1, AMC9) (Figure 7; Table 2).

Several genes homologous to PCD-associated nuclease genes, including those encoding to endonuclease 2 (ENDO2), endonuclease 4 (ENDO4) and bi-functional nuclease 1 (BFN1) were upregulated in the 1st scale. In parallel, degradation processes which were found to occur in the outer scale were associated by elevated expression of genes homologous to genes known to be involved with cell wall degradation, such as genes encoding pectinacetylesterase, involved in the degradation of plant cell wall pectin components. In addition genes homologous to genes encoding enzymes such as β-xylosidase 1 (BXL1), glucosyl hydrolase superfamily (GHF), pectin methylesterase inhibitor (PMEI), β-glucosidase 11 (BGLU11), and UDP-glucosyl transferase 88A1 (UGT88A1), were highly upregulated in the 1st scale (Figure 7; Table 2).

The sequential scales of the onion bulb, which have a common origin as leaf bases, differed in various characteristics. Their thickness increased gradually from the outer to inner scales, as also revealed by light microscopy (Figures 1B,C). The scales also differed in color, the skin having a typical brown color, whereas the other scales were light yellow to white (Figures 1A,D). The brown color of the skin has been associated with peroxidase-dependent oxidation of phenolics (Takahama, 2004). Onion scales contain quercetin 4′-glucoside and quercetin 3,4′-diglucoside as major phenolics, with their concentration increasing from the inner to outer scales (Tsushida and Suzuki, 1995; Hirota et al., 1998). Outer-scale browning is a result of loss of cellular compartmentalization due to cell death, since in living cells, enzymes that catalyze glucosidation and the respective substrates may be compartmentalized (Takahama and Oniki, 2000).

Structural analysis of the outer scale demonstrated cell death within that scale that initiated from the parenchyma cells and spread gradually toward the outer epidermis (Figure 2). Since cell death and tissue desiccation are unique to the outer scale and start from the inner tissue, we hypothesized that these processes are active and programmed. Similarly, it has been reported that cell death in leaf senescence also initiates from mesophyll cells and then proceeds to other cell types (Lim et al., 2007). Mochizuki-Kawai et al. (2013) demonstrated that PCD begins earlier in the mesophyll cells than in the epidermal cells of tulip petals during senescence.

A typical hallmark of PCD, DNA fragmentation, was detected by TUNEL assay in the outer scales and gradually disappeared toward the inner scales (Figure 3). This pattern of nuclear DNA fragmentation demonstrates, for the first time, that the mechanism of cell death in the outer scales during skin formation is programmed. PCD has been detected in naturally senescing leaves from a variety of plants (Yen and Yang, 1998; Simeonova et al., 2000; Cao et al., 2003; Uzelac et al., 2008). However, to the best of our knowledge, there are no studies on PCD during skin formation, although there are studies related to the role of PCD in protective dry tissue such as seed-coat development. DNA fragmentation has been found at early developmental stages of Arabidopsis seed-coat formation (Nakaune et al., 2005) and during the development of seed coats of cowpea (Lima et al., 2015). Radchuk et al. (2010) identified DNA fragmentation during barley pericarp development. Differences in the frequency of TUNEL positive nuclei were also detected between the inner and outer faces of the scales. As we observed in the structural analyses, cell death proceeds gradually from the inner to the outer face of the scale. This might be explained by the loss of the nuclei in the inner face dead cells which does not allow detection by the TUNEL assay.

Our TEM analyses of ultrastructural changes in the scales supported the results of PCD in the outer scales. The irregular features of the cells in the outer scale may be involved in the process of PCD, since no such abnormalities were found in cells of the inner scale that remained viable (Figure 4). Similar changes, including small vesicles in the cytoplasm and granular substances in the vacuoles, have been found during developmental PCD in other species (Teper-Bamnolker et al., 2012). Gunawardena et al. (2004) reported that dying cells contain many electron-dense knob-like bodies and electron-lucent blebs on the inner surface of the tonoplast in cells undergoing PCD during remodeling of leaf shape in the lace plant. Differentiation-induced PCD occurs in various processes, for example, in xylem tracheary elements (Bollhöner et al., 2012), leaf morphogenesis (Gunawardena, 2008), floral development (Rogers, 2006), root cap formation (Fendrych et al., 2014) and anther tapetum layer formation (Van Hautegem et al., 2015). Senescence-induced PCD is the last stage of organ senescence which occurs in all tissues of an organ (Thomas, 2013). It can be suggested that the PCD in the outer scale is a developmental PCD process, as it is associated with skin formation.

In this study, we found that the most overrepresented metabolic processes in the outer scale are associated with defense response, PCD processes, carbohydrate metabolism and flavonoid biosynthesis (Figure 6). It is well known that outer onion browning derives from loss of cellular compartmentalization due to cell death, which leads to enzymatic oxidation of phenolic compounds by peroxidase (Takahama, 2004). Upregulation of genes related to flavonoids and secondary metabolic processes may be relevant to this browning mechanism (Figure 6). In contrast to the outer scale, gene expression associated with increased metabolism and developmental growth processes was overrepresented in the inner scales (Figure 6). This difference in gene expression probably results in the inner scales having viable tissue that preserves its metabolic activity as compared to the outer scale, a distressed suicide tissue that is programmed to die and form skin.

We propose a conceptual model for PCD in the outer scale during skin formation, based on the PCD-related gene expression patterns (Figure 7). This hypothetical model suggests sequential processes from initiation to execution of PCD, within the outer scale (Figure 8). ROS-producing enzymes, including copper amine oxidases (CuAOs) and polyamine oxidases (PAOs), induce ROS accumulation (Gechev et al., 2006; Kärkönen and Kuchitsu, 2015). Apoplastic CuAOs and PAOs have been shown to play a key role as a source of ROS (such H₂O₂) in tissue differentiation organ development as well as PCD (Moore et al., 2003; Pourtau et al., 2006; Wingler and Roitsch, 2008; Tavladoraki et al., 2016). Sensing of ROS during senescence and PCD-related processes, activates two major signal cascades: mitogen-activated protein kinase (MAPK) and Ca²⁺/calmodulin, which leads to activation of key TFs to initiate these processes (Bieker et al., 2012; Rogers, 2012). Plant hormones have also a regulatory important role in PCD-related processes, thus genes involved in their biosynthesis and signal transduction are essential for these processes. Among these are ethylene and jasmonic acid (JA), which accelerate the processes (Altschul et al., 1990; Kim et al., 2015). We found that ethylene-response factor (ERF) genes and genes required for JA biosynthesis, such 12-oxophytodienoate reductase (OPR2) gene, and involved in controlling JA responses, such jasmonate ZIM-Domain 1 gene (JAZ1), were overrepresented in the 1st outer scale and may have crucial role in skin formation in the outer scales (Figure 7; Table 2). Ethylene is involved in developmental processes that employ PCD such as senescence (Graham et al., 2012), perforation formation in the lace plant (Dauphinee et al., 2012) and in maize endosperm development (Young et al., 1997).

In the outer senescing scale, there was overrepresentation of several genes from the MAPK cascade and the signaling cascade of calcium regulation, including genes encoding calcium and calmodulin-binding proteins. The MAPK signal cascade in plants has been shown to be involved in a wide range of cellular responses, such as biotic and abiotic stress responses, hormone responses, cell proliferation, cell death, and developmental processes (Nakagami et al., 2005). Calcium signaling may be a significant component in the regulation of senescence processes. Jones (2001) described the involvement of calcium in various types of cell death.

Various TF genes were overrepresented in the outer scale, mainly in the NAC, WRKY, MYB and zinc finger (C3HC4-type RING) families (Figure 7; Table 2). NAC and WRKY TFs have been associated with senescence in several tissues, including Arabidopsis leaves, petals and siliques (Miao et al., 2004; Wagstaff et al., 2009; Balazadeh et al., 2011; Besseau et al., 2012; Rogers, 2013). Several genes homologous to genes encoding different classes of proteases, including papain family cysteine protease (RD19C), cysteine proteinases (CYP), aspartic proteinase A1 (ASPA1), metacaspase 1 and 9 (AMC1, AMC9) were found to be highly expressed in the outer scale, and may have essential role as executors of PCD. In general, proteases have been shown to be involved in PCD, including cysteine proteases, serine proteases, aspartic proteases (Beers et al., 2004). Plant-specific papain-type KDEL-tailed cysteine endopeptidase (CEPs) gene has been found to be upregulated in the outer scale. This gene is similar to the Arabidopsis gene shown to be involved in PCD during vegetative development as executors for the last step in the process (Zhou et al., 2016). Aspartic proteinase PASPA3, from Arabidopsis was shown to be a marker of developmental PCD in root cap and seed endosperm (Fendrych et al., 2014; Fourquin et al., 2016) and other generally accepted cell death processes, e.g., the tapetum and the tracheary elements (Olvera-Carrillo et al., 2015). In plant PCD, special functions are described for metacaspases (Tsiatsiani et al., 2011). Metacaspases, which are cysteine-dependent proteinases, have been associated with various types of plant PCD, such as in controlled pathogen-induced PCD in Arabidopsis (Coll et al., 2010) and in embryogenesis in Picea (Suarez et al., 2004). In Arabidopsis, there are nine metacaspase genes, although only metacaspase 9 is highly expressed in senescence and PCD (Kwon and Hwang, 2013). Metacaspase 9 is also strongly upregulated in older Arabidopsis petals (Winter et al., 2007). At the final stages of the suggested model, PCD-associated endonuclease genes, mainly the bi-functional nuclease 1 (BFN1) gene were upregulated in the 1st outer scale. BFN1 gene, a senescence and PCD marker has been shown to be induced during both senescence and developmental PCD (Pérez-Amador et al., 2000; Fendrych et al., 2014; Sakamoto and Takami, 2014). This gene suggested to be involved in DNA degradation during senescence or PCD processes in plants (Aoyagi et al., 1998; Muramoto et al., 1999; Panavas et al., 1999). Farage-Barhom et al. (2008) demonstrated the intracellular localization of BFN1-GFP in Arabidopsis cells undergoing senescence and death, which initial in filamentous structures spread throughout the cytoplasm, which then clustered around the nuclei as the protoplasts senesced. This pattern of localization highlights BFN1’s function as a nucleic acid-degrading enzyme in senescence and PCD. Another recent study revealed that nuclease BFN1 is responsible for the rapid cell-autonomous corpse clearance and DNA fragmentation during Arabidopsis root cap cell death (Fendrych et al., 2014). Senescence-related degradation processes were reflected in the enrichment of genes for cell wall degradation in the outer scale. Several genes encode degrading-enzymes such as β-xylosidase 1 (BXL1), β-glucosidase 11 (BGLU11), pectin methylesterase inhibitor (PMEI) and UDP-glucosyl transferase 88A1 (UGT88A1), were found in the 1st outer scale, and they may have a role in scale senescence by controlling the degradation of cell wall components and releasing sugars for respiration (Lee et al., 2007). The transcriptome analysis in this work led to a putative model for senescence in the outer scale that supports the involvement of PCD. The identification of genes involved in senescence-associated PCD in the outer scale revealed a possible basis for cell-death processes in skin formation. Our results clearly demonstrate that cell death in the outer scales of onion bulbs is mediated by programmed processes accompanied by a unique set of morphological and transcriptional alterations, leading to formation of the protective dry skin.