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

doi:10.21203/rs.2.18308/v1

Download PDF

Research article

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

https://doi.org/10.21203/rs.2.18308/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

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.

Plant Physiology and Morphology

Plant Molecular Biology and Genetics

Flower form

RNA sequencing

Flower development

MADS-box genes

Flower forms of ornamental plants are typically controlled by multiple sets of genes. The variation of flower parts like petal, sepal, stamen and carpel are critical for the formation of the flower forms in ornamental plants. Natural mutations, human manipulations, and targeted selections have led to many desirable flower forms in ornamental plants for all humans to enjoy today.

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–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–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–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!

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). Hence, it is reasonable to assume that all assembled unigenes have a high degree of integrity.

Tabal 1. Summary of the sequencing results for four floral organs in tree peony
Samples	Number	Read Number	Base Number	GC Conten	%≥Q30
Petal	PE1	30,505,090	9,088,941,390	44.52%	89.99%
	PE2	34,668,124	10,331,991,108	44.57%	90.58%
	PE3	33,988,044	10,138,469,592	45.40%	91.04%
Sepal	SE1	32,785,787	9,774,077,598	44.89%	91.20%
	SE2	31,070,742	9,253,816,890	44.76%	91.13%
	SE3	29,309,927	8,742,359,852	44.83%	90.84%
Stamen	ST1	36,870,554	10,982,297,828	45.40%	92.12%
	ST2	32,823,379	9,781,335,460	45.41%	92.15%
	ST3	35,848,484	10,696,258,334	45.46%	91.78%
Carpel	CA1	34,168,656	10,192,298,092	44.54%	91.52%
	CA2	35,251,448	10,516,365,094	44.32%	91.63%
	CA3	31,763,113	9,479,948,186	44.49%	91.36%
Note: %≥Q30 represents the percentage of the base with the quantity be equal or greater than 30 in clean data.

Table 2. Summary of the assembly for four floral organs in the tree peony
Length Range	Transcript	Unigene
200–300	23,645(13.32%)	18,966(24.95%)
300–500	26,606(14.99%)	15,324(20.16%)
500–1000	45,226(25.48%)	18,128(23.85%)
1000–2000	47,716(26.88%)	13,952(18.36%)
2000+	34,315(19.33%)	9,637(12.68%)
Total Number	177,508	76,007
Total Length	223,165,477	73,357,029
N50 Length	1,903	1,639
Mean Length	1257.21	965.14
Note: N50 Length represents the length of the unigene of N50.

Functional annotation and clustering of gene expression profiles of tree peony floral parts transcriptome

BLAST was used to annotate unigenes by databases including Nr, Swiss-Prot, GO, COG, KOG, eggNOG4.5 and KEGG, resulting in 30,505 (40%) successfully annotated unigenes (Fig. 2A). Nighty eight percent of total annotated unigenes, or 29,872 of 30,505 were perfectly matched in Nr, while the smallest number of unigenes (8,321) were hit in KEGG. All unigenes were classified using a COG database search. In total, 8,321 unigenes were annotated in 25 COG classifications. Among them, the term ‘general function prediction only’ (2042, 24.5%) represented the highest number of unigenes, while ‘extracellular structures’ was represented by zero (Fig. 2B).

For the GO annotation, a total of 17,184 annotated unigenes were classified into three principal categories which were further subdivided into 52 functional terms (Fig. 2C). Among them, ‘cell’ (7113, 41.4%), ‘cell part’ (7113, 41.4%) and ‘organelle’ (5304, 30.9%) were the dominant function terms in ‘cellular component’ principal category. The term which was most represented in ‘molecular function’ principal category was ‘binding’ (9171, 53.4%), followed by ‘catalytic activity’ (8950, 52.1%). Within the ‘biological process’ principal category, unigenes (11857, 69%) belonged to ‘metabolic process’.

Identification of different expression unigenes and clustering of gene expression profiles of tree peony floral parts transcriptome

