Rosa hybrida cultivars, Akito, Fuchsiana, Furiosa, H30, Topsun, and Tropical Amazon, were grown on commercial flower farms in Naivasha, Kenya, and transported dry to the UK via air freight. Stems were collected within a day of arrival, transported dry to the laboratory, and rehydrated at 5°C overnight in tap water before being placed in vases. Stem ends were cut upon arrival to the laboratory and prior to being placed in vases, cut to a final stem length of 47 cm. Any foliage showing visible damage or on the lower 20 cm of the stem was removed prior to being placed in vases. All flower stems were kept at 21°C for the duration of vase life, with a 12 h light cycle of 15–20 μmol in m⁻² and s⁻¹, from cool white fluorescent tubes. Stems were placed in either 100 ml tap water, 2% sucrose (w/v), or commercial rose flower food (Chrysal) at the manufacturers' recommended concentration to replicate consumer conditions, in bunches, or individually in clear glass vases. Vases were arranged on benches in a completely randomized design and re-ordered every few days to reduce the effect of placement within the vase life room.

The vase-life was determined as the time from stems being placed in the vase (tap water 1 L in bunches of 8 stems with three replicate bunches) until a termination factor was reached. Flower quality was assessed daily, with the day of termination and terminating factor recorded for each stem. Termination factors included the following: bending of the peduncle (necking), visible color change (blueing or browning) of the petals, and wilting.

Individual stems of “H30” and “Fuchsiana” were held in 2% sucrose with the addition of Pseudomonas fluorescens (5 × 10⁴ cfu per mL) to induce necking. Stems were weighed to the nearest mg just prior to being placed in vases. Stems at each physiological necking stage (straight, <90°, and >90°; Figure 1E) were harvested and weighed as necking occurred. The percentage relative fresh weight (RFW %) for each necking stage was calculated using Equation 1, where W_t was the weight of the stem (g) at harvest and W_{t = 0} was the weight of the stem (g) at the start of vase life.

Equation 1:

Necking was induced in ‘H30', as above, using Pseudomonas fluorescens. Peduncle sections (1 cm) for straight, <90°, and >90° physiological necking stages were excised using a sterile razor blade and placed immediately in a standard formaldehyde alcohol acetic acid (FAA) fixative solution (Ruzin, 1999). Samples were fixed overnight at 5°C and stored at 5°C until used.

For light microscopy, peduncle samples were fixed in FAA, rinsed and dehydrated in graded ethanol at concentrations from 30% to absolute ethanol (1 h for each). After overnight in absolute ethanol, they were transferred to acetonitrile for 1 h, followed by infiltration with a 50:50 acetonitrile/ Spurr resin mixture for 24 h. Peduncles were then infiltrated with 100% Spurr resin and polymerised at 60°C. Transverse sections (500 nm) were cut using a microtome and stained with toluidine blue. Slides were scanned using an Olympus BX51 microscope, with an Olympus CC12 color microscope camera and an Olympus dotSlide automated scanner using a 40 × objective.

For scanning electron microscopy, peduncle samples were cut into 1-mm transverse sections and mounted uncoated onto specimen stubs. A FEI Quanta 250 scanning electron microscope was used (Biomedical imaging unit, Southampton), with 10 kV accelerating voltage, 8.7 mm working distance (WD), and low vacuum mode 60 Pa chamber pressure.

Bunches of “H30” stems were held in tap water with rose flower food (Chrysal) to replicate consumer conditions. As previously described (Lear et al., 2019), for each stage of necking (straight, <90°, and >90°) 2-cm stem sections were excised using a sterile razor blade and flash-frozen in liquid nitrogen. Necked and control (straight-stemmed) samples were taken as the necking stages appeared, and were stored at −80°C.

For each of the necking stages, three peduncle sections were ground to produce one pooled sample with three pooled samples per stage. Total RNA was extracted from each sample as previously described (Moazzam Jazi et al., 2015), and adapted for use with 1.5 mL microcentrifuge tubes. Genomic DNA (gDNA) contamination was removed from the RNA samples using the RapidOut DNA removal kit (Thermo Scientific). Samples were quality tested using a Qubit fluorometer and then sequenced using an Illumina NovaSeq 5000 to produce paired-end reads for each sample. One set of replicate samples was collected and sequenced in 2014 with a mean sequencing depth of 62.5 million, and the other two in 2016 with a depth of 25.3 million reads.

