A gynoecious cucumber (C. sativus L.) line 13-3B was used in this study. The seeds were germinated on wet filter paper in a Petri dish at 28°C in dark overnight. Then the resulting seedlings were grown in a growth chamber under 16 h/8 h with 25°C/18°C in day/night, respectively. Upon two true-leaf stage, plants were transferred to a greenhouse in the experimental field of the Northwest A&F University. Pest control and water management were carried out according to standard practices.

For male flowers induction in the gynoecious cucumber 13-3B, 1000 ppm GA₃ (dissolved in 0.1% ethanol) or deionized water with 0.1% ethanol (Control) were applied by foliar spray for three times at 7 day intervals, starting when the first true leaf was approximately 2.5 cm in diameter. The sex of the flowers on each node of the main stem was recorded until anthesis of flowers on node 25.

In addition, ethylene production in shoot apices was measured after 7 days of the third GA₃ treatment. And the RNA-Seq analyses were performed in shoot apices from the cucumber plants firstly treated with GA₃ for 6 h, 12 h and the Control, respectively. GA₃ was acquired from Sigma-Aldrich Chemical Co. (Shanghai, China).

The ethylene production was measured by gas chromatography as described previously with some modifications (Zhang et al., 2014b). In brief, the excised shoot apices from cucumber plants treated with exogenous GA₃ and the Control were enclosed in 10 mL vessels after weighing and sealed with rubber stoppers. After incubation at 25°C for 16 h, 1 mL of head gas was withdrawn from each vessel using a syringe and injected into a gas chromatograph (GC-9A, Shimadzu, Japan) equipped with a hydrogen FID and an activated alumina column for the measurement of ethylene production. Standard ethylene gas was used for calibrating the instrument. Amount of ethylene was calculated per 1 g fresh weight and per hour.

The shoot apices from the gynoecious cucumber plants were collected at 6 and 12 h after the first treatment with GA₃ and the Control. Samples were immediately frozen in liquid nitrogen and stored at -80°C for RNA-Seq analyses. Total RNA was isolated using the RNA extraction kit (Promega, USA). RNA was checked by RNase-free agarose gel electrophoresis to avoid possible degradation and contamination, and then examined using the NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA, USA) for RNA purity. RNA concentration and integrity were measured and assessed using the Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA) and RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), respectively.

Digital gene expression libraries were constructed using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, Ispawich, USA) following instructions of manufacturer and six index codes were added to attribute sequences to various samples (Wang et al., 2009). Briefly, poly (A) mRNA was isolated from 3 μg total RNA using oligo-dT magnetic beads (Life Technologies, Carlsbad, CA, USA), and then broken into short fragments by adding fragmentation buffer. First-strand cDNA was synthesized using random hexamer-primed reverse transcription, followed by the synthesis of the second-strand cDNA using RNase H and DNA polymerase I. After adenylation of the 3′ ends of cDNA fragments, NEBNext adapter oligonucleotides were ligated to prepare for hybridization, and then the cDNA fragments were purified using AMPure XP system (Beckman Coulter, Beverly, MA, USA) to select the fragments of preferentially 150–200 bp in length. The size-selected, adaptor-ligated cDNA fragments were enriched using Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index Primer in the PCR reaction. PCR products were purified with AMPure XP system and library quality was assessed using the Agilent Bioanalyzer 2100 system. At last, the cDNA libraries were sequenced on an Illumina HiSeq 4000 platform using the paired-end technology by Novogene Co. (Beijing, China).

Raw reads were pre-processed to remove low quality sequences (there were more than 50% bases with quality lower than 20 in one sequence), reads with more than 5% N bases (bases unknown) and reads containing adaptor sequences. Then the clean reads were mapped to the cucumber genome (Chinese long) v2¹ using TopHat (Huang et al., 2009; Trapnell et al., 2009), allowing up to one mismatch. Unigenes mapped by at least one read, in at least one sample, were identified for further analysis.

