Comparative transcriptome analysis reveals the regulatory networks of cytokinin in promoting the floral feminization in the oil plant Sapium sebiferum

Background Sapium sebiferum, whose seeds contain high level of fatty acids, has been considered as one of the most important oil plants. However, the high male to female flower ratio limited the seed yield improvement and its industrial potentials. Thus, the study of the sex determination in S. sebiferum is of significant importance in increasing the seed yield. Results In this study, we demonstrated that in S. sebiferum, cytokinin (CK) had strong feminization effects on the floral development. Exogenous application with 6-benzylaminopurine (6-BA) or thidiazuron (TDZ) significantly induced the development of female flowers and increased the fruit number. Interestingly, the feminization effects of cytokinin were also detected on the androecious genotype of S. sebiferum which only produce male flowers. To further investigate the mechanism underlying the role of cytokinin in the flower development and sex differentiation, we performed the comparative transcriptome analysis of the floral buds of the androecious plants subjected to 6-BA. The results showed that there were separately 129, 352 and 642 genes differentially expressed at 6 h, 12 h and 24 h after 6-BA treatment. Functional analysis of the differentially expressed genes (DEGs) showed that many genes are related to the hormonal biosynthesis and signaling, nutrients translocation and cell cycle. Moreover, there were twenty one flowering-related genes identified to be differentially regulated by 6-BA treatment. Specifically, the gynoecium development-related genes SPATULA (SPT), KANADI 2 (KAN2), JAGGED (JAG) and Cytochrome P450 78A9 (CYP79A9) were significantly up-regulated, whereas the expression of PISTILLATA (PI), TATA Box Associated Factor II 59 (TAFII59) and MYB Domain Protein 108 (MYB108) that were important for male organ development was down-regulated in response to 6-BA treatment, demonstrating that cytokinin could directly target the floral organ identity genes to regulate the flower sex. Conclusions Our work demonstrated that cytokinin is a potential regulator in female flower development in S. sebiferum. The transcriptome analysis of the floral sex transition from androecious to monoecious in response to cytokinin treatment on the androecious S. sebiferum provided valuable information related to the mechanism of sex determination in the perennial woody plants. Electronic supplementary material The online version of this article (10.1186/s12870-018-1314-5) contains supplementary material, which is available to authorized users.


Background
Sapium sebiferum, a member of the Euphorbiaceae family, is one of the most important perennial oil plants. It has attracted great attention due to its high oil production, pharmaceutical values, high ornamental values as landscape plant and its viability to grow in the marginal land (Peng et al. 2008;Yang et al., 2007). S. sebiferum has been cultivated as a valuable oil plant for more than 14 centuries, due to the high oil content of the seed coat and seed kernel, which can be used for the manufacture of industrial products, such as lubricant, soap and candles [1]. Aside from being a source for industrial purposes, the oil extracted from the S. sebiferum is also considered as a kind of renewable oil resources, which can be further exploited for the production of diesel fuel and fatty acid alkyl esters [2,3]. Nevertheless, low female to male flower ratio significantly limited the seed yield improvement and the industrial potentials of S. sebiferum. Most of the S. sebiferum plants are monoecious with separate male and female flowers on the same inflorescence, while some are androecious plants that only produce male flowers. There are over one hundred male flowers arrange in the narrow raceme-like inflorescence, while only with few female flowers growing at the base of the inflorescence. Thus, the manipulation of the flower sex differentiation by using the biochemical and genetic strategies was of significant importance in improving the seed yield in S. sebiferum.
Phytohormones are essential factors involved in the regulation of vegetative and reproductive growth in the plant kingdom. Phytohormones, such as cytokinin (CK), auxin, gibberellin (GA) and ethylene (ETH), were demonstrated to be involved in the regulation of flowering and sex determination in many species [4]. Generally, cytokinin, auxin and ethylene were considered as positive regulators in gynoecium development, whereas GA (gibberellin) promoted the development of androecious organs [5][6][7][8]. In Jatropha curcas and Plukenetia volubilis, exogenous application with cytokinin can increase the female flower number [9,10]. Nevertheless, the mechanisms how cytokinin regulated the sex determination remains largely unknown. In this study, we demonstrated that CK was a key regulator in sex determination in the oil plant S. sebiferum. CK treatment on the floral bud significantly induced the floral feminization and promoted the fruiting. As the seed yield improvement of S. sebiferum was greatly limited by the low female to male flower ratio, thus our findings will be of significant importance in the cultivation of S. sebiferum. Further, the feminization effect was also detected on the androecious genotype of S. sebiferum. To further investigate how CK regulated the transition of the floral bud from androecious to monecious, RNA-seq analysis was carried out on the floral buds of the androecious S. sebiferum at 6 h, 12 h and 24 h after 6-BA treatment. The present study not only provided an effective method in improving the seed yield during S. sebiferum cultivation, the transcriptome data also provided insights in understanding the mechanism of sex determination regulated by cytokinin in the woody plants.