The Galaxy web platform (usegalaxy.org) was used to process and analyse the raw sequencing data (Afgan et al., 2016). Raw reads were filtered and trimmed using Trim Galore! v0.4.3.1 with a phred quality score cut off of 20 (Krueger, 2017). The reads were then aligned to the Rosa chinensis Old Blush genome v2 using HISAT v2.1.0 (Kim et al., 2015; Raymond et al., 2018).

Read counts were produced per exon using DEXSeq count v1.20.1 (Anders et al., 2012), with a minimum alignment quality threshold of 10. Differential expression analysis was completed using DESeq2 v2.11.40.1 (Love et al., 2014). Two factors were considered in the statistical model, necking stage and sequencing year, with the necking stage used as the primary factor.

Differentially expressed genes were functionally annotated with gene ontology (GO) terms using the InterProScan annotation files for the Rosa chinensis genome (Jung et al., 2014). The GO terms were then grouped into broader Plant GO slim biological process groups for visual analysis. Homologous Arabidopsis thaliana (TAIR10) identifiers were used in enrichment analysis with AgriGO v2.0 and STRING (Tian et al., 2017; Szklarczyk et al., 2018), as identified for the Rosa chinensis genome using Blastp, with an e-value cut off of 1e⁻⁶ (Jung et al., 2014). For expression analysis of senescence-associated genes, all genes associated with the term “senescence” for Arabidopsis thaliana were downloaded from uniport and were analyzed using STRING. The expression of transcription factors was analyzed using the Plant transcription factor database v5.0 for Rosa chinensis (Jin et al., 2016).

Necking was induced in “H30” using Pseudomonas fluorescens as described above, and peduncle samples for each physiological necking stage (straight, <90°, and >90°) were collected as necking occurred. The gDNA removal and complementary DNA (cDNA) synthesis, using oligo-dT and random primers, were performed with a Quantitect Reverse Transcription Kit (Qiagen). The cDNA was 1/10 diluted with RNase-free water, to a working concentration of 5 ng/μL (assuming equal efficiency of reverse transcription) and stored at −20°C. Reactions (20 μL) were run using a Rotor-Gene Q thermocycler (Qiagen, Model: 5-Plex) with 1× Rotor-gene SYBR Green PCR master mix (Qiagen), 0.5 μM forward primer, 0.5 μM reverse primer, RNase-free water, and 10 ng cDNA (2 μl of 1/10 diluted template cDNA). All primers are listed in Supplementary Table 1. A two-step cycling program was used with an initial DNA polymerase activation at 95°C for 5 m, followed by 40 cycles of denaturation at 95°C for 5 s, and combined annealing and extension at 60°C for 10 s. All qPCR reactions were run with two technical replicates and four biological replicates for each stage. The Ct values were recorded at a 0.2 threshold. Expression levels were calculated (Pfaffl, 2001) relative to RhGACPC2, with the necking stage “Straight” as the control.

In a comparison across six commercial rose cultivars, mean vase life varied significantly between 5 and 14 days (Figure 1A). The Cv.s H30 and Akito were found to have significantly shorter vase lives than all the other cultivars tested, with mean vase lives of 7.6 and 5.3 days respectively (p < 0.05) (Figure 1A). When the cause of vase life termination was assessed, all of the “H30” and “Akito” stems were terminated due to water stress-associated factors, namely, necking and wilting (Figures 1B,C). In contrast, none of the cv. Fuchsiana, Topsun, or Tropical Amazon stems were terminated due to necking during the observed 14-day vase life period. The majority of cv. Topsun stems were also terminated due to water stress-associated factors. However, in this case, only petal wilting was seen, while the “Fuchsiana” stems were all terminated due to a change in petal color, including both blueing and browning (Figure 1C). In no cultivar were 100% of stems terminated by necking.

To investigate the effects of necking in more detail, the addition of Pseudomonas fluorescens cultures to the vase water was tested to induce necking in both “H30” and “Fuchsiana” flower stems. Necking was consistently induced by the P. fluorescens treatment and had a significant effect on the RFW of the flower head (p < 0.001; Figure 1D). There was no significant difference between the mean RFW of cv.s H30 and Fuchsiana at each necking stage and no significant interaction between necking stage and cultivar in terms of RFW (Supplementary Figure 1).

The necking process was divided into two phases: an early phase from straight peduncles to those bending by <90°, and the later stage of bending from <90° to ≥ 90°. Analysis of necking peduncles identified the mid-point of the bend in a necking peduncle to be on average at 19.7 (± 1.2) mm below the flower head, with a mean peduncle length of 97.8 (± 2.2) mm (Supplementary Figure 2). Transverse sections of the predicted mid-point in necking were, therefore, taken at ~20 mm below the flower head for straight peduncles, and <90° peduncles, where the bending point could not easily be identified. Light microscopy sections of necking peduncles showed shrinkage of the pith as necking progresses (Figures 2A–C), with an undulation of the epidermis, clearly seen between ≥90° and straight peduncles (Figures 2D,E). Shrinkage of the pith was also visible in scanning electron microscopy (SEM) images, with greater contrast between the straight and the ≥90° sections (Figures 2F,G).