In this study, samples from three treatments (GA 6 h, GA 12 h, and Control) were prepared for genome-wide expression analyses. Two biological replicates were performed for each treatment, and thus six DGE libraries were sequenced. 44.28–60.29 million raw reads from each library were generated. After removal of low-quality tags and adapter sequences, 42.31–57.63 million high-quality clean reads with a total of 6.35–8.64G bases were obtained. Among these clean reads, the percentage of Q20 (base quality more than 20) and GC was 97.28–97.52% and 43.46–43.63%, respectively (Table 1). Furthermore, we clustered these clean reads into unique tags, which were mapped to the cucumber genome using TopHat (Huang et al., 2009; Trapnell et al., 2009). About 36.28–49.30 million clean reads (85.21–85.74% of total clean reads) from RNA-Seq data in the six libraries were mapped uniquely to the cucumber genome (Table 1).

The R package edgeR used to identify the DGEs (Robinson et al., 2010). The expression level of each gene was calculated and normalized to FPKM. The FDR was used to determine the threshold of the P-value in multiple tests. In our study, the FDR < 0.05 and fold change > 2 were used as significance cut-offs of the expression differences.

Furthermore, GO enrichment analysis of DGEs was performed using the GOseq R package (Young et al., 2010). GO terms with corrected P-value < 0.05 were considered significantly enriched by differential expressed genes.

Quantitative Real-Time PCR analyses were performed using the independent cucumber shoot apices in the same time point of GA₃ application as those used for RNA-Seq. Total RNA was isolated using the RNA extraction kit (Promega, USA), and cDNA was synthesized using the PrimeScript RT reagent Kit (TaKaRa, China). qRT-PCR was carried out using SYBR Premix Ex Taq (TaKaRa, China) on an ABI 7500 Real-Time PCR System (Applied Biosystems, USA). The cucumber α-TUBULIN (TUA) gene was used as an internal control in analyzing gene expression. Three biological replicates were performed for each experiment. The gene specific primers for qRT-PCR are listed in Supplementary Table S5.

To verify the effect of GA on sex expression of cucumber, a gynoecious cucumber line 13-3B was treated with GA₃ and the sex of flowers on each node was recorded until anthesis of the flowers on node 25 of the main stems. As shown in Figure 1, there were no male flower nodes in the control plants, however, GA₃ treatment induced male flowers in the gynoecious cucumber line (Figure 1A), accounting for 51.2% of nodes with male flowers (Figure 1B). Interestingly, formation of male flowers occurred mainly at the lower node positions as compared with the location of female flowers on the main stems (Figure 1A). These observations suggested that exogenous GA can promote male flowers formation in cucumber.

It is well known that ethylene can also control sex determination of cucumber (Malepszy and Niemirowicz-Szczytt, 1991; Bai and Xu, 2013). To assess potential involvement of ethylene in GA-regulated sex expression of cucumber, ethylene production in the shoot apices was measured in GA₃-treated and control plants. As shown in Figure 1C, ethylene production was significantly decreased after GA₃ treatment, suggesting that ethylene might function as a negative factor in GA-regulated male flower formation in cucumber.

To identify the genes and gene networks that are involved in GA-regulated cucumber sex expression, the genome-wide expression analyses were performed to compare the transcriptome profiles of the shoot apices between the gynoecious cucumber plants treated with GA₃ for different time points (6 and 12 h) and the Control through the digital gene expression (DGE) approach (Eveland et al., 2010). Based on deep sequencing, 23,911 unigenes were collected in six libraries. Using fold change > 2 and FDR < 0.05 as the significance cut-offs, we identified 1073 DEGs including 727 up-regulated genes and 346 down-regulated genes after GA₃ treatment for 6 h compared with the Control. And we also found that 1590 genes were differentially expressed, in which 765 genes were up-regulated and 825 genes were down-regulated after GA₃ treatment for 12 h (Figure 2; Supplementary Tables S1 and S2). Moreover, 594 DEGs containing 303 up-regulated genes (Figure 2A) and 291 down-regulated genes (Figure 2B) were shared in the two sets of transcriptome comparisons.