Plant materials and growth conditions
In this experiment, three-year-old S. sebiferum trees, growing in the experimental field in Hefei Institutes of Physical Science (N 31°49′, E117°13′), Chinese Academy of Sciences, were used for the hormonal treatment. Trees were planted with 2.5 m × 2.5 m spacing, and fertilized every year before bud sprouting.

Phytohormone preparation and application
To make stock solutions (20 mM) of phytohormones, 6-BA, TDZ, trans-zeatin (tZ) or kinetin was dissolved in 0.2 M NaOH solution. The stock solutions were diluted with water to make different concentrations of working solutions containing 0.05% Tween-20 [11]. All working solutions for each treatment had the same solvent. The floral buds (one week after its formation) were treated with working solutions using a 100 ml plastic sprayer.
Sample collection and RNA isolation for transcriptome sequencing and quantitative real time PCR Floral buds of the androecious plants (Additional file 1: Figure S1) were treated with 6-BA (200 μM) and mock solution (containing the same solvent with 6-BA solution). The floral buds after removing the leaves (approximately 3 mm in length, Additional file 1: Figure S1) were collected at 6 h, 12 h and 24 h separately, and then immediately frozen in liquid nitrogen. Each treatment at each time point was prepared with six independent biological replicates (three for RNA seq, and three for quantitative real time PCR). The total RNA was extracted by using a Plant Total RNA Isolation Kit (Omega, China). RNA purity was checked by gel electrophoresis and NanoDrop 2000c spectrophotometer (Thermo Fisher, MA, USA). RNA concentration and integrity were separately determined using a Qubit 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA) and Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

Transcriptome sequencing and de novo transcriptome assembly
Sequencing libraries (each sample with three biological replicates) were prepared following the instructions of the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA). The literary quality was determined using the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed according to the instructions of cBot Cluster Generation System. Then the libraries were sequenced on the Illumina HiSeq 2000 platform at the Beijing Genomics Institute (BGI, Shenzhen, China) and the raw reads were generated in 90-bp paired-end format. The low-quality reads (< 20 Phred scores) were removed by Fastq_clean [12] and further assessed by FastQC [13]. The clean reads were assembled by using Trinity (version 2.0.6) [14]. The paired-end reads were then mapped to the de novo assemblies by using Bowtie (version 1.1.1) [15]. The abundance estimation was performed using Corset (version 1.03) [16].

Transcript annotation and gene expression analysis
All the assembled transcripts were searched against the NCBI non-redundant protein database (Nr), Swiss-Prot protein database and NCBI non-redundant nucleotide sequence (Nt) database. The GO functional classification was then carried out using Blast2GO program based on the Nr annotation [17]. The GO analysis was carried out by using the WEGO software [18]. The transcripts were aligned to the Clusters of Orthologous Groups of proteins (COG) database and Kyoto Encyclopedia of Genes and Genomes (KEGG) database [19,20]. To quantify the differentially expressed unigenes, the read counts were calculated by using the fragments per kilobase of transcript per million fragments mapped method (FPKM) [21].
Quantitative real time PCR (qPCR) analysis cDNA was synthesized according to the instructions of EasyScript One-Step gDNA removal and cDNA synthesis SuperMix (Transgen Biotechnology, Beijing, China). qPCR was performed using TaKaRa SYBR Premix Ex Taq™ II (TaKaRa Biotechnology, Dalian, China) on the Roche Light Cycler 96 System (Roche, Swiss). The PCR procedures were as follows: 95°C for 5 min; 45 cycles of 95°C for 10 s, 60°C for 10 s, and 72°C for 15 s; the melting curve analysis was from 65°C to 95°C. The information of the primers was listed as Additional file 2: Table S1.

