Pistils may be required for anthocyanin synthesis —— the whole-transcriptome analysis of mutant and normal capitula of Chrysanthemum morifolium


 Background Chrysanthemum morifolium is one of the most economically important and popular floricultural crops in Asteraceae. Chrysanthemums have many different flower colors and shapes. However, the molecular mechanism controlling the development of chrysanthemum floral colors and shapes is still an enigma. We obtained a cut chrysanthemum variety with mutant capitula in which the ray florets became green and the inside pistils became vegetative buds, while normal capitula have many rounds of purple ray florets and few disc florets. ﻿Results We conducted whole-transcriptome analysis of differentially expressed genes (DEGs) between the mutant and normal capitula using third-generation and second-generation sequencing techniques. We identified DEGs between the mutant and normal capitula to reveal important regulators underlying their differential development. Regulatory genes involved in the photoperiod pathway and the control of floral organ identification as well as important functional genes in the anthocyanin synthesis pathway were also identified. Therefore, a list of candidate genes for studying flower development and anthocyanin synthesis in chrysanthemums was generated. Qualitative analysis of pigments in the florets of normal and mutant capitula revealed anthocyanins were synthesized and accumulated in the florets of normal capitula, but not in the florets of mutant capitula. It was indicated that pistils may be required for anthocyanin synthesis in chrysanthemums. Conclusions These results will help to elucidate the molecular mechanisms of floral organ development and will contribute to the development of techniques for studying flower shape and color regulation to promote breeding in chrysanthemum.


Introduction
Chrysanthemum morifolium is one of the most economically important and popular floricultural crops in Asteraceae, and ranks second in the cut flower trade after rose [1].
The head-like inflorescence (capitulum) resembling a single large flower is the main ornamental part of C. morifolium, and is also regarded as the key innovation behind the evolutionary success of the Asteraceae [2]. The typical capitulum of chrysanthemum is formed by two morphologically distinct florets: the marginal ray florets and the central disc florets. Ray florets have ligulate and zygomorphic colorful corollas (petals) and aborted stamens, which function in attracting pollinators. The disc florets have radially symmetrical colorless corollas and their fertile pollens are hermaphroditic and are used for reproduction in chrysanthemum (Additional file 1). The colors and shapes of the flowers are the most visible and amazing products of evolution, and also connect humans to nature [3].
Flowering is a key developmental process in the plant life cycle that is very complex and is controlled by endogenous factors and environmental cues. As currently understood, flower development contains three phases: flowering determination, flower evocation, and floral organ development [4]. In Arabidopsis, tremendous progress has been achieved toward understanding the molecular mechanisms involved in flower development [5,6].
ABCE models have revealed that A-class together with E-class genes specify sepal identity, A-class, B-class and E-class genes specify petals, B-class, C-class and E-class genes determine stamens, and C-class and E-class genes determine carpel/gynoecium organ identity [7]. With the notable exception of A-class genes, all of these genes belong to the MADS-box family of transcription factors including the AP1, AP3, PI, AG and SEP genes.
The splendid colors presented mainly by flower petals have enabled plants to constantly develop new showy traits and prosper throughout millions of years of evolution.
Anthocyanins and carotenoids are the two major groups of pigments generated in plant petals. Anthocyanins are accumulated in the vacuoles of petal epidermal cells and confer orange-to-violet colors in flowers [8]. Beside attracting pollinators, anthocyanins also protect against damage from UV irradiation [9]. Anthocyanins provide chrysanthemum ray florets with colorful bright colors to attract pollinators, which improves the success rate of cross pollination between different species or varieties and promotes cultivar groups with large varieties of flower types in C. morifolium. Anthocyanins enhance the ornamental value of chrysanthemums, and many cut flower and pot flower varieties with bright colors are on sale every year to satisfy the demand of the market. Understanding the mechanism of anthocyanin biosynthesis and its regulation will contribute to the cultivation and improvement of new color varieties of chrysanthemums.
In chrysanthemums, a few regulatory genes involved in flower development have been isolated, such as MADS-box, TCP, and WUS-like genes [10][11][12]16]. Some important functional genes and transcription factors involved in the anthocyanin biosynthesis pathway have also been characterized including ANS, F3′H, F3H and MYB-like genes [13][14][15][16]. However, chrysanthemums have complex capitula containing two morphologically distinct florets and through a long period of breeding a variety of flower shapes and diverse colors have been created. The mechanism involved in chrysanthemum flower evolution and development is extremely complicated and little is known about it.
The development of RNA-seq technology has greatly improved transcriptomic studies in chrysanthemums [1]. However, the read-length offered by second-generation highthroughput sequencing platforms is much shorter than the typical length of a eukaryotic mRNA. Additionally, the differences in transcript abundance and the presence of different isoforms make the assembly of transcriptomes from short reads extremely challenging [17].
Despite these problems, Hideki Hirakawa et al. performed de novo whole-genome assembly in C. seticuspe using the Illumina sequencing platform and Chi Song et al.
sequenced the diploid C. nankingense genome using the Oxford Nanopore long-read technology; however, no more than 40% of the transcriptome sequencing reads from C. morifolium can be mapped to each of these two genome sequences, probably because of the extreme variation in chromosome ploidy and biological characteristics [18,19]. Thirdgeneration sequencing technology has dramatically increased the length of sequencing reads, which enables the precise location and sequencing of repetitive regions and isoforms within a single read.
Recently, we obtained a mutant plant of the cut chrysanthemum variety C. morifolium 'ZY' with both mutant and normal capitula. In the mutant capitula the ray florets became green and the inside pistils became vegetative buds, while the normal capitula had many rounds of purple ray florets and few disc florets. In this study, we analyzed a mixed sample of normal and mutant flowers, leaves, stems and roots of 'ZY' with the single-molecule long-read sequencing technology from Pacific Biosciences (PacBio). Based on the results, transcriptional sequencing and analysis of the mutant and normal capitula were performed using second-generation sequencing technology and RNA-Seq quantification analysis.
Thus, we combined third-generation and second-generation sequencing to generate a more complete/full-length C. morifolium transcriptome.
Based on the transcriptional sequencing and analysis, we identified differentially expressed genes (DEGs) between mutant and normal capitula to reveal important regulators controlling their differential development. Regulatory genes involved in the photoperiod pathway and the control of floral organ identification as well as important functional genes in the anthocyanin synthesis pathway were also identified to create a list of candidate genes for studying flower development and anthocyanin synthesis in chrysanthemums. These results will be helpful for elucidating the molecular mechanisms of floral organ development and will contribute to the development of techniques for studying flower shape and color regulation, as well as breeding and molecular biology in chrysanthemum.