To validate the DEGs identified by RNA-Seq, we performed quantitative real time RT-PCR (qRT-PCR) assays using the independent cucumber shoot apices in the same time point after GA₃ treatment as those used for RNA-Seq analysis. Twenty DEGs were randomly chosen for qRT-PCR analyses, in which 10 genes including five up-regulated genes and five down-regulated genes were from the set of GA 6 h vs. Control, and the other 10 genes containing five up-regulated genes and five down-regulated genes were from the set of GA 12 h vs. Control. As shown in Figure 3, all the 20 genes showed the similar expression patterns in the qRT-PCR analyses as those in the RNA-Seq data, although the particular values of fold-change were different. The Pearson correlation coefficient between qRT-PCR and RNA-Seq data was 0.975 (P = 3.5E-13), indicating that the RNA-Seq results were highly reliable.

Given that ethylene production was significantly decreased in the cucumber plants treated with GA₃ (Figure 1C), we screened the ethylene biosynthetic genes in the DGEs. We found that the transcript of M (CsACS2) gene encoding an ACC synthase was inhibited in both sets of transcriptome comparisons (GA 6 h vs. Control and GA 12 h vs. Control) (Table 2), and the qRT-PCR verification displayed the same expression pattern (Figure 3). Since the M gene is believed to inhibit stamen development in flower buds (Saito et al., 2007; Li et al., 2009, 2012), we speculated that GA might release the inhibitory effect of the M gene and allow male flowers to develop in cucumber.

In addition, another two ethylene biosynthetic genes, CsACO1 and CsACO3 which encode two ACC oxidases, were also differently expressed in the shoot apices after GA₃ treatment. CsACO3 expression was dramatically down-regulated in the set of GA 12 h vs. Control, but not changed in GA 6 h vs. Control. However, there was an increase in the transcript level of CsACO1 in both GA 6 h vs. Control and GA 12 h vs. Control, but the fold change was lower than that of M gene (Table 2). These results suggested that the decreased transcript levels of M and CsACO3 might inhibit ethylene biosynthesis in the cucumber plants treated with GA₃. Nonetheless, an increased expression of CsACO1 following GA₃ treatment was insufficient to rescue the effect of M and CsACO3 on ethylene production.

To further understand the potential functions of DEGs identified by DGE, GO term enrichment analyses (Corrected P-value < 0.05) were carried out in both sets of RNA-Seq data. We found that the DEGs were markedly enriched in biological process and molecular function (MF) groups. For the biological process category, the most significantly enriched GO terms were “cellular carbohydrate biosynthetic process” (P = 1.1E-02) and “regulation of cellular macromolecule biosynthetic process” (P = 5.8E-04) in GA 6 h vs. Control and GA 12 h vs. Control groups, respectively (Figures 4 and 5). While five GO terms including “oxidoreductase activity, acting on paired donors” (P = 1.1E-02), heme binding (P = 1.1E-02), tetrapyrrole binding (P = 1.5E-02), iron ion binding (P = 2.5E-02), and “sequence-specific DNA binding transcription factor activity” (P = 4.6E-02) in the set of GA 6 h vs. Control (Figure 4) and two terms containing “sequence-specific DNA binding transcription factor activity” (P = 2.2E-03) and “DNA binding” (P = 1.1E-02) in the GA 12 h vs. Control group (Figure 5) were detected in the MF category. Furthermore, the transcription factors were dramatically enriched in the DGEs in both sets of data. Accordingly, many ethylene-responsive transcription factors (ERFs) including four in GA 6 h vs. Control (Table 3 and Supplementary Table S3) and nine in GA 12 h vs. Control (Table 4 and Supplementary Table S4) were identified to be down-regulated, consistent with the reduced ethylene production (Figure 1C). Among them, three genes, ERF43 (Csa3M895680) and CRF2s (Csa5M139630 and Csa4M051360), were shared in the two sets. These observations indicated that ERFs may be implicated in cucumber sex expression. Because the ERFs act as positive regulators in ethylene signal transduction pathway (Wang et al., 2002; Guo and Ecker, 2004; Prescott et al., 2016; Zhang et al., 2016), we speculated that the decreased ethylene production inhibited the expression of ERFs, followed by modulated sexual development in cucumber plants treated with GA₃.