Statistical analyses
Data were analyzed using the Statistical Product and Service Solution (SPSS) version 13.0 software. The significance between the hormone treated and control group was determined using Student's test.

Results
6-BA and TDZ had strong feminization effects on the floral development and promoted fruiting in S. sebiferum Based on the floral structure, S. sebiferum can be mainly divided into two main genotypes, S. sebiferum var. conferticarpa and S. sebiferum var. laxiarpa [22]. Most of S. sebiferum var. conferticarpa plants are monoecious, with few female flowers growing at the floral base, while some are androecious plants that only produce male flowers (Additional file 1: Figure S1). The results showed that exogenous 6-BA or TDZ treatment significantly induced the female flowers, which appeared at the position where the male flowers located (Fig. 1). The feminization effect was increased with increasing concentrations of 6-BA (Fig. 1a, b), or TDZ (Fig. 1c, d), reaching an average of 57.7 and 78.5 female flowers per inflorescence, separately by 6-BA and TDZ treatment. TDZ treatment had better effects in promoting the floral development, resulting in larger floral size and higher female flower number than that by 6-BA treatment (Fig. 1). However, over-dosed 6-BA or TDZ (e.g. 1 mM) was actually toxic to the floral bud (Fig. 1e). Approximately 20-30% of the 6-BA-or TDZ-induced female flowers could further develop into fruits (Fig. 2). It was worthy to mention here that although trans-zeatin (tZ) or kinetin were two other types of cytokinins. Exogenous application with these tZ or kinetin had no effects on the flower development or fruiting in S. sebiferum (Figs. 1f and 2d, e).

6-BA or TDZ treatment feminized the androecious genotype of S. sebiferum
In S. sebiferum, there are a small population of androecious plants that only produces male flowers and bares no fruits (Additional file 3: Figure S2). The effects of cytokinin on the floral development on the androecious genotype of S. sebiferum was also investigated in this study. The results showed that 6-BA or TDZ treatment can induce the female flower development on the male inflorescence of the androecious S. sebiferum (Fig. 3a, b). The induced female flowers, which had similar structures with that of the normal flowers ( Fig. 3c), were located at the position of the male flowers all over the whole inflorescence, and can further develop into fruits (Fig. 3d).
Transcriptome sequencing of the floral buds of the androecious S. sebiferum in response to 6-BA treatment The androecious genotype provided great advantages in investigating the mechanism of sex differentiation regulated by cytokinin in S. sebiferum. After 6-BA or TDZ treatment, abundant female flowers were induced in the "male inflorescence", thus we presumed that cytokinin was not only the key factor in regulating the sex differentiation in S. sebiferum, but also possibly participated in the formation of the androecious plants. The transcriptome analysis of the floral buds at the early developing stages of the androecious plants would be helpful for identification of key genes involved in the female flower formation and development in response to cytokinin in S. sebiferum.
In S. sebiferum, the floral buds were located at the terminal of each branchlet. One week after the floral bud formation (approximately 3 mm in length, Additional file 1: Figure S1), 6-BA (200 μM) and mock solutions were directly applied to the floral buds. In total of six groups of samples (6-BA-and mock-treated floral buds at 6 h, 12 h and 24 h, each group contains six biological replications) were collected for the following transcriptome sequencing and quantitative real-time PCR, aiming at identification of the early responsive genes that were involved in the sex determination regulated by cytokinin. Eighteen cDNA libraries were constructed and sequenced on the Illumina high-throughput sequencing platform. The detailed information of the transcriptome assembly can be found in Table 1 (Table 2). Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analysis were further carried out to generally describe the biological functions of these differentially expressed genes (DEGs) ( Fig. 5 and Additional file 3: Figure S2). The results showed that the DEGs were mainly classified into 19 categories by KEGG analysis and the predominant categories were mostly related to "Metabolism" (58%) and "Genetic information processing" (21%) (Fig. 5). Further, the predominant categories by GO analysis were "metabolic process", "cellular process", "cell", and "catalytic activity" (Additional file 6: Figure S3), which was closed correlated to the basic biological functions of cytokinin in promoting cell growth, replication and metabolism.
Previous studies suggested that cytokinin signal transduction occurs though a two-component phosphorelay system. Arabidopsis response regulators (ARRs) can be classified into two groups (type A and type B), both of which are important components in cytokinin signal transduction [26]. Our results showed that the expression of four type A ARR genes, ARR3 (CL7852.Contig3_All), ARR4 (CL7852.Contig4_All), ARR5 (CL1788.Contig4_All) and ARR9 (Unigene12447), was up-regulated at 6 h, 12 h or 24 h after 6-BA treatment (Fig. 7). As the expression of type B ARRs was not significantly affected after 6-BA treatment, it suggested that in S. sebiferum, cytokinin signal transduction in the floral bud might be dependent on the regulation of type A ARRs.

