De novo transcriptome analysis reveals insights into dynamic homeostasis regulation of somatic embryogenesis in upland cotton (G. hirsutum L.)

Plant regeneration via somatic embryogenesis (SE) is the key step for genetic improvement of cotton (Gossypium hirsutum L.) through genetic engineering mediated by Agrobacteria, but the molecular mechanisms underlying SE in cotton is still unclear. Here, RNA-Sequencing was used to analyze the genes expressed during SE and their expression dynamics using RNAs isolated from non-embryogenic callus (NEC), embryogenic callus (EC) and somatic embryos (SEs). A total of 101, 670 unigenes were de novo assembled. The genes differentially expressed (DEGs) amongst NEC, EC and SEs were identified, annotated and classified. More DEGs were found between SEs and EC than between EC and NEC. A significant number of DEGs were related to hormone homeostasis, stress and ROS responses, and metabolism of polyamines. To confirm the expression dynamics of selected DEGs involved in various pathways, experiments were set up to investigate the effects of hormones (Indole-3-butytric acid, IBA; Kinetin, KT), polyamines, H2O2 and stresses on SE. Our results showed that exogenous application of IBA and KT positively regulated the development of EC and SEs, and that polyamines and H2O2 promoted the conversion of EC into SEs. Furthermore, we found that low and moderate stress is beneficial for proliferation of EC and SEs formation. Together, our global analysis of transcriptomic dynamics reveals that hormone homeostasis, polyamines, and stress response synergistically regulating SE in cotton. Electronic supplementary material The online version of this article (doi:10.1007/s11103-016-0511-6) contains supplementary material, which is available to authorized users.


Introduction
Crop losses due to multiple stresses have been predicted to be much greater than previously estimated in the coming future, especially due to climate change and industry development (Chen et al. 2015). Genetic improvement through breeding is the most efficient strategy for enhancing crop's stress tolerance. In addition to traditional breeding approaches, genetic transformation provides an alternative method to improve pest and disease resistance, as well as yield, of crops (Markram et al. 1997;Wu et al. 2008). For the majority of crops, Agrobacterium-mediated transformation is the method of choice for genetic modification. The prerequisite of this method is to have a highly efficient and (TEOSINTE B R ANCHED1-CYCLOIDEA-PCF15) and GhPIF4 (PHYTOCHOME INTE R ACTING FACTO R 4) (Min et al. 2015). Despite these reports, we still know little about the molecular mechanisms underlying cotton SE.
In this study, we used Xinluzao 33 (Gossypium hirsutum L.), one of the main cotton cultivars planted in Xinjiang Province, China, to investigate the transcriptome profiles of tissues from different stages during SE. Plant regeneration via SE using Xinluzao 33 is problematic because of the long time required for embryogenic callus (EC) induction and the low ratio of somatic embryo differentiation. Our transcriptome studies not only confirmed the importance of auxin and cytokinin homeostasis in cotton SE, but revealed the importance of genes involved in the polyamine metabolic pathway and stress response in cotton SE. Our results provided foundation and clues for establishing a reproducible and highly efficient plant regeneration system that can be used in a diverse of cotton genotypes for fast and precise genetic improvement by genetic engineering.

Plant materials and culture conditions
The method for cotton cultivar Xinluzao 33 SE used in this study has been previously described Sun et al. 2006). Non-embryogenic callus (NEC) (45 days), EC and SEs (a mixture of globular embryos, torpedo embryos and cotyledon embryos; SEs) (Supplementary Fig. 1) were sampled for RNA-Sequencing (RNA-seq), determinating the level of hormones, polyamines, H 2 O 2 , and qRT-PCR.

RNA extraction, library construction, and RNA-Sequencing
The RNA library construction, and RNA-Sequencing was performed by Illumina Genome Analyzer IIx at Beijing Genomics Institute (BGI)-Shenzhen, China. The NEC, EC and SEs samples were biologically replicated once, i.e. in total six libraries were used in RNA-Seq.