Interestingly, we also noticed that an ethylene receptor, CsETR1, was enriched in the GO term of “DNA binding” (Figure 5), and its expression was increased by 2.12-fold in shoot apices after GA₃ treatment for 12 h (Table 4 and Supplementary Table S4), but not changed in GA 6 h vs. Control group. While the qRT-PCR verification revealed the same expression pattern (Figure 3). ETR1 is an important member of ethylene receptors family that acts as negative regulator in the ethylene signaling pathway (Wang et al., 2002; Guo and Ecker, 2004; Light et al., 2016; Prescott et al., 2016). Previous studied have confirmed that CsETR1 plays a negative role in stamen arrest during development of flower buds in cucumber (Wang et al., 2010). In accordance with these findings, we speculated that up-regulated CsETR1 may promote male flowers formation by alleviating stamen arrest in cucumber after GA₃ treatment.

Through GO term enrichment analyses, we further identified an AG (C-class floral homeotic gene) homolog CAG2, which was enriched in the “sequence-specific DNA binding transcription factor activity” group. CAG2 is one of the three AG genes in cucumber that controls pistil development due to its specific expression in the carpel (Kater et al., 1998; Perl-Treves et al., 1998). We found that the transcript level of CAG2 was significantly decreased by 23.48-fold in shoot apices after 12 h of GA₃ treatment (Table 4 and Supplementary Table S4), and the qRT-PCR assay showed the same expression pattern (Figure 3). Our data implied that GA may restrain the femaleness via inhibiting the CAG2 expression in cucumber.

Upon exogenous GA₃ treatment in the gynoecious cucumber line 13-3B, the male flowers were markedly induced, meanwhile, ethylene production in the shoot apices was significantly decreased (Figure 1) due to collaborative regulation by three ethylene biosynthetic genes such as M (CsACS2), CsACO1, and CsACO3. Notably, M and CsACO1 were significantly down-regulated and up-regulated in both sets of transcriptome comparisons, respectively, however, CsACO3 transcript was only decreased in GA 12 h vs. Control group (Figure 3; Table 2). Since both CsACO1 and CsACO3 encode ACC oxidases, and they showed opposite expression patterns and similar fold changes in RNA-Seq data (Table 2), we speculated that interplay between up-regulation of CsACO1 and down-regulation of CsACO3 offset the ACC oxidase activity in shoot apices after 12 h of GA₃ treatment. This implied that the reduced M transcript might play major role in inhibiting ethylene biosynthesis. Given that the M gene can directly inhibit stamen development in cucumber flower buds (Saito et al., 2007; Li et al., 2009, 2012), our data revealed that GA might release the inhibitory effect of the M gene on stamen arrest and restrain ethylene production, by down-regulating the M gene expression (Figure 6, right panel).

In ethylene signal transduction pathway, the receptors such as ETRs function as negative regulators, while ERFs, downstream components of receptors, act as positive transcription factors (Wang et al., 2002; Guo and Ecker, 2004; Light et al., 2016; Prescott et al., 2016; Zhang et al., 2016). RNA-Seq data displayed that the ethylene receptor CsETR1 in cucumber, which is thought to be a negative regulatory factor in stamen arrest of flower buds (Wang et al., 2010), was markedly up-regulated after 12 h of GA₃ treatment (Figure 3; Table 4 and Supplementary Table S4). Up-regulation of CsETR1 occurred at relatively later time point than down-regulation of the M gene (Table 2), suggesting that the increased CsETR1 transcript was not directly caused by GA₃ rather by reduced ethylene production. Moreover, the expression levels of some ERFs were dramatically decreased after GA₃ treatment (Tables 3 and 4, Supplementary Tables S3 and S4), and they were probably involved in cucumber sex expression through inhibiting maleness or promoting femaleness. However, this understanding was based on bioinformatics analysis, the precise roles of ERFs on cucumber flower development remained unclear and should be verified in future studies using advanced physiological and molecular techniques. These observations indicated that increased CsETR1 expression may stimulate male tendency through direct inhibition in stamen arrest or down-regulation on the ERFs transcript in cucumber after GA₃ treatment (Figure 6, right panel).

