Comparative gene expression profile analysis of ovules provides insights into Jatropha curcas L. ovule development

Jatropha curcas, an economically important biofuel feedstock with oil-rich seeds, has attracted considerable attention among researchers in recent years. Nevertheless, valuable information on the yield component of this plant, particularly regarding ovule development, remains scarce. In this study, transcriptome profiles of anther and ovule development were established to investigate the ovule development mechanism of J. curcas. In total, 64,325 unigenes with annotation were obtained, and 1723 differentially expressed genes (DEGs) were identified between different stages. The DEG analysis showed the participation of five transcription factor families (bHLH, WRKY, MYB, NAC and ERF), five hormone signaling pathways (auxin, gibberellic acid (GA), cytokinin, brassinosteroids (BR) and jasmonic acid (JA)), five MADS-box genes (AGAMOUS-2, AGAMOUS-1, AGL1, AGL11, and AGL14), SUP and SLK3 in ovule development. The role of GA and JA in ovule development was evident with increases in flower buds during ovule development: GA was increased approximately twofold, and JA was increased approximately sevenfold. In addition, the expression pattern analysis using qRT-PCR revealed that CRABS CLAW and AGAMOUS-2 were also involved in ovule development. The upregulation of BR signaling genes during ovule development might have been regulated by other phytohormone signaling pathways through crosstalk. This study provides a valuable framework for investigating the regulatory networks of ovule development in J. curcas.

