Tingting Li · Fosheng Li · Lanju Mei · Na Li ·Min Yao · Lin Tang

Abstract Idesia polycarpa Maxim. var vestita Diels. is a dioecious tree species native to eastern Asia. There are diff iculties associated with distinguishing the sex of the plant at the seedling stage. In order to explore the mechanism of sex differentiation in f lower development, we conducted the transcriptome prof iles of male and female f lowers at early,metaphase and late developmental stages. Approximately 123,335 unigenes with a total length of 83,996 Mb and an average length of 168 bp were assembled. The unigenes were blasted into Nr, Nt, Pfam, KOG/COG, Swiss-prot, KEGG,GO databases. Homology analysis demonstrated that I.polycarpa and black cottonwood had the highest homology with the alignment of 92,871 sequences. This study identif ied 80 groups of transcription factor families with a total of 1475 unigenes, mainly including MYB, WRKY, AP2 and bHLH transcription factor families. KEGG pathway analysis showed that the expression of numerous plant hormones (cytokinin, gibberellin and ethylene) and f lavonoid biosynthesis pathway were different at various stages of female and male f lower development. In addition, a number of unigenes associated with f lowering were identif ied which were key genes associated with photoperiodic, vernalization, thermosensory, gibberellin, and autonomic pathways.The results show that I. polycarpa f loral organ development was in accordance with the ABCDE model, in which the down-regulation of the B gene family might affect stamen fertility in late stages of female f lower development. qRTPCR experiments validated that the expression patterns of 15 unigenes were consistent with those in RNA-seq results. The results highlight a central role for plant sex identif ication in seedling production and a sex-determining mechanism for dioecious plants. In addition, the transcriptome data provided a theoretical basis for I. polycarpa genetic diversity analysis and molecular- assisted breeding.

Keywords Dioecious p lants · Transcriptome s equencing ·Sex differentiation · Idesia p olycarpa · Floral o rgan development · Transcription f actor

Introduction

Many angiosperms have hermaphrodite f lowers with both female and male functions but about 30% produce unisexual f lowers (Richards 1997). However, the proportion of dioecism, where male and female f lowers are on separate individuals, is about 6% (Renner 1995), such as papaya, willow,poplar, and Jatropha. The evolution of dioecious plants may take two differentiation pathways: one where the female organs in f lowers with both sexes loose the female function to the male line, the other is a monoecious plant due to single sex evolution (Aryal and Ming 2014). In recent years, numerous efforts have been made to identify the sex determination in plants (Chen et al. 2017).

Idesia polycarpaMaxim. varvestitaDiels belongs to the Salicaceae family and is a dioecious tree species distributed mainly in eastern Asia (Jiang et al. 2014). Dried fruits of the species contain 26—47%I. polycarpaoil and are rich in oleic, linoleic, and linolenic acids and other unsaturated fatty acids which is more diverse than oil components of olive oil(Zhang 2011). The seed oil ofI. polycarpais highly nutritious and healthy, high quality edible oil (Kai-Lin et al. 2009;Guo et al. 2012). However, distinguishing the morphological characteristics of male and female f lowers has always been difficult. At present, seedlings are planted and male trees are removed after f lowering in three years, which results in high costs of management and cultivation. Research onI.polycarpaboth in China and in other countries is focused on the evaluation of cultivation techniques (Wu et al. 2010),and the extraction and processing of oil and fat components(Liu et al. 2011; Gong et al. 2012). Research has also been carried on transcriptome analysis ofI. polycarpamale and female f lowers, but the sampling time was late, the samples few in number, and did not show the variety of genes at different f lowering stages (Mei et al. 2017). Therefore, transcriptome information ofI. polycarpais still unknown at different stages of f lower development which seriously hinders an examination of the f lower development mechanism at a molecular level. Transcriptome analysis is the expression,at the cellular or tissue level of the complete set of RNA transcripts of the genome, and can identify low numbers of genes and unknown transcripts, and considerably facilitates research on the function of human genome species (Wang et al. 2009a, b; Varshney et al. 2009). Genomic expression of samples at different developmental stages can be obtained analyzed, and candidate genes associated with a developmental mechanism identif ied (Br?utigam and Gowik 2010).Advances in sequencing technology have greatly facilitated the development of transcriptome studies, especially high throughput sequencing based on second generation sequencing. In addition, transcriptome technologies do not require the genome information of a species to be known, which accelerates the study of species with unknown genomes(Bentley et al. 2008). Currently, transcriptome technologies have also been widely used in the area of plant f lower development. For example, Zhang et al. ( 2013) screened 120 genes associated with orchid f lowering, providing a theoretical basis for removing impediments through cloning f lowering orchids by transcriptome sequencing. In addition,transcriptome research can also compare female and male f lower expression differences, and analyze the mechanism of unisexual f lower formation, such as sex differentiation of the dioeciousSalix suchowensisW.C. Cheng ex. G. Zhu(Liu et al. 2013). Different f loral formation mechanisms have been investigated by transcriptome analysis in cucumber(Guo et al. 2010). The male and female f lower structures of wheat and the occurrence of male f lower mutations were examined by transcriptome analysis (Yang et al. 2015), as well as comparative transcriptome analysis of male and female f lowers of the monoeciousQuercus suberL. (Rocheta et al. 2014).