Sequencing and assembly
As shown in Figure 1, in the mutant capitula the ray florets became green and the inside pistils became vegetative buds, while the normal capitula had many rounds of purple ray florets and few disc florets. We analyzed a mixed sample of normal and mutant capitula, leaves, stems and roots of 'ZY' using PacBio sequencing and then analyzed the normal and mutant capitula separately using the Illumina paired-end sequencing technology. The resulting sequences were assembled into 130,097 isoforms with an N50 of 3013 bp and average length of 2510 bp (Table 1).

Gene annotation and functional classification
A total of 118,589 isoforms were annotated by BLAST in at least one of the four databases searched (nr, Swiss-Prot, KOG, and KEGG), leaving 11,508 (8.85%) isoforms without annotation. In total, 118,043, 101,048, 87,630, and 54,245 isoforms were annotated in the nr, Swiss-Prot, KOG, and KEGG databases, respectively. Gene Ontology (GO) has three ontologies-molecular functions, cellular components, and biological processes-that facilitate gene annotation and analysis. A total of 36,144 isoforms were classified into 47 functional categories, including 19, 17 and 11 GO categories in the biochemical process, cellular component and molecular function ontologies, respectively.
The dominant categories in the biochemical process, molecular function and cellular component ontologies were 'metabolic process' (20,871), 'catalytic activity' (22,818) and 'cell' (11,887), respectively, indicating that numerous metabolic activities were activated during the development of chrysanthemum capitula and that this process was regulated by a wide range of genes that interacted within cells. In addition, we observed a high percentage of genes from the 'cellular process', 'binding', and 'cell part' categories, but few from 'locomotion', 'transcription factor activity, protein binding', and 'extracellular matrix (2,089 members). In addition, 1,339 isoforms were assigned to 'Plant hormone signal transduction'.

Comparison of the transcriptomes of normal and mutant capitula
The set of isoforms common to normal and mutant capitula In total, 124,284 isoforms were shared by normal and mutant capitula ( Figure 5A). By contrast, 3269 and 955 isoforms showed specific expression in normal and mutant capitula, respectively. Therefore, more genes were expressed in normal capitula than in mutant capitula, because the pistils were mutated to vegetative buds in mutant capitula.