Necking was recorded under realistic consumer vase life conditions, including commercial flower food, but without induction of the necking. Libraries were prepared using RNA extracted from 20-mm segments of the peduncle centered on the mean mid-point of the necking bend. They were divided into those with a bend of <90° and ≥ 90°, as well as the equivalent region in non-necking stems. Combining data from three biological replicates, a total of 113, 126, and 100 million paired-end reads, respectively, were produced for straight, <90°, and ≥90° sample groups through next-generation sequencing (NGS) (Supplementary Table 2). Alignment of the trimmed reads to the Rosa chinensis' Old blush genome with HISAT2 achieved an overall average read alignment rate of 87.6% (Table 1).

A total of 38,956 unique Rosa chinensis genes were identified across all samples, with 34,252 (87.9 %) of the genes appearing in all three of the necking stages, and most genes identified in the straight and <90° stages (Table 1). Of the 38,956 identified genes, 76.2 % had homology to Arabidopsis thaliana proteins resulting in 14, 983 unique TAIR10 identifiers, representing 95.3 % of the total TAIR10 proteins with homology to the Rosa chinensis genome.

A total of 3,647 genes showed significant differential expression across all three comparison groups, based on DESeq2 (p adjust. FDR ≤ 0.05). Following differential expression analysis with DESeq2, there were most up- and downregulated genes (3,598 genes) in the comparison between ≥90° and straight necking stages, with a range of 2.2 to −1.79 log₂ fold change. In contrast, fewest differentially expressed genes (DEGs), only 56, were identified between necking stages <90° and straight. In two of the three comparisons (<90° vs. straight and >90° vs. <90°) more DEGs were down- rather than upregulated (Table 2). Comparing early and later bending, 43 genes were uniquely differentially expressed in the early part of the bending process (from straight to <90°), while many more (794) uniquely changed in expression in the later bending phase (<90° to ≥90°). Only 13 DEGs were shared between early and late bending. The highest number, however, 2,838 genes, only changed significantly in expression between straight peduncles and those bending by ≥90° (Supplementary Figure 3A).

Of the 3,598 Rosa chinensis genes whose expression changed in the >90 vs. straight peduncles, 92 % (3,312 genes) were assigned a TAIR10 identifier, equaling a total of 2,865 unique TAIR10 identifiers. Using the TAIR annotations, of the 761 DEGs whose expression changed in more than one of the comparisons (Venn intersections, Supplementary Figure 3A), 234 genes could be assigned a biological process GO term and could be arranged into a biological process heat map (Supplementary Figure 3B). Most changes were consistent in direction of change between the early (<90° compared to straight) and later (≥90° compared to <90°), response with more downregulation than upregulation of expression. Expression of most genes involved with protein phosphorylation, protein ubiquitination, and proteolysis were downregulated, and expression of genes involved with protein de-phosphorylation and protein glycosylation were upregulated in ≥90° necking peduncles. Expression of genes related to metabolic processes (including carbohydrate metabolism and oxidation-reduction), regulation of transcription, and transmembrane transport were both up- and downregulated across the necking stages. However, expression of genes involved with starch biosynthesis, photosynthesis, phosphorelay signal transduction, and metal ion transport was downregulated in ≥90° stage necking peduncles, compared to straight (control) and <90° peduncles. Similarly, expression of genes related to microtubule-based processes, lipid metabolism, and response to stimuli was downregulated in both <90° and ≥90° necking stages in comparison to straight (control) genes. In contrast, the expression of genes related to trehalose biosynthesis was upregulated.