Therefore, this study focused on male and female f lower buds ofI. polycarpaRNA sequencing at different developmental stages. It provides a theoretical basis and mechanism ofI. polycarpaf lower development to enrich the genetic information, and lays the foundation for subsequent research on reproductive biology ofI. polycarpa.

Materials a nd methods

Plant m aterial a nd RNA e xtraction

I. polycarpamaterials were collected at the Yanjiang District corn farm, Ziyang, Sichuan Province (30°08′N, 104°82′E).There were 10 male and 10 female plants, 5 m in height,selected for f ive periods of inf lorescence collection, starting from the end of March to the end of April, 2015. The f lower bud samples were f itted to a freezing tube, placed in liquid nitrogen for temporary preservation, and then transferred to a refrigerator at - 80 °C for RNA isolation. At the same time, a portion of the buds were f ixed in FAA solution for subsequent anatomical observation and identif ication of developmental stages (Fig. 1). Developmental stages were divided into f ive periods: early differentiation(F1, M1), pre-differentiation (F2, M2), differentiation (F3,M3), completion of differentiation (F4, M4), and late differentiation (F5, M5). Male f lowers (M1) formed earlier than female f lowers with stamens but female f lowers were still in a sepal primordia growth stage (Fig. 1). Stamen primordia in female and male f lowers continued to develop, anthers and f ilaments were inconspicuous, and ovaries and styles began to expand. Pistils and stamens were similar (Fig. 1,F2 and M2). Moreover, in F3 and M3 (Fig. 1), female and male f lowers continued to develop, f ilaments and anthers are seen, and stigmata are visible and ovules appear in the ovary.However, in F4 and M4 (Fig. 1), ovaries of the female f lower continued to expand, but stamen development stopped before the formation of tetrads so pollen grains could not be formed and a cavity developed in the anther chamber. Stamens of male f lowers continued to develop and the primordia of the pistils began to develop but with incomplete ovules in degenerated ovaries. Finally, differentiation between female and male f lowers was more obvious. Ovaries continued to grow in the female f lowers and stamens atrophied; male f lower anthers became larger and anther wall atrophy in F5 and M5 (Fig. 1). According to research (Chu et al. 2013),inf lorescence development ofI. polycarpawas divided into three stages: namely, an inf lorescence primordia differentiation stage, f loral organ primordia differentiation stage, and male and female structure differentiation stage. Therefore,three male and female mixed samples were used for transcriptome sequencing in order to explore the inf luence of male and female f lower differentiation factors. F1 and F2 are a mixed sample Fa, M1 and M2 a mixed sample Ma; the sample F3 is Fb and M3 Mb, F4 and F5 are a mixed sample Fc, and M4 and M5 a mixed sample Mc.

Fig. 1 Female (F1—F5) and male (M1—M5) f lowers at different stages of development

The total RNA ofI. polycarpaf lower buds was extracted using a Plant Total RNA Isolation Kit (FOREGENE). Degradation and contamination of RNA was monitored on 1%agarose gels and RNA purity was checked using the Nano-Photometer spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using a Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, CA, USA).Finally, RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system(Agilent Technologies, Santa Clara, CA, USA).

Library c onstruction a nd RNA-seq

These were performed by the Novogene Bioinformatics Technology Corporation, Beijing, China. An amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEB Next Ultra? RNA Library Prep Kit for Illumina(NEB, USA) following manufacturer’s recommendations,and index codes were added to attribute sequences to each sample. Library quality was assessed on the Agilent Bioanalyzer 2100 system. Clustering of the index-coded samples was then performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq? 2000 platform.

Sequence a ssembly a nd functional a nnotation of the unigenes

Raw data (raw reads) of fastq format were f irst processed through in-house perl scripts (Zhang et al. 2015). In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N, and low quality reads from raw data. At the same time, Q20, Q30,GC-content and sequence duplication levels of the clean data were calculated. All downstream analyses were based on clean data of high quality. In addition, the left f iles, (read1 f iles), from all libraries and samples were pooled into one large left.fq f ile, and right f iles (read2 f iles) into one large right.fq f ile. Transcriptome assembly was accomplished based on left.fq and right.fq using Trinity (Grabherr et al.2011) with a min_kmer_cov set to 2 by default and all other parameters set default.

To obtain theI. polycarpaunigenes annotation information, all unigene sequences were performed by BLASTx search (E-value ≤ 10 -5 ) against multiple nucleic acids and protein common databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences);Pfam (protein family); KOG/COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (a manually annotated and reviewed protein sequence database); KO (KEGG Ortholog database); and, GO (gene ontology). Unigenes and gene ontology databases were compared using blast2go software(Conesa et al. 2005). WEGO software (Ye et al. 2006) was then used to perform GO function classif ication statistics.Metabolic pathways involved in genes and their role in these pathways were analyzed using the KEGG database (Kanehisa et al. 2007).

Analysis of differentially expressed genes (DEGs)

