Elucidation of molecular mechanisms of flower form development in tree peony ( Paeonia suffriticosa ) through comparative transcriptome analysis of floral parts CURRENT STATUS: POSTED

Background: Tree peony (Paeonia suffruticasa) is an economically, medicinally ornamentally important woody flowering woody plants in East Asia and is a common also ornamental shrub in Europe and North America. It is well known and prized for their beautiful flowers in many different forms. Samen petalody has been shown to be the most effective way to modify flower forms. However, there is limited information on the molecular mechanisms of stamen petalody and flower form formation in tree peony. Results: In this study, RNA sequencing was used to assemble and annotate the unigenes in the tree peony to identify the critical genes related to flower parts formation and verify the key genes in different flower forms of tree peony cultivar. A total of 76,007 high quality unigenes were assembled and 30,505 were successfully annotated. A total of 1,833 TFs were identified in our study, among them 16 MADS-box genes were found and characterized. Six key genes were selected to verity their functions in stamen petalody. AG and SEP showed high expression level in carpals and sepals separately both in stamen petalody group and non-stamen petalody groups. PI and AP3 showed high expression levels in inter-petals in stamen petalody groups than in staments in non-stamen petalody. Conclusion: Sixteen MADS-box genes were identified for the first time in tree peony through RNA-seq method. We identified six key genes based on their differential expression levels in different flower parts. These six key genes represented all categories in the ABCDE model to verify the functions in stamen petalody. PI and AP3 were verified to likely play important roles in regulating stamen petalody in tree peony. Our study has helped establish the flower development model in tree peony, identified key molecular mechanisms in the development of different flower forms, and provided valuable information in improving genetic diversity of tree peony and many other woody plants. for 116 pathways at Calyx vs Stamen, for 111 pathways at Calyx vs Pistil and for 115 pathways at Stamen vs Pistil. Among the six paired comparisons, a total of 5,616 unigenes were mapped in KEGG pathways including 324, 505, 361, 821, 571 and 835 DEGs respectively in Petal Vs Calyx, Petal Vs Stamen, Petal Vs Pistil, Calyx Vs Stamen, Calyx Vs Pistil and Stamen Vs Pistil. The results of each comparison are c60425, c91686 and 75076 were nearly no expressed in all four floral parts. For stamen, c85161 (D class), c72042, c70349, and c80146 were of a high expression, meanwhile c70349 (B class), c80146 (B class), c87227, c58746 (A class) and c72216 (E class) were overexpressed in petal. Notably, c70164 and c69592 were specific expressed in sepal and c57631 (D class) in carpal. powerful and efficient approach for gene expression analysis at genome level, especially in the plant species without genomic sequence information. In our study, a total of 11,801 unigenes were detected in all four floral parts. Among them, 14.3% of unigenes were found to be specifically expressed in certain individual floral parts. Among all of the flower parts examined, the unigenes that were specifically expressed in carpels were found in the greatest abundance, followed by stamen, calyx, and petal, respectively. These results suggest that the unigenes specifically expressed in

Age pathway, Ambient Temperature Pathway Autonomous Pathway, Gibberellin Pathway, Photoperiod Pathway, and Vernalization pathway have been reported to be the key pathways in flower initiation in Arabidopsis thaliana [1] and numerous ornamental crops [2][3][4][5][6][7]. There are many regulatory genes involved in these six pathways of flower development. The ABC model of flower development based on exhaustive studies of Arabidopsis thaliana and Antirrhinum majus is widely regarded as the basic framework of the flower development of many higher plants [8]. The original ABC model, which had been greatly enhanced during the 1990s, was adopted as the most comprehensive model for flower development in 2000 and is now widely known as the ABCDE model [9][10][11][12][13]. Nearly all genes identified in the ABCED model (A class, B class, C class, D class and E class) are members of the MADS-box gene family [13][14][15][16][17][18]. Although there are anatomical differences between herbaceous annuals and perennial woody plants, the key regulatory genes involved in flower initiation and development have been shown to be highly conservative.
Tree peony (Paeonia suffruticasa), a deciduous woody shrub, has been cultivated for more than 2,000 years in China, is a well-known for its medicinal, nutraceutical and ornamental values in East Asia and highly treasured by many gardeners all over the world [19,20]. Tree peony is well known for its large and shower flowers with many different flower forms. The flower forms of tree peony are very diverse and have been divided into several classes including Simple, Lotus, Chrysanthemum, Rose, Anemone, Thousand Petal Crown, Hydrangea Globular and Hundred Proliferate [20]. These flower forms have been illustrated in several published papers [21,22] and are shown in Fig. 1 for clearer explanation.
The flower forms of tree peony vary a great deal in the arrangement, number, and shape of flower parts including petals, sepals, stamens and carpels [20]. Two main reasons for these highly distinctive flower forms are increases in the number of naturally occurring flower petals and stamen petalody.
The number of naturally flower occurring flower petals of tree peony ranges from 30 to 50, while stamen petalody is greater than 200 [20]. The later leads to more cultivars of different flower forms than the number of naturally flower petals of tree peony. In order to develop highly prized cultivars with even more attractive flower forms in tree peony, it is very critical to gain a deep understand of the molecular mechanism of flower development. A limited number of genes involved in flower developed have been cloned and real-time expressed. However, most regulator genes involved in flower development in tree peony have not been identified. In this study, we aim to identify the key genes involved in the development of flower parts using by RNA-Seq and verify their specific roles in controlling the presence or the absence of stamen petalody. The results from this study will help elucidate the molecular mechanisms of flower development of tree peony and could lead to more efficient development of new cultivars for the $14 million tree peony industry in China!

