BACKGROUND
flowering plants display remarkably variable floral architectures, which largely contribute to the diversity of sexual systems in angiosperms  <cit> .
andromonoecy is a rare sexual system, in which an individual bears both staminate and perfect flowers.
in plant kingdom, it has been found in nearly  <dig> plant species, comprising  <dig> % of flowering plants  <cit> .
as a transitional stage from hermaphroditism to monoecy, andromonoecy has independently evolved many times in different plant groups  <cit> .
given the important position in the evolution of plant sexual system, the developmental and genetic mechanisms underlying unisexual flower determination have attracted considerable attention in recent years .
taihangia rupestris yu & li , an endangered perennial herb belonging to the tribe dryadeae, is only disjunctively and sporadically distributed on cliff faces in the southern part of the taihang mountains of china  <cit> .
unlike other genera of this tribe, such as dryas and geum, taihangia always produces both bisexual flowers and unisexual male flowers in the same individual , forming the andromonoecious sexual system  <cit> .
taihangia unisexual male flowers are bisexual at initiation and become unisexual by arresting the pistil development in consequent developmental stages.
during the initiation of floral primordia, temperature is crucial for establishment of the uni- or bisexual flowers.
under low temperature conditions, the pistils of the bisexual flowers tend to abort, resulting in a high proportion of unisexual male flowers  <cit> .
several mads-box genes involved in the floral organ development in taihangia, including b-, c-, d-, and e-class genes, have been identified in previous studies .
of these identified mads genes, class b mads-box gene trpi was expressed in petals and stamens, while class e gene trsep <dig> was strongly expressed in the three inner whorls  <cit> .
for class c gene trag, it was initially expressed in the floral meristem domain that will initiate stamens and carpels, when the stamen primordia are firstly observed.
at the late stage of carpel development, trag is mainly expressed in the ovules, developing styles, and stigmas  <cit> .
despite the progresses mentioned above, underlying molecular mechanisms in the floral organ development and regulation of flower formation in taihangia remain poorly understood, due to limited genetic and sequence information for this non-model plant species.
currently, the next generation sequencing  technology provides opportunities for the efficient and comprehensive analysis of gene expression at whole-genome level in non-model plant without a reference genome  <cit> .
the powerful ngs technology has been widely applied to comparative transcriptome studies of male and female flowers on large-scale gene expression profiles in some plant species, generating considerable genomic data for the identification of functional genes involved in sex differentiation and flower development .
however, the knowledge about unisexual male flower specification and development in andromonoecious plant is still limited, to some extent, due to its infrequent occurrence.
up to now, only a few transcriptome studies have been performed for andromonocious plant species  <cit> .
to investigate the possible molecular mechanisms in the regulation of flower formation within andromonoecious taihangia, we performed de novo transcriptome sequencing of the male and hermaphroditic flowers at both early and late developmental stages by using the illumina platform.
based on the differences in unigene abundance among sequencing libraries, we identified the differentially expressed genes  between the male and hermaphroditic flowers at early and late developmental stages, respectively.
furthermore, go and kegg enrichment analyses were conducted to elucidate the main functions of degs, and to identify candidate genes involved in the regulation of flower formation for further functional analysis.
the present study would provide a genomic resource for gene discovery in the future and shed light on the underlying molecular mechanisms responsible for establishment of unisexual flowers and sex differentiation in andromonocious plant taihangia.
methods
plant material and rna extraction
the plant materials were transplanted from yidoushui  and zhuyufeng , henan, china in april  <dig>  and were grown in the greenhouse of henan polytechnic university.
in early spring  <dig>  the early and late developmental stages of male and hermaphroditic flowers, including male floral bud , hermaphroditic floral bud , male flower , and hermaphroditic flower , were harvested .
as shown in fig.
<dig>  flowers were classified into six developmental stages: young flower buds , elongated buds , pre-anthesis , anthesis , and fully opened flowers .
in this study, the flower samples at early developmental referred to flower buds including young and elongated buds .
the flower samples at late developmental stage were comprised of young flowers at pre-anthesis , flowers at anthesis , and mature flowers with fully opened petals .
the samples of the male and hermaphroditic flowers were collected from six different individuals at each of the developmental stage , respectively, and immediately frozen in liquid nitrogen and stored at −80 °c for rna extraction.fig.
<dig> male and hermaphroditic flowers within andromonoecious taihangia at early and late developmental stages used in rna-seq.
a hermaphroditic floral bud ; .
hermaphroditic flower ;  male floral bud ;  male flower .
early developmental stage: young and elongated buds; late developmental stage: young flowers at pre-anthesis, flowers at anthesis, and fully opened flowers
total rna was isolated from 100 mg homogenized plant material using an rneasy mini kit  for plant tissues according to the manufacturer’s manuals.
rna quality was preliminarily assessed by  <dig> % agarose gel electrophoresis, and then quantified with a uv-visible spectrophotometer .
rna samples for further analysis were selected according to the following criteria: a260/a <dig> ratio was  <dig> - <dig>  and a260/a <dig> ratio ranged from  <dig>  to  <dig> .
the equal amounts of rna from male flower buds at each of the early developmental stages  were mixed as one em sample, while one lm sample was prepared by combinations of the identical amounts of rna from male flowers at each of the late developmental stages .
similarly, eh and lh samples were obtained using the same method described above.
in total,  <dig> rna samples with three biological replicates were prepared for cdna synthesis and qpcr analysis.
after that, specific rna pools of em, eh, lm, and lh were obtained by combining identical quantities of total rna from the corresponding biological triplicates.
quality and quantity of the pooled rna were further analyzed using a nanodrop  <dig> instrument  and a chip rna  <dig> series ii bioanalyzer , respectively.
rna samples with 260/ <dig> ratios , 260/ <dig>  and rin  more than  <dig>  were used to cdna library construction and illumina sequencing.
library construction and rna-seq
construction of the four libraries and the rna-seq were performed by the biomarker biotechnology corporation .
five micrograms of total rna from each sample  were used to construct the sequencing libraries.
rna sequencing libraries were generated using the truseq rna sample prep kit  following the manufacturer’s recommendations.
briefly, the poly -containing mrna molecules were purified from the total rna by using the poly-t oligo-attached magnetic beads .
then, mrna was broken into short fragments, which were used as templates for cdna synthesis.
double-stranded cdna was synthesized using superscript ii, buffer, dntps, rnaseh, dna polymerase i, and random hexamer primers .
after that, the ‘a’ tail was added to the 3′ ends of the repaired cdna fragments and the illumina’s paired-end adapters were ligated to the cdna ends.
the products from the ligation reaction were amplified by pcr.
the pcr productions were purified by magnetic beads  and dissolved in eb solution to generate the sequencing libraries.
the quantity and quality of each sequencing library were tested by the agilent  <dig> bioanalyzer.
finally, the four libraries were sequenced separately on illumina hiseq™  <dig> platform.
sequence assembly and gene annotation
after sequencing, paired-end raw reads were firstly processed through in-house perl scripts  <cit> .
in this step, clean reads were acquired by removing the reads with adaptor contamination, low-quality sequences , and reads with more than 10% q <  <dig> bases.
at the same time, q <dig>  q <dig>  and gc-content of the clean data were calculated.
all the downstream analyses were based on clean data with high quality.
de novo assembly of the transcriptome was performed by trinity software   <cit> , with default parameters and no reference sequence.
in brief, reads were assembled into the contigs by the inchworm program at the first step.
the minimally overlapping contigs were clustered into connected components by the chrysalis program, and then the transcripts were constructed by the butterfly program  <cit> .
finally, the transcripts were further clustered based on nucleotide sequence identity, and the longest transcripts in the cluster units were regarded as unigenes to eliminate redundant sequences.
for detection of transcriptional profiles in the male and hermaphroditic flowers at different stages, clean reads from each library were initially assembled separately.
to obtain a uniform transcriptome reference across samples, all clean reads from four libraries were pooled together and de novo assembled to generate unigenes for assembly evaluation, gene annotation, and expression analysis.
the raw sequencing data were deposited in the ncbi short read archive  database  under the accession number srp <dig>
to identify putative functions of taihangia unigenes, we performed functional annotation by using a blastx search  against a series of protein databases: the national center for biotechnology information non-redundant , swissprot, pfam, clusters of orthologous groups , gene ontology , and kyoto encyclopedia of genes and genomes  databases.
the blast2go software package was used to compare and determine the unigene go annotations  <cit> .
finally, wego software was used to obtain the go functional classifications for all annotated unigenes  <cit> .
transcription factors  prediction
in order to identify the tfs represented in taihangia transcriptome, all assembled unigenes were searched against the plant transcription factor database planttfdb  <dig>  by using blastx with a cut-off e-value of  <dig> e− <dig>  <cit> .
to facilitate the transcription factor prediction, the identified tfs were further screened with similarity ≥  <dig>  and coverage ≥  <dig> .
analysis of differentially expressed genes
to calculate the amount of gene expression for each unigene, rsem  <cit>  was used to quantify the number of reads mapped to the assembled transcriptome, and read count for each gene was obtained from the mapping results.
we used the fpkm  <cit>   algorithm to normalize the gene expressional abundances in each library.
by pairwise comparisons of the four libraries, the degs were identified using the degseq  r package based on the read counts for each gene at different libraries  <cit> .
in addition, the false discovery rate  control method  <cit>  was used to identify the threshold of the p-value in the significance tests.
significance of differential expressed genes was defined at a false discovery rate <  <dig>  and an absolute value of log2-ratio ≥  <dig>  the classification of degs were performed with go and kegg analysis using the method described above.
go enrichment analysis of the differentially expressed genes  was implemented by using the cytoscape plugin bingo using a hypergeometric test with a corrected p-value <  <dig>   <cit> .
kegg pathway enrichment analysis of the degs was performed by using kobas with the hyper-geometric distribution model  <cit> .
the enrichment p-values were adjusted by using the benjamin and hochberg method.
validation of sequencing data by quantitative real-time pcr
sixteen degs were randomly selected and the expression profiles were investigated by qrt-pcr to confirm the transcriptome data.
first-strand cdna was synthesized from the same rna samples used for library construction by using the primescript rt reagent kit with gdna eraser   according to the instructions.
quantitative real-time pcr  was performed by using the ssofast evagreen supermix rt-pcr kit  and the miniopticon real-time detection system .
the pcr mix was composed of 10 μl evagreen supermix,  <dig>  μl of 1: <dig> diluted cdna,  <dig>  μl of each primer , and  <dig>  μl water in a final volume of 20 μl.
the reactions were incubated under following cycling conditions: 2 min at 50 °c, 2 min at 95 °c, 40 cycles of 95 °c for 15 s, 58 °c for 15 s, and 72 °c for 30s, and finally 72 °c for 2 min.
after the final cycle, a melting curve analysis was performed with a single cycle from 60 °c to 95 °c in 5 s intervals to verify the reaction specificity.
the efficiency of the primer sets was calculated by performing real-time pcr on serial dilutions of first-strand cdnas templates  <cit> .
by using the genorm software  <cit> , ef1α and ubq were identified as reference genes for stable expression across developmental stages.
the relative expression of each gene was calculated using the standard e-ΔΔcq method  <cit> , and normalized with the internal standard gene ef1α and ubq.
expression quantification and data analysis were performed according to miqe guideline suggested by bustin et al.
<cit> .
three independent biological replicates and three technical repetitions were performed for each of the quantitative pcr experiments.
the relative expression of each target gene was analyzed by one-way anova.
statistical analyses of the data were performed at p <  <dig>  significance level using spss  <dig>  software.
RESULTS
illumina sequencing and de novo assembly
to obtain a comprehensive flower transcriptome in andromonoecious taihangia, four cdna libraries  were constructed from male and hermaphroditic flowers covering early and late development stages.
paired-end sequencing of the four constructed libraries were performed on an illumina hiseq  <dig> platform, generating a total of  <dig> , <dig> raw reads.
after filtering for low quality region, adapters, and sequencing tags, we obtained  <dig> , <dig> qualified illumina reads , including  <dig> , <dig> ,  <dig> , <dig> ,  <dig> , <dig> , and  <dig> , <dig>  clean reads for the libraries of em, eh, lm, and lh, respectively .
overall, the library produced  <dig> g base pairs clean data with  <dig> % q <dig> bases , indicating that the throughput sequencing was accurate enough to warrant the further analysis.
the clean reads from each libraries were separately de novo assembled by trinity programs, generating  <dig>  ,  <dig>  ,  <dig>  , and  <dig>   unigenes, respectively, with mean lengths of  <dig> to 905 bp.
after that, we pooled all high quality reads from the four individual libraries to perform the de novo assembly for a uniform reference transcriptome.
by using the trinity program, a total of  <dig>  unigenes longer than 200 bp were assembled with an n <dig> length of  <dig>  bp, an average length of 840 bp, and a maximal transcript length of  <dig>  bp.
the length distribution of taihangia unigenes was shown in additional file 2: figure s <dig> table  <dig> summary of illumina transcriptome sequencing and assembly for taihangia rupestris flower
to evaluate the quality of de novo assembly without a reference genome, a common method is to calculate the percentage of reads that can be mapped back to unigene sets.
based on this metric, we found that more than 84% of clean reads were mapped back to the taihangia unigenes for each library .
in addition to in silico quality assessment, the standard pcr and sanger sequencing were conducted to evaluate the quality and accuracy of taihangia assembled unigenes.
of  <dig> designed primer pairs,  <dig> primer pairs for  <dig> unigenes could steadily generated expected pcr products.
all  <dig> pcr products were sequenced by sanger sequencer, and the obtained sequences were compared with the assembled unigenes.
the alignment analysis showed that the sequences from sanger sequencing met well with the assembled unigenes with an average sequence identity of  <dig> %, confirming the accuracy of assembled unigenes derived from transcriptome .
together, our results supported the high quality of de novo assembly in current sequencing data.
functional annotation and classification of the unigenes
all the assembled  <dig>  unigenes were subjected to blastx comparisons, setting a cut-off e-value of 10− <dig>  against the public databases including nr, pfam, swiss-prot, go, kegg and cog databases.
as a result, a total of  <dig>  unigene sequences showed significant similarity to known proteins in publicly available databases .
based on the annotation results of the nr database, the e-value distribution analysis showed that  <dig> % of the matched sequences had strong homology with the e-value <  <dig> e− <dig>  while  <dig> % of the matched sequences ranged between 1e− <dig> and 1e− <dig>  the similarity distribution analysis revealed that  <dig> % of the unigenes had a similarity higher than 75%, while  <dig> % of mapped sequences had a similarity ranging from 60% to 75%.
for species distribution, the nr database queries showed that  <dig> % of the t. rupestris annotated sequences were matched to the sequences of fragaria vesca, followed by prunus mume , prunus persica , malus domestica , and pyrus × bretschneideri .
characteristics of the homology search of taihangia unigenes against the nr database were shown in additional file 5: figure s <dig> table  <dig> summary for the annotation of unigenes of taihangia flower
based on the blastx results against the nr database, we added go terms to obtain go functional annotations and categorizations of these assembled unigenes.
of the  <dig>  unigenes with nr annotation, a total of  <dig>  unigenes with at least one go term were assigned into biological processes , cellular components , and molecular functions  categories.
the three main categories were further classified into  <dig> functional groups including  <dig> biological processes,  <dig> cellular components and  <dig> molecular functions.
in the category of molecular function, “catalytic activity” was the most highly represented group, followed by group of “binding”.
under the biological process category, “metabolic process” and “cellular process” were the most highly represented groups.
for the cellular component, the major groups were “cell” and “cell part”.
the go analysis indicated that a great number of identified unigenes were associated with various biological processes and molecular functions in taihangia floral tissues.
the annotated sequences were further subjected to a search against the cog database for functional prediction and classification.
as a result,  <dig>  unigenes were grouped into  <dig> cog classifications, among which the “general function prediction only” represented the largest group , followed by “replication, recombination and repair” , “transcription” , “signal transduction mechanisms” , and “posttranslational modification, protein turnover, chaperones” .
the category “nuclear structure” with  <dig> unigenes was the smallest group.
to obtain a better understanding of biological functions of the unigenes, the annotated sequences were searched against the kegg database.
among the  <dig>  annotated unigenes,  <dig>  had significant matches and were assigned to  <dig> kegg pathways.
the top five pathways were “ribosome” , “protein processing in endoplasmic reticulum” , “carbon metabolism” , “biosynthesis of amino acids”  and “spliceosome” .
functional go, cog, and kegg pathway annotations and classification of the unigenes were shown as additional file  <dig>  in summary, all these annotation and classification analyses would provide valuable information for gene discovery and functional genomic studies in the future.
identification of transcription factors
in the present study, we identified  <dig>  tfs, representing  <dig> % of taihangia unigenes and falling into  <dig> tf families classified by plant transcription factor database planttfdb  <dig>  .
among these detected tf gene families, bhlh  was the most abundant tf family, followed by myb_related , erf , nac , wrky , and c2h <dig> .
notably, a total of  <dig> mads-box genes, including  <dig> mikc-type and 17 m-type tfs, were found in taihangia based on transcriptome data.fig.
<dig> identification of transcription factors  and mads-box genes in taihangia.
a classification of tf families.
b phylogenetic tree of mikc mads-box proteins of taihangia, arabidopsis thaliana, strawberry, and peach.
red dots, taihangia genes; green dots, strawberry genes; blue dots, a. thaliana genes; purple dots, peach genes
to examine the phylogenetic relationships among the taihangia mads-box genes and to classify them into the established subfamilies, a neighbor–joining  phylogenetic tree was constructed based on a multiple sequence alignment of the predicted full-length mads-box protein sequences of f. vesca, p. persica, and arabidopsis .
the phylogenetic analysis revealed that  <dig> mads-box tfs, including eight mads-box genes reported in previous studies, were classified into classified into  <dig> mikcc and one mikc* subfamilies .
furthermore, several mads-box tfs, such as apetala <dig> , pistillata , agamous , shaterproof , and sepallata , described in abcde model were also included.
trumads <dig>  trumads <dig>  and trumads <dig> belonged to the ap1/ful subfamily in the a class, which included the ap1and ful clades.
trumads <dig> was classified into ap <dig> clade, while trumads <dig> and trumads <dig> were sorted into ful clade.
in the functional b-class  group, two genes  belonged to pi clade, while the other one gene  fell within tm <dig>  three mads-box tfs, trumads <dig>  trumads <dig>  and trumads <dig>  were assigned to the ag family, which were functionally classified as a c/d class group.
trumads <dig> belonged to the c-class group ag while trumads <dig> and trumads <dig> were most orthologous to the d-class functional group shp and stk, respectively.
in the e functional group, two mads-box tfs, trumads <dig> and trumads <dig>  were assigned to sep subfamily.
in general, the phylogenetic analysis showed that taihangia contains mads proteins in each putative functional group.
transcriptome profiles in male and hermaphroditic flowers
differences of gene expression profiles between male and hermaphroditic genders at different developmental stages were examined by pairwise comparisons of the four libraries .
a total of  <dig>  unigenes  were shared among all libraries, while approximate  <dig> % of unigenes  was detectable in individual libraries .
for hermaphroditic genders, there were  <dig>  and  <dig>  unigenes exclusively expressed in eh and lh, respectively.
meanwhile, we detected  <dig>  unigenes specific for em and  <dig>  unigenes only for lm.
a total of  <dig>  unigenes  were gender-specific, with  <dig>  unigenes detected only in the male flowers  and  <dig>  unigenes in hermaphroditic flowers , respectively.
in addition,  <dig> % of assembled unigenes  was only detected in floral buds , whereas  <dig> % of assembled unigenes  was observed in the mature flowers , but not in the floral bud samples.
overall, our results suggested that the unigenes specific expressed in the male or hermaphroditic flower might be potential candidates to be associated with the regulation of flower formation in andromonoecious taihangia.
differential expressed genes  analysis
based on the four comparisons of em vs eh, lm vs lh, em vs lm, and eh vs lh, we identified a total of  <dig>  differential expression genes .
among them,  <dig>  degs between gender types  and  <dig>  degs between developmental stages  were detected, indicating that majority of degs occurred among different developmental stages .
in order to identify related genes that involved in the regulation of flower formation,  <dig>  degs between male and hermaphroditic genders were further analyzed.
at the early developmental stage, we identified  <dig> up-regulated genes and  <dig> down-regulated genes expressed in eh.
on the other hand,  <dig> up-regulated genes and  <dig> down-regulated genes were detected in lh at the late developmental stage.
interestingly,  <dig> degs, including  <dig> hermaphroditic-biased and  <dig> male-biased unigenes, were differentially expressed in both early and late developmental stages, indicating that the gene expression pattern of the degs was maintained from the floral bud to the mature flower across developmental stages.fig.
<dig> identification of differentially expressed genes  in taihangia.
venn diagram of the number of degs based on four comparisons of em vs eh, lm vs lh, em vs lm, and eh vs lh.
total deg numbers in each of comparison were given outside the circles, and arrows indicate up- or down-regulation
among the  <dig>  tfs identified in this study,  <dig> tf genes were differentially expressed between male and hermaphroditic flowers.
these tf genes fall into  <dig> tf families classified by the planttfdb  <dig> , with the bhlh family  representing the largest amount of differentially expressed tfs, followed by erf , myb related , wrky , and nac .
tf genes showed different expression patterns in male and hermaphroditic flowers during flower development , suggesting possibly different roles of tfs in the early and late developmental stages.
at the early developmental stage,  <dig> putative tf genes belonging to  <dig> tf families were differentially expressed in comparison em vs eh.
the  <dig> tf genes belonging to  <dig> tf families, such as bhlh, myb related, zf-hd, nac, and erf, were significantly upregulated, while the remaining  <dig> tf genes belonging to  <dig> tf families  were downregulated in em.
in the late developmental stage, the identified  <dig> putative tf genes were classified into  <dig> tf families based on the comparison lm vs lh.
of them,  <dig> tf genes were detected with significantly upregulated expression in lh.
these tf genes might be important regulators contributing to the establishment and development of male and hermaphroditic flower within the andromonoecious system.
the expressed profiles of tfs derived from degs were shown as fig.
4b.fig.
<dig> the differential expressed transcription factors between the male and hermaphroditic flowers.
a the number of degs belonging to different transcription factor families detected in two comparisons of em vs eh and lm vs lh; .
expression pattern of differential expressed tfs between the male and hermaphroditic flowers
in summary, our results showed that the number of degs between male and hermaphroditic flowers remarkably increased alongside the flower development, and quite different deg profiles were observed between comparisons em vs eh and lm vs lh.
given that the morphogenesis of the taihangia uni- and bisexual flowers occurred in the early flower developmental stages, degs between em and eh were expected to be involved in sex differentiation, while degs between lm and lh were probably associated with floral organ development and reproduction.
go and pathway analyses of degs
go functional classifications of the degs between male and hermaphroditic samples were conducted by using the wego website.
of the identified  <dig> degs between gender types ,  <dig> degs were classified into  <dig> go subterms.
in the cellular components category, the cell , cell part , membrane , and organelle , were the most abundant subterms.
the catalytic activity  and the binding  subterms were dominant in the molecular function category.
for the biological process category, a high percentage of degs was assigned to the cellular process , metabolic process , and single-organism process  subterms .fig.
<dig> go classifications of degs between the male and hermaphroditic flowers in different developmental stages.
degs are grouped to the secondary classification of hierarchical go terms.
a go classification of degs in two comparisons of em vs eh and lm vs lh at early and late developmental stages; .
go classification of upregulated and downregulated degs in comparison em vs eh at early developmental stage; .
go classification of upregulated and downregulated degs in comparison lm vs lh at late developmental stage
fig.
<dig> top  <dig> enriched kegg pathways of degs in em vs eh and lm vs lh between male and hermaphroditic flowers.
number at the right of each column indicates the number of degs included in corresponding pathway.
the pathways were divided into five categories: “metabolism”, “cellular processes”, “environmental information processes”, “genetic information processes”, and “organismal systems”
to better understand the regulation mechanisms of flower formation across developmental processes in taihangia, we put emphasis on go functional classifications analysis of upregulated and downregulated degs in comparisons em vs eh and lm vs lh, respectively.
furthermore, go enrichment analysis was performed by using the hypergeometric test in bingo program with a corrected p value ≤  <dig>   <cit> .
at the early developmental stage,  <dig> degs were classified into  <dig> go subterms.
one hundred and seventy three upregulated genes and  <dig> downregulated genes in em were classified into  <dig> and  <dig> go subterms, respectively .
overrepresentations of go terms in em were found in “pollen tube development” , “regulation of pollen tube growth” , “hydrolase activity, acting on glycosyl bonds” , indicating that these genes involved in stamen development.
for up-regulated degs in eh, most of the overrepresented go terms were related to cell cycle, developmental process, and metabolic process, such as “regulation of macromolecule metabolic process” , “regulation of cell cycle” , “anatomical structure development” , “integral component of membrane” , and “regulation of transcription, dna-templated” .
in the late developmental stage,  <dig> go subterms were found for  <dig> degs in comparison lm vs lh.
one hundred and twenty six upregulated genes and  <dig> downregulated genes were classified into  <dig> and  <dig> go subterms, respectively .
the overrepresentations of go terms in lm were observed in “oxidoreductase activity”  while overrepresentations of go terms in lh mainly included “cell cycle process” , “dna methylation” , “flower development” , and “dna binding” .
in particular, a large number of go terms, such as “regulation of cell cycle”, “regulation of transcription, dna-templated”, and “regulation of developmental process”  overrepresented in early developmental stages were also enriched at late developmental stages.
the go enrichment analysis results are presented in additional file  <dig>
to further investigate the biochemical pathways of these degs, we mapped the  <dig>  degs to terms in the kegg database.
among them,  <dig> genes had a ko id and could be categorized into  <dig> pathways, of which “phenylpropanoid biosynthesis” , “plant-pathogen interaction” , “plant hormone signal transduction” , “pentose and glucuronate interconversions” , “fatty acid metabolism” , and “unsaturated fatty acid metabolism”  had the most deg numbers .
in addition, we conducted an enrichment analysis with the kegg pathway, which assigned the degs to  <dig> pathways such as biosynthesis of unsaturated fatty acids, phenylpropanoid biosynthesis, and pentose and glucuronate interconversions .
in particular, two pathways “unsaturated fatty acid metabolism”  and “dna replication”  were highly represented, respectively, which were consistent with the go enrichment analysis.fig.
<dig> expression pattern of functional degs involved in sex differentiation and flower organ development.
a flower development;  epigenetic processes;  cell proliferation.
the bar represents the scale of relative expression levels of degs.
red means high expression and blue means low expression
go and kegg pathway analysis revealed that degs between the two genders were associated with a wide range of functions, and participated in many biological processes, such as cell proliferation, chromatin modification, and flower development, suggesting that the formation of uni- and bisexual flowers might be controlled by many genes involved in numerous biological processes.
moreover, the identified go subterms and kegg pathways could provide further information on the molecular mechanism of sex differentiation in andromonoecious taihangia.
degs involved in regulation of taihangia male and hermaphroditic flower formation
flower development
based on deg analysis and gene expression profiles of male and hermaphroditic flowers, we identified a serial of degs involved in flower development .
among them, several genes related to floral organ development, such as mads-box, wuschel , yabby <dig> , and lrr-protein kinase bam <dig>  were up-regulated in hermaphroditic flowers.
particularly, four differentially expressed mads-box genes were found in the ag, b-sister , and mikc* subfamilies.
among them, both trumads <dig>  and trumads <dig>  were up-regulated in hermaphroditic flowers at late developmental stages.
trumads <dig> , belonging to d-class functional genes, showed dramatically higher expression in hermaphroditic than in male flowers throughout floral development.
our results indicated that these mads-box genes preferably expressed in lh were likely associated with gynoecium development in hermaphroditic flowers.
conversely, trumads <dig>  belonging to mikc* subfamilies, was upregulated in em at early developmental stages, suggesting that the mads-box tf probably involved in stamen development.
epigenetic process
epigenetic mechanisms, including dna methylation, histone modifications, chromatin conformation, and rna interference, play essential roles in modulating gene expression without altering the dna sequence.
in this study, a number of degs involved in dna methylation and histone modification were identified in taihangia flower at late developmental stage .
several genes, including three dna methyltransferase cmt <dig>  met, and drm <dig>  involved in dna methylation process showed differentially expression between male and hermaphroditic flowers.
on the other hand, histone modifications related genes, such as histone-lysine n-methyltransferase atxr <dig>  ashr <dig>  and suvr <dig>  were differentially expressed in comparison lm vs lh at late developmental stages.
no significant differences were observed in the expression levels of epigenetic related genes between male and hermaphroditic flowers at the early developmental stages, while all of them were significantly upregulated in lh at late developmental stages.
notably, one homologous of ago  <dig> gene involved in rna silencing and rna-directed dna methylation  pathway exhibited upregulated expression patterns in lh at late developmental stages.
our finding indicated that epigenetic processes may extensively occur in hermaphroditic flowers, especially at late developmental stages.
cell cycle process
a large number of degs involved in cell cycle process, such as dna replication, cell division, and regulation of cell cycle, were identified based on transcriptome data, and these degs showed significantly increased expression in hermaphroditic flowers at late developmental stage .
a number of degs, including mitotic spindle checkpoint protein mad2-like, atp-dependent dna helicase q-like, dna topoisomerase 3-alpha-like, and mitotic checkpoint serine/threonine-protein kinase bub <dig> beta-like were assigned into dna replication functional group.
meanwhile, it was found that majority of degs, belonging to cell division functional group, were down-regulated in lm at late developmental stage, including 65-kda microtubule-associated protein, g2/mitotic-specific cyclin-1-like, condensin complex subunit 2-like, dna replication licensing factor mcm2-like, and dna replication licensing factor mcm <dig>  similarly, several degs involved in cell cycle group also showed down regulated expression in male unisexual flowers, such as microtubule-associated protein rp/eb family member 1c, dna cross-link repair protein snm1-like, protein fizzy-related 3-like, and dna cross-link repair protein snm1-like.
in particular, many degs, including g2/mitotic-specific cyclin-1-like, protein fizzy-related 3-like, mitotic spindle checkpoint protein mad2-like, mitotic checkpoint serine/threonine-protein kinase bub <dig> beta-like, cyclin-a1-1-like and cyclin-d3-1-like, showed significantly upregulated expression in both eh and lh across the developmental stages.
overall, our results indicated that the cell cycle process was probably severely disturbed in male unisexual flowers.
biosynthesis of unsaturated fatty acids
unsaturated fatty acids are considered to be essential membrane components, and play key roles in many cellular events.
the synthesis of unsaturated fatty acid is a complicated process, requiring for involvement of a series of enzymes, such as palmitoyl-monogalactosyldiacylglycerol Δ <dig> desaturase  and acyl- desaturase  <dig> .
in the present study, several degs encoding several key enzymes invovled in the biosynthesis of unsaturated fatty acids pathway, including six fad <dig>  three very-long-chain 3-oxoacyl-coa reductase 1-like genes , and one acp <dig> genes, were found to be significantly downregulated in male flowers, indicating that biosynthesis of unsaturated fatty acids was likely restricted in male flowers.
validation of unigenes and gene expression profiling using rt-qpcr
to confirm the gene expression patterns identified by the rna-seq data, the transcript levels of  <dig> random selected degs were examined by quantitative pcr.
we performed qrt-pcr using these  <dig> genes with  <dig> biological replicates from each group  .
the detailed information of the selected degs, ids and primer pairs used in this study were shown in additional file 12: table s <dig>  it was found that all of  <dig> selected degs were successfully amplified with single bands with the expected sizes, and the patterns of gene expression detected by rt-qpcr were consistent with those from rna-seq data .
therefore, rt-qpcr experiments confirmed that the degs obtained from the assembled transcriptome were accurate and gene expression profiles from rna-seq data were reliable.fig.
<dig> validation of  <dig> randomly selected differentially expressed genes  derived from rna-seq using real-time quantitative rt-pcr.
quantitative gene expression data represented as the mean ± sd are shown on the left.
the ef1α and ubq were used as reference genes.
statistically significant differences of relative expression derived from one-way anova at the level of p <  <dig>  were marked as different lowercase letters.
the normalized fpkm, transformed by log <dig> , were indicated on the right
discussion
taihangia transcriptome profiles
as an andromonoecious plant species, taihangia provides an excellent model system within a uniform genetic background for the study of sex differentiation and establishment of unisexual flowers.
in this study, we conducted de novo transcriptome sequencing to compare transcriptional profiles of the male and hermaphroditic flowers by using the illumina rna-seq method.
to the best of our knowledge, this is the first report of transcriptome data in this rare cliff plant species with the andromonoecious sexual system.
a total of approximately  <dig> , <dig> qualified illumina reads  were obtained from four floral sequencing libraries and then assembled into  <dig>  unigenes.
of the assembled unigenes,  <dig>  unigenes  were annotated in the public database, while the function of other  <dig>  assembled unigenes  could not be identified.
given that the genomic and transcriptomic information is limited for taihangia, the lineage-specific genes associated with the unique floral phenotype may not be contained in function annotation, implying that potential novel genes were discovered from our transcriptome data.
by comparative analysis of male and hermaphroditic transcriptome data, we identified a large number of candidate degs involved in various biological processes, which is likely associated with the regulation of male and hermaphroditic flower formation.
the transcriptome data provide an important resource for the understanding of specific biological processes involved in sex differentiation and facilitating gene discovery in the future studies in andromonoecious taihangia.
floral organ identity-related mads-box genes may be crucial for flower organ formation
the mads-box transcription factors, initially identified as floral homeotic genes, are the major components in the well-known ‘abcde’ model that describes their roles in the determination of floral organ identity  <cit> .
based on taihangia transcriptome data, we identified  <dig> mads-box tfs, including several mads proteins in each putative functional group described in abcde model.
according to the “abcde” model  <cit> , similar levels of expression of class a-, b-, and e-class mads-box genes between male and hermaphroditic genders were expected, while class c- and d- mads genes should be upregulated expressed in the hermaphroditic flowers.
in the present study, most of mads-box tfs including a-, b- and e-class of genes did not show differential expression between male and hermaphroditic flowers.
as expected, two mads-box tfs, truag  and trustk , showed upregulated expression in hermaphroditic flowers.
however, it was noteworthy that upregulated expression of truag was not found in eh but in lh, suggesting that its expression pattern appears to follow the arrest of pistil development in male flowers.
a number of studies have found that the mads box genes did not play a part in sex determination in some dioecious and monoecious species that achieve unisexuality by organ abortion  <cit> .
in taihangia, pistils are initiated and then arrest in unisexual male flower, the molecular mechanism of sex determination must be acting downstream of organ identity.
thus, the transcriptional profile of several mads-box tfs more expressed in hermaphroditic flowers might be the consequence rather than the cause of sex differentiation in andromonoecious taihangia, and the mads-box tfs were related to pistil development in hermaphroditic flowers.
epigenetic-related genes are likely associated with reproduction in hermaphroditic flower
epigenetic processes such as dna methylation, histone modification, and small rna, play important roles in the regulation of gene expression during flower development  <cit> .
as key enzymes involved in dna methylation processes, dna methyltransferases and histone h3k <dig> methyltransferase are crucial for establishment and maintenance of dna methylation  <cit> .
in our study, three dna methyltransferases  and three histone h3k <dig> methyltransferase  showed preferable expression in lh, indicating their possible roles in taihangia reproduction at the late developmental stage.
in addition, several genes, such as argonaute  <dig> , suppressor of gene silencing  <dig> , rna-dependent rna polymerase  <dig> , dicer homolog 3a-like , and rpd <dig> , possibly involved in rna-induced silencing and rna-directed dna methylation  pathways were identified based on taihangia transcriptome data.
particularly, ago  <dig> showed approximately three folds higher expression in lh than in lm at late developmental stage.
as key regulators, ago proteins and dcl <dig>  along with rdr, are crucial for the epigenetic silencing  <cit> .
based on gene expression analysis, three genes, dcl, rdr, and rpd, showed more than two times higher expressions in lh than in lm, indicating that extensive small rna related processes possibly occurred in taihangia hermaphroditic flower at developmental stage.
currently, increasing evidences have supported that epigenetic processes played important roles in regulation of plant reproductive development .
for example, galla et al.
<cit>  reported that small rna generation related genes, including dcl <dig>  ago, sgs <dig>  and rdr, were involved in plant reproduction, such as cell identity and fate in the ovule, in hypericum perforatum.
the upregulated expressions of these epigenetic-related genes in taihangia hermaphroditic flowers, especially at the late developmental stages, may indicate that the epigenetic regulation played important roles in plant reproduction, such as female gametophyte formation.
the roles of epigenetic processes in t. rupestris are still unclear, and future studies are needed to investigate the roles of associated genes involved in epigenetic processes.
disturbance of cell cycle may contribute to pistil abortion in male flower
during flower development, arrest of pistil or stamen development is a normal phenomenon in monoecious and dioecious plant .
several biological processes, such as programmed cell death  and cessation of the cell division, are critical in establishment of unisexual flower  <cit> .
based on taihangia rna-seq data, a great number of degs involved in cell cycle, dna replication, and cell division were significantly downregulated in male flowers over the developmental stages, indicating that these biological processes played important roles in arrest of pistil development in taihangia male flowers.
for many plant species, pcd is essential to many aspects of plant morphogenesis, and it is a common process for the developmental arrest of reproductive organs in unisexual flowers  <cit> .
however, we did not found that pcd involved in pistil abortion in taihangia male flowers based on transcriptome data.
it has been reported that developmental arrest of sterile sex organs and the subsequent unisexuality of flowers could resulted from a cessation of cell division  <cit> .
similarly, down-regulated expressions of related genes involved in cell division probably led to restriction of cell proliferation, which may be closely associated with the arrest of pistil development in taihangia unisexual male flower.
as primary regulatory genes, cyclins play critical roles in controlling cell cycle progression, including dna replication, the g2/m transition, and mitosis  <cit> .
in this study, several cyclin genes  related to regulation of cell cycles were found to be uniformly down-regulated in unisexual male flowers.
particularly, cyclin-a1-1-like and cyclin-d3-1-like significantly decreased expression level in em and lm across developmental stages, indicating that there was a potential mechanistic link between the regulation of cell cycle and the establishment of unisexual flowers in taihangia.
during flower development, cyclin-a1and cyclin-d <dig> were considered as the core regulatory genes involved in regulation of flower primordia information and meristem activity  <cit> .
for example, cyclin-d <dig> could define distinct developmental zones, acting as a developmental boundary between lateral organs, in the shoot apical meristems by controlling the balance between cell division and cell differentiation  <cit> .
in taihangia unisexual male flowers, downregulated expression of cyca1and cycd <dig> probably resulted in the disturbance of cell cycle in gynoecium primordia, and subsequently determined arrest of pistil development.
these results suggested that the core regulatory genes involved in cell cycle progression probably contribute to pistil abortion in taihangia unisexual flowers by regulation of the cell cycle.
cold-responsive genes are likely to affect the formation of unisexual flowers
it is well known that the establishment of unisexual flower is influenced by environmental factors such as temperature, light intensity, photoperiod, and nutritional conditions via specialized gene expression .
for andromonoecious taihangia, temperature during the initiation of floral primordia is crucial for determining morphogenesis of gender type, and generation of unisexual male flowers is attributed to pistil abortion under low temperature  <cit> .
to deal with adverse cold stress, plants have developed versatile defense mechanisms to acquire cold tolerance by adjustments in molecular, and physiological, and developmental processes  <cit> .
for example, increase of membranes’ unsaturated fatty acid content could promote the stabilization of plasma and chloroplast membranes, and thus improve resistance to cold stress  <cit> .
as key enzymes involved in biosynthesis of unsaturated fatty acid, fad <dig>  acpr, and acp  <dig> were dramatically down-regulated in em, indicating that restriction of polyunsaturated fatty acids biosynthesis probably occurred in taihangia male flowers.
under low temperature condition, a number of cold-regulated genes such as lea family genes could be induced in higher plants  <cit> .
for instance, lea proteins are thought to play a role in membrane stabilization by interacting with membrane surfaces or with water molecules under cold stresses  <cit> .
based on transcriptome data, a number of cold-resistance genes, including lea dc <dig>  lea b <dig> , ascorbate oxidase, and peroxidase, were significantly down-regulated in male flowers.
these results agreed well with a recent study on comparative transcriptome analysis between male and hermaphroditic flowers in japanese apricot by shi et al.
<cit> , arguing that reduction of cold-resistance gene expression in male flowers is likely associated with pistil abortion under low temperature.
taken together, these results suggested that down-regulated expression of cold-responsive genes at the early stage is likely to affect the formation of unisexual male flowers in andromonoecious taihangia.
candidate tfs as important regulators in establishment of male and hermaphroditic flowers
in plant, tf genes play an important role in many aspects of plant growth and development  <cit> .
in the present study, we identified a large number of tf genes showing significantly differential expression between male and hermaphroditic flowers.
these tfs are expected to be involved in the regulation of formation of the male and hermaphroditic flowers in taihangia.
of these identified tfs, most stress-responsive tfs, such as nac, wrky, bhlh, gras, and erf showed increased expression in male flowers, implying that male flowers probably suffered more severe injuries than hermaphroditic flowers.
majority of upregulated tfs in bisexual flowers, including wox, yab, ap <dig>  and c2h <dig>  were associated with flower organ development.
for example, wuschel , a number of wox tf genes, is known to maintain the stem cell activity in the floral meristems  <cit> .
after floral induction, wus can induce ag expression and promote the up-regulation of ag in developing flowers  <cit> .
in our study, higher expression of one homolog of wus was not observed in hermaphroditic flowers at early developmental stages.
similarly, truag showed no difference in expression between em and eh, suggesting that a possible different regulatory mechanism occurred in taihangia flower development.
the yabby gene family is involved in adaxial–abaxial polarity in lateral organs, and participates in a diverse range of processes such as floral patterning, organ growth, and maintenance of meristem organization  <cit> .
one homologous of yabby <dig> genes identified in this study was upregulated in both eh and lh, indicating that yabby tfs were closely associated with the formation of bisexual flowers.
interestingly, five genes encoding growth-regulating factor  were detected and showed uniformly increased expression in hermaphroditic flowers.
grfs are plant-specific transcription factors, which play important roles in the control of cell proliferation via cell cycle regulation and central developmental processes including flower formation  <cit> .
recently, liang, et al.
reported that dramatically reduced expression levels of grf caused pistil abnormalities and suggested an important role of grf in controlling carpel number and pistil development  <cit> .
our result was consistent with the finding in arabidopsis, and confirmed that disturbance of cell proliferation in male flowers contributed to the pistil abortion.
CONCLUSIONS
de novo transcriptome sequencing of andromonoecious taihangia was performed by using the illumina paired-end sequencing technology in the present study.
in total,  <dig>  unigenes were obtained from the male and hermaphroditic flowers covering early and late developmental stages, and  <dig>  unique sequences were annotated in the public database.
focusing on the regulatory mechanism of flower formation in andromonoecious taihangia, we identified  <dig>  degs between male and hermaphroditic flowers.
go and pathway analysis showed that these degs involved in a wide range of biological processes, such as flower development and reproduction, cell cycle, epigenetic processes, and transcriptional regulation, and provided numerous candidate genes possible related to sex differentiation for future functional analysis.
the transcriptome data of taihangia could be helpful to improve understanding of the molecular mechanisms in sex differentiation and unisexual flower establishment, and would provide theoretical basis for the regulation of flower formation in this andromonoecious plant species.
additional files
additional file 1: figure s <dig>  male and hermaphroditic flowers produced within the same individual of taihangia on the cliff face.
additional file 2: figure s <dig>  length distribution of the unigene sequences.
additional file 3: table s <dig>  summary for the alignment of reads to unigene libraries.
additional file 4: table s <dig>  validation of assembled unigenes by using sanger sequencing.
additional file  <dig> figure s <dig>  characteristics of the homology search of taihangia unigenes against the nr database.
additional file 6: functional go, cog, and kegg pathway annotation and classification of the taihangia unigenes.
additional file 7: identification of transcription factors against plntfdb database.
additional file 8: table s <dig>
taihangia mads-box gene names and attributes.
additional file 9: figure s <dig>  venn diagram of the number of specific and shared unigenes in em, eh, lm, and lh samples.
additional file 10: over-representative go terms of degs in em vs eh, and lm vs lh, respectively.
additional file 11: pathway classification and enrichment analysis of the degs between the male and hermaphroditic flowers.
additional file 12: table s <dig>  the selected degs, ids, and primer pairs for rt-qpcr.
abbreviations
cogclusters of orthologous group
degdifferentially expressed gene
fpkmfragments per kilobase of exon per million fragments mapped
gogene ontology
keggkyoto encyclopedia of genes and genomes
qrt-pcrquantitative reverse transcription polymerase chain reaction
rddmrna direct dna methylation
rna-seqribonucleic acid sequencing
taihangia rupestris
taihangia
tftranscription factor

 flowering plants display remarkably variable floral architectures, which largely contribute to the diversity of sexual systems in angiosperms. andromonoecy is a rare sexual system, in which an individual bears both staminate and perfect flowers. in plant kingdom, it has been found in nearly plant species, comprising % of flowering plants. as a transitional stage from hermaphroditism to monoecy, andromonoecy has independently evolved many times in different plant groups. given the important position in the evolution of plant sexual system, the developmental and genetic mechanisms underlying unisexual flower determination have attracted considerable attention in recent years. taihangia rupestris yu & li, an endangered perennial herb belonging to the tribe dryadeae, is only disjunctively and sporadically distributed on cliff faces in the southern part of the taihang mountains of china. unlike other genera of this tribe, such as dryas and geum, taihangia always produces both bisexual flowers and unisexual male flowers in the same individual, forming the andromonoecious sexual system. taihangia unisexual male flowers are bisexual at initiation and become unisexual by arresting the pistil development in consequent developmental stages. during the initiation of floral primordia, temperature is crucial for establishment of the uni- or bisexual flowers. under low temperature conditions, the pistils of the bisexual flowers tend to abort, resulting in a high proportion of unisexual male flowers. several mads-box genes involved in the floral organ development in taihangia, including b-, c-, d-, and e-class genes, have been identified in previous studies. of these identified mads genes, class b mads-box gene trpi was expressed in petals and stamens, while class e gene trsep was strongly expressed in the three inner whorls. for class c gene trag, it was initially expressed in the floral meristem domain that will initiate stamens and carpels, when the stamen primordia are firstly observed. at the late stage of carpel development, trag is mainly expressed in the ovules, developing styles, and stigmas. despite the progresses mentioned above, underlying molecular mechanisms in the floral organ development and regulation of flower formation in taihangia remain poorly understood, due to limited genetic and sequence information for this non-model plant species. currently, the next generation sequencing technology provides opportunities for the efficient and comprehensive analysis of gene expression at whole-genome level in non-model plant without a reference genome. the powerful ngs technology has been widely applied to comparative transcriptome studies of male and female flowers on large-scale gene expression profiles in some plant species, generating considerable genomic data for the identification of functional genes involved in sex differentiation and flower development. however, the knowledge about unisexual male flower specification and development in andromonoecious plant is still limited, to some extent, due to its infrequent occurrence. up to now, only a few transcriptome studies have been performed for andromonocious plant species. to investigate the possible molecular mechanisms in the regulation of flower formation within andromonoecious taihangia, we performed de novo transcriptome sequencing of the male and hermaphroditic flowers at both early and late developmental stages by using the illumina platform. based on the differences in unigene abundance among sequencing libraries, we identified the differentially expressed genes between the male and hermaphroditic flowers at early and late developmental stages, respectively. furthermore, go and kegg enrichment analyses were conducted to elucidate the main functions of degs, and to identify candidate genes involved in the regulation of flower formation for further functional analysis. the present study would provide a genomic resource for gene discovery in the future and shed light on the underlying molecular mechanisms responsible for establishment of unisexual flowers and sex differentiation in andromonocious plant taihangia. methods plant material and rna extraction the plant materials were transplanted from yidoushui and zhuyufeng, henan, china in april and were grown in the greenhouse of henan polytechnic university. in early spring the early and late developmental stages of male and hermaphroditic flowers, including male floral bud, hermaphroditic floral bud, male flower, and hermaphroditic flower, were harvested. as shown in fig.  flowers were classified into six developmental stages: young flower buds, elongated buds, pre-anthesis, anthesis, and fully opened flowers. in this study, the flower samples at early developmental referred to flower buds including young and elongated buds. the flower samples at late developmental stage were comprised of young flowers at pre-anthesis, flowers at anthesis, and mature flowers with fully opened petals. the samples of the male and hermaphroditic flowers were collected from six different individuals at each of the developmental stage, respectively, and immediately frozen in liquid nitrogen and stored at −80 °c for rna extraction.fig. male and hermaphroditic flowers within andromonoecious taihangia at early and late developmental stages used in rna-seq. a hermaphroditic floral bud;. hermaphroditic flower; male floral bud; male flower. early developmental stage: young and elongated buds; late developmental stage: young flowers at pre-anthesis, flowers at anthesis, and fully opened flowers total rna was isolated from 100 mg homogenized plant material using an rneasy mini kit for plant tissues according to the manufacturer’s manuals. rna quality was preliminarily assessed by % agarose gel electrophoresis, and then quantified with a uv-visible spectrophotometer. rna samples for further analysis were selected according to the following criteria: a260/a ratio was - and a260/a ratio ranged from to. the equal amounts of rna from male flower buds at each of the early developmental stages were mixed as one em sample, while one lm sample was prepared by combinations of the identical amounts of rna from male flowers at each of the late developmental stages. similarly, eh and lh samples were obtained using the same method described above. in total, rna samples with three biological replicates were prepared for cdna synthesis and qpcr analysis. after that, specific rna pools of em, eh, lm, and lh were obtained by combining identical quantities of total rna from the corresponding biological triplicates. quality and quantity of the pooled rna were further analyzed using a nanodrop instrument and a chip rna series ii bioanalyzer, respectively. rna samples with 260/ ratios, 260/ and rin more than were used to cdna library construction and illumina sequencing. library construction and rna-seq construction of the four libraries and the rna-seq were performed by the biomarker biotechnology corporation. five micrograms of total rna from each sample were used to construct the sequencing libraries. rna sequencing libraries were generated using the truseq rna sample prep kit following the manufacturer’s recommendations. briefly, the poly -containing mrna molecules were purified from the total rna by using the poly-t oligo-attached magnetic beads. then, mrna was broken into short fragments, which were used as templates for cdna synthesis. double-stranded cdna was synthesized using superscript ii, buffer, dntps, rnaseh, dna polymerase i, and random hexamer primers. after that, the ‘a’ tail was added to the 3′ ends of the repaired cdna fragments and the illumina’s paired-end adapters were ligated to the cdna ends. the products from the ligation reaction were amplified by pcr. the pcr productions were purified by magnetic beads and dissolved in eb solution to generate the sequencing libraries. the quantity and quality of each sequencing library were tested by the agilent bioanalyzer. finally, the four libraries were sequenced separately on illumina hiseq™ platform. sequence assembly and gene annotation after sequencing, paired-end raw reads were firstly processed through in-house perl scripts. in this step, clean reads were acquired by removing the reads with adaptor contamination, low-quality sequences, and reads with more than 10% q <  bases. at the same time, q q and gc-content of the clean data were calculated. all the downstream analyses were based on clean data with high quality. de novo assembly of the transcriptome was performed by trinity software, with default parameters and no reference sequence. in brief, reads were assembled into the contigs by the inchworm program at the first step. the minimally overlapping contigs were clustered into connected components by the chrysalis program, and then the transcripts were constructed by the butterfly program. finally, the transcripts were further clustered based on nucleotide sequence identity, and the longest transcripts in the cluster units were regarded as unigenes to eliminate redundant sequences. for detection of transcriptional profiles in the male and hermaphroditic flowers at different stages, clean reads from each library were initially assembled separately. to obtain a uniform transcriptome reference across samples, all clean reads from four libraries were pooled together and de novo assembled to generate unigenes for assembly evaluation, gene annotation, and expression analysis. the raw sequencing data were deposited in the ncbi short read archive database under the accession number srp to identify putative functions of taihangia unigenes, we performed functional annotation by using a blastx search against a series of protein databases: the national center for biotechnology information non-redundant, swissprot, pfam, clusters of orthologous groups, gene ontology, and kyoto encyclopedia of genes and genomes databases. the blast2go software package was used to compare and determine the unigene go annotations. finally, wego software was used to obtain the go functional classifications for all annotated unigenes. transcription factors prediction in order to identify the tfs represented in taihangia transcriptome, all assembled unigenes were searched against the plant transcription factor database planttfdb by using blastx with a cut-off e-value of e−. to facilitate the transcription factor prediction, the identified tfs were further screened with similarity ≥  and coverage ≥ . analysis of differentially expressed genes to calculate the amount of gene expression for each unigene, rsem was used to quantify the number of reads mapped to the assembled transcriptome, and read count for each gene was obtained from the mapping results. we used the fpkm algorithm to normalize the gene expressional abundances in each library. by pairwise comparisons of the four libraries, the degs were identified using the degseq r package based on the read counts for each gene at different libraries. in addition, the false discovery rate control method was used to identify the threshold of the p-value in the significance tests. significance of differential expressed genes was defined at a false discovery rate <  and an absolute value of log2-ratio ≥  the classification of degs were performed with go and kegg analysis using the method described above. go enrichment analysis of the differentially expressed genes was implemented by using the cytoscape plugin bingo using a hypergeometric test with a corrected p-value < . kegg pathway enrichment analysis of the degs was performed by using kobas with the hyper-geometric distribution model. the enrichment p-values were adjusted by using the benjamin and hochberg method. validation of sequencing data by quantitative real-time pcr sixteen degs were randomly selected and the expression profiles were investigated by qrt-pcr to confirm the transcriptome data. first-strand cdna was synthesized from the same rna samples used for library construction by using the primescript rt reagent kit with gdna eraser according to the instructions. quantitative real-time pcr was performed by using the ssofast evagreen supermix rt-pcr kit and the miniopticon real-time detection system. the pcr mix was composed of 10 μl evagreen supermix,  μl of 1: diluted cdna,  μl of each primer, and  μl water in a final volume of 20 μl. the reactions were incubated under following cycling conditions: 2 min at 50 °c, 2 min at 95 °c, 40 cycles of 95 °c for 15 s, 58 °c for 15 s, and 72 °c for 30s, and finally 72 °c for 2 min. after the final cycle, a melting curve analysis was performed with a single cycle from 60 °c to 95 °c in 5 s intervals to verify the reaction specificity. the efficiency of the primer sets was calculated by performing real-time pcr on serial dilutions of first-strand cdnas templates. by using the genorm software, ef1α and ubq were identified as reference genes for stable expression across developmental stages. the relative expression of each gene was calculated using the standard e-ΔΔcq method, and normalized with the internal standard gene ef1α and ubq. expression quantification and data analysis were performed according to miqe guideline suggested by bustin et al.. three independent biological replicates and three technical repetitions were performed for each of the quantitative pcr experiments. the relative expression of each target gene was analyzed by one-way anova. statistical analyses of the data were performed at p <  significance level using spss software.  illumina sequencing and de novo assembly to obtain a comprehensive flower transcriptome in andromonoecious taihangia, four cdna libraries were constructed from male and hermaphroditic flowers covering early and late development stages. paired-end sequencing of the four constructed libraries were performed on an illumina hiseq platform, generating a total of, raw reads. after filtering for low quality region, adapters, and sequencing tags, we obtained, qualified illumina reads, including,,,,,, and, clean reads for the libraries of em, eh, lm, and lh, respectively. overall, the library produced g base pairs clean data with % q bases, indicating that the throughput sequencing was accurate enough to warrant the further analysis. the clean reads from each libraries were separately de novo assembled by trinity programs, generating,,, and unigenes, respectively, with mean lengths of to 905 bp. after that, we pooled all high quality reads from the four individual libraries to perform the de novo assembly for a uniform reference transcriptome. by using the trinity program, a total of unigenes longer than 200 bp were assembled with an n length of  bp, an average length of 840 bp, and a maximal transcript length of  bp. the length distribution of taihangia unigenes was shown in additional file 2: figure s table summary of illumina transcriptome sequencing and assembly for taihangia rupestris flower to evaluate the quality of de novo assembly without a reference genome, a common method is to calculate the percentage of reads that can be mapped back to unigene sets. based on this metric, we found that more than 84% of clean reads were mapped back to the taihangia unigenes for each library. in addition to in silico quality assessment, the standard pcr and sanger sequencing were conducted to evaluate the quality and accuracy of taihangia assembled unigenes. of designed primer pairs, primer pairs for unigenes could steadily generated expected pcr products. all pcr products were sequenced by sanger sequencer, and the obtained sequences were compared with the assembled unigenes. the alignment analysis showed that the sequences from sanger sequencing met well with the assembled unigenes with an average sequence identity of %, confirming the accuracy of assembled unigenes derived from transcriptome. together, our results supported the high quality of de novo assembly in current sequencing data. functional annotation and classification of the unigenes all the assembled unigenes were subjected to blastx comparisons, setting a cut-off e-value of 10− against the public databases including nr, pfam, swiss-prot, go, kegg and cog databases. as a result, a total of unigene sequences showed significant similarity to known proteins in publicly available databases. based on the annotation results of the nr database, the e-value distribution analysis showed that % of the matched sequences had strong homology with the e-value <  e− while % of the matched sequences ranged between 1e− and 1e− the similarity distribution analysis revealed that % of the unigenes had a similarity higher than 75%, while % of mapped sequences had a similarity ranging from 60% to 75%. for species distribution, the nr database queries showed that % of the t. rupestris annotated sequences were matched to the sequences of fragaria vesca, followed by prunus mume, prunus persica, malus domestica, and pyrus × bretschneideri. characteristics of the homology search of taihangia unigenes against the nr database were shown in additional file 5: figure s table summary for the annotation of unigenes of taihangia flower based on the blastx results against the nr database, we added go terms to obtain go functional annotations and categorizations of these assembled unigenes. of the unigenes with nr annotation, a total of unigenes with at least one go term were assigned into biological processes, cellular components, and molecular functions categories. the three main categories were further classified into functional groups including biological processes, cellular components and molecular functions. in the category of molecular function, “catalytic activity” was the most highly represented group, followed by group of “binding”. under the biological process category, “metabolic process” and “cellular process” were the most highly represented groups. for the cellular component, the major groups were “cell” and “cell part”. the go analysis indicated that a great number of identified unigenes were associated with various biological processes and molecular functions in taihangia floral tissues. the annotated sequences were further subjected to a search against the cog database for functional prediction and classification. as a result, unigenes were grouped into cog classifications, among which the “general function prediction only” represented the largest group, followed by “replication, recombination and repair”, “transcription”, “signal transduction mechanisms”, and “posttranslational modification, protein turnover, chaperones”. the category “nuclear structure” with unigenes was the smallest group. to obtain a better understanding of biological functions of the unigenes, the annotated sequences were searched against the kegg database. among the annotated unigenes, had significant matches and were assigned to kegg pathways. the top five pathways were “ribosome”, “protein processing in endoplasmic reticulum”, “carbon metabolism”, “biosynthesis of amino acids” and “spliceosome”. functional go, cog, and kegg pathway annotations and classification of the unigenes were shown as additional file in summary, all these annotation and classification analyses would provide valuable information for gene discovery and functional genomic studies in the future. identification of transcription factors in the present study, we identified tfs, representing % of taihangia unigenes and falling into tf families classified by plant transcription factor database planttfdb. among these detected tf gene families, bhlh was the most abundant tf family, followed by myb_related, erf, nac, wrky, and c2h. notably, a total of mads-box genes, including mikc-type and 17 m-type tfs, were found in taihangia based on transcriptome data.fig. identification of transcription factors and mads-box genes in taihangia. a classification of tf families. b phylogenetic tree of mikc mads-box proteins of taihangia, arabidopsis thaliana, strawberry, and peach. red dots, taihangia genes; green dots, strawberry genes; blue dots, a. thaliana genes; purple dots, peach genes to examine the phylogenetic relationships among the taihangia mads-box genes and to classify them into the established subfamilies, a neighbor–joining phylogenetic tree was constructed based on a multiple sequence alignment of the predicted full-length mads-box protein sequences of f. vesca, p. persica, and arabidopsis. the phylogenetic analysis revealed that mads-box tfs, including eight mads-box genes reported in previous studies, were classified into classified into mikcc and one mikc* subfamilies. furthermore, several mads-box tfs, such as apetala, pistillata, agamous, shaterproof, and sepallata, described in abcde model were also included. trumads trumads and trumads belonged to the ap1/ful subfamily in the a class, which included the ap1and ful clades. trumads was classified into ap clade, while trumads and trumads were sorted into ful clade. in the functional b-class group, two genes belonged to pi clade, while the other one gene fell within tm three mads-box tfs, trumads trumads and trumads were assigned to the ag family, which were functionally classified as a c/d class group. trumads belonged to the c-class group ag while trumads and trumads were most orthologous to the d-class functional group shp and stk, respectively. in the e functional group, two mads-box tfs, trumads and trumads were assigned to sep subfamily. in general, the phylogenetic analysis showed that taihangia contains mads proteins in each putative functional group. transcriptome profiles in male and hermaphroditic flowers differences of gene expression profiles between male and hermaphroditic genders at different developmental stages were examined by pairwise comparisons of the four libraries. a total of unigenes were shared among all libraries, while approximate % of unigenes was detectable in individual libraries. for hermaphroditic genders, there were and unigenes exclusively expressed in eh and lh, respectively. meanwhile, we detected unigenes specific for em and unigenes only for lm. a total of unigenes were gender-specific, with unigenes detected only in the male flowers and unigenes in hermaphroditic flowers, respectively. in addition, % of assembled unigenes was only detected in floral buds, whereas % of assembled unigenes was observed in the mature flowers, but not in the floral bud samples. overall, our results suggested that the unigenes specific expressed in the male or hermaphroditic flower might be potential candidates to be associated with the regulation of flower formation in andromonoecious taihangia. differential expressed genes analysis based on the four comparisons of em vs eh, lm vs lh, em vs lm, and eh vs lh, we identified a total of differential expression genes. among them, degs between gender types and degs between developmental stages were detected, indicating that majority of degs occurred among different developmental stages. in order to identify related genes that involved in the regulation of flower formation, degs between male and hermaphroditic genders were further analyzed. at the early developmental stage, we identified up-regulated genes and down-regulated genes expressed in eh. on the other hand, up-regulated genes and down-regulated genes were detected in lh at the late developmental stage. interestingly, degs, including hermaphroditic-biased and male-biased unigenes, were differentially expressed in both early and late developmental stages, indicating that the gene expression pattern of the degs was maintained from the floral bud to the mature flower across developmental stages.fig. identification of differentially expressed genes in taihangia. venn diagram of the number of degs based on four comparisons of em vs eh, lm vs lh, em vs lm, and eh vs lh. total deg numbers in each of comparison were given outside the circles, and arrows indicate up- or down-regulation among the tfs identified in this study, tf genes were differentially expressed between male and hermaphroditic flowers. these tf genes fall into tf families classified by the planttfdb, with the bhlh family representing the largest amount of differentially expressed tfs, followed by erf, myb related, wrky, and nac. tf genes showed different expression patterns in male and hermaphroditic flowers during flower development, suggesting possibly different roles of tfs in the early and late developmental stages. at the early developmental stage, putative tf genes belonging to tf families were differentially expressed in comparison em vs eh. the tf genes belonging to tf families, such as bhlh, myb related, zf-hd, nac, and erf, were significantly upregulated, while the remaining tf genes belonging to tf families were downregulated in em. in the late developmental stage, the identified putative tf genes were classified into tf families based on the comparison lm vs lh. of them, tf genes were detected with significantly upregulated expression in lh. these tf genes might be important regulators contributing to the establishment and development of male and hermaphroditic flower within the andromonoecious system. the expressed profiles of tfs derived from degs were shown as fig. 4b.fig. the differential expressed transcription factors between the male and hermaphroditic flowers. a the number of degs belonging to different transcription factor families detected in two comparisons of em vs eh and lm vs lh;. expression pattern of differential expressed tfs between the male and hermaphroditic flowers in summary, our results showed that the number of degs between male and hermaphroditic flowers remarkably increased alongside the flower development, and quite different deg profiles were observed between comparisons em vs eh and lm vs lh. given that the morphogenesis of the taihangia uni- and bisexual flowers occurred in the early flower developmental stages, degs between em and eh were expected to be involved in sex differentiation, while degs between lm and lh were probably associated with floral organ development and reproduction. go and pathway analyses of degs go functional classifications of the degs between male and hermaphroditic samples were conducted by using the wego website. of the identified degs between gender types, degs were classified into go subterms. in the cellular components category, the cell, cell part, membrane, and organelle, were the most abundant subterms. the catalytic activity and the binding subterms were dominant in the molecular function category. for the biological process category, a high percentage of degs was assigned to the cellular process, metabolic process, and single-organism process subterms.fig. go classifications of degs between the male and hermaphroditic flowers in different developmental stages. degs are grouped to the secondary classification of hierarchical go terms. a go classification of degs in two comparisons of em vs eh and lm vs lh at early and late developmental stages;. go classification of upregulated and downregulated degs in comparison em vs eh at early developmental stage;. go classification of upregulated and downregulated degs in comparison lm vs lh at late developmental stage fig. top enriched kegg pathways of degs in em vs eh and lm vs lh between male and hermaphroditic flowers. number at the right of each column indicates the number of degs included in corresponding pathway. the pathways were divided into five categories: “metabolism”, “cellular processes”, “environmental information processes”, “genetic information processes”, and “organismal systems” to better understand the regulation mechanisms of flower formation across developmental processes in taihangia, we put emphasis on go functional classifications analysis of upregulated and downregulated degs in comparisons em vs eh and lm vs lh, respectively. furthermore, go enrichment analysis was performed by using the hypergeometric test in bingo program with a corrected p value ≤ . at the early developmental stage, degs were classified into go subterms. one hundred and seventy three upregulated genes and downregulated genes in em were classified into and go subterms, respectively. overrepresentations of go terms in em were found in “pollen tube development”, “regulation of pollen tube growth”, “hydrolase activity, acting on glycosyl bonds”, indicating that these genes involved in stamen development. for up-regulated degs in eh, most of the overrepresented go terms were related to cell cycle, developmental process, and metabolic process, such as “regulation of macromolecule metabolic process”, “regulation of cell cycle”, “anatomical structure development”, “integral component of membrane”, and “regulation of transcription, dna-templated”. in the late developmental stage, go subterms were found for degs in comparison lm vs lh. one hundred and twenty six upregulated genes and downregulated genes were classified into and go subterms, respectively. the overrepresentations of go terms in lm were observed in “oxidoreductase activity” while overrepresentations of go terms in lh mainly included “cell cycle process”, “dna methylation”, “flower development”, and “dna binding”. in particular, a large number of go terms, such as “regulation of cell cycle”, “regulation of transcription, dna-templated”, and “regulation of developmental process” overrepresented in early developmental stages were also enriched at late developmental stages. the go enrichment analysis results are presented in additional file to further investigate the biochemical pathways of these degs, we mapped the degs to terms in the kegg database. among them, genes had a ko id and could be categorized into pathways, of which “phenylpropanoid biosynthesis”, “plant-pathogen interaction”, “plant hormone signal transduction”, “pentose and glucuronate interconversions”, “fatty acid metabolism”, and “unsaturated fatty acid metabolism” had the most deg numbers. in addition, we conducted an enrichment analysis with the kegg pathway, which assigned the degs to pathways such as biosynthesis of unsaturated fatty acids, phenylpropanoid biosynthesis, and pentose and glucuronate interconversions. in particular, two pathways “unsaturated fatty acid metabolism” and “dna replication” were highly represented, respectively, which were consistent with the go enrichment analysis.fig. expression pattern of functional degs involved in sex differentiation and flower organ development. a flower development; epigenetic processes; cell proliferation. the bar represents the scale of relative expression levels of degs. red means high expression and blue means low expression go and kegg pathway analysis revealed that degs between the two genders were associated with a wide range of functions, and participated in many biological processes, such as cell proliferation, chromatin modification, and flower development, suggesting that the formation of uni- and bisexual flowers might be controlled by many genes involved in numerous biological processes. moreover, the identified go subterms and kegg pathways could provide further information on the molecular mechanism of sex differentiation in andromonoecious taihangia. degs involved in regulation of taihangia male and hermaphroditic flower formation flower development based on deg analysis and gene expression profiles of male and hermaphroditic flowers, we identified a serial of degs involved in flower development. among them, several genes related to floral organ development, such as mads-box, wuschel, yabby, and lrr-protein kinase bam were up-regulated in hermaphroditic flowers. particularly, four differentially expressed mads-box genes were found in the ag, b-sister, and mikc* subfamilies. among them, both trumads and trumads were up-regulated in hermaphroditic flowers at late developmental stages. trumads, belonging to d-class functional genes, showed dramatically higher expression in hermaphroditic than in male flowers throughout floral development. our results indicated that these mads-box genes preferably expressed in lh were likely associated with gynoecium development in hermaphroditic flowers. conversely, trumads belonging to mikc* subfamilies, was upregulated in em at early developmental stages, suggesting that the mads-box tf probably involved in stamen development. epigenetic process epigenetic mechanisms, including dna methylation, histone modifications, chromatin conformation, and rna interference, play essential roles in modulating gene expression without altering the dna sequence. in this study, a number of degs involved in dna methylation and histone modification were identified in taihangia flower at late developmental stage. several genes, including three dna methyltransferase cmt met, and drm involved in dna methylation process showed differentially expression between male and hermaphroditic flowers. on the other hand, histone modifications related genes, such as histone-lysine n-methyltransferase atxr ashr and suvr were differentially expressed in comparison lm vs lh at late developmental stages. no significant differences were observed in the expression levels of epigenetic related genes between male and hermaphroditic flowers at the early developmental stages, while all of them were significantly upregulated in lh at late developmental stages. notably, one homologous of ago gene involved in rna silencing and rna-directed dna methylation pathway exhibited upregulated expression patterns in lh at late developmental stages. our finding indicated that epigenetic processes may extensively occur in hermaphroditic flowers, especially at late developmental stages. cell cycle process a large number of degs involved in cell cycle process, such as dna replication, cell division, and regulation of cell cycle, were identified based on transcriptome data, and these degs showed significantly increased expression in hermaphroditic flowers at late developmental stage. a number of degs, including mitotic spindle checkpoint protein mad2-like, atp-dependent dna helicase q-like, dna topoisomerase 3-alpha-like, and mitotic checkpoint serine/threonine-protein kinase bub beta-like were assigned into dna replication functional group. meanwhile, it was found that majority of degs, belonging to cell division functional group, were down-regulated in lm at late developmental stage, including 65-kda microtubule-associated protein, g2/mitotic-specific cyclin-1-like, condensin complex subunit 2-like, dna replication licensing factor mcm2-like, and dna replication licensing factor mcm similarly, several degs involved in cell cycle group also showed down regulated expression in male unisexual flowers, such as microtubule-associated protein rp/eb family member 1c, dna cross-link repair protein snm1-like, protein fizzy-related 3-like, and dna cross-link repair protein snm1-like. in particular, many degs, including g2/mitotic-specific cyclin-1-like, protein fizzy-related 3-like, mitotic spindle checkpoint protein mad2-like, mitotic checkpoint serine/threonine-protein kinase bub beta-like, cyclin-a1-1-like and cyclin-d3-1-like, showed significantly upregulated expression in both eh and lh across the developmental stages. overall, our results indicated that the cell cycle process was probably severely disturbed in male unisexual flowers. biosynthesis of unsaturated fatty acids unsaturated fatty acids are considered to be essential membrane components, and play key roles in many cellular events. the synthesis of unsaturated fatty acid is a complicated process, requiring for involvement of a series of enzymes, such as palmitoyl-monogalactosyldiacylglycerol Δ desaturase and acyl- desaturase. in the present study, several degs encoding several key enzymes invovled in the biosynthesis of unsaturated fatty acids pathway, including six fad three very-long-chain 3-oxoacyl-coa reductase 1-like genes, and one acp genes, were found to be significantly downregulated in male flowers, indicating that biosynthesis of unsaturated fatty acids was likely restricted in male flowers. validation of unigenes and gene expression profiling using rt-qpcr to confirm the gene expression patterns identified by the rna-seq data, the transcript levels of random selected degs were examined by quantitative pcr. we performed qrt-pcr using these genes with biological replicates from each group. the detailed information of the selected degs, ids and primer pairs used in this study were shown in additional file 12: table s it was found that all of selected degs were successfully amplified with single bands with the expected sizes, and the patterns of gene expression detected by rt-qpcr were consistent with those from rna-seq data. therefore, rt-qpcr experiments confirmed that the degs obtained from the assembled transcriptome were accurate and gene expression profiles from rna-seq data were reliable.fig. validation of randomly selected differentially expressed genes derived from rna-seq using real-time quantitative rt-pcr. quantitative gene expression data represented as the mean ± sd are shown on the left. the ef1α and ubq were used as reference genes. statistically significant differences of relative expression derived from one-way anova at the level of p <  were marked as different lowercase letters. the normalized fpkm, transformed by log, were indicated on the right discussion taihangia transcriptome profiles as an andromonoecious plant species, taihangia provides an excellent model system within a uniform genetic background for the study of sex differentiation and establishment of unisexual flowers. in this study, we conducted de novo transcriptome sequencing to compare transcriptional profiles of the male and hermaphroditic flowers by using the illumina rna-seq method. to the best of our knowledge, this is the first report of transcriptome data in this rare cliff plant species with the andromonoecious sexual system. a total of approximately, qualified illumina reads were obtained from four floral sequencing libraries and then assembled into unigenes. of the assembled unigenes, unigenes were annotated in the public database, while the function of other assembled unigenes could not be identified. given that the genomic and transcriptomic information is limited for taihangia, the lineage-specific genes associated with the unique floral phenotype may not be contained in function annotation, implying that potential novel genes were discovered from our transcriptome data. by comparative analysis of male and hermaphroditic transcriptome data, we identified a large number of candidate degs involved in various biological processes, which is likely associated with the regulation of male and hermaphroditic flower formation. the transcriptome data provide an important resource for the understanding of specific biological processes involved in sex differentiation and facilitating gene discovery in the future studies in andromonoecious taihangia. floral organ identity-related mads-box genes may be crucial for flower organ formation the mads-box transcription factors, initially identified as floral homeotic genes, are the major components in the well-known ‘abcde’ model that describes their roles in the determination of floral organ identity. based on taihangia transcriptome data, we identified mads-box tfs, including several mads proteins in each putative functional group described in abcde model. according to the “abcde” model, similar levels of expression of class a-, b-, and e-class mads-box genes between male and hermaphroditic genders were expected, while class c- and d- mads genes should be upregulated expressed in the hermaphroditic flowers. in the present study, most of mads-box tfs including a-, b- and e-class of genes did not show differential expression between male and hermaphroditic flowers. as expected, two mads-box tfs, truag and trustk, showed upregulated expression in hermaphroditic flowers. however, it was noteworthy that upregulated expression of truag was not found in eh but in lh, suggesting that its expression pattern appears to follow the arrest of pistil development in male flowers. a number of studies have found that the mads box genes did not play a part in sex determination in some dioecious and monoecious species that achieve unisexuality by organ abortion. in taihangia, pistils are initiated and then arrest in unisexual male flower, the molecular mechanism of sex determination must be acting downstream of organ identity. thus, the transcriptional profile of several mads-box tfs more expressed in hermaphroditic flowers might be the consequence rather than the cause of sex differentiation in andromonoecious taihangia, and the mads-box tfs were related to pistil development in hermaphroditic flowers. epigenetic-related genes are likely associated with reproduction in hermaphroditic flower epigenetic processes such as dna methylation, histone modification, and small rna, play important roles in the regulation of gene expression during flower development. as key enzymes involved in dna methylation processes, dna methyltransferases and histone h3k methyltransferase are crucial for establishment and maintenance of dna methylation. in our study, three dna methyltransferases and three histone h3k methyltransferase showed preferable expression in lh, indicating their possible roles in taihangia reproduction at the late developmental stage. in addition, several genes, such as argonaute, suppressor of gene silencing, rna-dependent rna polymerase, dicer homolog 3a-like, and rpd, possibly involved in rna-induced silencing and rna-directed dna methylation pathways were identified based on taihangia transcriptome data. particularly, ago showed approximately three folds higher expression in lh than in lm at late developmental stage. as key regulators, ago proteins and dcl along with rdr, are crucial for the epigenetic silencing. based on gene expression analysis, three genes, dcl, rdr, and rpd, showed more than two times higher expressions in lh than in lm, indicating that extensive small rna related processes possibly occurred in taihangia hermaphroditic flower at developmental stage. currently, increasing evidences have supported that epigenetic processes played important roles in regulation of plant reproductive development. for example, galla et al. reported that small rna generation related genes, including dcl ago, sgs and rdr, were involved in plant reproduction, such as cell identity and fate in the ovule, in hypericum perforatum. the upregulated expressions of these epigenetic-related genes in taihangia hermaphroditic flowers, especially at the late developmental stages, may indicate that the epigenetic regulation played important roles in plant reproduction, such as female gametophyte formation. the roles of epigenetic processes in t. rupestris are still unclear, and future studies are needed to investigate the roles of associated genes involved in epigenetic processes. disturbance of cell cycle may contribute to pistil abortion in male flower during flower development, arrest of pistil or stamen development is a normal phenomenon in monoecious and dioecious plant. several biological processes, such as programmed cell death and cessation of the cell division, are critical in establishment of unisexual flower. based on taihangia rna-seq data, a great number of degs involved in cell cycle, dna replication, and cell division were significantly downregulated in male flowers over the developmental stages, indicating that these biological processes played important roles in arrest of pistil development in taihangia male flowers. for many plant species, pcd is essential to many aspects of plant morphogenesis, and it is a common process for the developmental arrest of reproductive organs in unisexual flowers. however, we did not found that pcd involved in pistil abortion in taihangia male flowers based on transcriptome data. it has been reported that developmental arrest of sterile sex organs and the subsequent unisexuality of flowers could resulted from a cessation of cell division. similarly, down-regulated expressions of related genes involved in cell division probably led to restriction of cell proliferation, which may be closely associated with the arrest of pistil development in taihangia unisexual male flower. as primary regulatory genes, cyclins play critical roles in controlling cell cycle progression, including dna replication, the g2/m transition, and mitosis. in this study, several cyclin genes related to regulation of cell cycles were found to be uniformly down-regulated in unisexual male flowers. particularly, cyclin-a1-1-like and cyclin-d3-1-like significantly decreased expression level in em and lm across developmental stages, indicating that there was a potential mechanistic link between the regulation of cell cycle and the establishment of unisexual flowers in taihangia. during flower development, cyclin-a1and cyclin-d were considered as the core regulatory genes involved in regulation of flower primordia information and meristem activity. for example, cyclin-d could define distinct developmental zones, acting as a developmental boundary between lateral organs, in the shoot apical meristems by controlling the balance between cell division and cell differentiation. in taihangia unisexual male flowers, downregulated expression of cyca1and cycd probably resulted in the disturbance of cell cycle in gynoecium primordia, and subsequently determined arrest of pistil development. these results suggested that the core regulatory genes involved in cell cycle progression probably contribute to pistil abortion in taihangia unisexual flowers by regulation of the cell cycle. cold-responsive genes are likely to affect the formation of unisexual flowers it is well known that the establishment of unisexual flower is influenced by environmental factors such as temperature, light intensity, photoperiod, and nutritional conditions via specialized gene expression. for andromonoecious taihangia, temperature during the initiation of floral primordia is crucial for determining morphogenesis of gender type, and generation of unisexual male flowers is attributed to pistil abortion under low temperature. to deal with adverse cold stress, plants have developed versatile defense mechanisms to acquire cold tolerance by adjustments in molecular, and physiological, and developmental processes. for example, increase of membranes’ unsaturated fatty acid content could promote the stabilization of plasma and chloroplast membranes, and thus improve resistance to cold stress. as key enzymes involved in biosynthesis of unsaturated fatty acid, fad acpr, and acp were dramatically down-regulated in em, indicating that restriction of polyunsaturated fatty acids biosynthesis probably occurred in taihangia male flowers. under low temperature condition, a number of cold-regulated genes such as lea family genes could be induced in higher plants. for instance, lea proteins are thought to play a role in membrane stabilization by interacting with membrane surfaces or with water molecules under cold stresses. based on transcriptome data, a number of cold-resistance genes, including lea dc lea b, ascorbate oxidase, and peroxidase, were significantly down-regulated in male flowers. these results agreed well with a recent study on comparative transcriptome analysis between male and hermaphroditic flowers in japanese apricot by shi et al., arguing that reduction of cold-resistance gene expression in male flowers is likely associated with pistil abortion under low temperature. taken together, these results suggested that down-regulated expression of cold-responsive genes at the early stage is likely to affect the formation of unisexual male flowers in andromonoecious taihangia. candidate tfs as important regulators in establishment of male and hermaphroditic flowers in plant, tf genes play an important role in many aspects of plant growth and development. in the present study, we identified a large number of tf genes showing significantly differential expression between male and hermaphroditic flowers. these tfs are expected to be involved in the regulation of formation of the male and hermaphroditic flowers in taihangia. of these identified tfs, most stress-responsive tfs, such as nac, wrky, bhlh, gras, and erf showed increased expression in male flowers, implying that male flowers probably suffered more severe injuries than hermaphroditic flowers. majority of upregulated tfs in bisexual flowers, including wox, yab, ap and c2h were associated with flower organ development. for example, wuschel, a number of wox tf genes, is known to maintain the stem cell activity in the floral meristems. after floral induction, wus can induce ag expression and promote the up-regulation of ag in developing flowers. in our study, higher expression of one homolog of wus was not observed in hermaphroditic flowers at early developmental stages. similarly, truag showed no difference in expression between em and eh, suggesting that a possible different regulatory mechanism occurred in taihangia flower development. the yabby gene family is involved in adaxial–abaxial polarity in lateral organs, and participates in a diverse range of processes such as floral patterning, organ growth, and maintenance of meristem organization. one homologous of yabby genes identified in this study was upregulated in both eh and lh, indicating that yabby tfs were closely associated with the formation of bisexual flowers. interestingly, five genes encoding growth-regulating factor were detected and showed uniformly increased expression in hermaphroditic flowers. grfs are plant-specific transcription factors, which play important roles in the control of cell proliferation via cell cycle regulation and central developmental processes including flower formation. recently, liang, et al. reported that dramatically reduced expression levels of grf caused pistil abnormalities and suggested an important role of grf in controlling carpel number and pistil development. our result was consistent with the finding in arabidopsis, and confirmed that disturbance of cell proliferation in male flowers contributed to the pistil abortion.  de novo transcriptome sequencing of andromonoecious taihangia was performed by using the illumina paired-end sequencing technology in the present study. in total, unigenes were obtained from the male and hermaphroditic flowers covering early and late developmental stages, and unique sequences were annotated in the public database. focusing on the regulatory mechanism of flower formation in andromonoecious taihangia, we identified degs between male and hermaphroditic flowers. go and pathway analysis showed that these degs involved in a wide range of biological processes, such as flower development and reproduction, cell cycle, epigenetic processes, and transcriptional regulation, and provided numerous candidate genes possible related to sex differentiation for future functional analysis. the transcriptome data of taihangia could be helpful to improve understanding of the molecular mechanisms in sex differentiation and unisexual flower establishment, and would provide theoretical basis for the regulation of flower formation in this andromonoecious plant species. additional files additional file 1: figure s male and hermaphroditic flowers produced within the same individual of taihangia on the cliff face. additional file 2: figure s length distribution of the unigene sequences. additional file 3: table s summary for the alignment of reads to unigene libraries. additional file 4: table s validation of assembled unigenes by using sanger sequencing. additional file figure s characteristics of the homology search of taihangia unigenes against the nr database. additional file 6: functional go, cog, and kegg pathway annotation and classification of the taihangia unigenes. additional file 7: identification of transcription factors against plntfdb database. additional file 8: table s taihangia mads-box gene names and attributes. additional file 9: figure s venn diagram of the number of specific and shared unigenes in em, eh, lm, and lh samples. additional file 10: over-representative go terms of degs in em vs eh, and lm vs lh, respectively. additional file 11: pathway classification and enrichment analysis of the degs between the male and hermaphroditic flowers. additional file 12: table s the selected degs, ids, and primer pairs for rt-qpcr. abbreviations cogclusters of orthologous group degdifferentially expressed gene fpkmfragments per kilobase of exon per million fragments mapped gogene ontology keggkyoto encyclopedia of genes and genomes qrt-pcrquantitative reverse transcription polymerase chain reaction rddmrna direct dna methylation rna-seqribonucleic acid sequencing taihangia rupestris taihangia tftranscription factor 