The expression level of unigenes was normalized by FPKM,fragments per kilo base of gene per million mapped reads(Trapnell et al. 2010) in order to eliminate any errors caused by the amount of the sequenced data and the length of the sequence. FDR (False Discovery Rate) was used to perform multiple hypothesis testing correction. For each sequenced library, the read counts were adjusted by the edgeR program package through one scaling normalized factor. Differential expression analysis of two samples was performed using the DEGseq (Wang et al. 2009a, b) R package (Wang et al.2009a). Thep- value was adjusted using aq-value (Paul et al. 2003). Theq-value ＜ 0.005 and |log2 foldchange| ＞ 1 was set as the threshold for signif icantly differential expressions. The classif ication of DEGs was performed on GO and KEGG analysis with correctedpvalue ≤ 0.05 andq-value ≤ 0.05. Gene ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GOseq R packages based on Wallenius non-central hyper-geometric distribution (Young et al. 2010), which can adjust for gene length bias in DEGs. KOBAS software (Mao et al. 2005) was used to test the statistical enrichment of differential expression genes in KEGG pathways.

Validation of s equencing dat a b y qPCR

Fifteen differentially expressed unigenes were randomly selected from theI. polycarpatranscriptome data. The actin homologous gene was chosen as a reference gene qRT-PCR from transcriptome database, which was used to design primers (Table S1). Total RNA was extracted by using Plant Total RNA Isolation Kit (FOREGENE). Subsequently, the RNA of these samples was transcribed into cDNAs using the Takara reverse transcription kit (PrimeScript? RT reagent kit, Dalian, China) with a 10 μL reverse transcription reaction system. Real-time assays were performed with SYBR Green Dye using CFX96 Touch? Real-Time PCR (Bio-Rad, Hercules, CA, USA) detecting platform. The qRT-PCR amplif ication program had the following cycle conditions:95 °C for 3 min, followed by 40 cycles of 95 °C for 10 s,55 °C for 10 s, and 72 °C for 20 s. The qRT-PCR 10 μL reaction system was constructed using the reference kit SsoFast? EvaGreen?Supermix (Bio-Rad, Hercules, CA,USA). Three biological replicates and technical duplication were used for each gene. The primer amplif ication effectiveness was analyzed using CFX Manager? (v3.0) software(Vandesompele et al. 2002), and the expression of unigenes was calculated by the 2 -ΔΔCt method (Pfaffl2001).

Results

Illumina sequencing and de novo assembly

For an overview of theI. polycarpaf lower transcriptome,six cDNA libraries were constructed of male and female f lowers covering early, mid and late development stages and were sequenced on the Illumina Hiseq?2000 platform(Table S2). The total output of 268,177,572 raw reads was submitted to the NCBI (Nation Center for Biotechnology Information) Sequence Read Archive, the SRA accession:SRP138341. After the f iltering of raw data, 259,662,844 high-quality clean reads were obtained and the total data of the six libraries’ clean reads were approximately 38.95 G (Table S3). The clean reads were de novo assembled to obtain a reference sequence for subsequent analysis by Trinity programs (transcriptome splicing software), resulting in 123,335 unigenes (Fig. S1), with a total length of 83,996 Mb, an average length of 681 bp and a N50 (length of spliced transcripts not less than 50% of total length) of 1074 bp. The length of the unigenes ranged from 200 to 16,715 bp. Unigenes less than 1000 bp were 102,542, 83.1%of the total; lengths between 1001 and 2000 bp were 12,330 and accounted for 10.0%. Unigenes more than 2000 bp in length were only 8463 and accounted for only 6.9% of the total volume. These results suggest that high-quality transcription data ofIdesia polycarpahad been produced, providing a good foundation for subsequent analysis.

Functional a nnotation a nd classif ication of the unigenes

All 123,335 unigenes were matched by BlastX comparisons against seven public databases, including Nr, Nt, Pfam,KOG/COG, Swiss-prot, KEGG and GO with unigene numbers of 38,172 (30.9%), 48,419 (39.3%), 25,935 (21.0%),11,761 (9.5%), 25,077 (20.3%), 8493 (6.9%), and 27,894(22.6%) (Table S4). A total of 59,471 (48.2% of all unigenes)were annotated sequences in at least one of the databases.In addition, 51.8% were without annotation information. It is possible that these unigenes were shorter in length and reduced the accuracy of the alignment because it has been reported that the validity of data alignment is related to the length of the sequence being compared, and the longer the length, the higher the accuracy (Novaes et al. 2008).

Fig. 2 Characteristics of sequence homology of I. Polycarpa unigenes against Nr databases. a The E-value distribution of the NR annotation; b similarity distributions of NR annotations; c NR annotated species distribution

