Analysis of differentially expressed genes during embryogenesis in ovary culture of cucumber (Cucumis sativus L.)

Background: Ovary culture has been a useful way to generate double haploid (DH) plant in cucumber (Cucumis sativus L.). However, the rate of embryo induction is low, and the probability for the induced embryo to grow into normal embryo is low as well. This is largely due to unknown of the mechanism of embryogenesis in cucumber. In this study, the differentially expressed genes during embryogenesis, including the early stages of embryo formation, embryo maturation and shoot formation, was investigated with transcriptomic technique to set up basis for a more efficient ovary culture technology Results: Cytological observations led to suggestions that cell enlargement is the symbol that gametophytes had switched to the sporophyte development pathway during the early embryogenesis stage. In this stage, RNA-seq revealed 3468 up-regulated genes, including hormone signal transduction genes, hormone response genes and stress-induced genes. The reported embryogenesis-related genes BBM, HSP90 and AGL were also actively expressed during this stage. The total of 480 genes that function in protein complex binding, microtubule binding, tetrapyrrole binding, tubulin binding and other microtubule activities were continuously up-regulated during the embryo maturation stage, indicating that the cytoskeleton structure was continuously being built and maintained by the action of microtubule-binding proteins and enzyme modification during embryo development. In the shoot formation stage, 1383 genes were up-regulated, which were mainly enriched in phenylpropanoid biosynthesis, plant hormone signal transduction, phenylalanine metabolism, and starch and sucrose metabolism. The shoot formation stage might be regulated by 6 transcription factors that contained a B3 domain, 9 genes in the AP2/ERF family and 2 genes encoded WUS homologous domain proteins. Conclusions: Findings from this study offer a valuable framework for explaining the transcriptional regulatory mechanism underlying embryogenesis in cucumber ovary culture.


Background
Gametophyte culture involved a stress-induced reprogramming of male or female gametophytes to develop into embryo-like structures, which can be directly regenerated into completely homozygous, doubled haploid (DH). The study of male gametophyte (microspore) embryogenesis began in the 1960s [1]. Microspore embryogenesis in vitro can be easily monitored and provides a convenient experimental platform for large-scale physiological and biochemical analyses [2]. However, the study of gynogenesis is not simple; the embryo sac is small and embedded in surrounding tissues, making early embryos difficult to observe and isolate. Despite these difficulties, to inducing gynogenesis in vitro, especially when the haploid induced by the male gametophyte is not successful or the induction rate is too low and when self-incompatibility, male-sterile and dioecious plants occurs, is an effective way to obtain haploid or double haploid plants [3].
In vitro gynogenesis is mainly used in the Cucurbitaceae, and previous studies have revealed that this technique is still imperfect and that the embryogenesis rate is low [4][5][6][7][8][9][10][11]. The study of the mechanism of embryogenesis has therefore been delayed, and the understanding of the mechanisms underlying gynogenesis induction is extremely limited.
In general, embryogenesis is a stress-induced phenomenon. Micropores are cultured and develop into embryos in vitro, and stress treatments involving cold or heat shock and hormones are used as trigger factors to induce gametocytes to follow the sporophyte development pathway [12]. The genes associated with the reprogramming phase and the early stage of embryogenesis were slowly characterized, such as AGL15-related proteins and the AGL15 gene, the BBM gene and HSPgene [13][14][15]. More recently, the use of functional genomics tools has made it possible to identify more genes associated with is a cross-pollinated plant with obvious heterosis. Therefore, almost all the cucumber varieties in production are hybrids; however, the purification process of the hybrids is slow and long. In the breeding of parents, it takes 6-8 years of artificial self-breeding to develop a stable inbred line. To accelerate the purification of cucumber parents and improve breeding efficiency, haploid gametophyte cultures can be used to induce embryoids, and homozygous double haploids can be obtained in 1-2 years. DH technology is a powerful tool to speed up plant breeding. However, additional research is needed to improve our understanding of the genes or the roles of genes involved in embryogenesis and the mechanism of haploid induction in embryo sacs [25].
This study was based on the highly efficient in vitro ovary culture technology system of cucumber (unpublished). We divided the process from acquiring embryogenic potential to plant regeneration into three stages: early embryo development, embryo maturation (from pre-embryos to cotyledon embryos) and the shoot formation stage. We evaluated the embryo morphology and transcriptomes at different developmental stages during embryogenesis in cucumbers via ovary culture. The metabolism and biological process of embryogenesis in the gynogenesis of cucumber are discussed, and several key genes regulating embryogenesis were identified.