Differences of gene expression profiles among four floral parts were examined by pairwise comparisons of the 12 libraries. A total of 11,801 (36.7%) unigenes were shared among all libraries, while approximate 14.3% of unigenes (4355) was detectable in individual floral parts. For petal, there were 328 unique unigenes which was the smallest number of unigenes of all four parts. Pistil had the largest number of unique unigenes, followed by stamen (1129) and calyx (935) (Fig. 3A). Those genes expressed in certain floral parts, probably played important roles in development of the coincident floral parts.

Based on the six comparisons of Petal Vs Sepal, Petal Vs Stamen, Petal Vs Carpel, Sepal Vs Stamen, Sepal Vs Carpel and Stamen Vs Carpel, we identified 9592 different expression unigenes (DEGs) in total (Fig. 3C). Among them, comparison of Petal Vs Sepal had the smallest number of 1706 (17.8%) DEGs, while comparison of Sepal Vs Stamen had the largest number of 5517 (57.5%) DEGs followed closely by Stamen Vs Carpal 5310 (55.4%) and Petal Vs Stamen 3505 (36.5%).

To further investigate the gene expression profiles, we performed K-Means clustering of gene expression levels of DEGs. We identified 9 gene clusters with distinctive expression patterns (Fig. 3B). Cluster 7 and Cluster9 containing 499 and 452 unigenes respectively were showing opposite patterns. While Cluster 7 showed gradual increase from petal to pistil, cluster 9 showed gradual decrease, indicating different roles regarding the floral organ development. Cluster 1, 5 and 6 including 1342, 894 and 550 unigenes were showing similar patterns with up-regulated unigenes in stamen, on the contrary, Cluster 2, 3, 4 and 8 containing 1320, 1311, 1096 and 466 respectively with down-regulated unigenes in stamen. These facts indicted that stamen might be the floral part contained more complex molecular process than others.

KEGG pathway enrichment of DEGs

In order to identify active pathways in the four floral organs in tree peony, the DEGs were annotated in the KEGG pathways. KEGG pathway enrichment analysis was performed to search the dominant pathways in the four floral organs. Specific enrichment of genes was obtained for 102 pathways at Petal vs Calyx, for 113 pathways at Petal vs Stamen, for 102 pathways at Petal vs Pistil, 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 listed in the 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 KEGG pathway analysis revealed a general view of the involvement of many biological processes among the floral parts development. To explore the detail expression patterns of the key genes involving in flowering time pathways, BlastT was used to search the orthologous of these key genes in Arabidopsis thaliana. A total of 24,163 pairs of orthologs were identified. Among them, key genes related to flowering time were selected and performed expression analysis (Fig. 5). We found that genes such as red/far-red light signaling pathways (PHYA, PHYB, CRY1, ZTL) were lower expressed in stamen, which indicated that light signaling was less involved in stamen development. The expression level of genes related to floral integrators and flora parts development (AP1, SOCI, FT, FLC, LHY, SVP) were showed expected patterns which were high expressed in petal or in sepal. The blast search also displayed that genes involved in autonomous pathways (FCA, FY, FLD, FVE, FPA) were all shared similar expression patterns which reached peaks in carpal and throughs in stamen, which indicated that these genes more involved in carpal development.

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 c87227 were high expressed in all four floral parts, while 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.

Closer Examination of Six Key Genes Involved in Controlling the Development of Flower Forms Using Stamen Petalody Types

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-stamen petalody cultivars. Meanwhile, AG showed similar expression profiles as SEP, but in carpals instead of sepals. AP3 showed high level expressions in petals and much lower level expressions in stamens in the non-stamen petalody group. However, the expression profiles of AP3 showed high level in both the ow-petals and the iw-petals in the stamen petalody group. The expression profiles of PI are similar to those of AP3. PI showed lower expression in stamen in the non-stamen petalody group and very higher expression in the stamen petalody cultivars. This indicates that PI likely has the similar function as AP3 in regulating stamen petalody in tree peony. The expression profiles of AP1-1 and AP1-2 were quite inconsistent in both stamen petalody types and non-stamen petalody types.