Furthermore, despite we have shown that GA can regulate cucumber sex expression in cooperation with ethylene, how GA modulated ethylene and which genes participated in this process were unknown. Given that the DELLA proteins are central repressors of GA responses (Sun, 2010, 2011; Plackett et al., 2014), and accumulating evidence suggested that DELLA proteins play important roles in ethylene-mediated plant growth and development processes through interactions with some regulatory factors in ethylene signaling pathway, such as CTR1 (CONSTITUTIVE TRIPLE RESPONSE1), EIN3/EIL1 (ETHYLENE INSENSITIVE 3/EIN3-LIKE 1), RAP 2.3 (RELATED TO APETALA 2.3), and ERF11 (ETHYLENE RESPONSE FACTOR 11) (Achard et al., 2003, 2007; Pierik et al., 2009; An et al., 2012; Luo et al., 2013; Marín-de la Rosa et al., 2014; Zhou et al., 2016). Therefore, we proposed that DELLA proteins may be involved in the sex differentiation of cucumber coupled with GA and ethylene in a collaborative regulation at the protein level, since the expressions of four DELLA homologs in cucumber, CsGAIP, CsGAI1, CsGAI2, and CsGAI3 (Zhang et al., 2014a), had no change after GA₃ treatment (Supplementary Tables S1 and S2).

In addition, GA and ethylene can also cooperatively regulate other aspects of plant growth and development. For example, GA promoted apical hook development of Arabidopsis, in part through transcriptional regulation of several genes in ethylene biosynthetic pathway mediated by DELLA proteins (Gallego-Bartolome et al., 2011). In this process, GA induced ethylene production, that was opposite to our results in cucumber sex expression. We speculated that this distinction may be due to different roles of hormones in various plant developmental processes. In fact, GA and ethylene showed similar functions in hook development, but they may play opposite roles in sex determination, thus, the regulatory mechanisms were different.

AGAMOUS, the C-class floral homeotic gene, belongs to the MADS-box family. In Arabidopsis, GA can induce the expression of AG gene (Yu et al., 2004). However, our data provided a novel point that an AG homolog CAG2 in cucumber was down-regulated upon GA₃ treatment (Figure 3; Table 4 and Supplementary Table S4). This distinction may be due to different abilities of these two genes to induce reproductive organ fate, in which, AG in Arabidopsis controls stamen and carpel development (Lohmann and Weigel, 2002), but CAG2 is particularly restricted to the carpel (Kater et al., 1998; Perl-Treves et al., 1998). Previous studies showed that CAG2 transcripts were not mediated by ethylene (Perl-Treves et al., 1998), hence, we speculated that GA probably suppressed pistil development through inhibiting the CAG2 expression, thereby allowing male flowers to develop, and ethylene was not involved in this process. Moreover, in our previous study, we demonstrated that transcript of CsGAMYB1 was upregulated by GA₃ treatment in male flower buds and silencing of CsGAMYB1 could suppress masculinization of cucumber, but the ethylene production and expression of F and M genes were not changed in the CsGAMYB1-RNAi lines (Zhang et al., 2014b). These observations suggested that GA can also regulate sex expression of cucumber via an ethylene-independent pathway (Figure 6, left panel). However, the relationship between CAG2 and CsGAMYB1 regulations remains obscure and that should be verified in further studies.