DEGs between mutant and normal capitula
The transcriptomes of the normal and mutant capitula were compared, and the resulting reads were mapped to the reference transcriptome. A total of 35,419 DEGs (8,232 upregulated and 27,187 down-regulated in the mutant capitula relative to the normal capitula) were identified between the normal and mutant capitula ( Figure 5B). The correlation coefficient of gene expression between the normal and mutant capitula was 0.8897, which was determined by an algorithm developed from the correlation scatter plot.
A total of 131 DEGs were specifically expressed in the mutant capitula relative to the normal capitula, including TCP1 and AP2/ERF domain-containing genes. Conversely, 2,132 genes were not expressed in the mutant capitula but were specifically expressed in the normal capitula, including some important transcription factor genes (MYB, GRAS, and BTF3 genes), ubiquitin-conjugating enzyme genes, zinc finger protein genes and many genes without annotations. Annotation information on the DEGs specifically expressed in the mutant and normal capitula is provided in Additional files 2 and 3, respectively. These genes may play important roles during the flower development process especially during pistil determination and development in chrysanthemums. Additionally, these genes should be important candidate regulators in chrysanthemum flower development. Therefore, functional research should be conducted on these genes in the future.
GO and KEGG pathway enrichment analyses were conducted on the DEGs to identify differences in biological processes and pathways between the mutant and normal capitula.
It was found that 256 genes enriched in the 'reproduction' term (GO:0000003) in biochemical processes were all remarkably down-regulated in the mutant capitula relative to the normal capitula, among which 11 genes were specifically expressed in the normal capitula including WD40 and UBA1-like protein-encoding genes. These genes may play important roles in the regulatory pathways related to reproduction in chrysanthemums (Additional file 4).
In total, 6,733, 7,216, and 3,879 DEGs were enriched in the biological process, molecular function and cellular component categories, respectively (Additional files 5-7).
In the biological process category, the dominant terms were 'metabolic process'

Identification of genes involved in the photoperiod pathway in chrysanthemum
As a typical short-day plant, chrysanthemum can flower in response to a single short day. identified in this study, which is an upstream regulatory gene of LFY. Interestingly, its expression level was significantly up-regulated in the mutant capitula relative to the normal capitula. As an 'A-class'-like gene, the expression of AP1 is directly activated by LFY [20,21]. Homologs of AP1 were identified, which were all remarkably up-regulated in the mutant capitula relative to the normal capitula.
As another 'A-class' gene, AP2 is not a MADS-box transcription factor. Two homologs of AP2 were identified, both of which were up-regulated in the normal capitula relative to the mutant capitula. In most core eudicot species, B-class genes include three different lineages: PI, euAP3 and TM6; however, TM6-like genes seem to have been lost in Arabidopsis and Antirrhinum [22]. In this study, homologs of PI and AP3 were identified, but the expression of TM6-like genes was not detected. In contrast, the expression of TM6like genes was detected in chrysanthemums in earlier studies [23,24]. We also identified homologs of the C-class gene AG and E-like MADS-box genes in this study. We found that Interestingly, all four DFR homologs and one CHI homolog showed extremely low expression levels (close to 0) in the mutant capitula relative to the normal capitula. Realtime PCR showed that the CHI homolog was not expressed in the mutant capitula but was expressed in the normal capitula. Annotation information for these genes in the anthocyanin biosynthetic pathway is listed in Additional file 12.
The pigment types and contents in flowers determine the variety of flower colors.
Qualitative analysis of pigments in the florets of normal and mutant capitula was performed by HPLC. Anthocyanins were detected in the florets of normal capitula, but not in the florets of mutant capitula ( Figure 8B). As shown in Figure 8A, the detectable anthocyanins mainly included delphinidin, cyanidin, petunidin, pelargonidin, peonidin, and malvidin.  Figure 9).

Discussion
After the pistils mutated to become vegetative buds, anthocyanins were not synthesized, which indicated pistils may be required for anthocyanin synthesis in chrysanthemums As one of the major flower pigment groups in higher plants, anthocyanin synthesis and accumulation is an integral part of flower development in most plant species and is tightly linked with petal cell expansion. It is thought that the activation of the anthocyanin synthesis pathway during petal development requires both environmental and endogenous signals.
In Petunia hybrida, gibberellins (GAs), sugars and light were revealed to be required to induce the transcription of anthocyanin synthesis genes and accumulate pigment in the developing corolla. Research results have also indicated that GAs, sugars and light are involved in the regulation of various pathways to complete the entire flower development process [25]. David Weiss presented evidence that in the early stages of Petunia hybrida flower development, GAs produced by the anthers controlled anthocyanin synthesis and accumulation in the corolla by activating the transcription of related genes in the anthocyanin synthesis pathway, and the removal of stamens at the early stage of flower development inhibited anthocyanin synthesis in corollas [26]. However, in C. morifolium, the capitula of many varieties, such as those with the 'Pingpang' shape, have only ray florets, because of the long period of double flower breeding. Interestingly, in C. morifolium 'ZY', most capitula have ray florets without disc florets. Therefore, in chrysanthemums, anthocyanin synthesis may be associated with pistils.
In the mutant capitula, the pistils of ray florets mutated to become vegetative buds and the corollas of ray florets became green because of a lack of anthocyanins. When the pistils of ray florets mutated into vegetative buds, anthocyanin synthesis in the corolla may have been blocked. Therefore, it was hypothesized that pistils may be required for anthocyanin synthesis in chrysanthemums.

