Rose flowers (Rosa chinensis Jacq. cv. Gold Medal) were grown at a greenhouse on the campus of China Agricultural University, Beijing, China. Rose flower opening was divided into six stages: stage 1, partially opened bud; stage 2, completely opened bud; stages 3 and 4, partially opened flower; stage 5, fully opened flower with visible anthers; stage 6, fully opened flower with abscised petals (Figure 1). Rose flowers at different opening stages were harvested. The flower stems were placed immediately in water, and transported to the laboratory within 15 min. The flower stems were re-cut to 20 cm in length under water and placed in deionized water until further processing. The petal was shed at separation layer of abscission zone. Both distal and proximal sides of separation layer belong to abscission zone. Distal side attaches to petal organ, and proximal side attaches to receptacle (Figure 1). Therefore, we took sample of petal abscission zone by excising the base of petal (less than 1 mm in length) and the receptacle where petals locate (less than 1 mm in length). Since petals at stage 6 started to abscise, we only focused on the stages prior to petal shedding. Therefore, AZ samples at stages 1, 3, and 5 were collected and used for the transcriptome profiling. Three biological samples were collected for each stage.

Our preliminary tests demonstrated that R. hybrida cv. Samantha showed better responses to virus-induced gene silencing (VIGS) and much higher silencing efficiency than R. chinensis Jacq. cv. Gold Medal that was used for the transcriptome profiling (data not shown). In addition, the plantlets of R. hybrida cv. Samantha bloom as quickly as 40 days after rooting under our growth conditions. Therefore, rose plantlets of R. hybrida cv. Samantha were selected for VIGS. Rose plantlets were propagated by tissue culture. Rose shoots with at least 1 node and approximately 2 cm in length were used as explants and cultured on Murashige and Skoog (MS) medium supplemented with 1.0 mg/L 6-Benzyl aminopurine, 3 mg/L Gibberellic Acid, and 0.05 mg/L α-Naphthalene acetic acid for 30 days, then transferred to 1/2-strength MS medium supplemented with 0.1 mg/L NAA for 30 days for rooting.

Total RNA was extracted using the hot borate method according to previously described (Wan and Wilkins, 1994), and treated with RNase-free DNase I (Promega) to remove any contaminating genomic DNA. Three biological replicates were performed for each developmental stage (stages 1, 3, and 5). Strand-specific RNA libraries were constructed using the protocol described previously (Zhong et al., 2011), and sequenced on a HiSeq 2500 system according to the manufacturer’s instructions. The raw reads were deposited into NCBI SRA database under accession no. PRJNA325324.

RNA-Seq reads were first processed to remove low quality and adaptor sequences using Trimmomatic (Bolger et al., 2014). Reads shorter than 40 bp were removed. The resulting RNA-Seq reads were then aligned to the ribosomal RNA (rRNA) database (Quast et al., 2013) using Bowtie with default parameters (Li and Durbin, 2010). Reads mapped to the rRNA database were discarded. The high-quality cleaned reads were assembled de novo into contigs using the Trinity program (Grabherr et al., 2011). To remove the redundancy of Trinity-assembled contigs, the contigs were further assembled de novo using iAssembler (Zheng et al., 2011). The final assembled rose contigs were blasted against the UniProt (Swiss-Prot and TrEMBL; The UniProt Consortium, 2014) and Arabidopsis protein (version TAIR10) databases (Lamesch et al., 2012) with a cutoff e-value of 1e^-5. Based on the UniProt and Arabidopsis protein blast results, functional descriptions (human readable descriptions) were assigned to each unigene using AHRD¹. Gene ontology (GO) terms were assigned to the rose assembled transcripts based on the GO terms annotated to their corresponding homologues in the UniProt database. Biochemical pathways were predicted from the rose transcripts using the Pathway Tools (Karp et al., 2002).

To identify differentially expressed genes, high-quality cleaned reads were aligned to rose contigs using Bowtie allowing up to two mismatches. Only the best alignments for each read were retained. Following alignments, raw read count for each rose contig in each sample was derived and normalized to reads per kilo base exon model per million mapped reads (RPKM). The significance of differential gene expression between different samples was determined using DESeq (Anders and Huber, 2010), and raw p-values of multiple tests were corrected using false discovery rate (Benjamini and Hochberg, 1995).

To analyze the transcript abundance of selected genes for confirmation of RNA-Seq data, quantitative RT-PCR (qRT-PCR) reactions were performed, as previously described (Ma et al., 2015b). Briefly, total RNAs were isolated from AZ samples with three biological repeats. 1 μl of the first strand cDNA was used as template with the Step One PlusTM real-time PCR system (Applied Biosystems) using KAPA^TM SYBR^® FAST quantitative PCR kits (Kapa Biosystems). RhActin5 was used as a reference gene (Pei et al., 2013). The primers used for determining transcript abundance were listed in Supplementary Table S1.