Sixteen MADS-box genes have been successfully identified for the first time in tree peony through our systematic comparison of unigenes using a powerful and efficient RNA-seq method. Six of them exhibited unique expression in specific part(s) of tree peony flower and matched those of the ABCDE model as established in Arabidopsis. Our results have laid a solid foundation toward building a comprehensive model of flower parts formation in tree peony and possibly in other woody ornamental plants. Construction of a flower development model could ultimately lead to more efficient creation of new tree peony cultivars with new floral forms.

According the ABCDE flower development model, A class genes determined the development of sepals; C class genes and D class genes determined the development of carpels and pistil independently; B class genes together with A class genes determined the development of petals; B together with C class genes determined the development of stamens, E class genes were interacted with all the other four class genes to function together. [9–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 that comparisons of Calyx Vs Stamen and Stamen Vs Carpel contained more DEGs than other comparisons. These facts indicated that the molecular progress regulating the development in stamen and carpel involved more complex development evens than in calyx and petal, which partially coincided with the specified unigenes results. K-Means clustering of gene expression levels of DEGs found 9 gene clusters with distinctive expression patterns, in which 7 of them showed up-regulated or down-regulated unigenes in stamen. It is too early to predict which sets of unigenes are more critical for the development of specific flower parts especially. Further investigations will reveal more information about the specific roles of these unigenes in regulating flower form development. Our results indicated that the formation of every flower part is quite complex and may be regulated by at least hundreds of different unigenes. The large numbers of flower part specific regulating unigenes may be linked to the evolution of higher plants due to the unique functions of flowers as the organ for pollination and seed production, and ultimately the survival of plant species. Flower colors, forms, and sizes may play critical roles in attracting pollinators for seed production. Humans enjoy different flower forms for our ‘selfish’ reasons. However, successful production of different flower forms may involve the manipulation of several set of unigenes from

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. AG, AP1-1, AP1-2 and SEP did not show any important role in regulation stamen petalody. However, they may still be important in regulating the formation of other flower parts. Effective manipulation of PI and AG should lead to much faster development of many exciting tree peony flower forms. Our results provided a solid foundation for the successful establishment the flower development model in tree peony based the ABCDE model in arabdopsis and will others develop similar models in highly complex woody plants.

Plant material and floral organs collection

Simple form tree peony cultivar ‘Fengdanbai’ was selected for this study. Flower samples were collected at the early bloom stage in April 2017 from the ten-year old specimens planted at the experimental farm of Henan Agricultural University in Zhengzhou, Henan (34.785259°N, 113.663036°E). Three flowers were taken from three different tree peonies and were defined as biological replicates in this study. Freshly collected flowers were immediately separated into four parts which included carpels, petals, sepals, and stamens. These flower parts were then frozen in liquid nitrogen and stored at -80℃ for subsequent sequencing (Fig. 8). Floral parts were labeled as SE1 (sepal replication one), SE2 (sepal replication two), and SE3 (sepal replication three); PE1 (petal replication one), PE2 ((petal replication two), PE3 (petal replication three); CA1 (carpel replication one), CA2 (carpel replication two), and CA3 (carpel replication three); ST1 (stamen replication one), ST2 (stamen replication two), and ST3 (stamen replication three) for subsequent RNA extraction and sequencing.

RNA extraction and cDNA library construction

The total RNA was extracted from all 12 samples with an Easy-spin^TM kit (Aidlab, Beijing, China) following the instructions of the manufacturer. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was measured with NANOPhotometer Spectrophotometer (Implen, CA, USA). RNA concentration and RNA integrity were assessed using Qubit RNA Assay Kit in Qubit 2.0 Fluometer (Life technologies, CA, USA) and RNA Nano 6000 Assay Kit of the agilent bioanalyzer 2100 system (Agilent technologies, CA, USA), separately.

A total amount of 3 μg of purified RNA were used to construct the cDNA library using NEBNext Ultra RNA library prep kit for Illumine (NEB, USA) following the manufacture instruction. The brief steps are as follow: 1) mRNA was purified from total RNA using poly-T oligo-attached magnetic beads; 2) Fragmentation was carried out using divalent cations under elevated temperature in NEBNext first strand synthesis reaction buffer; 3) First strand cDNA was synthesized using random hexamer primer and M-MuLV reverse transcriptase; 4) Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H; 5) Remaining overhangs were converted into blunt ends via polymerase activities and nebnextadaptor with hairpin loop structure were ligated to prepare for hybridization; 6) PCR was performed with phusion high-fidelity DNA polymerase, universal PCR primers and index primer. The library was checked by Qubit 2.0 and Agilent 2100 systems to examine the concentration and insert size. Sequencing was performed by Illumina HiSeq with pair-end reads of 150 bp.

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.