In the NR database, 38,172 (30.9%) unigenes were annotated results, in which the matching degree, similarity and distribution of species are shown (Fig. 2). In this study, E-value ＜ 1.0e -60 unigenes were considered as a higher homology and E-values between 1.0e -5 and 1.0e -60 had moderate homology. The E-value distribution analysis revealed that homologous unigenes accounted for 25.2%; the other 50.8% moderately homologous. The ratio of identical to other plant genomes (E-value = 0) accounted for 24.0%inI. polycarpa(Fig. 2 a). The similarity analysis indicated that the similarity degree, (above 80.0%), accounted for 76.2% and unigenes, (between 60.0 and 80.0%) accounted for 19.8%, while the similarity degree (less than 60.0%) was only 4.0% (Fig. 2 b). Comparative analysis of the E-value distribution and the similarity distribution revealed that the unigenes matching degree was higher in the Nr database.A homology comparison ofPopulus tomentosaannotated sequences found that 75.3% were matched to the unigenes ofI. polycarpa, followed byVitis vinifera(3.4%),Jatropha curcasL. (2.2%),Ricinus communis(2.0%) andTheobroma cacao(1.5%) (Fig. 2 c). There are higher homologies or similar relations betweenI. polycarpaand tree species compared with vines and shrubs becauseIdesiais a macrophanerophyte, a species with buds between 30 and 50 m above the ground. Most of the annotated unigenes exhibited high hits with the sequence from the family Salicaceae.Therefore, the analysis suggests thatI. polycarpais closely related to the generaSalixandPopulus, which is in accord with previous research (Yang et al. 2016; Mei et al. 2017).

Fig. 3 function classify of All-Unigene. a KOG function classif ication of All-Unigene; b GO classif ication of All-Unigene; c KEGG pathway classif ication of All-Unigene; only the f irst 23 metabolic pathways involving unigenes are shown

Some 13,238 unigenes (10.75) of the transcripts were annotated to the KOG database, where one unigene may correspond to one or more KOG subcategories (Fig. 3 a).The functional category with the largest group of unigenes was ‘general function prediction’ (2156 unigenes), followed by ‘posttranslational modif ication, protein turnover,chaperones’ (1554 unigenes), ‘signal transduction mechanisms’ (1123 unigenes), ‘translation, ribosomal structure and biogenesis’ (905 unigenes), and ‘secondary metabolites biosynthesis, and transport and catabolism’ (704 unigenes).However, the least functional categories were ‘cell motility’(4 unigenes) and unnamed proteins (2 unigenes). According to GO functional classif ication statistics, the 27,894 unigenes, (22.6%) of the GO annotation, were normalized to 56 functional groups, including 23 biological processes (72,929 unigenes), 14 cellular components (51,239 unigenes) and 19 molecular functions (33,303 unigenes). In the biological processes category, the groups with the largest number of unigenes were “cellular process”, “metabolic process”, “and single organism process”. In addition, “cell”, “cell part”and “organelle” were the most highly represented groups in the cellular components category. Under the molecular functions category, “binding” and “catalytic activity” were the most highly represented groups (Fig. 3 b). In this study,8493 unigenes (6.9%) were classif ied to the KEGG database,involving a total of 274 metabolic pathways. Unigenes were divided into f ive branches based on the KEGG metabolic pathway involved, namely: “metabolism”, “genetic information processing”, “environmental information processing”,“cellular processes”, and “organism systems”. As shown in(Fig. 3 c), these branches were classif ied into 33 categories,in which most representations by the unigenes were “signal transduction” (1076 members; 12.7%), “translation”(863 members; 10.2%), “carbohydrate metabolism” (834 members; 9.8%), and “energy metabolism” (740 members;8.7%). All annotations and categorizations suggest that f lower development involves a large number of metabolic pathways and play an extremely important role in future genomic function studies.

Identif ication of transcription factors,f lowering-associated genes and MADS-box genes

Transcription factors are important regulators of plant growth and development and are divided into numerous families according to their domains. Different transcription factors play different roles in the process of ontogeny. In this study, we identif ied 80 transcription factor families for a total of 1475 unigenes. These families mainly included MYB (121), AP2 (71), bHLH (64) and WRKY (55) (Fig. 4).

Flower development is a complex process and induction is regulated by numerous pathways such as the photoperiod pathway, the vernalization pathway, the thermosensory pathway (Blazquez and Ahn 2003), the gibberellin pathway and the autonomic pathway. In addition, f lowering induction signals also require a key integration link, for example,f lowering integron genesFTandSOC1, f lowering inhibitory genesTFL1andSVP, and f loral meristem determining genesLEAFY(LFY). Flowering related genes and unigenes of theI. polycarpadatabase were identif ied (Table S5). These included photoperiod pathway genes such asCONSTANS(CO),FRIGIDA(FRI),Constitutively Photomorphogenic(COP1),Gibberellic Acid Insensitive(GAI), andGibberellic Insensitive Dwarf1(GID1). This study has shown that some related homeobox (homologous DNA sequences affecting development) genes exist in the database, and that some of the f lower development- related genes in the male and female f lower structures were differentially expressed.

The ABCDE model describes in detail a series of orthologous genes involved in fl oral organ determination and their mode of action. The f loral organ homeobox genes inI. polycarpaare listed in Table S6. A family AP1 gene was up-regulated in female f lowers at the mid stage of differentiation, and a B gene,PMADS, similar to one found inPetunia hybrida, was signif icantly down- regulated in female f lowers at the late development stage. However, the homologous genes ofAGAMOUS(AG) (C gene) inI. polycarpawas expressed in similar amounts in male and female f lower buds. In addition, a gene similar toAGL11(D gene)was identif ied in the database which was similar in male and female f lowers. Finally, three unigenes related to the E functional genesSEPALLATA 2(SEP2) andSEPALLATA3(SEP3) (Cui et al. 2010) were found and an E homologue was up-regulated in male f lowers at the early and mid-stages of development. The research showed that theI. polycarpaf lower organs differentiation model generally agreed with the ABCDE model.