As the ≥90° vs. straight comparison of gene expression yielded the largest number of differentially expressed genes (3,598 genes), this group of genes was further analyzed with the enrichment analysis tools AgriGO v2 (Tian et al., 2017) and Search Tool for Retrieval of Interacting Genes/Proteins (STRING; Szklarczyk et al., 2018) using Arabidopsis thaliana (TAIR10) identifiers to detect associated processes and pathways. Of the 2,865 unique TAIR10 annotated genes in this comparison, 762 were upregulated and 1,261 were downregulated (Log₂ fold change >0.5 and <-0.5, p adjust. < 0.05). Parametric analysis of gene set enrichment (PAGE) for biological process with AgriGO v2 identified more significantly down- than upregulated GO terms (Supplementary Figure 4). However, transport (GO:0006810) was found to be a significantly upregulated GO term, with a z score of 4.21 and included a total of 315 transport-associated genes (p adjust. < 0.05). Under the parent terms of biological regulation and developmental process: growth (GO:0040007), regulation of cell size (GO:0008361), and cell differentiation (GO:0030154) all showed significant down-regulation, with 86, 61, and 107 associated genes, respectively (p adjust < 0.05). Interestingly, the child terms of response to a stimulus: tropism (GO:0009606) and response to abiotic stimulus (GO:0009628) were also found to be significantly downregulated with 23 and 310 identified genes, respectively (Supplementary Figure 4). Of the 23 tropism-related genes identified, 18 were found to be associated with gravitropism and three with phototropism. All three phototropism genes (NPH3, PHYB, and PKS1) were downregulated, as well as all the tropism genes associated with auxin-activated signaling pathways (ACB1, ABCB19, AUX1, PIN1, and RAC3) (Supplementary Table 3).

Using singular enrichment analysis (SEA) with STRING (Szklarczyk et al., 2018), and considering the up- and downregulated genes separately, galactose metabolism (GO:ath00052) was identified as a significantly upregulated pathway (FDR < 0.05; Supplementary Table 4A). In contrast, photosynthetic pathways (ath00196 and ath00195) and phenylpropanoid biosynthesis (ath00940) were significantly downregulated (FDR < 0.05; Supplementary Table 4B). However, starch and sucrose metabolism (ath00500), and the more general biosynthesis of secondary metabolites (ath01110) and metabolic pathways (ath01100) included genes that were both significantly up- and downregulated (Supplementary Table 4B).

Parametric analysis of gene set enrichment (PAGE) was also carried out with STRING, using the entire >90° vs. straight dataset and corresponding Log₂ FC values (p adjust, < 0.05). This approach revealed water channel activity to be significantly enriched, with an enrichment score of 4.10. Out of the nine aquaporin genes identified, eight showed downregulation (Supplementary Table 5). TIP1-3 was the only water channel gene to be upregulated. However, TIP1-3 also had the lowest mean count of all the aquaporin genes, with an average count of just 34 compared to an average count of over 9,500 for PIP2-7 and almost 1,400 for TIP4-1.

Annotation of the ≥90° vs. straight and <90° vs. straight rose peduncle differentially expressed genes against the transcription factor gene list identified 225 genes from 43 transcription factor families in at least one of the two contrast groups (p adjust. <0.05; Supplementary Figure 5). The majority of the families (27) were on average downregulated during necking, while 15 were upregulated. Of particular interest are transcription factor families related to stress and growth responses (Figure 3A). All the genes belonging to the growth-regulating factor (GRF), squamosa promoter binding proteins (SBPs), TCP, YABBY, and Zinc-finger homeodomain (ZF-HD) transcription factor families were downregulated throughout the necking process, and the majority of ARF, B3, and MYB transcription factors were also downregulated. In contrast, all of the differentially expressed genes in the WRKY family, and the majority of genes in the MYB-related, NAC, and WUS homeobox-containing (WOX) transcription factor families showed significant upregulation. Other transcription factor families including ethylene-responsive factor (ERF) and GRAS comprised a mixture of significantly up- and downregulated genes in response to necking.

Based on a comparison to the Uniprot database, 64 of the differentially expressed genes between ≥90° and straight peduncles were identified as senescence-related. Using STRING's multiple protein search, 44 of these genes had connected nodes (Figure 3B). A cluster of eleven autophagy-related genes (ATG) members were significantly upregulated and co-expressed, with many of the connections having been experimentally determined in other systems. A rose gene with homology to histone deacetylase 9 (HDA9) was also significantly upregulated, and co-expressed with APG9 (ATG9), linking the ATG cluster to the four histone protein genes H2AXA, HTB9, HTA12, and AT1G09200, all of which were significantly downregulated. A cluster of three genes with homology to mitogen-activated protein kinases (MPKs), MPK1, 4, and 13, and an MPK kinase (MKK2) were also differentially expressed, with the majority being upregulated. Histone 3.1 (AT1G09200/ AT5G65360), the bidirectional sugar transporter (SAG29), and ribulose carboxylase/oxygenase activase (RCA) were the most downregulated genes with connected nodes, (Log₂ FC −1.16, −1.09, and −0.91 respectively), while transcription factor TCP9 and mechanosensitive ion channel protein 10 (MSL10) (Log₂ FC −1.02 and −0.99) were also strongly downregulated, but with no known connections to the other senescence genes.