A 290 bp fragment in the 3′end of RhIAA16 gene was amplified by PCR from rose cDNAs to specifically silence RhIAA16. The fragment was then inserted into pTRV2 vector to generate the pTRV2-RhIAA16 construct. The primers used for amplifying RhIAA16 are listed in Supplementary Table S1.

Constructs of pTRV1, pTRV2, and pTRV2-RhIAA16 were transformed into Agrobacterium tumefaciens GV3101. A. tumefaciens were cultured in Luria-Bertani medium supplemented with 10 mM MES, 20 μM acetosyringone, 50 μg/ml gentamicin sulfate, and 50 μg/ml kanamycin. The cultured bacterium was collected by centrifuge at 4,000 rpm for 10 min, and re-suspended in infiltration buffer (10 mM MgCl₂, 200 μM acetosyringone, 10 mM MES, pH 5.6) to OD₆₀₀ of ~1.5 (Yin et al., 2015). Mixtures of cultures containing an equal ratio (v/v) of pTRV1 and pTRV2 or pTRV2-RhIAA16 were used for inoculation. VIGS of rose plantlets was performed as previously described (Tian et al., 2014). Rose plantlets, as shown in Supplementary Figure S1, were immersed in bacterial suspension solution and infiltrated under a vacuum at 0.7 MPa for 2 min. After infiltration, plantlets were washed in deionized water, and transplanted into pots containing a mixture of 1:1 (v/v) of peat and vermiculite. The plantlets were immediately placed in dark at 8°C for 3 days in a low temperature incubator (MIR-253, SANYO), and then grown in a culture room at 22 ± 1°C, 40% relative humidity. Three independent experiments were performed. 30 plantlets were used for each experiment. Prior to the petal abscission test, we PCR-screened the plants and determined the transcript abundance of RhIAA16 in petals among the 30 plantlets. We found that the transcript levels of RhIAA16 in more than 30% plantlets were reduced, compared to un-inoculated and empty vector controls. These plants with down-regulated RhIAA16 expression were used for petal abscission assay.

To perform transcriptomic analysis, petal AZs of rose flowers (Figure 1) with three biological replicates at stage 1, 3, and 5 were collected to construct a total of nine RNA-Seq libraries. A total of 75,752,884 paired-end raw reads with length of 100 nucleotides (nt) were obtained. After further filtering and cleaning, a total of 57,312,389 clean read pairs were obtained. De novo assembly of these high-quality cleaned reads generated 80,226 unique transcripts with an average length of 743 bp (Table 1). The size distribution indicated that the lengths of the 19,414 transcripts were more than 1000 bp (Supplementary Figure S2). Correlation coefficients of transcriptome profiles among the nine libraries and between the biological replicates were calculated (Supplementary Table S2). High correlation coefficients were obtained, suggesting the robustness of our RNA-Seq dataset.

To further validate the expression profiles of RNA-Seq data, four selected transcripts were analyzed by qRT-PCR. The results from qRT-PCR analysis were generally in agreement with the expression profiles obtained by RNA-Seq data (Figure 2).

Differentially transcribed genes (DTGs) were determined using a cutoff ratio of >2 or <0.5 when comparing their transcript abundance in stage 3 to that in stage 1 (S3 vs. S1), and/or in stage 5 to that in stage 3 (S5 vs. S3). A total of 2592 DTGs were obtained (Supplementary Table S3). Based on the change in ratio of DTG transcript abundance, the number of DTGs at stage 3/1 and stage 5/3 was counted (Figure 3). Compared with stage 1, 782 DTGs were up-regulated and 300 DTGs were down-regulated in stage 3. Compared with stage 3, 1179 transcripts were increased and 1408 transcripts were decreased in stage 5 (Figure 3), suggesting that major transcriptional dynamic for petal abscission occurs just prior to petal shedding (stage 5).

To evaluate the potential functions of genes that showed transcriptional changes during petal abscission, we identified the GO terms of the DTGs in the biological process category (Figure 4). At stages 3 and 5, the biological processes were enriched in the metabolic process and defense responses including response to abiotic stimulus, external stimulus, organic substance (Figures 4A,B).