Unigenes were annotated in seven databases including Nr, COG, GO, KEGG, KOG, Pfam, Swissprot and eggNOG by BLAST with E value equal or less than 1e-5 [33-35]. KOBAS2.0 was used to obtain KEGG Orthology results [36, 37]. Predicted amina acid sequences of unigenes were blasted against the Protein family (Pfam) database by HMMER with E value equal or less than 1e-10 [38].

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.

Funding

This research project was supported by a grant from the Key Science Research Projects (KSRP) of Colleges and Universities of Henan Province (19A220003) and another grant from the Science Creation Projects (SCP) of Henan Agricultural University (KJCX2018A05). KSRP’s support included sample collection, RNA sequencing and analysis, and lab materials and supplies while SCP provided financial support for real time PCR in Figure 7, personnel, travel, the publishing fee, and all other miscellaneous expenses.

Authors’ contributions

JL initiated the hypothesis, designed the experiment, coordinated the sample collection and analysis, and wrote the first draft. YZ and HW collected samples and performed real time PCR; ZW provided a part of the literature review, YL conducted a part of the data analysis; YHL wrote a part of the materials and methods, and provided critical input in the experimental design. YGG wrote a part of the abstract, results and discussion. All authors have read and approved the manuscript.

Acknowledgments

We thank Biomarker Technologies Corporation (Beijing, China) for providing technical support.

AG: Agamous; AP1: Apetala1; AP3: Apetala3; COG: Cluster of orthologous groups of proteins; CRY1: Cryptochrome 1 gene; DEG: Different expression unigene; eggNOG: Nonsupervised orthologous groups; FCA: Arabidopsis flowering time gene; FLC: Flower locus C; FLD: Flower locus D; FPA: Autonomous pathway independent flowering time gene; FT: Flower locus T gene; FPKM: Fragments per kilobase of gene per million mapped reads; FUL: Fruitfull; FVE: A classical flowering time locus of the autonomous pathway; GO: Gene ontology; KEGG: Kyoto encydopedia of genes and genomes; KOG: Eukaryotic ortholog groups; LHY: Late elongated hypocotyl; Nr: Non-redundan protein; PFam: Protein family; PHYA: Phytochrome-directed repression A gene; PHYB: Phytochrome-directed repression B gene; PI: Pistilata; SHP: Shaterproof; SEP: Sepallata; SOC1: Suppressor of overexpression of constans 1; STK: Seedstick; SVP: Short vegetative phase; Swiss-Prot: The manually annotated and reviewed protein sequence database; TFs: Transcription factors; ZTL: Zeitlupe.