Differentially expressed genes (DEGs) analysis

An experiment with both male and female f lower organs was carried out with f ive groups of differentially expressed genes over the three stages of female f lower development(Fa—Fb—Fc), among the three stages of male f lower development (Ma—Mb—Mc), and between the three periods of female and male flower formation (Fa—Ma, Fb—Mb,Fc—Mc) (Table S7; Fig. 5). There were 494 differentially expressed genes identif ied in the early and mid- stages of female f lower development (Fig. 5). There were 119 differentially expressed genes in the mid- and late stages, and 537 genes in early and late stages. Compared with the early to mid- stages, the differentially expressed genes were more developed in the medium to late stages, indicating that many differential genes are for later differentiation than in the earlystages of development. The DEGs between early and midstages were signif icantly more than those between the midand late stages. The difference of male f lower development expression was similar to that of female f lowers. At the early stage of development, the expression of differential genes was more active, and the up regulation and down regulation of entire genes were more signif icant than those of late stages.

Fig. 4 Statistical analysis of transcription factor expression differences of male and female f lowers

Fig. 5 Venn diagrams of differentially expressed genes in different samples. a (1) The total female DEGs in early,middle and late stage of f lower development; a (2) The total male DEGs in early, middle and late stage of f lower development; b (1) Female f lowers up- regulated differentially expressed genes in early, middle and late stages of development;b (2) male f lowers up- regulated differentially expressed genes in early, middle and late stages of development; c (1) Female f lowers down- regulated differentially expressed genes in early, middle and late stages of development; c (2) male f lowers down- regulated differentially expressed genes in early, middle and late stages of development

GO and pathw ay a nalysis of D EGs

GO enrichment of DEGs between male and female f lowers was analyzed at early, mid- and late stages. At early differentiation, 938 DEGs were identif ied in female and male f lowers. There were 773 genes (82.4%) for GO annotation, of which 25 GO subterms or subordinate terms were signif icantly different (correctedp-value ＜ 0.05) in annotation (Table S8). In the biological processes category, metabolic process (GO: 0008152), single-organism metabolic process (GO: 0044710) and oxidation—reduction process(GO: 0055114) were the most abundant subterms. The most signif icant enrichment of a subterm was cell periphery (GO:0071944) in the cellular components category. The oxidoreductase activity (GO: 0016491) subterm was dominant in the molecular function category (Fig. 6 a). Up-regulated DEGs were mainly concentrated in the “metabolic process”(GO: 0008152) and the “single organism metabolic process”(GO: 0044710); down-regulated DEGs were enriched in the“single organism metabolic process” (GO: 0044710), the“oxidation—reduction process” (GO: 0055114) and the “oxidoreductase activity” (GO: 0016491) due to the metabolic activities of male and female f lowers (Fig. S2 and S3). Themost signif icant 15 subterms are listed in Tables S9 and S10. Many of the up-regulated genes in female f lowers are related to photosynthesis, such as “response to red light”,“red light”, “chloroplast”, “chloroplast part”, and “response to blue light”. Reports in the literature indicate that, in early stages of f lower spike formation more chloroplasts are accompanied by strong photosynthesis, with f lower development gradually weakening, but in fruit development,

photosynthesis then increased (Hui et al. 2007). This suggests that, in early female f lower development, inf lorescence photosynthesis is more active for subsequent f lower development by the accumulation of prime metabolites.

Fig. 6 GO classif ication and enrichment analysis of differentially expressed genes between female and male f lowers. a GO classif ication and enrichment analysis of differentially expressed genes between female and male f lowers in the early stage of differentiation; b GO classif ication and enrichment analysis of differentially expressed genes between female and male f lowers in the middle stage of differentiation; c GO classif ication and enrichment analysis of differentially expressed genes between female and male f lowers in the later stage of differentiate

At the mid- stage of differentiation, 512 differentially expressed genes (DEGs) of male and female f lower structures were annotated by GO, accounting for 81.4% of the total 629 unigenes. (629).There were 20 GO signif icantly different subterms (correctedp-value ＜ 0.05), the top f irst 15 in Table S11. In the biological processes category, a high percentage of DEGs was assigned to the metabolic process(GO: 0008152), the single-organism metabolic process(GO: 0044710) and to the organic cyclic compound biosynthetic process (GO: 1901362). The fatty acid synthase complex (GO: 0005835) was the most abundant subterm in the cellular components category. The catalytic activity(GO: 0003824) and oxidoreductase activity (GO: 0016491)represented the most abundant subterms in the molecular function category (Fig. 6 b). As shown in Fig. S4 and S5,the up-regulated DEGs over representations of GO terms were in the “organic cyclic compound biosynthetic process” (GO: 1901362), followed by the “aromatic compound biosynthetic process” (GO: 0019438), the “heterocycle biosynthetic process” (GO: 0018130), the “regulation of nucleobase-containing compound metabolic process” (GO:0019219) and the “DNA binding” process (GO: 0003677).The “metabolic process” (GO: 0008152), the “single-organism metabolic process” (GO: 0044710), the “carbohydrate metabolic process” (GO: 0005975) and the “carboxylic ester hydrolase activity” (GO: 0052689) were enriched in the down-regulated DEGs. The most signif icant 15 subterms are listed in Tables S12 and S13. At the mid- stage, male pollen continued to develop and ovules began to cease to develop; however, pollen still developed in the female f lower and the ovaries began to dilate. Therefore, a large number of transcripts and translations of genes in the female f lower were involved in its growth. The results suggest that a large number of genes were expressed to provide the energy base for female f lowers.