Analysis of sequencing data and identification of differentially expressed genes
The trimmed and filtered of raw reads and identification of differentially expressed genes in databases were based on a previous published paper (Xiang et al. 2010). Annotation of the differentially expressed unigenes was based on the Gossypium raimondii genome and A. thaliala genome (Kaul et al. 2000;Mayer et al. 1999;Wang et al. 2012). The genes mentioned in this study were also BLAST in other two recently reported G. hirsutum genome (TM-1) databases (Li repeatable system for tissue culture and plant regeneration via somatic embryogenesis (SE) or organogenesis.
Many studies have investigated physiological and biochemical changes during SE in various plant species with a focus on understanding the mechanisms of gene regulation related to SE. These efforts identified genes differentially expressed in somatic embryos (SEs), highlighted the pathways likely to be involved in SE and discovered molecular or protein markers for SE (Mantiri et al. 2008). Some of the identified genes have been experimentally demonstrated to play an important role in SE. For instance, Arabidopsis thaliana plants transformed with the SOMATIC EMBRYO-GENESIS R ECEPTOR -LIKE KINASE 1(AtSERK1) gene showed a marked increase in SE compared to wild-type cultures (Hecht et al. 2001). In Arabidopsis, transcription factors LEAFY COTYLEDON 1 (LEC1) is required for embryo development and LEC2 induces somatic embryo development in vegetative cells. Ectopic expression of LEC1 (Lotan et al. 1998) or LEC2 (Stone et al. 2001) caused spontaneous formation of SEs on intact plants or explants. AGL15 is another transcription factor that promotes SE in Arabidopsis (Harding et al. 2003). The BABY BOOM (BBM) gene can promote SEs in transgenic Arabidopsis and Brassica when ectopically overexpressed (Boutilier et al. 2002). Medicago truncatula SOMATIC EMBRYO RELATED FACTOR1 (MtSERF1) is essential for somatic embryo development and is induced by ethylene (Mantiri et al. 2008).
Cotton is one of the most important economic crops in the world, providing excellent natural fiber, and is also a source of oil. Cotton plantlets can be regenerated via SE using various combinations of plant growth-regulators, such as 2,4-Dichlorophenoxyacetic acid (2,4-D), indole-3-butyric acid (IBA), and naphthalene acetic acid (NAA) in combination with Kinetin (KT) (Sun et al. 2006). But cotton remains to be one of the notoriously recalcitrant plant species for plant regeneration via SE. So far, only few genotypes have been successfully used in gene transformation and genetic engineering. SE in cotton usually requires a very long culture time, and has a high frequency of abnormal embryos. To understand the molecular mechanisms underlying cotton SE and identify genes and critical pathways for cotton SE, suppression subtractive hybridization, macroarray and transcriptome during SE has been recently investigated (Xu et al. 2013;Yang et al. 2012;Zeng et al. 2006). Global transcriptome analyses suggested that auxin and cytokinin signaling pathways are critical in dedifferentiation of somatic cells and redifferentiation of SEs in cotton. In addition, stress-responsive genes and pathways were also found to be involved in SE (Yang et al. 2012). A more recent study has shown that the function of GhCKI (CASEIN KINASE I), a unique key regulatory factor strongly affecting cotton SE, is achieved by regulating auxin homeostasis through the network including other three genes, i.e. GhLEC1, GhTCP15 as the expression level of the genes analyzed. Primers used in qRT-PCR are listed in Supplementary Table 1.

Statistical analysis
Analysis of variance was performed using the SPSS16.0 statistical analysis package. Differences between means were compared by Fisher's least-significant-difference test at the 5 and 1 % probability levels.

RNA-sequencing and transcriptome de novo assembly
Generally, the process of cotton SE includes NEC induction, EC induction, differentiation of SEs, and plant regeneration ( Supplementary Fig. 1). To investigate the transcriptome dynamics during cotton SE, we performed RNA-Seq using mRNA extracted from NEC, EC and SEs.
In total, 168,196,388 75-bp paired-end raw reads were generated. For all three samples, the percentage of nucleotides with a quality score above 20 (Q20) was over 98.00 %, the percentage of reads containing N was lower than 0.02 % and the GC percentage was ~44 % (Supplementary Table 2). After filtering, 160,407,732 clean reads were retained and used in transcriptome de novo assembly using the Trinity program. As a result, 97,857, 87,997 and 84,176 unigenes were assembled in NEC, EC and SEs, respectively (Supplementary Table 3). Combining unigenes from all three samples, in total, 101,669 non-redundant unigenes were assembled. The majority of unigenes were 200-500 bp in length, and the numbers of unigenes and contigs decreased with the increase of the unigene length (Supplementary Table 3).

Functional annotation of unigenes
Unigenes were functionally annotated by searching against the NT, Swiss-Prot, KEGG, COG and GO databases using the Blast2GO algorithm with an E value ≤10 −5 . Of the total 101,669 non-redundant unigenes, 61,389, 43,333, 39,721, 25,923 and 50,884 had a hit in the NR, NT, Swiss-Prot, KEGG, COG and GO database, respectively. The predicted amino acid sequences of the unigenes were further aligned against the Pfam database (E value ≤10 −6 ) using the HMMER software. Altogether, 74,306 (73 %) unigenes had an annotated function. According to the COG (http://www. ncbi.nlm.nih.gov/COG) functional classification, the unigenes assembled were categorized into 25 groups, and were mostly covered by the inherent function, such as RNA processing, transcription, replication, recombination, repair and nucleotide or amino acid metabolism (Fig. 1a). Some unigenes were annotated to be involved in signal transduction, et al. 2015;Zhang et al. 2015), in order to confirm the gene sequences and ID.

Establishment of suspension cultures of uniform embryogenic callus
Suspension culture of uniform EC that was used in various physiological experiments outlined below was established following our previously described method .

Exogenous treatments of polyamines, H 2 O 2 , IBA, KT and simulated stresses
To determine the effects of polyamines, IBA (Indole-3-butytric acid) and KT on the conversion of EC into SEs, 60 µL of uniform embryogenic calli was cultured on somatic embryo induction medium supplemented with 1 mM of putrecine (Put), 1 mM H 2 O 2 , or 0.15 mg/L IBA + 0.15 mg/L KT. All reagents were purchased from Sigma Chemical Co. (St. Louis, MO). The simulated stress was performed with Phytagel (2.5, 4 g/L) or NaCl (0, 25, 50, 100 mM) in the MS medium. Each treatment was repeated three times biologically.

Determination of polyamine concentration
Polyamine (PA) concentrations were determined using the high performance liquid chromatography (HPLC) method following the procedure previously described . The assays were biologically repeated three times.

3′-diaminobenzidine (DAB)
Detection of hydrogen peroxide by 3, 3′-diaminobenzidine (DAB) staining in each sample was achieved using our previously described method . The assays were biologically repeated three times.

Determination of IAA, KT, GA and ABA
IAA, KT, GA and ABA were determined using an ELISA Assay Kit and corresponding antibodies produced by the Jiancheng Biochemistry Company (Nanjing, China) according to the manufacturer's instructions. The assays were biologically repeated three times.

RNA extraction and quantitative real time PCR (qRT-PCR)
RNA extraction, qRT-PCR and expression assay methods were described in our previous paper . The means of the three biological experiments were calculated Fig. 1 Function classification of all unigenes in COG and GO. a COG function classification. b GO classification of unigenes in embryogenic callus compared with non-embryogenic callus at three levels. c GO classification of unigenes in somatic embryos compared with embryogenic callus at three levels unigenes were only differentially expressed in NEC, 786 differentially expressed unigenes were only switched on/off EC, and 545 genes changed their expression level in SEs, which suggested that distinct spatial transcriptional profiles were present (Fig. 2a). Compared with NEC, EC had 5905 and 5983 unigenes down-and up-regulated, respectively. Compared with EC, SEs had 23,172 down-regulated unigenes and 19,361 up-regulated unigenes (Fig. 2b). These results suggest that the cotton transcriptome undergoes significantly dynamic changes during SE, particularly during the period from EC to SEs.

Pathways of differentially expressed unigenes based on KEGG
In total, there were 101,669 differentially expressed unigenes amongst NEC, EC and SEs, which were mainly involved in 127 pathways based on KEGG annotation. The pathways related to metabolites (e.g.purine, glycerophospholipid, starch and sucrose, ether lipid, amino sugar, nucleotide sugar, pyruvate) and biosynthesis of secondary metabolites (e.g. phenylpropanoid, stilbenoid, diarylheptanoid, gingerol and flavonoid), covered almost half of the differentially expressed unigenes. A considerable number of unigenes involved in DNA repair and replication, transcription, mRNA translation and protein folding, sorting and degradation were differentially regulated during the development of EC and SEs. Additionally, some stimulus pathways (e.g. plant hormone signal transduction and plant-pathogen interaction), ubiquitin mediated proteolysis, oxidative phosphorylation, ABC transporters and circadian rhythm were also found to be covered by differentially expressed unigenes. The top 30 pathways with the most number of differentially expressed unigenes are shown in Table 1.
cell cycle control and cell division. Genes related to plant growth, developmental progress and responses to stimuli were also well covered by the de novo assembly (Fig. 1a).
Based on GO (http://www.geneontology.org) classification, the majority of the unigenes under cellular component were involved in cell, cell part, organelle, membrane, organelle part, macromolecular complex, membrane part, cell junction and extracellular region. Molecular function categories were mostly clustered in binding, catalytic activity and transporter activity. The annotated unigenes clustered under biological process mainly involved in cellular process, metabolic process, single-organism process, response to stimulus and regulation of biological process. Comparing EC with SEs, there were more unigenes involved in locomotion at biological process level, more unigenes involved in macromolecular complex and nucleoid at cellular component level, and more unigenes involved in structural molecule activity at molecular function level (Fig. 1b, c).

Comparative analysis of differentially expressed unigenes in NEC, EC and SEs
Growth of NEC is in the process of dedifferentiation of vegetative growth, while formation of EC and development of SEs are regarded as in the process of dedifferentiation of reproductive growth. Spatial analysis was also performed on differentially expressed unigenes to ascertain the degree of overlap existing between the three different developmental processes during cotton SE. There were 94,020, 84,001 and 86,358 differentially expressed unigenes in NEC, EC and SEs (Fig. 2a). Among these, more than half (75.7 %) of the differentially expressed genes were present in all three developmental processes. Significant numbers of genes were present in one developmental process only: 13,505

Auxin and cytokinin played an important role in the transition from NEC to EC and from EC to SEs
A number of genes involved in auxin and cytokinin synthesis and signal transduction pathways were differentially expressed from NEC to EC and from EC to SEs ( Fig. 3; Table 2). For example, from NEC to SEs, the expression levels of AAO1-2 and CYP71A13, genes involved in IAA biosynthesis, were first down-regulated in EC and then slightly up-regulated in SEs. The gene encoding myrosinase, an IAA biosynthesis related enzyme, was continuously down-regulated from NEC to SEs (Fig. 3). For AUX1, AUX/IAA and ARF, genes involved in auxin signal transduction, the highest expression levels were observed in

Transcription factors (TFs) involved in cotton SE
We identified 302 and 112 differentially expressed TFs in EC (EC vs NEC) and SEs (SEs vs EC), respectively. Among these TFs, members of the following families were overrepresented: AP2/EREBP, C3H/C2H2/C2C2-Dof zinc finger proteins, MYB domain-containing proteins, and WRKY domain transcription factors ( Fig. 3 and Supplementary  Table 4). In addition, some TFs, such as SERK, LEC, WUS, BBM, CKI and AGL15, have been previously demonstrated to play a role in SE were also found to be differentially expressed during cotton SE ( Fig. 3; Table 2). These TFs are multifunctional regulators in both zygotic and SE with a role in hormone signaling and stress responses (Mahdavi-Darvari et al. 2015). Some of these TFs have been used as markers of totipotency in plant species. These results suggest that dynamic changes of the expression levels of plant  For homologous genes, the unigene with the highest absolute value of the fold change was selected. The positive and negative number represent up-and down-regulation, respectively. -Indicates unchanged the number of SEs (Fig. 4e, f). These results demonstrated that polyamines played a critical role in promoting the conversion of EC to SEs, Put and Spd might function similarly and are mutually substitutable.

Cotton SE is regulated by reactive oxygen species (ROS)
According to our transcriptomic data, some oxidative phosphorylation related genes, including NADH dehydrogenase, succinate dehydrogenase, cytochrome C oxidase and A/F/Vtype ATPase, were differentially expressed amongst NEC, EC and SEs ( Fig. 3; Table 2). The NADH dehydrogenase related genes, such as the ND, Ndufs and Ndh family genes were up-regulated in EC and SEs compared to NEC. The genes in succinate dehydrogenase pathway, such as SDHC, SDHA and SDHB were also differentially expressed during cotton SE. The cytochrome C oxidase encoding genes, including COXs, ISP, Cyt1/b, COR1 and QCRs showed an up-regulated pattern in EC and SEs. The ATPase (includes A, F and V-type) genes, such as beta, alpha, delta and OSCP had a relatively high expression level in EC ( Fig. 3; Table 2 and Supplementary Fig. 2). Differential expression of genes related to oxidative phosphorylation indicated that energy metabolism is very active in the transition from NEC to EC and development of SEs. Additionally, genes encoding the protective system enzymes, including CAT, SOD, POD and APX were upregulated in EC and SEs ( Fig. 3; Table 2 and Supplementary  Fig. 2). These genes are involved in the antioxidant system and might play a role in maintaining ROS homeostasis during the differential and development of EC and SEs. The concentration of endogenous H 2 O 2 was slightly higher in EC and SEs than in NEC as revealed by DAB staining (Fig. 4). Applying 1 mM H 2 O 2 in the medium significantly increased the number of SEs (Fig. 4j, k), suggesting that ROS played a crucial role in the development of EC and SEs during cotton SE.

Stress induced cotton SE
We found that many stress response related genes, such as CDPK, HSP90, SGT and genes encoding the protective system enzymes like NOX, SOD, CAT, POD, APX, SAMDC, PAO were differentially expressed during cotton SE ( Fig. 3; Table 2). To know the role of stresses in differential and development of SEs, the effects of different concentrations of NaCl (0-100 mM) and simulated drought stress on cotton SE were investigated by adding Phytagel (4 g/L) into the MS medium. After 4-week of NaCl treatment, EC growth was promoted by relatively lower concentrations of NaCl (≤50mM) but inhibited by the high concentration of NaCl (100 mM). Compared to the control, the 25 and the 50 mM NEC. From EC to SEs, the expression level of AUX1 was further decreased whereas the expression level of ARF was increased (Table 2; Supplementary Fig. 2). TRIT1, a gene related to cytokinin (zeatin) biosynthesis, was up-regulated in EC but down-regulated in SEs. CRE1 involved in cytokinin signal transduction was up-regulated in EC and remained high in SEs ( Fig. 3; Table 2; Supplementary Fig. 2). Another two genes involved in cytokinin signal transduction, A-ARR and B-ARR, showed a similar expression pattern as CRE1, although compared to EC, SEs had a decreased level of A-ARR ( Supplementary Fig. 2).
To monitor changes of auxin and cytokinin during SE, we measured the concentrations of endogenous IAA and KT in NEC, EC and SEs by ELLISA assay. The IAA concentration was the highest in NEC and decreased in EC and then slightly increased in SEs. The concentration of KT increased in EC and then significantly decreased in SEs. These results were consistent with the expression levels of AAO1-2, CYP71A13 and TRIT1 in these samples (Figs. 3,4a). SE can be controlled using different combinations of 2,4-D and KT, or IBA and KT (Fig. 3). With the reduction of 2,4-D concentration in the media, the callus became loose and granular in texture and yellow-green in color, a sign for efficient EC induction. The combination of IBA and KT significantly increased the fresh weight and the number of total embryos (Fig. 4b, c). Together, these results confirmed that auxin and cytokinin are critical regulators during cotton SE and that the role of IBA and KT in transition from NEC to SEs was supported by their dynamic changes of expression levels during the transition.

Polyamines promoted cotton SE
Based on KEGG analysis and GO annotation, the genes involved in polyamine synthesis or metabolism, including SAMDC, ADC, SPDS, SPMS and PAO, were found to be differentially expressed during cotton SE. All of these genes, particularly SAMDC and PAO, showed a higher expression level in EC compared with NEC. The expression levels of SAMDC, ADC and SPDS were decreased in SEs compared with EC, but the expression level of SPMS was further increased in SEs while that of PAO remained unchanged in SEs ( Fig. 3; Table 2 and Supplementary Fig. 2).
We measured the concentrations of the three common PAs (Put, Spd and Spm) in NEC, EC and SEs. The concentrations of Spd and total PA showed significant increase in EC and SEs compared to NEC. Although the levels of all three types of PAs and the total PAs did not show significant difference between NEC and EC, the levels of Spd as well as the total PAs were significantly increased in SEs whereas the levels of Put and Spm largely remained unchanged in SEs (Fig. 4d). Nevertheless, application of exogenous Put in the medium significantly increased the growth of EC and  Stone et al. 2001). BBM was similar to the AP2/ERF family genes and expressed preferentially in developing embryos (Boutilier et al. 2002). Over-expression of AGL15 enhanced production of secondary embryos from cultured zygotic embryos in A. thaliana (Harding et al. 2003). In our study, these referenced TF genes were found to be differentially expressed during cotton SE (Fig. 3; Table 2).
We found that 302 and 112 TFs were differentially expressed during the transition from NEC to EC and the development of SEs, respectively. According to Yang et al. (2012), of the total 466 differentially expressed TFs during cotton SE, 338 were found during the transition from NEC to EC and 342 were found during development of SEs. More differentially expressed TFs identified by Yang et al. (2012) during development of SEs could be due to the difference in the samples used. Yang et al. (2012) used globular embryos, torpedo embryos and cotyledon embryos individually, while we used a mixture of three types of embryos, which could have compromised the sensitivity of DEG identification.

Auxin and cytokinin are important regulators for cotton SE
Auxin and cytokinin are critical plant growth regulators (PGRs) for the induction of SE. Their regulatory role is probably achieved by establishing auxin and cytokinin gradients during the induction phase of SE, and is essential for initiating dedifferentiation and cell division of already differentiated cells before they can express embryogenic competence (Sharma et al. 2008;Su et al. 2011). Despite the absolute requirement for exogenous auxin to sustain growth of plant cells cultured in vitro, cultured plant cells produce substantial amounts of the native auxin, IAA. Comprehensive studies have been conducted on callogenesis from several cotton cultivars using various explants, and different combinations of growth regulators. The concentrations of IAA and KT were different in NEC, EC and SEs during cotton SE. Concentration of IAA decreased but concentration of KT increased in the EC stage compared to the NEC stage, consistent with the expression levels of genes involved in IAA and cytokinin biosynthesis (Table 2). This dynamics of auxin and cytokinin played a very important role in the transition from NEC to EC and the development of SEs (Yang et al. 2012;Xu et al. 2013). Previous studies had also shown that sharp changes in endogenous auxin and cytokinin levels might be one of the first steps leading to SE (Thomas et al. 2002). Our experimental results on the role of exogenously applied IBA and cytokinin on SE (Fig. 3) were consistent with previous studies, i.e. induction of NEC from hypocotyl requires addition of 2,4-D into the medium, but to induce EC and SEs, the concentration of 2,4-D in the medium was decreased to a lower level or replaced by combination of IBA and KT (Kumar et al. 2015; Sun treatments slightly and significantly increased the fresh weight of EC, respectively, but the 100 mM treatment had very significant adverse effect on EC growth (Supplementary Fig. 3), and no EC finally survived in the MS medium containing 100 mM of NaCl. Consequently, after 4-weeks of culture, significantly more SEs were observed in the medium containing 50 mM of NaCl than in the control medium (Fig. 3). Similarly, in the drought stress treatment, more SEs in the medium with 4 g/L Phytagel were observed than in the control medium (2.5 g/L) after 4-weeks of culture ( Fig. 3 and Supplementary Fig. 4). These results suggest that low and moderate stresses were beneficial for EC proliferation and differential of SEs.

Fatty acid biosynthesis and metabolism during cotton SE
Fatty acids, as a type of protective chemical in reproduction and seed preservation, play a very significant role in the evolution and developmental processes. Fatty acids in plants, as in all other organisms, are the major structural components of membrane phospholipids and triacylglycerol storage oils. During cotton SE, genes related to fatty acid biosynthesis, including FabD, FabH, accC and FatB, were up-regulated in EC compared to NEC, but down-regulated in SEs compared to EC (Table 2 and Supplementary Fig. 2). The genes involved in metabolism of fatty acids, such as ACOX1, ACOX3, fadA, fadB and fadI showed an expression pattern opposite to those in EC and SEs. These results indicated that fatty acids accumulated during the development of EC, a situation similar to initialization of plant seeds.

Transcription factors regulate SE
During the past three decades, numerous transcription factors and protein kinases, and a wide range of hormones have been shown to play a role in SE in different plant species (El Ouakfaoui et al. 2010;Mahdavi-Darvari et al. 2015). The SERK family genes, LEC, WUS and genes encoding germin-like proteins were implicated in the signal transduction pathway and transcriptional regulation during SE (aan den Toorn et al. 2015;Gaj et al. 2005;Podio et al. 2014;Zuo et al. 2002). LEC1, LEC2 and FUS3 are key genes that control SE process (Fambrini et al. 2006;Gaj et al. 2005). The capacity of SE is completely repressed in double (lec1 lec2, lec1 fus3, lec2 fus3) or triple (lec1 lec2 fus3) mutants in Arabidopsis thaliana (Gaj et al. 2005). The expression level of LEC2/FUS changes rapidly in response to auxin treatment (Stone et al. 2008), suggesting that LEC/FUS may be the upstream genes in the auxin signaling pathway (Gaj et al. that the key genes for biosynthesis of PAs were expressed at a higher level in EC (e.g. ADC and SAMDC), or in both EC and SEs (e.g. SPDS and SPMS; Table 2; Fig. 3). Both ADC mRNA and protein were localized in dividing cells of embryo meristems, probably required for mitosis (Vuosku et al. 2006). Exogenous application of Put and Spm enhanced the growth of embryogenic cultures of Araucaria angustifolia and significantly affected the endogenous concentrations of PAs, IAA and ABA in embryogenic tissues (Steiner et al. 2007). Our results also showed that exogenous application of Put enhanced the growth and the development of cotton SEs (Fig. 3), which could be achieved by PAs-mediated regulation of specific physiological processes important for cellular differentiation.

Possible roles of other hormones in cotton SE
In addition to auxin and cytokinin, genes related to biosynthesis and signaling of other hormones, including gibberellin, abscisic acid, ethylene, brassinosteroid, jasmonic acid and salicylic acid, were also differentially expressed in NEC, EC and SEs (Table 2 and Supplementary Fig. 5). The gibberellin synthesis and signal transduction genes, GA3/CYP701 and KAO, were down-regulated in EC but upregulated in SEs. Endogenous GA was at higher levels in EC and SEs than NEC (Fig. 4a). DELLA, a negative regulatory factor (Ueguchi-Tanaka et al. 2008), was up-regulated in EC and down-regulated in SEs. ABA synthesis and signal et al. 2006). However, the role of auxin and cytokinin is not achieved alone, but by interacting with other pathways and components involved in differential and development of SEs, such as their upstream TF regulators, hormone signaling pathways and ROS. Further studies are required to understand the interactive relationship of these components involved in the complex network regulating cotton SE.

Relationship between polyamines and SE
Polyamines have been previously linked to plant stress response and SE (Yoda et al. 2009). Many studies have demonstrated the importance of PAs for SE in several plant species (De-la-Peña et al. 2015;Feirer et al. 1984), probably by promoting cellular differentiation during SE (Gemperlová et al. 2009;Montague et al. 1978Montague et al. , 1979Niemenak et al. 2012). The concentrations of polyamines increased during the early stages of SE in conifers but decreased during the late stages (Gemperlová et al. 2009;Paul et al. 2009). Application of putrescine, spermidine and spermine significantly improved SE in Theobroma cacao L, Citrus sinensis and Hurst Ecotype (Malá et al. 2012;Silva et al. 2009;Wu et al. 2009). In our study, the levels of Spd significantly increased in SEs and the levels of Put and Spm almost did not change, a result reflected by the expression levels of the genes, such as SAMDC, ADC and SPDS, encoding enzymes catalyzing the production of PAs (Table 2; Fig. 3 and Supplementary  Fig. 2). Both our transcriptomic and qRT-PCR data showed  Fig. 5 A proposed model for a role of TFs, hormones, stresses, polyamines and ROS in cotton somatic embryogenesis. Italic words represent genes with down-and up-regulated genes highlighted in green and red color, respectively transduction related genes were up-regulated in EC and SEs, including ABA2, NCED, PP2C, ABF and PYR/PYL. ABA levels increased in EC and SEs, which was consistent with the expression levels of the genes for ABA synthesis (Fig. 4a). SE of carrot (Daucus carota L.) and Hevea brasiliensis required ABA as the source of growth regulator in medium (Dunstan et al. 1988;Etienne et al. 1993). Ethylene and polyamines shared the common precursor SAM, always played a competitive role in SE (Bai et al. 2013;Linkies and Leubner-Metzger 2012). Genes related to ethylene synthesis and signal transduction, CTH, cysK, GOT1, EIN2 and EIN3, were down-regulated in EC and SEs compared with NEC, and our data showed a competitive role between PAs and ethylene in cotton SE. JA synthesis and signal transduction related genes, such as JAR1, COI-1 and JAZ were up-regulated in SEs, but these genes were down-regulated in EC. SA related genes had an opposite expression pattern with JA in such process. Genes of CYP90C1/ROT3, CYP85A1/BR6OX1 and CYP734A1/BAS1, involved in BRs synthesis and signaling, were down-regulated in EC but upregulated in SEs. These hormones may play different roles in cotton SE, and the exact regulatory mechanisms will be discussed in our future work.
In conclusion, the transcriptome analysis unveiled the complex and dynamic nature of the gene network involved in cell dedifferentiation and redifferentiation during cotton SE, and identified the key regulators important for cotton SE. Our results further confirmed the main findings previously reported for cotton SE, mainly revealed the importance of genes involved in the polyamine metabolic pathways and stress response in cotton SE. Based on our results and previous reports, we proposed a working model for the potential gene network essential for successful SE (Fig. 5). In this model, we propose that TFs are the top layer players regulating reprogramming of the transcriptome required for cell dedifferentiation and redifferentiation, and that the hormone synthesis and signaling pathways, stress-response pathways and PA metabolic pathways function side-by-side and synergistically to promote embryogenic cell initiation and plantlet regeneration.