. The depth, coverage and homogenization of sequencing were high, which suggested that these results could reflect the actual expression of genes, and the data were thus suitable for further analysis (    (Table S4). In total, 1723 genes that were differentially expressed during ovule development in female flower buds were screened out and are shown in Fig. 4. Of these differentially expressed genes (DEGs), 9 genes were downregulated and 12 genes were upregulated throughout the entire process of ovule development; 579 genes were downregulated and 1036 genes were upregulated in one stage; 42 genes were downregulated and 100 genes were upregulated at two stages (Fig. 5).
Gene expression patterns of selected genes. In total, 16 DEGs were selected for qRT-PCR analysis to validate the expression profiles obtained by digital gene expression analysis, and the results showed that the expression patterns of these genes were consistent with the results obtained by digital gene expression analysis, except AUX22 (c11413_g1) (Fig. 6), indicating that transcriptome data in this study were reliable.
KEGG pathway enrichment analysis of differentially expressed genes. After KEGG analysis, the significantly enriched pathways were selected for further analysis. The results indicated that sugar metabolism and protein biosynthesis/processing were respectively enhanced during the maturation of ES and the formation of MES, and the metabolisms of auxin and JA were also upregulated during the formation of MMC and the maturation of EC, but the metabolism of SA declined during ovule development. Among the significantly enriched pathways (Fig. 7), plant hormone signal transduction was downregulated in pairwise JCFI vs JCFII and JCFII vs JCFIII but showed a trend to upregulate in pairwise JCFIII vs JCFIV; starch and sucrose metabolism was upregulated in pairwise JCFIII vs JCFIV; protein processing in endoplasmic reticulum was upregulated in pairwise JCFIII vs JCFII; tryptophan metabolism associated with auxin metabolism was upregulated in pairwise JCFI vs JCFII; alpha-linolenic acid metabolism associated with JA metabolism was upregulated in pairwise JCFIII vs JCFIV; phenylalanine metabolism associated with SA metabolism was downregulated in pairwise JCFI vs JCFII.
Differentially expressed transcription factor genes during ovule development. A total of 156 transcription factor (TF) genes were differentially expressed during ovule development, of which the top six TF families were bHLH, WRKY, MYB, ethylene-responsive transcription factors (ERF), NAC and TCP: 26 from bHLH, 21 from WRKY, 20 from MYB, 18 from ERF, 13 from NAC, and 8 from TCP. Among the DEGs of these six TF families, the upregulated TF genes included 17 bHLH, 11 WRKY, 12 ERF, 10 NAC, 10 MYB and 2 TCP (Fig. 8). Of these 10 MYB genes, a R2R3-type MYB gene (c28570_g1) was confirmed to be upregulated with the development of the ovule by qRT-PCR analysis (Fig. 9). These results suggested that bHLHs, WRKYs, ERFs, NACs and MYBs would extensively participate in ovule development.
Prediction of gene regulating ovule development. In total, 9 genes annotated as ovule development-related genes were upregulated during ovule formation as suggested by transcriptome analysis, including 5 MADS-box protein genes (AGAMOUS-2 (AG-2), AG-1, AGL1/SHP1, AGL11/STK and AGL14), YABBY 5, UPERMAN (SUP), SLK3 and Argonaute 2 (Table 1). These five MADS-box proteins were then subjected to sequence alignment with their homologs in the protein database UniProt. The results showed that these five MADS-box proteins have a MADS-MEF2-LIKE and a K-box domain and belong to the MADS-MEF2-LIKE subfamily of MADS (Fig. 9). The MADS-MEF2-LIKE domain showed high homology between these five MADS-box proteins and their homologs (Fig. 9). Of these 9 genes, AG-2 was upregulated throughout ovule development, as suggested by both transcriptome analysis and qRT-PCR analysis ( Fig. 10; Table 1); AGL1, AGL14, AG-1 and SUP were upregulated in pairwise JCFI vs JCFII; SLK3 was upregulated in pairwise JCFIV vs JCFIII; AGL11/STK was continuously upregulated from JCFI to JCFIII; YABBY 5 was continuously upregulated from JCFII to JCFIV; argonaute 2 was upregulated in pairwise JCFII vs JCFIV and JCFI vs JCFIV. On the other hand, CRABS CLAW (CRC), annotated as an ovule development-related gene, was upregulated from JCFI to JCFIV, as suggested by qRT-PCR analysis (Fig. 10).
Phytohormone signaling responsive genes associated with ovule development. In all, 71 phytohormone signaling-responsive DEGs were identified, 61 of which were screened out by transcriptome analysis, 9 by qRT-PCR analysis, and IAA1 by both qRT-PCR and transcriptome analysis (Figs 10, 11 www.nature.com/scientificreports www.nature.com/scientificreports/ of IAA1 was first downregulated at JCFII and then significantly upregulated during the following development stages (Fig. 10B,L).
Differentially expressed genes related to phytohormone metabolism. The DEGs related to phytohormone metabolism were screened out as follows: auxin metabolism-related DEGs were screened out from the DEGs enriched in the tryptophan metabolism pathway, CK metabolism-related DEGs from the DEGs enriched in the enriched zeatin biosynthesis pathway, gibberellin metabolism-related DEGs from the DEGs enriched in the diterpenoid biosynthesis pathway, BR metabolism-related DEGs from the DEGs enriched in metabolism pathway, JA metabolism-related DEGs from the DEGs enriched in the alpha-linolenic acid metabolism pathway, and SA biosynthesis-related DEGs from the DEGs enriched in phenylalanine biosynthesis ( Table 2).

Discussion
The embryo sac of J. curcas belongs to the polygonum type. A macrospore mother cell (MMC) only produces one functional macrospore that then develops into an embryo sac with 8 nuclei after three consecutive mitoses. Of these 8 nuclei, two move to the center of the embryo sac and form a polar, and the remaining six form an egg, two synergids, and three antipodal cells, respectively (Figs 1 and 2). A transcriptome covering the entire development process of sporogenesis and gametogenesis in both male and female flowers was constructed using Illumina sequencing.
Overall, plant hormone signal transduction declined during ovule development but showed a trend toward enhancement during the maturation of ES. Protein biosynthesis and processing were enhanced during the formation of MES. Sugar metabolism seemed to be enhanced during the development process of the ovule. A strong relationship exists between flower development and carbohydrates [12][13][14][15] . In particular, the concentration of starch in the ovule at specific steps of development is closely correlated with fertility 16 . Transcription factors participating in ovule development. Five TF families, such as bHLH, WRKY, MYB, NAC and ERF, would extensively contribute to ovule development, as suggested by the upregulation of 17 bHLH, 11 WRKY, 12 ERF, 10 NAC and 10 MYB during ovule development. The members of these five TF families are thought to be involved in ovule development. bHLHs play important roles in carpel, stigma and anther development and phytochrome signaling [17][18][19] . MYBs can regulate the development of anther, petal, and embryogenesis [20][21][22] . R2R3-MYBs are the largest group of plant MYB factors, and several R2R3-MYBs have been documented to function in ovule development. FOURLIPS and MYB88 are thought to regulate megasporogenesis 23 , and MYB98 is specifically expressed in synergid cells 24 . For WRKY, WRKY71 can accelerate flowering via the direct www.nature.com/scientificreports www.nature.com/scientificreports/ activation of FT and LEAFY 25 , and WRKY34 and WRKY2 are required for male gametogenesis in Arabidopsis 26 . Several members of the AP2/ERF family are thought to be involved in microspore, somatic and zygotic embryogenesis 27 . Members of the NAC family, such as CUC genes, contribute to ovule primordial development 28,29 .

MADS-BOX, SUP
and SLK3 genes participate in ovule development. In the present study, five MADS-BOX genes (AG-1, AG-2, AGL1, AGL11 and AGL14), SUP and SLK3 are expected to function in the development of the ovule in J. curcas. AG has been reported to determine stamen, carpel and ovule identity [30][31][32] . In this study, AG-1 functioned during the stage of occurrence of MMC and AG-2 throughout the entire development www.nature.com/scientificreports www.nature.com/scientificreports/ of the ovule. SHP1 (AGL1) and STK (AGL11) are known as ovule identity genes that control ovule cell fate and regulate sporophyte and gametophyte development [33][34][35] . SHP1 (AGL1) is also reported to function in regulating cell divisions in the ovule and in promoting stigma, style and medial tissue development 36,37 . Our results show that AGL1 and AGL11 play roles during the stage of MMC, controlling ovule cell fate, and regulating sporophyte development and cell divisions in the ovule. However, AGL11 is also predicted to function during the formation of MES and regulate gametophyte development. AGL14 can control auxin transport via PIN transcriptional regulation 38 , and the PIN protein-dependent auxin gradient is responsible for pattern formation during early embryogenesis, organogenesis and ovule development 39 . Our results suggest that AGL14 functions during the    www.nature.com/scientificreports www.nature.com/scientificreports/ GA plays important roles at the stage from MMC in ovule development. GA is known to play various roles in plant reproductive development, and different species respond differently to GA 48 . The GA signaling pathway is involved in the development of stamens and pollen in Arabidopsis and tomato [49][50][51] . GA is a negative regulator in the development of stamen in maize 52,53 but a positive regulator in Arabidopsis and rice 54 . GA is thought to play essential and complicated roles in floral organ development in J. curcas. In J. curcas, GA treatment increased both female and male flowers in the studies of Pan et al. 10 and Gayakvad et al. 55 promoted the development of stamens to produce bisexual flowers in research of Pi et al. 56 and increased the number of female flowers in the research of Makwana et al. 57 . In gynoecious plants, GA treatment represses pistil development in female flowers to produce neutral flowers but does not resume the development of stamens 58 . In the present study, GA was predicted to promote the development of macrospores and embryo sacs in monoecious female flowers. The GA concentration in flower buds was increased at the stages from the occurrence of MMC to the maturation of ES, which was due to the upregulation of KO and KAO1. KO and KAO1 catalyzed the reaction to convert ent kaurene to GA 12 during GA biosynthesis 59 . Furthermore, the upregulation of several GA-responsive genes (such as GASL7 GASA3 GID1b and GAI) during the stages from the formation of MMC to the maturation of SE may be the result of the increase in GA. In cotton, GASL7 is predominantly expressed in cotyledons, and GASA3 is expressed in fiber 60 . In Arabidopsis thaliana, GASA3 accumulates in root, meristems, and flower seeds 61 . In Arabidopsis thaliana, 3 GA hormone receptors, such as GID1a, GID1b, and GID1c, have partially specialized functions in both proteolytic and nonproteolytic GA signaling. GID1b plays a stronger role in nonproteolytic GA signaling 62 . As a member of DELLA, GAI is a negative regulator of GA signaling 63 .

JA plays important roles in the stage of maturation of EC.
JA is thought to play critical roles in regulating reproductive development in plants, such as sex determination and stamen development in maize, and the development of inflorescences and flowers in rice and Arabidopsis 64 . In J. curcas, JA biosynthesis is lower in gynoecious than in monoecious florescence, indicating that the decrease of JA biosynthesis would likely contribute to the abortion of male flowers in gynoecious plants 60 . However, our results suggest that JA contributes to the maturation of EC in females; JA levels are increased during EC maturation, which is the result of the enhancement of JA biosynthesis, suggesting the upregulation of four JA biosynthesis genes, such as 3-ketoacyl-CoA thiolase 2, linoleate 13S-lipoxygenase 3-1, allene oxide synthase, and 12-oxophytodienoate reductase 3 59 . On the other hand, several JA-responsive genes, such as MYC2, JAZ, TIFY 6B and TIFY 10B, were upregulated at the maturation of EC, supporting our results. TIFY 10B and TIFY 6B belong to the JAZ subfamilies of the TIFY family and could be induced by JA 65,66 . In Arabidopsis, JAZ1/TIFY10b can work as a transcriptional repressor in the models of JA signaling 67,68 . Activation of JA response genes is mediated in part by MYC2 [69][70][71] ; in the absence of JA, MYC2's function as a transcriptional activator is repressed by members of the JAZ family proteins, and elevated levels of JA induce the release of MYC2 from this repression by increasing the rate of degradation of JAZs 70,72 . Furthermore, MYC2 is a "master switch" in the crosstalk between JA and the other hormone signaling pathways 73 . Auxin may regulate the occurrence of MMC in ovule development. Auxin is thought to regulate ovule development 74,75 , possibly by determining the cell identity of the female gametophyte cell 39 , and the gradient of auxin is responsible for pattern formation during ovule development 39,76 . Similarly, our results also suggest that auxin contributes to ovule development. The auxin biosynthesis in flower buds is expected to be enhanced during the occurrence of MMC, as suggested by the upregulation of YUCCA8, ALDH 3H1 and ALDH 2B4-1. In Arabidopsis thaliana, the YUCCA gene encodes a flavin monooxygenase-like enzyme that appears to oxidize tryptamine to N-hydroxytryptamine 77,78 . The ALDH family is thought to function in IAA biosynthesis by catalyzing the NAD-dependent formation of IAA from indole-3-acetaldehyde 79 . The expression levels of SAURs (22D, 10A5 and two 15A), two GH3s (GH3.17 and GH3.1), two AUX/IAA (IAA1, IAA4), AFR3, ARG7-1, AIT12, ATLP2 and 5NG4 were upregulated during the development of the ovule, which would be associated with the elevated level of JA. Also, 10A5 and 15A are members of SAURs and could be rapidly induced by auxin treatment within 2.5 min 80,81 . The AUX/IAA family is short-lived transcriptional factors and repressors of early auxin response genes at low auxin concentrations, and increased auxin could reduce the level of AUX/IAA proteins by accelerating their degradation. Different GH3s respond differentially to auxin treatment. In Arabidopsis, under IAA treatment, GH3.1 showed strong transcriptional induction, while GH3.5 and GH3.6 showed weaker responses, and GH3.17 showed little or no induction 82,83 . AFR3 is a member of group II AFRs and is thought to be involved in floral meristem, gynoecium, stamen and perianth patterning in Arabidopsis 84 . Exogenous auxin treatment upregulates ARF3 85 .

Cytokinin may regulate the occurrence of MMC and the maturation of ES in ovule development.
In Arabidopsis thaliana, cytokinin is reported to regulate ovule formation and ovule number [86][87][88][89][90][91] and to specify the functional megaspore in the female gametophyte 92 . In J. curcas, the application of CK increased the number of female flower 10,55 . Similarly, our results suggest that CK contributes to female development. CK biosynthesis in female flower buds would be enhanced at stages of MMC occurrence and ES maturation, as suggested by the upregulation of IPT and the downregulation of cytokinin dehydrogenase at these two stages. IPT catalyzes the initial step in the de novo biosynthesis of cytokinin in higher plants, whereas cytokinin dehydrogenase catalyzes the degradation of cytokinin 59 . In response to the increase in CK, several CK response genes, such as HP4-2, HP5, HP1-2, and five A-type ARRs (ARR9-2, ARR5, ARR9-1 and two ARR3) were upregulated. The cytokine signaling cascade typically consists of three functional modules: HK, HP, and response regulator (ARR) 93 . HP proteins serve as phosphorelay carriers between cytokinin receptors and downstream nuclear responses, activating B-type ARR proteins, which in turn, activate A-type ARRs 94 . Additionally, A-type ARRs are reported to participate in ovule development. Megagametophyte defects could be observed in the Arabidopsis mutants arr7 and arr15 95 . In ovules lacking a functional embryo sac, ARR7 and other A-type ARRs (ARR4, ARR5 and ARR6) are

Brassinosteroids may function in the stages from meiosis of MMC to ES maturation.
Brassinosteroids (BRs) suppress the development of the female flower in cucumber 101 , whereas they positively regulate the ovule number in Arabidopsis by regulating related gene expression by BZR1 102 . BRs influence ovule development by regulating the transcription of genes, such as HUELLENLOS (HLL) and AINTEGUMENTA (ANT), which are redundant in the control of the ovule primordial growth 103 . ANT is a direct target of BRZ1, while HLL is regulated in an indirect way 102 . However, in this research, BR biosynthesis in flower buds is expected to decline during ovule development, as suggested by the downregulation of CYP90D1 and CYP90B1. In Arabidopsis, CYP90D1 plays an important role in the early C-22 oxidation of BR biosynthesis, and CYP90B1 catalyzes multiple 22 a-hydroxylation steps in BR biosynthesis 104,105 . In contrast to the decline of BR biosynthesis, BR signaling-responsive DEGs were all upregulated during ovule development, such as XTH22, BZR1, BAK1 and BES1/BZR1, indicating that cross-talk might exist between BRs and other phytohormone signaling pathways. Several reports proposed that BRs and GAs could interact at the signaling leve [106][107][108] . The DELLA protein GAI is thought to inactivate BZR1 by inhibiting the ability of BZR1 to bind to target promoters and negatively regulate BR signaling [106][107][108] . Additionally, auxin can also crosstalk with BRs through G-protein signaling 109 . XTHs, which can be upregulated by BRs, are thought to modify the length of xyloglucans, enabling the cell wall to expand 110,111 . BAK1 is a coreceptor of BRI1, and the binding of BRs to BRI1 relieves the repression of BKI and causes the association of BRI1 with BAK1, as well as a series of phosphorylation events 112,113 . BZR1 and BES1 are two major transcription factors regulated by BIN2 and mediate BR-regulated gene expression 77,114,115 .

Conclusion
To the best of our knowledge, this report is the first to perform a transcriptome analysis of ovule development in J. curcas. A set of ovule development-related genes has been identified in female flower buds through expression profiling analysis, which provided comprehensive information on the ovule development of J. curcas, revealing that five TF families (bHLH, WRKY, MYB, NAC and ERF), five hormones signaling (auxin, GA, CK, BR and JA), five MADS-box genes (AG-1, AG-2, AGL1, AGL11 and AGL14), SUP and SLK3 participate in ovule development. GA and JA are involved in ovule development, as confirmed by their accumulation in flower buds during ovule development. Additionally, CRC and SUP are demonstrated to be involved in ovule development by qRT-PCR analysis. Transcriptome analysis is only an initial step in the exploration of the molecular mechanism of ovule development in J. curcas, but the ovule development-related genes identified in this study will establish a foundation for investigations into the molecular mechanisms of ovule development in J. curcas, and further analyses of these genes are needed to elucidate their roles in ovule development in J. curcas. Grouping on the flower bud samples and total RNA extraction. The flower bud samples stored in RNAlocker were then divided into eight development phases according to the morphology of flower buds at different development phases using an anatomical lens as follows: JCMI, the formation of ten stamen primordia; JCMII, from the occurrence of sporogenous cells to the occurrence of microspore mother cell; JCMIII, from the meiosis of microspore mother cell to the formation of meiotic tetrads; JCMIV, from free microspore stage to the maturation of pollen grain; JCFI, the growth of three carpels; JCFII, from the occurrence of sporogenous cells (SC) to the occurrence of macrospore mother cell (MMC); JCFIII, from the meiosis of MMC to the formation of mononuclear embryo sac (MES); JCFIV, from the formation of MES to the maturation of embryo sac (ES) (Figs 1  and 2).

Materials and Methods
For total RNA extraction, the flower buds from those eight development stages (each phase has three replicates) were ground in liquid nitrogen, and total RNA was isolated using the E.A.N. A Plant RNA Kit (Omega, USA) according to the manufacturer's manual. The obtained RNA samples were then used for library construction.
Library construction and Illumina sequencing. The RNA of flower buds from those eight flower development stages was applied for library construction. A total amount of 3 μg RNA per sample was used for the RNA sample preparations. RNA purity was checked using a NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using a Qubit ® RNA Assay Kit in Qubit ® 2.0 Flourometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Sequencing libraries were generated using the (2019) 9:15973 | https://doi.org/10.1038/s41598-019-52421-0 www.nature.com/scientificreports www.nature.com/scientificreports/ NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB, USA) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. The quality of the library was assessed on the Agilent Bioanalyzer 2100 system and ABI Step OnePlus Real-Time PCR System in Novogene (Beijing, China). The clustering of the index-coded samples was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform, and paired-end reads were generated (Novogene Bioinformatics Technology, Beijing, China).
De Novo Transcriptome Assembly and Abundance Estimation. The raw reads produced from sequencing machines were cleaned by removing reads with adaptors or unknown nucleotides larger than 10% and low quality with the percentage of low-quality bases (base quality ≤10) more than 50%. The left files (read1 files) from all libraries/samples were pooled into one large left. fq file, and right files (read2 files) into one big right. fq file. Transcriptome assembly was accomplished based on the left. fq and right. fq using Trinity 116 with min-kmer-cov set to 2 by default and all other parameters set default. Transcriptome de novo assembly was carried out with the short read assembly program Trinity. The resulting sequences of trinity were called transcripts. When a gene has several transcripts, the longest one was selected as a unigene of the gene. Multiple samples from the same species are sequenced, and unigenes from each sample's assembly can be taken into a further process of sequence splicing and redundancy removal with sequence clustering software to acquire nonredundant unigenes as long as possible. For protein coding region sequence (CDS) prediction, unigenes were first aligned by blastx (E-value < 1.0E −5 ) to protein databases in the priority order of NR, Swiss-Prot. Proteins with the highest ranks in blast results were taken to decide the CDS of unigenes, and then both the nucleotide sequences (5′-3′) and amino sequences of the unigene CDS were acquired by translating the CDS into amino sequences with the standard codon table. Unigenes that cannot be aligned to any database were scanned by ESTScan (3.0.3) to obtain nucleotide sequences (5′-3′) and amino sequences of the CDS.

Different expression genes (DGEs) determination.
Clean read data of each sample were mapped back onto the assembled transcriptome. The resulting read count for each gene was obtained and then applied for gene expression level estimation in each sample using RSEM 117 . Differential expression analysis between two groups was performed using the DESeq R package (1.10.1). The resulting P values were adjusted using Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 were classified as differentially expressed genes (DEGs). The obtained DEGs were then subjected to KEGG enrichment analysis, and KOBAS software was used to test the statistical enrichment of differential expression genes in KEGG pathways 118 . Expression profile of genes involved in ovule development. Total RNA samples from flower buds at four different development phases (JCFI, JCFII, JCFIII and JCFIV) (Fig. 2) were used for expression pattern analysis of 16 selected DEGs (Table S5). The total RNA (1 µg) of each sample was used for first-strand cDNA synthesis using AMV RNA PCR Kit 3.0 (Takara). qPCR was performed on an ABI StepOnePlus Real-Time PCR System (Applied Biosystems, Inc. USA) using 2 × SYBR green PCR mix (QIAGEN, Shanghai, China). The β-tubulin gene was chosen as the endogenous reference gene for the qPCR analysis. The primers used for qPCR analysis are listed in Supplementary Table S1. Each sample has three biological replicates. The average threshold cycle (Ct) was calculated, and the relative expression level of each gene was then calculated according to the 2 −ΔΔCt method. A one-way ANOVA was performed on the gene expression level of the samples at different development stages using the software Statistical Package for the Social Science (SPSS) version 11.5 (SPSS Inc., Chicago, IL, USA). The individual treatment means were compared using the LSD (least significance difference) test.
Quantification of phytohormone. The female flower bud samples at four development phases, JCFI, JCFII, JCFIII and JCFIV, were used for quantification of phytohormones. A 0.1-g sample frozen in liquid nitrogen was ground in liquid nitrogen and then dipped in 1 mL mixture composed of cold acetone: deionized water: hydrochloric acid (37%) = 2:1:0.002 (V:V:V). After incubation at 4 °C for 30 min, 1 mL dichloromethane was added. After incubation at 4 °C for 30 min, the sample was then centrifuged for 15 min at 13,000 g and 4 °C, and the resulting subnatant was collected for drying by liquid nitrogen. Phytohormones in the dried residue were extracted with 0.1 mL methanol. Quantification of phytohormones was performed as described previously using a liquid chromatography-mass chromatography system (Shiseido SP HPLC-Thermo TSQ Quantum Ultra MS/ Ms) with an ODS column (SHISEIDO C18, 5 μm, 2.0 × 150 mm) 119,120,121 .
Availability of Supporting Data All RNA-Seq data for this project has been deposited in NCBI under the SRA accession SRP115141. Summary description of library Illumina sequencing, De novo transcriptome assembly Table S1 and annotation are listed in Supplementary Data Tables S2-S5, Figs S1-S7.