At the late stage of differentiation, among the 538 DEGs of male and female f lowers, 425 were annotated with GO(80.0%), and 30 GO subterms were signif icantly different(correctedp-value ＜ 0.05) (see Table S14 for the f irst 15 subterms). The oxidation—reduction process (GO: 0055114),lipid metabolic process (GO: 0006629) and cellular lipid metabolic process (GO: 0044255) were the most highly represented groups in the biological processes category. In the cellular components category, the highest percentage of genes came from the fatty acid synthase activity (GO:0005835). The oxidoreductase activity (GO: 0016491), followed by hydrolase activity, hydrolyzing O-glycosyl compounds (GO: 0004553), coenzyme binding (GO: 0050662)and hydrolase activity acting on glycosyl bonds (GO:0016798), were dominant groups in the molecular function category (Fig. 6 c). The corrected_p-value ＞ 0.05 of the upregulated DEGs are shown in Fig. S6. In addition to the“l(fā)ipid metabolic process” (GO: 0006629), the “cellular lipid metabolic process” (GO: 0044255) and the “l(fā)ipid biosynthetic process” (GO: 0008610) were the highest represented groups in the down-regulated DEGs. The most signif icant 15 terms are shown in Tables S15 and S16. At this stage of development, male f lower pollen had gradually matured and ovaries had degraded, female f lower pollen had aborted and disappeared, and the ovaries had enlarged. The results suggest that the large amount of lipid metabolism in male f lowers might provide an energy source for another development. However, down-regulation of development-related genes in female f lowers was due to pollen abortion and ovary formation. GO functional enrichment analysis showed that gene functions of male and female f lowers at various developmental stages were different and one of the main reasons for sex differentiation in the development of male and female f lowers.

Pathway enrichment of differentially expressed genes between male and female f lowers at different periods is shown in Table S17. At the early stage among the 938 DEGs, 140 were annotated into the KEGG pathway, and 10 had highly signif icant differences. The pathways most represented were “protein processing in endoplasmic reticulum”(ko04141), “plant hormone signal transduction” (ko04075),“photosynthesis-antenna proteins” (ko00196) and “estrogen signaling” (ko04915) (Table S17[1]). However, the up-regulation of photosynthesis-antenna proteins might indicate that female f lowers photosynthesized at earlier stages of development due to the accumulation of resources for later f lower development. Protein processing in the endoplasmic reticulum was up-regulated in female fl owers as a large amount was necessary for early stages of f lower development, and prepared for massive gene expression during the medium term. At the mid- stage, 629 DEGs were assigned to 160 KEGG pathways, in which 14 are signif icantly different. The dominant pathways were “glycolysis/gluconeogenesis” (ko00010) and “plant hormone signal transduction”(ko04075) (Table S17 [2]). It was speculated that endogenous hormones play a signif icant role in the development of male and female f lowers. Female hormones were signif icantly up-regulated and enriched while the metabolism of these hormones in male f lowers was signif icantly lower. The hormones involved in regulation include cytokinin, gibberellin, ethylene, and salicylic acid. At the late stage of development, 538 DEGs were assigned to 90 KEGG pathways of which 13 were signif icantly different. The “phenylpropanoid biosynthesis” (ko00940) pathway and the “f lavonoid biosynthesis” (ko00941) pathway were the most highly represented(Table S17 [3]). Flavonoids, a large class of natural pigments are widely distributed in plants and accumulate in leaves and f lowers. Previous studies have reported that f lavonoids not only attract insects based on color and specif ic odors but also attract insects pollinated plants (Harborne and Williams 1973). Therefore, female f lowers at a late development stage had a high expression of f lavonoids, making the ovary appear yellow- green and increasing a specif ic odor to attract insects prepare for subsequent f lowering and pollination.

Validation of unig enes and g ene e xpression pr of iling

To verify that gene expression patterns were consistent with RNA-seq results, 15 differentially expressed genes were selected from the early, mid, and late stages of f lower development (Table 1), each of which was three biological repetitions (Fig. 7). The RT-qPCR results were consistent with the expected results and similar to the RNA-seq results. These results conf irm that the transcriptome differences between the different libraries and the transcriptome database ofI.polycarpaare valid.

Discussion

Effect of hormones on f lower organ development