Results
Transcriptome sequencing and assembly A total of 118.98 Gb clean data from four floral parts with three biological replications each was acquired following a systematic quality check and data filtering. The base number of each sample was greater than 8.5 Gb (Table 1). The percentages of the GC bases of all 12 samples were greater than 44%, and the percentages of the reads with an average quality value greater than 30 were over 89%.
These parameters confirmed that sequencing data were reliable and appropriate for further analysis.
A total of 76,007 unigenes were obtained with median contig (N50) length of 1639 bp and the mean length of 965.14 bp ( Table 2). The distribution of the unigene numbers in different length ranges were highly uniform (Additional Fig. 1 Fig. 4. These pathways were divided into five groups, which include "Cellular Process," "Environmental Information Processing," "Genetic Information Processing," "Metabolism," and "Organizational systems." A common theme was identified through these six paired comparisons and the dominant pathways were annotated with either a circle or a triangle (Fig. 4). In the "Metabolism" group, the total numbers of DEGs regulating carbon metabolism, photosynthesis, biosynthesis of amino acids, starch and sucrose metabolism and phenylpropanoid biosynthesis were shown in greater abundance and these DEGs may be critical in the formation of flower parts in tree peony. In the "Environmental Information Processing" group, the total number of DEGs for "plant hormone signal transduction" was also high.

Key genes involved in regulating flowering time in tree peony
The

Analysis of MADS-box genes in tree peony
A total of 1833 TFs were identified in our study, which represented 6% of annotated unigenes and fall into 205 TF families classified by plant transcription factor database PlantTFDB. Among these TF gene families, AP2/ERF-ERF (68) was the most abundant TF family, followed by C2H2 (63) and bHLH (61).
Especially, a total of 29 MADS-box genes including 16 MADS-MIKC type were found in tree peony transcriptome data (Fig. 6A).
The MADS-box genes have been involved in many aspects of floral development and flower organs differentiation. To further examine the phylogenetic relationships among the tree peony MADS-box genes and to classified them into the established subfamilies, a phylogenetic Neighbor-joining tree containing MADS-box genes in tree peony, grape, kiwifruit and arabidopsis was made (Fig. 6B). The phylogenetic analysis revealed that 16 MADS-box genes were classified into 10 subfamilies including a new clade which was not reported in arabidopsis. In addition, a few key MADS-box genes described in ABCDE models were also included, such as AP1 (APETALA1), FUL (FRUITFULL), PI (PISTILATA), AP3 (APETALA3), AG (AGAMOUS), STK (SEEDSTICK), SHP (SHATERPROOF) and SEP (SEPALLATA). In the AP1/FUL subfamily (A class), c58746 was belonged to AP1 clade, while c87227 fell within AGL79. Two genes, c70349 and c80146, were assigned to PI/AP3 subfamily, which were functionally classified as B class. c57631, c77234 and c85161 belonged to AG/SHP subfamily in C/D class, which included AG, SHP and STK three subclade. c77234 was classified into AG clade in C class, while c57631 and c85161 were sorted into STK and SHP subclade respectively in D class. Only one gene c72216 was assigned to SEP subfamily in the E functional group. In conclusion, the phylogenetic relationships showed that the tree peony contained all five classes flower parts differentiation MADS-box genes.
To further explore how these regulators related to the flower developmental process, we examined the expression patterns of the 16 MADS-box genes (Fig. 6C). Genes including c77234 (C class) and Out of the 16 MADS-Box genes examined is this study, six of them were specifically selected for determining their specific roles in controlling stamen petalody based on two criteria. They include significant gene expression in the specific flower parts and their fit for the ABCDE model. These six keys genes were AG, AP1-1, AP1-2, AP3, PI, and SEP, respectively. Eight tree peony cultivars classified to stamen petalody and non-stamen petalody groups were used to verify the functions of the six candidate genes by real-time PCR (Fig. 7).
The results showed that SEP had a high-level expression of sepals in both stamen petalody and non- with all the other four class genes to function together. [9][10][11][12][13]23]. The interactions of these five class transcription factors forming tetrame have led to the creation of diverse flower forms in quite a few plant species. In Matthiola incana, the mutated AG alleles corresponded to the double-flower cultivars [24]. In Magnolia stellate, the expression changes in an AG orthologours gene resulted in double flowers [25]. In Phalaenopsis orchid, the tepal became a leaf-like organ when SEP-like gene was silenced by virus-induced silencing. SEP played important roles in floral form formation throughout the developmental process by formation of various multiple protein complexes [26]. In Prunus lannesiana, overexpressed AP1 transgenic plants showed carpelloid stuctures and of petal to stamen-like stuctures [27]. In our study, AG and SEP showed significantly high expression levels in carpel and sepal in both the stamen petalody and non-stamen petalody groups. This indicates that both AG and SEP play important roles in carpels and sepals, respectively, in tree peony regardless of the presence or absence of stamen petalody. A highly significant finding in our study is that both AP3 and PI showed elevated expression levels in stamen petalody cultivars compared to non-stamen cultivars. It is safe to assume that the increased expression levels of AP3 or PI in stamen trigger the stamen petalody, which is an important way to modify the flower forms in tree peony. The inconsistent expressions of AP1-1 and AP1-2 indicate that they may have multiple functions in regulating flower forms that may not be tied stamen petalody. In Arobdoipsis, AP1 was reported to regulate the development of sepals and petals with other MADS-box genes [9]. It is still quite possible that AP1-1 and AP1-2 may play important roles in regulating the formation of petals and sepals in tree peony.
Addition research will be needed to characterize the specific functions of AG, AP1-1, AP1-2 AP3, PI, and SEP during flower formation for more targeted manipulation of flower forms.
RNA-seq is a powerful and efficient approach for gene expression analysis at genome level, especially in the plant species without genomic sequence information. In our study, a total of 11,801 unigenes were detected in all four floral parts. Among them, 14.3% of unigenes were found to be specifically expressed in certain individual floral parts. Among all of the flower parts examined, the unigenes that were specifically expressed in carpels were found in the greatest abundance, followed by stamen, calyx, and petal, respectively. These results suggest that the unigenes specifically expressed in certain individual floral parts might be potential candidates for the regulation of flower part formation in tree peony. Further investigations of each set of these unigenes could lead to targeted manipulation of individual flower part development. A recent study showed that 28,199 differentially expressed genes were identified in Paeonia rockii by comparative transcriptome analysis. Among them, 15 candidate genes regulating the carpel quantitative variation were selected [28]. Another study showed that a set of 4876 DEGs were identified in Asparagus officinalis, including 43 DEGs in hormone response and biosynthesis associated with sex expression and reproduction [29]. A set of 3,7000 DEGs were identified in blueberry, among them 89% of flowering pathway genes, 86% of MADS-box genes, and 84% of cold-regulated genes were up or down regulated in chilled flower buds compared to nonchilled flower buds [30]. A total of 1502 DEGs were detected in sterile flower buds [31].The application of this method has yielded valuable data in tree peony. It is interesting to note