Many transcription factors are implicated in the developmental regulation of florets in chrysanthemums
In this study, 963 of 3,921 detected important transcription factors genes were dramatically differentially expressed between the normal and mutant capitula, including members of the ERF, C2H2, MYB, bHLH and WRKY transcription factor families. This C2H2 zinc finger protein genes functioning as transcription factors play important roles in many biological processes related to plant growth and development, hormone signaling and stress responses [27]. In many plants, C2H2 zinc finger protein genes are involved in salt, cold, drought and oxidative tolerance, light stress and pathogen defense [28].
Additionally, some members participate in the developmental regulation of flowers. For example, SIZF2 controls flower and leaf shape in A. thaliana. SUPERMAN (SUP) determines the boundary between the stamen and carpel whorls, suppresses class B gene expression and promotes stem cell termination in the fourth whorl of A. thaliana flowers [28][29][30]. In this study, 35 C2H2 genes showed no expression in the mutant capitula, which indicated that these C2H2 zinc finger protein genes may have important functions in the determination and growth regulation of pistils or anthocyanin synthesis and accumulation in the ray floret corolla in chrysanthemums.
As one of the largest transcription factor families, basic leucine zipper (bZIP) family transcription factors play key roles in controlling plant development and stress responses [31]. bZIP genes have been reported to be involved in various flower developmental processes in plants, including pollen development, the floral transition and flower initiation [32,33]. In this study, the expression deficiency of six bZIP genes in mutant capitula suggested that these bZIPs are also important regulators of ray floret development in chrysanthemums.
MYB transcription factors comprise one of the largest and most diverse transcription factor families in the plant kingdom, exist widely in eukaryotes, and play essential roles in a wide range of physiological and biochemical processes controlling plant growth and development [34,35]. bHLH transcription factors possess the highly conserved bHLH domain including a basic region and a HLH region, and play important roles in plant growth and development, metabolic regulation, and responses to environmental changes. The regulatory functions of bHLH transcription factors in active secondary metabolism especially anthocyanin synthesis have been a focus of research [36].
It was revealed that R2R3-MYB, bHLH, and WD40 proteins form a ternary complex called the MBW complex that regulates the anthocyanin biosynthesis pathway in model plants [37]. We found that five bHLH genes and seven MYB genes showed no expression in the mutant capitula, which indicated that these genes are essential regulators of normal

Carpels may originate from vegetative buds
The female reproductive organ, named the carpel, is specific to the angiosperms, or flowering plants. The carpel is probably the major factor contributing to the success of angiosperms, because it allowed angiosperms to diversify from an unconfirmed, possibly gymnosperm-like ancestor to form the more than 300,000 species alive today [39]. As in most species, the chrysanthemum carpel is differentiated into stigma, style, and ovary tissues. Flower developmental mechanism studies show that flowers develop from a floral apex. On the flanks of the floral apex, floral organ development starts in a centripetal sequence: the outermost organs are initiated first and the innermost organs are initiated last. However, the evolutionary sequence of floral organs is reversed. The ovules in the center of the flower are the evolutionarily oldest organs and the perianth organs at the periphery are the youngest [40]. Each organ has its own evolutionary history.
In this study, the pistils in the mutant capitula were mutated and became vegetative buds.
In the early stages, the vegetative buds mutated from the pistils looked very much like normal pistils ( Figure 1). As the mutant florets grew, the two seemingly bifurcate stigma became two tender leaves and the center part between the two seemingly bifurcate stigma grew into a shoot apex. Here we present a bold hypothesis; that the pistil is derived from a vegetative bud. The predecessor of the stigma may be tender leaves developed from the phyllopodium of the vegetative bud; thus, the ovule would have evolved out of the apical meristem of the vegetative bud. This hypothesis supports the inference that the ovule did not originate from a phyllome, and each ovule represents a short shoot [41,42]. Additionally, the dorsiventral gene expression in the outer integument of the ovule indicates that the outer integument of ovules has a leaf-derived origin, which is consistent with the inference that the stigma is derived from leaves [43,44].