taihangia rupestris, an andromonoecious plant species, bears both male and hermaphroditic flowers within the same individual. however, the establishment and development of male and hermaphroditic flowers in andromonoecious taihangia remain poorly understood, due to the limited genetic and sequence information. to investigate the potential molecular mechanism in the regulation of taihangia flower formation, we used de novo rna sequencing to compare the transcriptome profiles of male and hermaphroditic flowers at early and late developmental stages. by pairwise comparison of unigene abundance between the samples, we identified differential expressed genes, including transcription factors between the male and hermaphroditic flowers. at the early developmental stage, we found up-regulated genes and down-regulated genes expressed in hermaphroditic floral buds, while up-regulated genes and down-regulated genes were detected in hermaphroditic flowers at the late developmental stage. comparative transcriptome analysis revealed the significant differences in gene expression profiles between male and hermaphroditic flowers at early and late developmental stages. four cdna libraries, including male floral bud, hermaphroditic floral bud, male flower, and hermaphroditic flower, were constructed and sequenced by using the illumina rna-seq method. male flowerhermaphroditic flowerflower formationsex differentiationtaihangiatranscriptomehttp://dx.doi.org/ /501100001809national natural science foundation of china3110018231170354gong hongmei zhang yanwen scientific innovation fund from henan polytechnic universityb2008-25t2013-2li weiguo issue-copyright-statement© the author 2017 go and kegg enrichment analyses showed that a large number of degs were associated with a wide range of functions, including cell cycle, epigenetic processes, flower development, and biosynthesis of unsaturated fatty acid pathway. the transcriptome data of taihangia would be helpful to improve the understanding of the underlying molecular mechanisms in regulation of flower formation and unisexual flower establishment in andromonoecious plants. in addition,,, and unigenes were assigned into gene ontology functional groups, clusters of orthologous group categories, and kyoto encyclopedia of genes and genomes pathways, respectively. in this study, transcriptome data of this rare andromonoecious taihangia were reported for the first time. finally, real-time quantitative pcr was conducted to validate the degs identified in the present study.
 taihangia rupestris, an andromonoecious plant species, bears both male and hermaphroditic flowers within the same individual. however, the establishment and development of male and hermaphroditic flowers in andromonoecious taihangia remain poorly understood, due to the limited genetic and sequence information. to investigate the potential molecular mechanism in the regulation of taihangia flower formation, we used de novo rna sequencing to compare the transcriptome profiles of male and hermaphroditic flowers at early and late developmental stages.  four cdna libraries, including male floral bud, hermaphroditic floral bud, male flower, and hermaphroditic flower, were constructed and sequenced by using the illumina rna-seq method. totally,, qualified illumina reads were obtained and then assembled into unigenes, of which unigenes were annotated in the ncbi non-redundant protein database. in addition,,, and unigenes were assigned into gene ontology functional groups, clusters of orthologous group categories, and kyoto encyclopedia of genes and genomes pathways, respectively. by pairwise comparison of unigene abundance between the samples, we identified differential expressed genes, including transcription factors between the male and hermaphroditic flowers. at the early developmental stage, we found up-regulated genes and down-regulated genes expressed in hermaphroditic floral buds, while up-regulated genes and down-regulated genes were detected in hermaphroditic flowers at the late developmental stage. go and kegg enrichment analyses showed that a large number of degs were associated with a wide range of functions, including cell cycle, epigenetic processes, flower development, and biosynthesis of unsaturated fatty acid pathway. finally, real-time quantitative pcr was conducted to validate the degs identified in the present study.  in this study, transcriptome data of this rare andromonoecious taihangia were reported for the first time. comparative transcriptome analysis revealed the significant differences in gene expression profiles between male and hermaphroditic flowers at early and late developmental stages. the transcriptome data of taihangia would be helpful to improve the understanding of the underlying molecular mechanisms in regulation of flower formation and unisexual flower establishment in andromonoecious plants. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords male flowerhermaphroditic flowerflower formationsex differentiationtaihangiatranscriptomehttp://dx.doi.org/ /501100001809national natural science foundation of china3110018231170354gong hongmei zhang yanwen scientific innovation fund from henan polytechnic universityb2008-25t2013-2li weiguo issue-copyright-statement© the author 2017 