Result
Morphological and cytological characterization of embryogenesis in cucumber The embryogenesis early stage, embryo maturation stage and shoot formation stage were observed to discern the cytological changes during embryogenesis by ovary culture in cucumber. Ovaries at 6 h before anthesis were collected and pretreated at 4 o C. The ovaries were sterilized, sliced, inoculated into the medium, incubated for 2 d at 33 o C in the dark and then transferred to 25 o C until plants formed. Interestingly, based on the morphological and cytological observations, critical developmental transitions in embryos were discovered at 2 (T1), 10 (T2), 20 (T3), 30 (T4) and 60 (T5) d of culture. Fresh unpollinated ovaries before culture were named T0 ( Fig. 1). At T0, fresh unpollinated ovaries were selected, in which the ovules were obvious (Fig. 1a) and the embryo sacs were mature (Fig. 1a1).At T1, heat shock stress induced the embryogenesis of cucumber (Fig. 1b). The most obvious sign that referred to the development of the embryo was the enlargement of the cells in the embryo sac (Fig. 1b1) respectively. At T2, the explants were cultured for 10 d, the ovules that successfully protruded from the surrounding tissues, and a pro-embryo had formed (Fig. 1c). At T3, after 20 d of culture, the total number of embryos increased continuously in one ovary slice (Fig. 1d). From T3 to T4 (culture for 30 d), the whole embryo had formed, including globular embryos, heart-shaped embryos, torpedo-shaped embryos and cotyledon-shaped embryos ( Fig. 1 g). At the same time, the color of the embryonic cells appeared green, indicating that the development of plants from the embryo would begin (Fig. 1e). At T5 (60 d after culture), shoots had formed (Fig. 1f).
At T1, ovule enlargement and cell enlargement could be clearly observed, suggesting that the transformation of the development pathway from the gametophyte to the sporophyte was induced by stress and that the stage between T0 and T1 was the key to inducing embryogenesis. After the switching of cell development, embryogenic potential was acquired. Mitosis then occurred continuously to form pro-embryos, globular embryos, heart-shaped embryos, torpedo-shaped embryos and cotyledon-shaped embryos, similar to zygote development, and this process continued until T4. After 6 d of culture, the mature embryos developed into shoots. Therefore, we divided the six time points of ovary culture into three stages: T0 to T1 (embryogenesis early stage), T1 to T4 (embryo maturation stage), and T4 to T5 (shoot formation stage). Transcriptome data were generated via Illumina II HiSeq TM 2000 sequencing, and the statistics of RNA-seq alignment were shown in Table S1.