Conclusions
In our study, comparative transcriptome analysis revealed significant differences in gene expression and signaling pathways between the mutant and normal capitula. We identified DEGs between the mutant and normal capitula to reveal important regulators underlying their differential development. The transcription factors specifically expressed in the normal capitula in this study are important candidate genes to explore the regulatory molecular mechanisms in chrysanthemum flower development. Regulatory genes involved in the photoperiod pathway and the control of floral organ identification as well as important functional genes in the anthocyanin synthesis pathway were also identified. Qualitative analysis of pigments in the florets of normal and mutant capitula revealed anthocyanins were synthesized and accumulated in the florets of normal capitula, but not in the florets of mutant capitula. It was indicated that pistils may be required for anthocyanin synthesis in chrysanthemums. This study represents the first step in exploring the molecular mechanism of floral organ development and will contribute to the development of techniques for studying flower shape and color regulation to promote breeding in chrysanthemum.

Plant materials and RNA extraction
The tissues (normal and mutant capitula) used in this study were obtained from a cut- RNA quantity and quality assessment was performed using a NanoDropND2000 instrument (Thermo Scientific).

Library construction, PacBio sequencing and Data processing
Library construction and PacBio sequencing were conducted with the PacBio Sequel system (PacBio, CA, USA). Briefly, 1 μg of total RNA extracted from each of the five tissues was equally pooled together and the mRNA was enriched by Oligo (dt) magnetic beads. Second, the low-quality isoforms were further corrected using Illumina short reads obtained from the same samples with the LoRDEC tool (version 0.8) [45]. Then, the final transcriptome isoform sequences were filtered by removing the redundant sequences with the CD-HIT-v4.6.7 software using a threshold of 0.99 identities.

Illumina sequencing
The total RNA isolated from the normal and mutant capitula was used for Illumina sequencing on an Illumina HiSeq™ 2000 system (Illumina, San Diego, CA, USA). We purified the poly (A) mRNAs, fragmented them into small pieces, and then synthesized the first-and second-strand cDNAs.
After the double-stranded cDNAs were purified and resolved for end reparation and poly (A) tail addition, the short fragments were connected with sequencing adapters. Briefly, a cDNA library with average insert sizes of 300-500 bp was created and cDNA sequencing was performed using the Illumina HiSeq™ 2000 system, with paired end 2×100 nt multiplexes.  [38]. GO annotation was also performed with the Blast2GO software [46]. Isoforms ranking in the 20 highest scores and no shorter than 33 HSP (High-scoring Segment Pair) hits were selected for Blast2GO analysis. Then, functional classification of the isoforms was performed using the WEGO software [47].

Analysis of chrysanthemum transcriptome sequencing results
RNA-Seq quantification analysis using the number of reads per kilobase of exon model per million mapped reads (RPKM) was performed to calculate the expression level of each isoform [39]. The chrysanthemum transcriptome of the samples was used as a reference for the screening and analysis of DEGs. A rigorous algorithm was created based on the method of Audic et al. to screen the DEGs [48]. The false discovery rate (FDR) was used to affirm the threshold of the P-value in multiple tests and analyses [41]. The absolute value of log2 (ratio) ≥ 2 and FDR < 0.05 were used as the thresholds to determine significant differences in gene expression [49]. Only the DEGs with a minimum of a twofold change in expression level were used in the differential gene expression analysis.

Alternative splicing detection
To analyze alternative splicing events in the transcript isoforms, the Coding Genome reconstruction Tool (Cogent) was employed to partition the transcripts into gene families based on k-mer similarity and then each family was reconstructed into a coding reference genome with De Bruijn graph methods [50]. The alternative splicing events of the transcript isoforms were analyzed using the SUPPA tool [51].

Gene expression analysis based on qRT-PCR
Total RNA was extracted from the normal and mutant capitula as described above. First, we treated the total RNA with DNase (Promega, USA), and then subjected it to reverse transcription into cDNA using a reverse transcription system (Tiangen, China  [52].

Qualitative analysis of pigments in chrysanthemum flowers
The analysis of anthocyanin profiles in the normal and mutant capitula was performed using high pressure liquid chromatography. A 1.0-10.0 g sample was ground into fine powder in liquid N2 and homogenized into 50 ml anthocyanin extracts [ethyl alcohol: distilled water: hydrochloric acid (2:1:1, v/v/v)] assisted by sonication at 20 °C for 30 min [53]. Then, the mixture was heated in boiling water for 1 h using a water bath and centrifuged at 16000 ×g for 10 min at 20 °C. min. The injection volume was 20 μL and the photodiode array detector was set at 530 nm for anthocyanins [54]. Three biological replicates of each sample type were analyzed.