KEGG pathway analysis of DEGs determined that plant hormones were involved in the development ofI. polycarpamale and female f lowers. At the early stage, male f lowers had a large number of hormone- regulated pathways (auxin,cytokinin, abscisic acid, ethylene), which were up-expressed,conf irming that these hormones play an important regulatory role in the development of male f lowers. Male f lower development in Japanese white birch (Betula platyphyllaSukaczev) maintain high auxin levels which is benef icial to the differentiation of stamen primordia and the formation of sporulation cells (Yang et al. 2002). Research on cytokinin and sex regulation is mainly focused on the mechanism of estrogen promotion. For example, exogenous cytokinin can signif icantly promote female structures inJatropha curcasL. (Ai et al. 2002; Pan and Xu 2011). However, little is known about the mechanism of androgenic cytokinin and further studies are needed to determine how to promoteI.polycarpamale f lower development. In this study, the upregulation of abscisic acid is benef icial in promoting the differentiation of f lower buds, participating in the adjustment of the photoperiod, and regulating gibberellic acid inhibition (Xu et al. 2007). Ethylene is an estrogen- stimulating hormone (Bai 2010) but it also has been reported to have effects on male stimulation, such as the application of ethephon toLuffa cylindricaL. which helps male f lower differentiation (Mahida and Valia 2015). However, little work has been done on how to contribute to male f lower development. The transcriptome analysis in this study found the ethylene metabolism- related gene increase signif icantly in early male f lower formation, which is different from other reports. Therefore, the effect of ethylene onI. polycarpamale f lower development needs to be further explored. At the- mid stage, salicylic acid is a very low organic compound synthesized in plants. A previous report indicated that salicylic acid accelerated pollen germination and pollen tube growth (Zhao et al. 2012). Cytokinin, gibberellin and ethylene are common estrogenic hormones which play important roles in promoting pistil development at the metaphase stage of female f lower differentiation. These genes were downregulated at the mid- stage of male f lower differentiation,indicating that they might be an important reason for the abortion of pistil development in male f lowers.

Table 1 The 15 differentially expressed genes data in RT-qPCR experiments

Fig. 7 Results of qPCR validation

Candidate t ranscription f actors a s regulators o f f lower organ development

Previous studies have indicated that transcription factors regulate the expression of many downstream genes and play an important role in growth and development,metabolism and responses to the external environment(Hobert 2008). In this study, we mapped 80 transcriptome families, the most representative of which were MYB(121), AP2 (71), bHLH (64) and WRKY (55). Research has shown that MYB transcription factors are associated with the development of pistils and stamens (Dubos et al. 2010). For example,MYB124andMYB88can dedicate to the development of female flowers (Makkena et al. 2012), while inArabidopsis,AtMYB21,AtMYB24,AtMYB57,AtMYB108,AtMYB35/TDF1,AtMYB80 (formerlyAtMYB103),AtMYB99,AtMYB33andAtMYB65)play signif icant roles in the development of anthers (Millar and Gubler 2005; Cheng et al. 2009; Mandaokar and Browse 2009; Phan et al. 2011). In the early and late stages of differentiation, a gene similar toMYB44was highly expressed in female f lowers, and aMYB122family gene was signif icantly up-regulated in male f lowers at the mid and late stages of differentiation. According to related research,MYB44was involved in plant responses to abiotic and biotic stress regulation (Lü et al. 2011; Persak and Pitzschke 2013), and regulation ofArabidopsis thalianaseed germination (Nguyen et al. 2012). However, the impact on the development of female and male f lowers has not been thoroughly investigated. In addition, it has been reported thatMYB122involved in the biosynthesis of indole glucosinolates (Frerigmann and Gigolashvili 2014),and indole glucosinolate myrosinase hydrolyzed to form a precursor of IAN, suggests thatMYB122may regulate indole acetic acid pathway, thereby affecting the development of male f lowers.DYSFUNCTIONAL TAPETUM 1(DYT1) (Zhang et al. 2006) andABORTED MICROSPORES(AMS) (Sorensen et al. 2003), two kinds of bHLH transcription factor family genes, are important for the normal development of anther and tapetum(the innermost cell layer of the young sporangia of vascular plants). In this study, anAMShomolog gene (c76147_g2) was identif ied, weakly expressed in female f lowers but strongly in mid- and late stages of male f lower development, inferring that this gene may be a key in male f lower development. Stamen development in female f lowers stagnated before the fourth division with disintegration in each layer of the anther wall by the mid stage. Microspore mother cells were decomposed and could not form normal anthers.Therefore,AMSmay be an important gene for mid and late development of male f lower differentiation. Moreover, a unigene homologous withbHLH14was up-regulated in the early stages of differentiation of male f lowers and has been reported to regulate pollen development through the jasmonic acid pathway (Song et al. 2014). AP2 transcription factors are associated with a variety of biological processes such as growth and development, pathogen defense,environmental stress response (drought, high salinity), and hormone regulation (Kang et al. 2011). An AP2 homolog unigene (c57799_g1) was up-regulated in female f lowers at the early stage. Three unigenes, (c57799_g1, c65425_g1,c70825_g1), were down-regulated in female f lowers and signif icantly up-regulated in male f lowers in late developmental stages, suggesting that the genes may participate in f lower development by hormone regulation. WRKY is a specif ic transcription factor of higher plants and is often regulated by hormones (Eulgem et al. 2000). Our research identif ied 58 differentially expressed transcription factors in male and female f lowers which may be involved in the regulation of hormones in f lower development.

Flowering- r elated u nigenes o f Idesia polycarpa