Plant material and floral organs collection
Simple form tree peony cultivar 'Fengdanbai' was selected for this study.

Transcriptome assembly and annotation
Trinity software was used to assembly the clean reads. Firstly, K-mers were obtained by breaking clean reads into short fragments [32]. Secondly, high frequency K-mers were extended from two directions with an overlap of k-1 bp until all K-mers in the library was used out. Thirdly, components was constructed by the linked K-mers (contigs) with an overlap of k-1 bp. Fourthly, De Brujin graphs were simplified and real reads were used to solve the graphs. Finally, transcripts were assembled and the longest transcripts were unigenes.

Differential expression and unigene expression level analysis
The threshold values of false discovery rate and fold change were less than 0.01 and more than 2, respectively [39].

Unigene expression level analysis
To calculate the expression of unigenes, 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 [40,41]. We used the FPKM (fragments per kilobase of gene per million mapped reads) algorithm to normalize the gene expressional abundances in each library [42].

Real-time quantitative PCR
Six specific expressed MADS-box genes were selected to perform Real-time quantitative PCR. Eight tree peony cultivars with different flower forms (classified to stamen petalody and non-stamen petalody groups) were used as the test plant materials. Flower parts in non-stamen petalody groups were separated into carpels, sepals, petals, and stamens, while stamen petalody groups were separated into carpels, sepals, outside whorl petals, and inside-whorl petals (which are produced by stamen petalody). All tree peony cultivars in this study were provided and identified by Zhanying Wang, the Chief Scientist and the curator of the Tree Penoy Germplasm Center at Luoyang Agricultural Science Academy, Luoyang, Henan, China. These specimens can only obtained with a special permit from the Chinese government and are not available to the general public.
A total of 2 μg RNA was extracted from each of these samples. Qualified RNA was tested using the APPlied Biosystems 7500 Real-Time PCR system with a 20 μl reaction volume containing 10 μl of SYBR Premix Ex Taq II (Takara, China), 300 pM final concentration of each primer, and 1 μl template. The reaction conditions were as follows: 30s at 95 ℃ 35 cycles of 5 s at 95 ℃ and 34s 60 ℃. The primers were list in Additional table 1. Ubiquitin was selected as the reference gene [43]. Each sample was performed with three replicates and the relative expressions were calculated using the 2 -ΔΔCt method.

Ethics approval and consent to participate
Not applicable

Consent for publication
Not applicable

Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.   The tree peony cultivar 'Fengdanbai' and the separated flower parts for RNA-sequencing.