Discussion
Cytokinin has strong feminization effects on the flower development in both monoecious and androecious genotypes of S. sebiferum The regulation of the female to male flower ratio is critical for improvement of the seed yield in the monoecious plants. Various phytohormones were demonstrated to be important regulators not only in vegetative growth, but also in the flower initiation and sex differentiation. In plants, phytohormones such as ethylene [29], gibberellins [30] and auxins [5,31], have feminization effects in various species. Cytokinins are essential for numerous decisions throughout the plant developmental processes and adaptation to the biotic and abiotic environment [32]. In Jatropha curcas, Plukenetia volubilis and Luffa cylindrica, exogenous CK treatment can increase the number of female flowers [9,10,33]. Our present study in S. sebiferum were consistent with these findings that 6-BA or TDZ treatment has strong feminization effects on the floral development (Fig. 1). Although tZ and kinetin were two other types of cytokinins that can induce the division of plant cells [34], exogenous application with tZ and kinetin had no effects on the floral development in S. sebiferum (Fig. 1f ). We postulated that in S. sebiferum, the regulation of floral sex might be controlled by specific cytokinin types.
Our results further demonstrated that 6-BA or TDZ application can also effectively induce the female flowers on the androecious S. sebiferum (Fig. 3). These findings suggested that cytokinin was not only the key regulator in controlling the female flower development, but also possibly was the potential determinant that involved in the formation of the androecious plants in S. sebiferum. The transcriptome sequencing of the male floral buds subjected to 6-BA treatment was carried out and abundant cytokinin-responsive genes that were possibly involved in female flower induction were further identified. The floral buds of the androecious S. sebiferum seemed to be less sensitive to 6-BA treatment in the transcriptomic level, since the number of the differentially expressed genes in response to 6-BA treatment (especially at 6 h) was fewer (Fig. 4), compared with the transcriptome results from the other plant species [7,[35][36][37]. We postulated that the formation of the androecious plants in S. sebiferum was possibly due to the genetic mutations that caused the less-sensitivity to cytokinins. Nevertheless, exogenous 6-BA treatment can still effectively induce the female flower development on the androcious S. sebiferum (Fig. 3). Transcriptome analysis of the transition of the inflorescences from androecious to monoecious on the androecious plants will provide valuable information of the regulative gene network of the cytokinin-regulated sex differentiation in S. sebiferum.

Cytokinin regulated the cell cycle and nutrient translocation
The originally identified function of cytokinin was to promote the cell division [38]. In this study, ten cell cycle-related genes were identified in response to 6-BA treatment in the transcriptome database (Fig. 6a). Three D-type CYCD6 (CL3342.Contig16_All and CL3342.Con-tig11_All), CYCD3 (CL6841.Contig3_All) and one H-type CYCH1 (CL384.Contig6_All) cyclin genes that played important roles in the regulation of cell division [39], were up-regulated by 6-BA treatment (Fig. 6a). E2F Transcription Factor 3 (E2F3) was a key component in the cyclin D/retinoblastoma/E2F pathway, was also up-regulated at 24 h after 6-BA treatment [40]. Origin of Replication Complex 1B (ORC1B) that regulates the initiation of DNA replication [41], was up-regulated at 24 h after 6-BA treatment (Fig. 6a). NEDD1 encoding a WD40 repeat protein that plays a critical role in the cell mitosis in interactions with tubulin complex [42], was significantly up-regulated by 6-BA treatment (Fig. 6a). Two other transcription factors Block of Cell Proliferation 1 (BOP1) and NAC Domain Containing Protein 68 (NAC068) that involved in the cell cycle regulation [43,44], were also up-regulated in response to 6-BA (Fig.  6a). Additionally, several genes encoding the DNA/RNA polymerases, such as Nuclear RNA Polymerase C1 (NRPC1, CL2915.Contig1_All), Poly(ADP-Ribose) Polymerase 1 (PARP1, CL8049.Contig5_All), NRPB4 (CL4312.Contig3_All) and Poly(A) Polymerase 2 (PAPS2, CL14256.Contig6_All), were also up-regulated by 6-BA treatment (Fig. 6a). These results suggested that exogenous application with cytokinin could cause immediate cell replication and organ enlargement via regulating the cell-cycle related genes in the floral buds in S. sebiferum.
Further, a few genes that involved in the translocation of the nutrients, such as sugar, peptide, magnesium, ammonium, nucleotide, calcium, zinc and phosphate, were also differentially expressed in response to 6-BA treatment (Fig. 6b). The altered nutrients transportation to the floral bud after 6-BA treatment could be associated with the vigorous female flower induction in the floral bud. It has been demonstrated that sugars, such as sucrose, glucose and trehalose-6-phosphate were the initiation signals involved in the flower initiation in Arabidopsis [23]. Trehalose-6-Phosphate Synthases (TPS) were the key enzymes involved in trehalose-6-phosphate biosynthesis. The homologs of TPS1 in S. sebiferum were up-regulated, whereas TPS9 and TPS10 were down-regulated in response to 6-BA treatment (Fig. 6b). Several genes such as Early Response to Dehydration 6 (ERD6), Glucose-6-Phosphate Translocator 2 (GPT2), Nucleotide/Sugar Transporter (NST) that involved in the sugar translocation [45][46][47], were also differentially regulated by 6-BA treatment (Fig. 6b). These results suggested that the cytokinin-responsive cell cycle and cell metabolism-related genes could be involved in the regulation of flowering.
Cytokinin regulated the biosynthesis, transportation and signaling of other phytohormones in the regulation of sex determination in S. sebiferum Phytohormones, such as auxin, cytokinin (CK), abscisic acid (ABA), ethylene (ETH), brassinosteroid (BR), gibberellin (GA), jasmonic acid (JA) and salicylic acid (SA), coordinate plant growth and development by modulating various cellular changes in response to the environmental or intrinsic signals [4]. The transcriptome data showed that many genes that were directly related to the biosynthesis, transportation and signaling of other phytohormones, were differentially regulated by 6-BA treatment in the floral buds of the androecious S. sebiferum. Specifically, ten genes were identified in auxin biosynthesis and signaling, eight genes in ABA biosynthesis and signaling, four genes in ethylene signaling, three gene in GA signaling, and three genes in JA signaling (Fig. 7). Previous studies demonstrated that cytokinin and auxin mutually regulated their biosynthesis and signaling in the regulation of plant growth and development [48]. Exogenous application with cytokinin in young shoot up-regulated the expression of auxin biosynthesis genes, such as YUCCA (YUC) [49]. In this study, two YUC genes YUC2 and YUC4 were found to be significantly up-regulated at 6 h, 12 h and 24 h in response to 6-BA treatment (Fig. 7). Ring Domain Ligase 2 (RGLG2) was previously reported to be involved in the regulation of auxin metabolism. Mutation of RGLG2 could cause significant decrease of auxin content in both leaf and root in Arabidopsis [50]. The expression of RGLG2 was increased in response to 6-BA treatment in S. sebiferum (Fig. 7). The results also showed that the expression of Dioxygenase for Auxin Oxidation 1 (DAO1), which encoded a major IAA oxidase in plants [28], was found to be down-regulated at 12 h, whereas promoted at 24 h (Fig. 7). These results suggested that the auxin biosynthesis was transiently promoted, at least in a short time after 6-BA treatment at the floral buds. Additionally, several auxin signaling and transportation-related genes, such as Auxin Response Factor 2 (ARF2), Indole-3-Acetic Acid Inducible 13 (IAA13), Indole-3-Acetic Acid Inducible 26 (IAA26), GH3.9, Auxin Resistant 3 (AXR3) and Protein Phosphatase 2A (PP2A) were found to be mostly up-regulated in response to 6-BA treatment (Fig. 7). It has been demonstrated that auxin was important for the female reproductive organ development in Arabidopsis [51], and also exhibited a strong feminization effect on Cannabis sativus and Opuntia stenopetala [5,31]. The expression of auxin signaling genes was also positively correlated with the transition to female flowers on the monoecious Ricinus communis [52]. Our results were in accord with these findings that auxin could also play important role in the female flower induction regulated by cytokinin in S. sebiferum.
In cucumber, ABA did not change the sex expression, nevertheless it promoted the female tendency of gynoecious plants [53]. XERICO and Nine-cis-epoxycarotenoid Dioxygenase 3 (NCED3) are two important regulators in ABA biosynthesis [54,55]. The transcriptome data showed that these two genes were significantly down-regulated by 6-BA treatment (Fig. 7). The ABA signaling-related genes, such as Ring Domain Ligase 2 (RGLG2), ATMYB44 (MYB44), Growth Regulating Factor 3 (GRF3) and Abscisic Acid Responsive Elements-binding Factor 3 (ABF3) were down-regulated by 6-BA treatment (Fig. 7). Previous studies showed that ABA and cytokinin antagonistically regulated many aspects of plant growth and development [56,57]. In S. sebiferum, the cytokinin responsiveness of the ABA biosynthesis and signaling-related genes suggested an antagonistic interactions between CK and ABA in the regulation of the floral development.
Cytokinin also interacted with ETH, GA and JA in the regulation of flowering [58]. PUCHI, a key component in ethylene-activated signaling pathway, was required for morphogenesis in the early floral meristem growth [59]. The transcript abundance of PUCHI in the floral bud was up-regulated in response to 6-BA treatment (Fig. 7). The S. sebiferum orthologs of three additional ethylene signaling-related genes, Ethylene Response Factor 74 (ERF74), Ethylene and Salt Inducible 3 (ESE3) and Transcription Factor IID (TFIID) were also differentially regulated by 6-BA treatment (Fig. 7), suggesting a potential interactions between CK and ETH in the regulation of sex determination in S. sebiferum. Among these differentially expressed genes after 6-BA treatment, three JA signaling genes ATMYC2 (MYC2), Jasmonate Insensitive 1 (JAI1), BHLH14 and three GA signaling genes GA-Stimulated Arabidopsis 1 (GASA1-1 and GASA1-2) and Glabrous Inflorescence Stems 2 (GIS2) were differentially regulated in response to 6-BA treatment (Fig. 7). Conclusively, these results implied that cytokinin interacted with auxin, ABA, ethylene, JA and GA in sex determination in S. sebiferum (Fig. 9).
Cytokinin regulated the expression of gynoecium and androecium identity-related genes 6-BA treatment profoundly induced the formation of female flowers in the androecious S. sebiferum (Fig. 3). The differential expression of ABCE floral organ-identity genes has been implicated in the regulation of male and female flower development [60,61]. Thus it suggested that cytokinin could directly target the floral organ-identity genes in the regulation of sex determination. A total of 21 flowering-related genes were identified in response to 6-BA treatment at the floral bud of the androecious plants (Fig. 8). Several early floral meristem growth-related genes, such as PUCHI, Target of Early Activation Tagged 1 (TOE1), Ethylene Response DNA Binding Factor 4 (EDF4) and WIGGUM (WIG), were differentially regulated in response to 6-BA treatment (Fig. 8), suggesting cytokinin was directly involved in the regulation of the early floral meristem development in S. sebiferum. Among the DEGs, a few genes that were closely related to the gynoecium and androecium development, were found to be differentially regulated in response to 6-BA treatment. SPATULA (SPT), a bHLH transcription factor, was involved in gynoecium development [62], was up-regulated at 12 h after 6-BA treatment (Fig. 8). KANADI 2 (KAN2), which encodes a Fig. 9 Role of cytokinins in the regulation of sex determination in S. sebiferum