Validation of differentially expressed genes
Validation of the Illumina sequencing data and the expression patterns of the DEGs revealed by RNA-seq was performed to examine the expression patterns of 10 DEGs, including 7 genes involved in plant hormone signal transduction and 3 plant-pathogen interaction genes. The results showed that the relative expression levels revealed by RNAseq and qRT-PCR were closely correlated (Pearson's r = 0.53, 0.77, 0.57, 0.74, 0.79, 0.69, 0.81, 0.81, 0.86, 0.62, 0.81). The fold changes in the qRT-PCR analysis were different from those in the RNA-seq analysis, which might be due to the difference in sensitivity between the qRT-PCR analysis and the RNA-seq technique. However, the qRT-PCR analysis showed that the up-and down-regulation trends of the differential gene expression were consistent with those of the RNA-seq analysis (Fig. 2).

Expression of embryogenesis-related genes in cucumber ovary culture
At the early stage of embryogenesis, after heat shock stress, ovule enlargement was obvious (Fig. 1b). Cytological observations also showed that one of the cells in the embryo sac expanded (Fig. 1b1), suggesting acquisition of embryogenic potential.
The differentially expressed genes (DEGs) at this stage (T0 vs T1) were identified as comprising 3468 up-regulated genes and 3065 down-regulated genes (Fig. 3). We performed a GO enrichment analysis of these genes, and the results showed that the proteins encoded by these genes were assigned to 4 biological processes, 3 molecular functions and 9 cellular components (Table S2). Different genes were involved in different biological processes, and the same gene was involved in a variety of biological processes.
The majority of DEGs were involved in the 'single-organism process' (GO:0044699) and the 'oxidation-reduction process' (GO:0055114). They were mainly distributed among the

Stress response proteins related to plant defense mechanisms
DEGs were significantly enriched in the plant-pathogen interaction pathway, which was related to plant defense mechanisms and complex physiological responses to heat shockinduced embryogenesis. In this regulatory pathway, the expression of Pti1 (Csa7M420160.1), HSP90 (Csa3M183950.1), EDS1 ( Csa1M006320.1) and other related protein genes were up-regulated. Heat shock protein 90 ( HSP90) was a protein that was up-regulated under heat stress. This gene had been highly conserved during evolution and had plant defense functions.
Response of genes involved in microtubule organization in cucumber ovary culture In the stage of embryo maturation, the gametophyte development pathway switched to the sporophyte development pathway, in which the cells enlarged and initiating mitosis, forming 2-cell, 4-cell and multicellular pro-embryos after heat stress (Fig. 1b1~b4); afterward, globular embryos, heart-shaped embryos, torpedo-shaped embryos and cotyledon embryos were formed (Fig. 1g). The DEGs analyses of T1 vs T2, T2 vs T3 and T3 vs T4 revealed 480 continuously expressed genes. Additionally, the numbers of uniquely expressed genes between T1 and T2, T2 and T3 and T3 and T4 were 2236, 371 and 1558, respectively (Fig. 6).  (Table S3). This findings indicated that, during certain stages, the cytoskeleton structure was continuously built and maintained by the action of microtubule binding proteins and enzyme modifications, providing the basis for positioning the various organelles and implementation functions, which ensured the orderly activities in various cells in time and space.
To determine the involvement of these differentially expressed genes in embryo development, we performed a pathway analysis to identify the potential target genes, and 17 significant pathways were obtained by mapping to the KEGG database ( Fig. 7): protein processing in the endoplasmic reticulum, phenylpropanoid biosynthesis, photosynthesisantenna proteins, limonene and pinene degradation, meiosis-yeast, the estrogen signaling pathway, cell cycle, cell cycle-yeast, diterpenoid biosynthesis, chloroalkane and chloroalkene degradation, carotenoid biosynthesis, progesterone-mediated oocyte maturation, glycerolipid metabolism, histidine metabolism, fatty acid degradation, antigen processing and presentation, and ascorbate and aldarate metabolism. Among them, the pathways of protein processing in the endoplasmic reticulum and phenylpropanoid biosynthesis were the most significant.
Expression of main oxidation-reduction and metabolic process-related genes in cucumber ovary culture In the stage of shoot formation, the embryos further differentiated into shoots (Fig. 1f). In total, 3320 genes were differentially expressed between T4 and T5, among which 1383 genes were up-regulated and 1837 were down-regulated in T5 (Fig. 3). Next, we performed an enrichment analysis of the genes using Gene Ontology (GO), which revealed 18 biological processes, 15 molecular functions and 5 cellular components (Table S4) The DEGs were subsequently annotated using the KEGG database to identify pathway enrichments. A variety of pathways were found to be significantly enriched (Fig. 8 In addition, the DEGs involved in starch and sugar metabolism were also significantly enriched and mainly encoding 3-βD glucosidase, T6Ps, GlgB, endoglucanase and alphatrehalase, which satisfied the requirements of embryo growth and development.

Discussion
Many previous experiments have shown that early embryogenesis of microspore cultures can be divided into three main phases: acquisition of embryogenic potential, initiation of cell division, and pattern formation [27]. In the present study, the whole process was investigated from acquisition of embryogenic potential to plant regeneration. We divided the ovary culture into three stages. In the early stage of embryogenesis, the acquisition of embryogenic potential by stress (e.g., low or high temperature and hormone induction) was observed together with the repression of gametophytic development, leading to the dedifferentiation of cells. In the embryo maturation stage, cell divisions gave rise to the formation of pro-embryos (cell clumps), globular embryos, heart-shaped embryos, torpedo-shaped embryos and cotyledon embryos (mature embryos). In the stage of shoot formation, mature embryos developed into shoots. There are active molecular events that regulate embryo development at different stages of development.
Embryogenesis-related genes are expressed in cucumber ovary culture The embryogenic ability and transformation of somatic cells were regulated by various plant hormones, such as auxin, abscisic acid, cytokinin and ethylene [28]. Moreover, the synthesis of auxin and ethylene increased significantly after ovary pollination and fertilization [29]. In normal growth, mitotic asymmetry of fertilized eggs was found in the gametophyte development pathway. Previous work using molecular and genetic approaches has identified auxin as a critical signal for the proper development of the asymmetric structure of the female gametophyte in Arabidopsis [30]. In in vitro culture, a certain stress condition is needed to block the development direction of the original gametophyte, followed by turning into the direction of the sporophyte and carrying on the symmetry splitting to eventually lead to embryogenesis [31]. In ovary culture of cucumber, heat shock pretreatment, silver nitrate, genotype and hormone combination factors could play key roles in embryo and callus production independently and simultaneously [32]. In the ovary culture of linseed, cultivar, combination of growth regulators, type of carbohydrates and their interaction significantly influenced callus induction and shoot formation frequency, and a relatively high shoot regeneration frequency was obtained when the medium was supplemented with TDZ and NAA [33].
These results implied that the expression of specific genes in embryogenesis was activated by exogenous hormone regulation, which laid the foundation for DNA replication in the cell division stage. Heat shock protein 90 (HSP90) is an evolutionarily conserved molecular chaperone induced by abiotic stress. HSP90 plays an important role in cell cycle control, genomic silencing, protein transduction and signal transduction [20]. The HSP90 chaperone is involved in maintaining phenotypic plasticity and developmental stability [34][35].
Previous studies have shown that HSP genes were specifically expressed in the spore stress-induction process of many crop species [36][37][38]. HSPs had been described as genes associated with the reprogramming phase and the early stages of embryogenesis [39].
After heat shock stress treatment, we found that 26 up-regulated genes encoded proteins and transcription factors related to heat shock, and HSP90 (Csa3M183950.1) was involved in the regulation of the hypersensitive response in the plant and pathogen interaction pathway. The expression of these genes might play a regulatory role in the process of embryogenesis.
The BBM gene was the first key gene to be isolated in the process of spore cell division; which was first expressed in zygotic embryogenesis and microspore embryogenesis [14,40]. We found that two BBM genes (Csa3M827310.1, Csa3M827320.1) were up-regulated at the beginning of embryogenesis and were subsequently down-regulated. In addition, AGL15 a member of the MADS family, was also considered a regulatory protein that acts at the start of cell division [13]. We found that the AGL gene (Csa3M258140.1) was upregulated in our study, indicating that all three of these genes play important roles in embryogenesis. We considered that BBM, HSP90 and AGL might be the critical genes involved in the induction of embryogenesis by ovary culture in cucumber.

Microtubule organization genes play an important role in the embryogenesis of cucumber ovary culture
There is a similar process between gametophyte embryogenesis and zygotic embryogenesis, and the only difference is that the former occurs in sex cells, while the latter occurs in fertilized egg cells. The process of gametophyte embryogenesis development is similar to that of zygote differentiation into embryos. Microtubules play a transport role in cells, participating in the construction of cell walls and promoting cell differentiation and division. Electron microscope-based observations of carrot culture showed that the appearance of microtubules was accompanied with the formation of somatic embryos [41]. The different patterns of microtubule organization in the cells of the mature embryo sac reflect their structural adaptations for their future function [42].
In our study, after embryogenic potential was acquired, the different embryo shapes At the stage of shoot formation, a large number of transcription factors were expressed.
The first constitutes transcription factors that contain a B3 domain [43]. These transcription factors can encode regulatory proteins involved in the embryonic development process, maintaining embryo development during late embryonic development [43,44]. In our study, total of 6 transcription factors(Csa2M359980.1, that contain a B3 domain were up-regulated.The second was the AIL gene from the AP2/ERF family, which is involved in key developmental processes throughout the whole plant life cycle. Some genes in the AP2/ERF gene family are expressed in many tissues and participate in many plant development processes, such as embryogenesis and shoot development [14,45). In addition, we found that the expression levels of 9 genes (Csa1M423190.1, Csa3M114480.1 , Csa3M652380.1, Csa3M114470.1 , Csa4M644740.1 ,   Csa5M175970.1, Csa5M608380.1, Csa6M496390.1 ,and CsaUNG031640) in theAP2/ERF family were up-regulated. Furthermore, WUS homologous domain proteins not only alter the cell fate of the shoot and flower meristem but also promote the development of somatic embryos into seedlings. The function of WUS proteins has no direct connection with the characteristics of the embryo but will alter the development state of the tissue by maintaining cells in an undifferentiated state in response to different stimuli [46]. The ability to transform the vegetative growth phase to the embryonic stage by WUS, as well as to eventually form somatic cells, indicates that this homologous domain protein also plays an important role in embryo maturation in addition to its role in the development of the meristem [47]. During this stage, two genes encoding WUS homologous domain proteins (Csa6M301060.1) and (Csa6M505860.1) were up-regulated. We considered that those homologous domain proteins might directly or indirectly regulate shoot formation.

Conclusion
Several studies have reported on the embryogenesis of ovary culture in cucumbers.
However, the mechanism that drives the process from embryogenic acquisition to the formation of embryos and plant regeneration is not well understood. In this study, we explored embryogenesis mechanism of ovary culture in cucumber. Inducing cucumber embryogenesis could be divided into three stages: early embryogenesis, embryo maturation (from proembryos to cotyledon embryos) and shoot formation (Fig. 9a).
The early stage of embryogenesis was the turning point for the formation of embryos, which had experienced dedifferentiation and loss of photosynthetic capacity, requiring the provision of exogenous nutrients and carbon sources such as plant hormones and sucrose in the media. Some physical factors, such as temperature, light quality, photoperiod and the presence of specific hormones, could affect the ability of the embryo sac to adapt to these conditions and survive the developmental transition. These dynamic changes were helpful for cell physiology reprogramming, metabolic alterations, revival of cell division and differentiation, morphogenesis and so on. Therefore, the cells of the embryo sac began to divide by dedifferentiation and started to form cell clumps (Fig. 1b1~b4).
Although there was no obvious change in appearance, the metabolism of some macromolecules in the cells changed markedly. A large number of hormone-related genes, cell protection-related genes, and some unique protein kinases and transcription factors (Csa6M147590.1 and Csa3M389850.1 ) were expressed only at the stage of early embryo development. The reported embryogenesis-related genes BBM , HSP90 and AGL were also actively expressed during this stage (Fig. 9b).
In the stage of embryo maturation, GO analysis showed continuously expressed genes participated in microtubule-based movement, movement of the cell or subcellular component processes, giving a good function of protein complex binding, microtubule binding, tetrapyrrole binding, tubulin binding and other microtubule activities, which are involved in protein processing in the endoplasmic reticulum and phenylpropanoid biosynthesis (Fig. 9b).

RNA isolation
The materials were harvested separately at the following time points: T0 (the ovules were cultured for 0 d, fresh unpollinated ovaries), T1 (the ovules were cultured for 2 d), T2 (the embryos were cultured for 10 d), T3 (the embryos were cultured for 20 d), T4 (the embryos were cultured for 30 d), and T5 (the shoots after culture for 6 d). These materials were stored at -80 o C for RNA extraction, and the remaining materials were maintained in culture to observe plantlet regeneration. Furthermore, all experimental procedures, such as culture medium replacement and sample collection, were performed under similar conditions to minimize possible circadian effects. Ovules at T0 and T1 were extracted by hand under a stereoscopic microscope, and embryo-like structures with integuments were selected directly at the time points from T2 to T5.
As described above, samples were collected from six time points of ovary culture for RNAseq analysis. For each sample, the ovules or embryo-like structures were immersed in liquid nitrogen in a mortar and ground into powder (three biological replicates per sample). Total RNA was isolated using TRIzol (Invitrogen), and DNase I (Fermentas) digestion was performed for 30 min at 25 o C to remove DNA, according to the manufacturer's instructions. The integrity and quality of the total RNA were checked using a NanoDrop 1000 spectrophotometer and formaldehyde-agarose gel electrophoresis. RNA was used only when the OD 260 : OD 280 ratio was above 1.8.

Transcriptome analysis
Transcriptomes were analyzed by an Illumina HiSeq TM 2000. The reads with more than 10% unknown bases and those of low quality were removed from the raw sequencing data using the software Illumina GA Pipeline (v1.6). The filtered reads were aligned to the cucumber Chinese long v2 genome using TopHat2 (http://tophat.cbcb.umd.edu/) software.
The reference genome and gene database information were obtained from the public website: http://cmb.bnu.edu.cn/Cucumis_sativus_v20/ [42]. Differential expression analysis was determined using DESeq (version 1.18.0). The resulting P values were adjusted by using Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P value of <0.05 and an absolute value of |log2fold change≥1| according to DESeq were considered differentially expressed. GO annotation was performed using Blast2GO software (GO association was performed by a BLASTX search against the NCBI NR database). GO enrichment analysis of differentially expressed genes (DEGs) was then performed via the BiNGO plugin for Cytoscape. Overrepresented GO terms were identified using a hypergeometric test with a significance threshold of 0.05 after the Benjamini-Hochberg FDR correction. KEGG enrichment analysis of the differentially expressed genes was performed using the KOBAS (2.0) [48] software. A 'P value ≤ 0.05' was used as the threshold to judge the significantly enriched pathways of the differentially expressed genes.
Real-time quantitative PCR Ten DEGs were selected from all the DEGs for analysis using real-time quantitative PCR (qRT-PCR). For this step, the total RNA was treated with the DNase I enzyme and subsequently converted to single-strand cDNA by applying the GoScript TM Reverse Transcription System (Promega) according to the manufacturer's protocol. The specific gene primers were designed using Primer 5.0 according to the cDNA sequences in Table   S5. Actin was used as the internal control to normalize small differences in template quantities. Quantitative real-time PCR was performed with GoTaq qPCR Master Mix

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The data sets supporting the results of this article are included within the article and its additional files.

Competing Interests
The authors declare that they have no conflict of interest.  Table S1. Statistics of RNA-Seq alignment Table S2. GO classification of common expressed genes in gynogenesis stage. Table S3. GO classification of common expressed genes from embryogenesis to the mature stage. Table S4. GO classification of common expressed genes in organogenesis stage.   The total number of up-regulated and down-regulated genes Red color represents up-regulated, and green represents down-regulated   Venn diagram showed differentially expressed genes in the stage from embryogenesis to the mature stage The DEG sets (T1 vs T2, T2 vs T3 and T3 vs T4) described in Fig. 3 were analyzed by using the Venn method. The numbers marked in the diagram indicated the number of continue express genes significantly expressed among the three DEGs sets.