Besides, DELLA proteins can down-regulate the expression of floral homeotic genes, AP3, PI and AG, subsequently inhibit flower development in Arabidopsis (Yu et al., 2004). CsGAIP, a DELLA homolog in cucumber, may restrain staminate development through transcriptional repression of AP3 and PI in Arabidopsis (Zhang et al., 2014a). However, the expression levels of AP3 and PI were not changed after GA₃ treatment (Supplementary Tables S1 and S2), indicating that they did not participate in GA-regulated male tendency. Accordingly, there may be a regulatory relationship between DELLA proteins and floral homeotic gene CAG2 during GA-modulated sexual development in cucumber, however, the mechanism would be possibly different from that of Arabidopsis.

In summary, our data revealed a novel viewpoint that GA might control sex differentiation of cucumber via both ethylene-dependent and ethylene-independent pathways, and DELLA proteins were likely to be involved in both processes. However, this model was proposed by bioinformatics data, therefore, elucidation of the critical roles of DELLA proteins in flower development by cucumber transformation, and identification of the relationships among DELLAs and ethylene regulatory factors, GA-DELLA-CsGAMYB1 signaling and CAG2 gene, will shed new light on the molecular details of GA-regulated sex expression in cucumber.

Generally, typical unisexual flowers have two morphological types. The type I is unisexual by abortion. Initiation of stamen and pistil occurs in all flowers, followed by the developmental arrest in one or another organ. The type II is unisexual from inception. Only stamen or pistil is initiated and it does not go through a hermaphroditic stage (Lebel-Hardenack and Grant, 1997; Ainsworth, 2000; Mitchell and Diggle, 2005). Now, it is believed that the morphology of cucumber flowers belongs to the type I (Bai et al., 2004), but its evolutionary mechanism is largely unknown. Bai and Xu proposed a “miR initiative” hypothesis (Bai and Xu, 2013), where they speculated that unisexual cucumber flowers are evolved from a hermaphrodite ancestor. The first step in the evolutionary process might be the miRNA-regulated arrest of ovary development, and this predication is based on the altered expression of miRNAs, such as miR396a, 156b, 159a, 171b, and 166a, in male flowers (Bai et al., 2004; Sun et al., 2010). And this event leads to environment-dependent andromonoecy which has no progeny. Then, the M gene is recruited. On the one hand, the M gene promotes ethylene biosynthesis, resulting in the rescue of ovary development for seed set, because the ethylene might regulate the miRNA production. On the other hand, the M gene inhibits stamen development to avoid self-pollination and maintain cross-pollination. So, the monoecious genotype is generated through the cooption of the M gene. The andromonoecious genotype is produced by the loss-of-function m gene, which is regarded as a reverted point mutation. Further, the F gene is coopted and generate the gynoecious genotype (Sun et al., 2010).

Until now, a potential role of GA in unisexual flower evolution of cucumber has not been reported. But based on the possible function of M gene on evolutionary development of cucumber flower and the effect of GA on the transcript of M gene (Table 2), we speculated that GA might be involved in the process of cucumber flower evolution through interaction with ethylene. In addition, previous studies showed that GA signaling system can regulate anther development by modulation of miR159/GAMYB (Achard et al., 2004). In this pathway, miR159 acts as a post-transcriptional regulator of GAMYB transcript levels. GA relieves the DELLA repression of GAMYB, which is mediated by the GA activation of miR159. As mentioned above, miR159 is likely to participate in the arrest of ovary development in the evolutionary process of cucumber flower. And the GAMYB homolog CsGAMYB1 can regulate cucumber sex expression via an ethylene-independent pathway (Zhang et al., 2014b). These observations further revealed the possible involvement of GA in unisexual flower evolution of cucumber, but this process might be dependent on miR159 and GAMYB, and have no relationship with ethylene. Finally, it is worth noting that these viewpoints are built on the basis of the “miR initiative” hypothesis and needed to be tested in further work.