To identify the biochemical pathways involved in petal abscission, we analyzed DTGs using the Pathway Tools (Karp et al., 2002). The 108 DTGs at stage 3 and 261 DTGs at stage 5 were classified into 42 and 92 biochemical pathways, respectively (Supplementary Table S4). The major pathways both at stage 3 and 5 included ethylene biosynthesis, starch degradation, superpathway of cytosolic glycolysis, pyruvate dehydrogenase and TCA cycle, and photorespiration (Supplementary Table S4). In addition, one of the major pathways at stage 5 is related to the lactose degradation III pathway (Supplementary Table S4). These results suggested that alterations in carbon metabolism play an important part in rose petal abscission.

Of 2592 DTGs, 150 encoded putative transcription factors (TFs; Figure 5A). More specifically, 15.3% (23/150) belonged to zinc finger family, 13.3% (20/150) to WRKYs, 12.7% (19/150) to ERFs (ethylene responsive factors), and 9.3% (14/150) to Aux/IAAs (Figure 5A). As the largest group of abscission-responsive TFs, the transcript abundance of most of the zinc finger family members (18/23) was increased in stage 3 or stage 5 (Supplementary Table S5). Furthermore, the transcript abundance of all the WRKY family members was induced during rose flower opening (Supplementary Table S5). The results suggested a complex transcriptional reprogramming of petal abscission.

Hormones act as internal cues to initiate abscission process (Addicott, 1982; Estornell et al., 2013). We identified 108 DTGs related to hormone pathways. Among them, DTGs related to auxin and ethylene pathways including 52 DTGs and 38 DTGs, respectively, were the largest group (Figure 5B; Table 2; Supplementary Table S6). In addition, 12 DTGs related to gibberellin were obtained, and 6 DTGs in jasmonic acid pathway (Figure 5B; Supplementary Table S6). No DTGs involved in the biosynthesis/signaling of abscisic acid, cytokinin, brassinosteroid, salicylic acid pathways were detected (Figure 5B).

Among DTGs related to auxin pathway, six auxin transporter genes were identified including five auxin efflux carrier genes (Table 2). In addition, 14 Aux/IAA family members were obtained, of which six members were up-regulated in stage 5, and eight members were down-regulated in stage 5 compared to stage 3 (Table 2). Among DTGs related to ethylene pathway, 15 DTGs encoded ethylene biosynthesis related 1-Aminocyclopropane-1-carboxylic acid oxidase (ACO). Transcript abundance of 10 ACO genes was accumulated in stage 5 compared to stage 3. Furthermore, 19 DTGs encoded ERFs were detected (Supplementary Table S6). The transcript abundance of 16 ERFs was increased in stage 5 compared to stage 3 (Supplementary Table S6). Overall, our results suggested that auxin and ethylene may play central roles in petal abscission of rose.

To identify key regulators governing petal abscission, we initiated functional analysis of DTGs using VIGS. We primarily focused on the up-regulated DTGs, hypothesizing that VIGS-down regulation of these genes might lead to changes in the petal abscission processes. Given the potential important roles that auxin plays in the regulation of abscission, we first examined the functions of several Aux/IAA genes (RSA33069, RSA04500, RSA45030) that were up-regulated in the rose abscission zone using rose cut flowers. We found that VIGS-silencing of the contig of RSA33069 exhibited accelerated petal abscission phenotype. Analysis of the RSA33069 cDNA sequence revealed that it encoded a deduced protein of 253 amino acids with a 762 bp predicted open reading frame (Figure 6A). The predicted amino acid sequence of RSA33069 showed that it belongs to the Aux/IAA family, and has the four canonical conserved domains known for this family (Figure 6A) (Overvoorde et al., 2005). In addition, phylogenetic tree analysis suggested that the protein has high degree of sequence homology to FvIAA16 from Fragaria vesca, therefore was designated as RhIAA16 (Figure 6B). RT-PCR analysis demonstrated that the transcript abundance of RhIAA16 in petal AZ was significantly induced during flower development to peak at stage 5 (Figure 7A). These results were consistent with the transcript abundance changes of RhIAA16 in the RNA-Seq data.

To further confirm the potential role of RhIAA16 in petal abscission, we chose a fragment from RhIAA16-specific 3′ un-translated region (UTR) to silence RhIAA16 in rose plantlets. qRT-PCR results showed that transcript abundance of RhIAA16 in RhIAA16-silenced (TRV2-RhIAA16) petal was significantly reduced compared to TRV2 control (Figure 7B). Petal abscission was detected at 5 days after full opening in the RhIAA16-silenced plantlets in contrast to 9 days after full opening in TRV2 control plantlets (Figures 7C,D), suggesting that silencing of RhIAA16 accelerated the timing of initial petal abscission. At 9 days and 11 days after full opening, petals in ~40.5 and 69.0% of flowers in RhIAA16-silenced plantlets had abscised whereas petals in only 17.8 and 37.9% of TRV2 control flowers were abscised, respectively (Figure 7C).