These were detected by transcriptome database, for example, key genes involved in the photoperiod, vernalization,autonomous, and gibberellin(GA) pathways were transmitted to three f lowering pathway integrators, f lowering integron genesFT,SOC1, and the f loral meristem geneLEAFY(LFY) (Moon et al. 2005). In the photoperiod control pathway, three genes ofCONSTANS(CO) in male f lowers were signif icantly up-expressed at an early stage. Two genes were signif icantly up-expressed in female f lowers and two others signif icantly increased in male f lowers by the mid stage. The complex ofFKF1andGIcould inhibit the expression ofCO, moreover,COP1, andSPAdegradedCOby the ubiquitination proteasome pathway (Laubinger et al. 2006; Liu et al. 2008). In early stages, theFKF1,COP1andSPAgenes in female f lowers increased signif icantly, which led to the decrease ofCO. In addition, threeGID1(GA cytoplasmic and nuclear localization receptors) homologous genes were identif ied and signif icantly increased in early male f lower differentiation, indicating that gibberellin may have a critical effect on f lower induction in male f lowers. It has been shown that required for proper f loral organ growth and elongation(Yu et al. 2004) and believed to be strongly devoted to regulate f lower bud formation. The expression ofFLOWING,LOCUS,C(FLC) f lowering inhibitors was up-regulated in early female f lowers.FLCis a f lowering inhibitor and represses the expression of the f loral activators,SOC1andFT(Amasino 2010). In addition, it was observed that male f lowers generally developed earlier than female f lowers.Therefore, it was concluded that photoperiod- related genes were down-regulated, and f lowering repressors up-regulated to delay female f lower development. Photoperiod- related genes were up-regulated, and gibberellin receptor up-regulation and down-regulation of the f loral repressor promoted the development of early molecular male formation in the f loral induction pathway.

MADS-box g enes in fl oral organs Idesia polycarpa

In this study, a number of f loral homeotic genes have been identif ied inIdesia polycarpa.Some were signif icantly upregulated in the mid-stage of female f lower development,andAP1was involved in the formation of the calyx as well as the regulation of other genes. For example,AP1andAP2can promoteSEP3expression either directly or through theLFYpathway.SEP3is a coactivator of B and C function genes withLFY, which is suppressed bySOC1,AGL24, andSVP(Moyroud et al. 2011). Therefore, during mid- f lower development,SEP3was increased by the up- regulation ofAP2.A decrease ofSEP3in female f lowers would affect their development but there was little difference inSEP3expression between male and female f lowers at the late stage of development.PMADSis the B gene ofPetunia hybrid, of which the homologous expression was signif icantly reduced in female f lowers (Krol et al. 1993).PISTILLATA(PI) andAPETALA3(AP3) that had obvious expression in f lowers,were recruited for petal identity in Arabidopsis (Winter et al. 2002). B class gene was less expressed in female f lowers as well as for single f lower ofIdesia polycarpa.This is in accordance with the fact that B genes participate in the development of stamens with C genes (Jack et al. 1992;Wang et al. 2011). Moreover, the expression ofAGAMOUS(AG) did not differ signif icantly between male and female f lower buds. Therefore, C genes not only participate in stamen development with B class genes, but also in the formation of pistils. In addition, we mapped anAGL11homeotic D gene which is reported to be related to the development of carpels (Pinyopich et al. 2003). Finally, class E functional genes,SEPALLATA 2(SEP2) andSEPALLATA3(SEP3),(Cui et al. 2010), were up-regulated at early and mid- stages of male development for petal, stamen and carpel development and prevented indeterminate growth of the f lower meristem (Pelaz et al. 2000). However, further studies are needed to determine the effects of these genes on the development of male and female f lowers. Therefore, as the f lower differentiation model ofI. polycarpagenerally conforms to the ABCE model, we can further study various family gene expression at different developmental stages and explore the similarities and differences between the classical ABCDE model and f lower development ofIdesia polycarpa.

Conclusions

In this study, a large number of genes were identif ied by transcriptome sequencing analysis at different stages ofI.polycarpamale and female f lower development. A total of 123335 unigenes were identif ied, and annotated to Nr, Nt,Pfam, KOG/COG, Swiss-prot, KEGG, and GO databases,accounting for 100% of the total sequenced data. A large number of differentially expressed genes were found in the development of male and female f lowers. GO and KEGG pathway analysis of differentially expressed genes indicated that these were associated with photosynthesis- related protein pathways and endoplasmic reticulum protein processing pathway, f lavonoid biosynthesis, and plant hormones.In addition, a number of unigenes related to f lowering and the ABCDE model of f lower development were identif ied in the study. The transcriptome data ofI. polycarpaprovided numerous candidate genes related to the development mechanism of male and female f lowers and laid a foundation for understandingI. polycarpareproductive biology. The transcriptome resources provide a theoretical basis for the regulation of f lower structure in dioecious plants and clearer differences need to be further explored.

Author’s contributionsLin Tang designed the experiment, Fosheng Li, Lanju Mei, Na Li and Min Yao were responsible for the collection of plant materials and the extraction of RNA; Tingting Li and Lanju Mei analyzed the data and wrote the manuscript and Lin Tang critically reviewed the manuscript. All authors read and approved the f inal manuscript.