Formar F, de Montaigu A, Coupland G. Snapshot: control of flowering in arabidopsis. Cell 2010; 3:550.
Chen J, Källman T, Ma X, Gyllenstrand N, Zaina G, Morgante M et al. Disentangling the roles of history and local selection in shaping clinal variation of allele frequencies and gene expression in norway spruce (Picea abies). Genetics 2012; 191:865-881.
Burgarella C, Chantret N, Gay L, Prosperi JM, Bonhomme M, Tiffin P et al. Adaptation to climate through flowering phenology: a case study in Medicago truncatula. Mol Ecol 2016; 25:3397-3415.
Jung WY, Lee A, Moon JS, Kim Y, Cho HS. Genome-wide identification of flowering time genes associated with vernalization and the regulatory flowering networks in chinese cabbage. Plant Biotechnol Rep 2018; 12:347-363.
Li YX, Li C, Bradbury PJ, Liu X, Lu F, Romay CM et al. Identification of genetic variants associated with maize flowering time using an extremely large multi-genetic background population. The Plant Journal 2016; 86:391-402.
Wang J, Qiu Y, Cheng F, Chen X, Zhang X, Wang H et al. Genome-wide identification, characterization, and evolutionary analysis of flowering genes in radish (Raphanus sativus). Bmc Genomics 2017; 18:981.
Peng Z, Lu Y, Li L, Zhao Q, Feng Q, Gao Z et al. The draft genome of the fast-growing non-timber forest species moso bamboo (Phyllostachys heterocycla). Nat Genet 2013; 45:456.
Meyerowitz EM, Bowman JL, Brockman LL, Drews GN, Jack T, Sieburth LE et al. A genetic and molecular model for flower development in Arabidopsis thaliana. Development 1991; 113:157-167.
Pelaz S, Gustafson-Brown C, Kohalmi SE, Crosby WL, Yanofsky MF. Apetala1 and sepallata3 interact to promote flower development. Plant Journal for Cell & Molecular Biology 2001; 26:385-394.
Angenent GC, Franken J, Busscher M, Van D, Went JV. A novel class of mads box genes is involved in ovule development in petunia. Plant Cell 1995; 7:1569-1582.
Bowman JL, Smyth DR, Meyerowitz EM. Genetic interactions among floral homeotic genes of arabidopsis. Development 1991; 112:1-20.
Pelaz S, Ditta GS, Baumann E, Wisman E, Yanofsky MF. B and c floral organ identity functions require sepallata MADS-box genes. Nature 2000; 405:200-203.
Colombo L, Franken J, Koetje E, Went JV, Dons HJ, Angenent GC et al. The petunia mads box gene fbp11 determines ovule identity. Plant Cell 1995; 7:1859-1868.
Ditta G, Pinyopich A, Robles P, Pelaz S, Yanofsky MF. The sep4 gene of Arabidopsis thaliana functions in floral organ and meristem identity. Curr Biol 2004; 14:1935-1940.
Jofuku KD, den Boer BG, Van MM, Okamuro JK. Control of arabidopsis flower and seed development by the homeotic gene apetala2. Plant Cell 1994; 6:1211-1225.
Qin Q, Yin T, Chen J, Xie M, Zhang S. Apetala3/deficiens and pistillata/globosa genes with floral development of plant. Chinese Journal of Cell Biology 2006; 28:571-576.
Yanofsky MF, Ma H, Bowman JL, Drews GN, Feldmann KA, Meyerowitz EM. The protein encoded by the arabidopsis homeotic gene agamous resembles transcription factors. Nature 1990; 346:35-39.
Honma T, Goto K. Complexes of mads-box proteins are sufficient to convert leaves into floral organs. Nature 2001; 409:525-529.
Cheng F, Li J, Chen D: Chinese flare mudan, 1 edn. Beijing: Forestry publishing house of China; 2005.
Li J, Zhang X, Zhao X: Tree penoy of china. Beijing: Encyclopedia of China publishing house; 2011.
Yan Han X, Wang L, Shu Q, An Liu Z, Xia Xu S, Tetsumura T. Molecular characterization of tree peony germplasm using sequence-related amplified polymorphism markers. Biochem Genet 2008; 46:162-179.
Han XY, Wang LS, Shu QY, Liu ZA, Xu SX, Tetsumura T. Molecular characterization of tree peony germplasm using sequence-related amplified polymorphism markers. Biochem Genet 2008; 46:162-179.
Coen ES, Meyerowitz EM. The war of the whorls: genetic interactions controlling flower development. Nature 1991; 353:31-37.
Nakatsuka T, Koishi K. Molecular characterization of a double-flower mutation in Matthiola incana. Plant Sci 2018; 268:39-46.
Zhang B, Liu Z, Ma J, Song Y, Chen F. Alternative splicing of the agamous orthologous gene in double flower of Magnolia stellata (Magnoliaceae). Plant Sci 2015; 241:277-285.
Pan ZJ, Chen YY, Du JS, Chen YY, Chung MC, Tsai WC et al. Flower development of phalaenopsis orchid involves functionally divergent sepallata-like genes. New Phytol 2014; 202:1024-1042.
Chen M, Lin I, Yang C. Functional analysis of three lily (Lilium longiflorum) apetala1-like mads box genes in regulating floral transition and formation. Plant Cell Physiol 2008; 49:704-717.
Liu N, Cheng F, Zhong Y, Guo X. Comparative transcriptome and co-expression network analysis of carpel quantitative variation in Paeonia rockii. Bmc Genomics 2019; 20:1-18.
Li S, Zhang G, Zhang X, Yuan J, Deng C, Gao W. Comparative transcriptome analysis reveals differentially expressed genes associated with sex expression in garden asparagus (Asparagus officinalis). Bmc Plant Biol 2017; 17:143.
Song G, Chen Q. Comparative transcriptome analysis of nonchilled, chilled, and late-pink bud reveals flowering pathway genes involved in chilling-mediated flowering in blueberry. Bmc Plant Biol 2018; 18:98.
Liu H, Tan M, Yu H, Li L, Zhou F, Yang M et al. Comparative transcriptome profiling of the fertile and sterile flower buds of a dominant genic male sterile line in sesame (Sesamum indicum). Bmc Plant Biol 2016; 16:250.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I et al. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol 2011; 29:644.
Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W et al. Gapped blast and psi-blast: a new generation of protein database search programs. Nucleic Acids Res 1997; 25:3389-3402.
Apweiler R, Bairoch A, Wu C, Barker W, Boeckmann B, Ferro B et al. Uniprot: the universal protein knowledgebase. Nucleic Acids Res 2004; 32:115-119.
Koonin E, Fedorova N, Jackson J, Jacobs A, Krylov D, Makarova K et al. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biol 2003; 5:60.
Kanehisa BM, Goto S, Kawashima S. Et al: the KEGG resource for deciphering the genome. Nucleic Acids Res 2013; 32:D277-D280.
Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S et al. Kobas 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res 2011; 39:W316-W322.
Finn RD. Pfam: the protein families database. Nucleic Acids Res 2014; 42:D222-D230.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol 2010; 11:R106.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short dna sequences to the human genome. Genome Biol 2009; 10:1-10.
Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BM et al. Ebseq: an empirical bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics 2013; 29:1035-1043.
Trapnell C, Williams B, Pertea G, Mortazavi A, Kwan G, van Baren M et al. Transcript assembly and quantification by RNA-seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol 2010; 28:511-515.
Liu CJ, Wang SL, Xue JQ, Zhu F, Ren X, Li M et al. Molecular cloning of ubiquitin protein gene and study on this gene as reference gene in tree peony. Acta Hortic Sin 2015; 42:1983-1992.

Download PDF

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

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

Status:

Version 1

Abstract

Figures

Background

Results

Transcriptome sequencing and assembly

Functional annotation and clustering of gene expression profiles of tree peony floral parts transcriptome

Identification of different expression unigenes and clustering of gene expression profiles of tree peony floral parts transcriptome

KEGG pathway enrichment of DEGs

Key genes involved in regulating flowering time in tree peony

Analysis of MADS-box genes in tree peony

Closer Examination of Six Key Genes Involved in Controlling the Development of Flower Forms Using Stamen Petalody Types

Discussions

Conclusions

Methods

Declarations

Abbreviations

References

Status:

Version 1