A gene, with homology to a nitrate transporter (NRT1.5), was the most upregulated of all the senescence-related genes identified, (Log₂ FC 1.71); NDR1/HIN1-lie protein 10 (YLS9) was the second most upregulated (Log₂ FC 1.64), however, also had no known connections to the other senescence-associated genes and was, therefore, also not shown.

Several changes in expression were noted within the sucrose metabolism pathway (ath00500), when all the up- and downregulated genes with TAIR10 identifiers from the ≥90° vs. straight gene set (p adjust. <0.05) were included (Supplementary Figure 6, Supplementary Table 6). Genes with homology to APL1 and ADG1 encoding the large and small subunits for ADP glucose pyrophosphorylase, (shown as EC: 2.7.7.27 on Supplementary Figure 6) were both downregulated, with an average Log₂ FC of −0.99. A rose gene with homology to the glucose-starch glycosyltransferase, GBSS1 (EC: 2.4.1.242) was also downregulated (Log₂ FC of −0.75). In contrast, genes involved in starch catabolism, α-amylase AMY2 (EC: 3.2.1.1), and β-amylase BAM1 (EC: 3.2.1.2) were both upregulated. The β-amylase 2 (BAM2), also included in EC: 3.2.1.2, was downregulated, however, there was a mean upregulation of 0.93 (Log₂ FC) of BAM1 and BAM2. Moreover, the enzymatic activity of BAM2 is much weaker than BAM1 (Fulton et al., 2008).

Glycosidase genes glucan endo-1,3-β-glucosidase genes 1, 2, and 4 (EC: 3.2.1.39) and glycosyl hydrolase 9B1, 9B3, and 9C2 (EC: 3.2.1.4) were downregulated in ≥90° necking peduncles indicating a reduction in 1,3-β-Glucan and cellulose degradation. However, genes involved with the conversion of uridine diphosphate (UDP)-glucose to sucrose and D-fructose (EC: 2.4.1.13), sucrose to D-glucose and D-fructose (EC: 3.2.1.26), D-fructose to D-fructose-6P, and D-glucose to D-glucose-6P (EC: 2.7.1.1) were all upregulated. Together, these lead to a potential increase in both D-fructose-6P and D-glucose-6P for use in glycolysis. Trehalose biosynthetic genes were also upregulated (EC: 2.4.1.15 and EC: 3.1.3.12), including trehalose-phosphatase/synthase 9 and 11, and trehalose-6-phosphate phosphatase G and J. None of the genes encoding trehalase activity (EC: 3.2.1.28) were differentially expressed (p adjust. < 0.05) in the ≥90° necking peduncles, indicating a potential increase in trehalose accumulation.

Nine galactose metabolism-related genes were significantly upregulated (> 0.5 Log₂ FC, p adjust < 0.05) in the ≥90° necking peduncles, compared to straight peduncles (Figure 4A, Supplementary Table 7). The expression pattern of three of these genes, RhUGE5 (EC: 5.1.3.2), RhSIP2 (EC: 2.4.1.82), and RhDIN10 (EC: 2.4.1.82) was confirmed by real-time PCR (Figure 4B), inducing necking using Pseudomonas fluorescens. This verified that consumer vase life conditions and the induced necking were consistent in expression changes elicited. All three genes showed a progressive increase in expression with peduncle necking. The mean expression of RhATBF1, RhPGM2, and RhAGAL3 was also higher in ≥90° necking than in straight or <90° necked peduncles, however, the differences were not significant (p < 0.05).

Seventeen phenylpropanoid biosynthesis genes were significantly downregulated (>0.5 Log₂ FC, p adjust <0.05) in ≥90° necked, compared to straight peduncles (Figure 5A), associated with five EC codes. To verify the RNAseq data expression of four genes: RhOMT1, RhPXPX42, RhCAD9, and RhHCT was also analyzed by real-time PCR (Supplementary Table 7). Where several rose genes were associated with a single EC point, the representative was selected as the gene with the highest expression (based on the mean count).

Expression of RhPXPX42, RhCAD9, and RhHCT genes was significantly downregulated between ≥90° necked and straight peduncles in both the RNA seq. and qPCR datasets, and, indeed, real-time PCR also revealed a significant downregulation between straight and <90° necked peduncles. Mean RhOMT1 expression was also lower in <90° compared to straight peduncles, however, the change was not significant (p < 0.05) in the rt-PCR analysis. Thus, three out of the four genes analyzed by rt-PCR mirrored the downregulation seen with RNAseq as necking progressed (Figure 5B).