Floral Development Stage-Specific Transcriptomic Analysis Reveals the Formation Mechanism of Different Shapes of Ray Florets in Chrysanthemum

The formation mechanism of different ray floret shapes of chrysanthemum (Chrysanthemum × morifolium) remains elusive due to its complex genetic background. C. vestitum, with the basic ray floret shapes of the flat, spoon, and tubular types, is considered a model material for studying ray floret morphogenesis. In this study, the flat and tubular type lines of C. vestitum at specific stages were used to investigate the key genes that regulate morphological differences in ray florets. We found that the expression levels of genes related to auxin synthesis, transport, and response were generally higher in the tubular type than in the flat type. CvARF3 was highly expressed in the flat type, while CvARF5 and CvARF6 were highly expressed in the tubular type. Additionally, the transcription levels of Class B and E genes closely related to petal development, including CvPI, CvAP3, Cvdefh21, CvSEP3, and CvCDM77, were expressed at higher levels in the tubular type than the flat type. Based on the results, it is proposed that auxin plays a key role in the development of ray florets, and auxin-related genes, especially CvARFs, may be key genes to control the morphological difference of ray florets. Simultaneously, MADS-box genes are involved in the co-regulation of ray floret morphogenesis. The results provide novel insights into the molecular mechanism of different petal type formation and lay a theoretical foundation for the directional breeding of petal type in chrysanthemums.


Introduction
Chrysanthemum (Chrysanthemum × morifolium) has a wide variety of flower types and is a valuable ornamental and commercial crop. However, the genetic mechanism of chrysanthemum flower pattern formation remains elusive because of the complex genetic background. The most attractive part of the chrysanthemum is the unique capitulum composed of central disc florets and peripheral ray florets, whose morphology is susceptible to change due to internal and external factors [1]. The morphology and number of the ray florets principally determine the ornamental traits of the chrysanthemum [2]. The ray floret shapes are defined as the petal types, which are divided into flat, spoon, and tubular types according to the corolla tube merged degree (CTMD) [3,4]. Other complex petals, such as marginal elaboration and appendages of the corolla, evolved on the basis of these three basic petal types and were combined with their number, orientation, and location to form a rich variety of flower patterns. Therefore, the regulation mechanism of CTMD in ray florets could lay the foundation for the analysis of complex petal types and flower pattern formation. Petals experience four main processes during development: initiation, identity determination, morphogenesis, and maturation [5]. Plants construct different morphological petal structures by controlling the expression patterns of genes [6][7][8][9]. Therefore, elucidating the development and evolution of petals is the key to deciphering the diversification of petals and flowers [10]. The complex mechanism of single floral petals has been explored deeply, but the genetic mechanism of morphological variation of ray florets in chrysanthemums is still a mystery.
In order to find the key genes controlling ray floret types, a lot of forward genetics studies have been conducted. CTMD was an important morphological index for defining petal types [4] and could be described by a B-2 genetic model via two additive-dominance major genes [11]. Three major quantitative trait loci (QTLs) controlling CTMD were detected using a high-density genetic linkage map [12]. In genetic analysis of gerbera (Gerbera hybrida), it was found that the laciniated outer corolla lips of the ray floret were controlled by a dominant gene [13], while it was suggested that the CTMD in sunflower (Helianthus annuus) may be controlled by a pair of recessive genes [14]. Therefore, the genetic laws of the ray floret type in Asteraceae vary substantially due to the complexity of the genetic background and species differences.
In addition, studies on various ray floret types in Asteraceae have focused on the genetic regulation of petal symmetry, with the CYCLOIDEA (CYC) genes of the TCP family being research hotspots. Different ray floret types were closely related to the expression modifications of CYC2s. Studies on Senecio vulgaris [15], H. annuus [14,16,17], and C. lavandulifolium [18][19][20] have found that expression level changes or mutations of CYC2s could affect the ray floret types [16,21,22]. However, the function of CYC2 genes in the ray florets of various Asteraceae plants was different, which could not completely explain the formation of different ray floret types.
According to numerous studies, ray floret morphology and flower patterns were influenced by plant hormones, such as auxin, ethylene, cytokinin (CTK), gibberellin (GA), abscisic acid (ABA), jasmonic acid (JA), and brassinosteroid (BR). Endogenous auxin derived the successive and centripetal initiation of floret primordia in an approximately circular pattern with a Fibonacci number, which influenced the number of florets [23]. Exogenous application of auxin led to homeotic conversions of florets and phyllaries [24]. The BR-related transcription factor BRI1-EMS-SUPPRESSOR 1 (BES1) had a repressive effect on organ boundary identity genes. In chrysanthemums, overexpressing CmBES1 resulted in an increased fusion of the peripheral ray florets [25]. Additionally, JASMONATE ZIM DOMAIN (CmJAZ1), a JAZ repressor, could repress petal cell expansion and decrease the size of ray florets [26]. The above studies showed that plant hormones significantly regulate the morphology of ray florets.
In order to explore the key genes regulating the morphological differences of ray florets, C. vestitum, the closest hexaploid plant relative to chrysanthemum in this genus [27][28][29], was used for the transcriptome analysis. C. vestitum can be used as a model material to explore the formation mechanism of various ornamental characters of chrysanthemum [30,31], whose ray florets could be divided into the same basic types according to CTMD as chrysanthemum. In the previous study, stable C. vestitum lines of CVW with all ray florets being flat type and CVT with all ray florets being tubular type were obtained. Phenotypic observation revealed that stages 5-6 were the ray floret primordia formation period, stages 7-8 were the petal primordia development period, and stages 9-10 were the key stages for the formation of the difference between flat and tubular ray florets [30]. In order to explore the key genes regulating the different types of ray floret at early developmental stages, the stage-specific materials CVW and CVT were used for transcriptome sequencing. This study not only extends the knowledge of the regulation mechanisms of different ray floret types but also lays the theoretical foundation for directional breeding of flower types in chrysanthemums.

Plant Materials and Growth Conditions
The C. vestitum lines with different ray floret types distributed in the central China and were collected non-destructively through cuttings for ex-situ conservation [30]. The flat-type line CVW (CTMD = 0, Figure 1A) and the tubular-type line CVT (CTMD = 1, Figure 1B) of C. vestitum were grown in an artificial climatic chamber. During the vegetative growth period, plant materials were in long-day (LD) conditions (16 h light/8 h dark). And after completion of this period, plant materials were transferred to short-day (SD) conditions (12 h light/12 h dark). The temperature was maintained at 22 ± 1 • C, the air humidity was controlled at 40-50%, and the light intensity was 100-110 µmol·m −2 ·s −1 .

Plant Materials and Growth Conditions
The C. vestitum lines with different ray floret types distributed in the central Chin and were collected non-destructively through cuttings for ex-situ conservation [30]. Th flat-type line CVW (CTMD = 0, Figure 1A) and the tubular-type line CVT (CTMD = 1 Figure 1B) of C. vestitum were grown in an artificial climatic chamber. During the vegeta tive growth period, plant materials were in long-day (LD) conditions (16 h light/8 h dark And after completion of this period, plant materials were transferred to short-day (SD conditions (12 h light/12 h dark). The temperature was maintained at 22 ± 1 °C, the ai humidity was controlled at 40-50%, and the light intensity was 100-110 μmol·m −2 ·s −1 .

RNA-Seq, Functional Annotation, and Data Processing
The materials were taken from the apical flower buds of CVW and CVT at specifi stages, including stages 5-6 when the ray floret primordia were initiating and developing stages 7-8 when the petal primordia of ray florets began to develop, and stages 9-10 whe the morphological difference of ray florets was formed ( Figure S1) [30]. A total of 18 l braries (W5-6, W7-8, W9-10, T5-6, T7-8, T9-10) were constructed for RNA-seq. In brie total RNAs were extracted from the samples and cDNA libraries were sequenced on th Illumina NovaSeq 6000 sequencing platform (Illumina, San Diego, CA, USA) by B omarker Technologies Corporation (Beijing, China). After connectors of the raw reads an low-quality sequences were removed, clean reads were obtained and then assembled us ing Trinity (version: v2.5.1, major parameter: --min_contig_length 200 --group_pairs_dis tance 500 --min_kmer_cov 1) [32]. The unigene sequences were aligned to the NR, Swiss Prot, COG, KOG, eggNOG4.5, and KEGG databases using DIAMOND software (version v2.0.4) [33] to obtain annotations. Raw sequence data were submitted to the National Cen ter for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database unde accession number PRJNA934692.

RNA-Seq, Functional Annotation, and Data Processing
The materials were taken from the apical flower buds of CVW and CVT at specific stages, including stages 5-6 when the ray floret primordia were initiating and developing, stages 7-8 when the petal primordia of ray florets began to develop, and stages 9-10 when the morphological difference of ray florets was formed ( Figure S1) [30]. A total of 18 libraries (W5-6, W7-8, W9-10, T5-6, T7-8, T9-10) were constructed for RNAseq. In brief, total RNAs were extracted from the samples and cDNA libraries were sequenced on the Illumina NovaSeq 6000 sequencing platform (Illumina, San Diego, CA, USA) by Biomarker Technologies Corporation (Beijing, China). After connectors of the raw reads and low-quality sequences were removed, clean reads were obtained and then assembled using Trinity (version: v2.5.1, major parameter: -min_contig_length 200 -group_pairs_distance 500 -min_kmer_cov 1) [32]. The unigene sequences were aligned to the NR, Swiss-Prot, COG, KOG, eggNOG4.5, and KEGG databases using DIAMOND software (version: v2.0.4) [33] to obtain annotations. Raw sequence data were submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under accession number PRJNA934692.

Weighted Gene Co-Expression Correlation Network Analysis (WGCNA)
In order to screen genes involved in different ray floret types of C. vestitum, WGCNA analysis was performed with the procedure described by Lu et al. [38]. Subsequently, the modules with the highest correlation with the key stages of CVW and CVT were identified for further analysis. Genes in the module were subjected to KEGG analysis, and a hub gene co-expression network was constructed using Cytoscape software (version: v3.5.1) [39].

Real-Time Quantitative Polymerase Chain Reaction (qRT-PCR) Analysis
The expression patterns of DEGs were verified by the qRT-PCR. The differences in gene expression were verified in CVW and CVT capitula at different developmental stages (stages 5-6, 7-8, and 9-10). According to the SYBR Premix Ex Taq kit (Takara, Kyoto, Japan), qRT-PCR analysis was performed on a CFX96™ real-time system (Bio-Rad Laboratories, Hercules, USA). Each qRT-PCR data point was derived from three technical replicates. The specific primer sequences for qRT-PCR were listed in Table S1. CvSAND was used as a reference gene to normalize the qRT-PCR data [40]. Relative expression levels were calculated using the 2 −∆∆CT method [41], and the data are presented as the mean ± SD.

Sequencing, Assembly, and Functional Annotation of the Transcriptome
To investigate the molecular mechanism of different ray floret types formation, stagespecific transcriptomes of CVW and CVT were carried out. The capitula of CVW and CVT at stages 5-6, 7-8, and 9-10 provided the templates for RNA-Seq analysis, and each group of samples contained three biological replicates. Summary statistics for 18 samples of sequencing data evaluations are shown in Table S2. A total of 125.09 Gb of clean data with a Q30 (percentage of sequences with sequencing error rates lower than 0.1%) of 94.88-95.78% was generated from 18 samples. Following assembly, 71,980 unigenes were obtained, of which 45,202 were above 1 kb in length, and the value of N50 was 2128 ( Figure S2). A total of 47,348 unigenes were annotated, of which 23,980 (50.65%) were annotated in KOG, 31,299 (66.10%) were annotated in Pfam, 29,613 (62.54%) were annotated in Swissprot, and 46,385 (97.97%) were annotated in the Nr database ( Figure 2A).
Functional classification of unigenes was performed using GO, COG, and KEGG assignments. 38,338 (80.97%) unigenes were divided into three functional groups (cell component, molecular function, and biological process) with a total of 44 categories by GO annotation ( Figure 2B). The cellular component functional group had the most proteins associated with cellular anatomical entities and intracellular. A larger number of unigenes related to binding and catalytic activity were annotated in the molecular function group. The majority of proteins were related to cellular processes, metabolic processes, and biological regulation in the biological processes functional group. In addition, 12,071 (25.49%) unigenes were annotated in COG and classified into 25 functional groups ( Figure 2C), among which the number of unigenes associated with translation, ribosomal structure and biogenesis, and signal transduction mechanisms were higher than those in other functional groups. A total of 29,780 (41.37%) unigenes were mapped into 136 KEGG pathways, with the most represented pathway being "plant-pathogen interaction (ko04626)", followed by "plant hormone signal transduction (ko04075)" ( Figure 2D).

Screening for the Genes Related to the Development of Different Ray Floret Types through WGCNA
To understand the overall transcriptional profile during the development of ray florets with different CTMDs, WGCNA analysis was performed on the genes obtained by sequencing. According to the expression trends, the obtained genes were divided into 15 modules ( Figure 3A), and the number of genes contained in each module is shown in Table S3. Eigengenes (the first major component in the module) represented the gene expression profile of the whole module, and 14 different expression patterns were obtained except for the gray module, the invalid module, with genes not involved in clustering ( Figure S3). In previous studies, it was found that stages 9-10 are key stages for the formation of differences between flat and tubular ray florets. The results of the present study showed that the modules with the highest correlation with stages 9-10 of CVW and CVT were lavenderblush2 (correlation coefficient: 0.99) and tan (correlation coefficient: 0.98) ( Figure 3B). It was speculated that the genes in these two modules may be closely related to the formation of ray floret morphological differences.
Further analysis of the lavenderblush2 module revealed that this module contained 2,281 genes and was abundantly expressed in stages 9-10 of the CVW ( Figure 4A). A KEGG function analysis was performed on these genes. It was found that the most genes were enriched in the plant hormone signal transduction pathway (ko04075) ( Figure 4B). The genes with KME values in the top 50 were selected to construct a co-expression network, and the hub genes in this network were screened according to their connectivity ( Figure 4C). Furthermore, some transcription factors belonging to nine gene families were found, including basic helix-loop-helix (bHLH), MYBs, APETALA2/ETHYLENE

Screening for the Genes Related to the Development of Different Ray Floret Types through WGCNA
To understand the overall transcriptional profile during the development of ray florets with different CTMDs, WGCNA analysis was performed on the genes obtained by sequencing. According to the expression trends, the obtained genes were divided into 15 modules ( Figure 3A), and the number of genes contained in each module is shown in Table S3. Eigengenes (the first major component in the module) represented the gene expression profile of the whole module, and 14 different expression patterns were obtained except for the gray module, the invalid module, with genes not involved in clustering ( Figure S3). In previous studies, it was found that stages 9-10 are key stages for the formation of differences between flat and tubular ray florets. The results of the present study showed that the modules with the highest correlation with stages 9-10 of CVW and CVT were lavenderblush2 (correlation coefficient: 0.99) and tan (correlation coefficient: 0.98) ( Figure 3B). It was speculated that the genes in these two modules may be closely related to the formation of ray floret morphological differences.
Further analysis of the lavenderblush2 module revealed that this module contained 2,281 genes and was abundantly expressed in stages 9-10 of the CVW ( Figure 4A). A KEGG function analysis was performed on these genes. It was found that the most genes were enriched in the plant hormone signal transduction pathway (ko04075) ( Figure 4B). The genes with KME values in the top 50 were selected to construct a co-expression network, and the hub genes in this network were screened according to their connectivity ( Figure 4C). Furthermore, some transcription factors belonging to nine gene families were found, in-       The tan module contained 1602 genes, most of which were significantly expressed at stages 9-10 of the CVT ( Figure 5A). KEGG analysis showed that these genes were enriched in the highest number of the protein processing pathways in the endoplasmic reticulum (ko04141) ( Figure 5B). A co-expression network was constructed using the genes with KME values in the top 50, and the top 20 Hub genes with high connectivity to other genes were obtained ( Figure 5C). Meanwhile, the statistics of transcription factors revealed that the tan module contained 5 bHLH, 5 MYB, 3 TCP, 2 AP2/ERF, 2 WRKY, 2 ARF, 2 AUXIN/INDOLE-3-ACETIC ACID (AUX/IAA), 1 MADS-box, and 1 NAC ( Figure 5D). The tan module contained 1602 genes, most of which were significantly expressed stages 9-10 of the CVT ( Figure 5A). KEGG analysis showed that these genes were enrich in the highest number of the protein processing pathways in the endoplasmic reticulu (ko04141) ( Figure 5B). A co-expression network was constructed using the genes wi KME values in the top 50, and the top 20 Hub genes with high connectivity to other gen were obtained ( Figure 5C). Meanwhile, the statistics of transcription factors revealed th the tan module contained 5 bHLH, 5 MYB, 3 TCP, 2 AP2/ERF, 2 WRKY, 2 ARF, AUXIN/INDOLE-3-ACETIC ACID (AUX/IAA), 1 MADS-box, and 1 NAC ( Figure 5D).

DEGs Identified by K-Means Cluster Analysis
The overall expression pattern of DEGs was shown on the clustering map with means cluster analysis, and all DEGs were classified into 11 clusters ( Figure 6A). T genes in cluster 1 were expressed at high levels at stages 5-6 and 7-8 of CVT samples. T DEGs in clusters 2, 4, 8, and 11 showed high expression in W5-6, T5-6, T9-10, and W9samples, respectively. The expression pattern of cluster 5 showed that there was no si nificant change in transcription level at different stages of CVT, while the expression lev increased gradually in CVW with the continuous development of capitula. In additio CVW and CVT had similar expression patterns in cluster 9, while the expression levels CVT at different stages were significantly higher than those of CVW. Since the key stag

DEGs Identified by K-Means Cluster Analysis
The overall expression pattern of DEGs was shown on the clustering map with Kmeans cluster analysis, and all DEGs were classified into 11 clusters ( Figure 6A). The genes in cluster 1 were expressed at high levels at stages 5-6 and 7-8 of CVT samples. The DEGs in clusters 2, 4, 8, and 11 showed high expression in W5-6, T5-6, T9-10, and W9-10 samples, respectively. The expression pattern of cluster 5 showed that there was no significant change in transcription level at different stages of CVT, while the expression level increased gradually in CVW with the continuous development of capitula. In addition, CVW and CVT had similar expression patterns in cluster 9, while the expression levels of CVT at different stages were significantly higher than those of CVW. Since the key stages of morphological differences in ray florets are stages 9-10 and genes involved in regulating ray floret development may function at earlier stages, it was suggested that the key candidate genes may be concentrated in clusters 5, 8, 9, and 11. Through K-means cluster analysis, 66 DEGs encoding transcription factors, including ARF, AP2/ERF, MADS-box, NAC, CYC2, bHLH, etc., were detected ( Figure 6B). of morphological differences in ray florets are stages 9-10 and genes involved in regulating ray floret development may function at earlier stages, it was suggested that the key candidate genes may be concentrated in clusters 5, 8, 9, and 11. Through K-means cluster analysis, 66 DEGs encoding transcription factors, including ARF, AP2/ERF, MADS-box, NAC, CYC2, bHLH, etc., were detected ( Figure 6B).

Screening for Candidate DEGs Regulating Morphological Differences in Ray Florets by Pairwise Comparison
To identify more important genes regulating CTMD differences in ray florets, a pairwise comparison was performed between CVW and CVT. It was shown that there were a large number of DEGs between CVW and CVT at the same developmental stage in the Venn diagram ( Figure 7A)

Screening for Candidate DEGs Regulating Morphological Differences in Ray Florets by Pairwise Comparison
To identify more important genes regulating CTMD differences in ray florets, a pairwise comparison was performed between CVW and CVT. It was shown that there were a large number of DEGs between CVW and CVT at the same developmental stage in the Venn diagram ( Figure 7A). A total of 25,141 DEGs were detected in the three comparisons, with 13,899 (6768 up-regulated and 7131 down-regulated), 16,763 (8312 up-regulated and 8451 down-regulated) and 12,266 (4648 up-regulated and 7618 down-regulated) DEGs in W5-6_vs._T5-6, W7-8_vs._T7-8, and W9-10_vs._T9-10, respectively ( Figure 7B). Based on the developmental characteristics of different ray floret types, three gene sets are considered to be very important, including (i) genes differentially expressed at stages 5-6, 7-8, and 9-10; (ii) genes differentially expressed at stages 7-8 and 9-10; and (iii) genes differentially expressed at stages 9-10 only. The three gene sets contained 478 transcription factors, among which were 25 AP2/ERF family genes, 15 MADS-box family genes, 7 NACs, 6 ARFs, 5 AUX/IAA family genes, and 4 TCPs possibly involved in the regulation of morphological differences ( Figure 7C). Furthermore, DEGs in the plant hormone signal transduction pathways of the three sets were screened, and STRING 11.5 (https://cn.string-db.org/, 6 June 2022) and Cytoscape software (version: v3.5.1) were used for interaction network analysis. It was found that DEGs were mainly auxin-related genes, including TRANSPORT INHIBITOR RESPONSE (TIR1, Unigene_033158), IAA14 (Unigene_009970), IAA3/SHORT HYPOCOTYL 2 (IAA3/SHY2, Unigene_085389), MONOPTEROS/ARF5 (MP/ARF5, Uni-gene_174210), ARF19 (Unigene_093663), IAA18 (Unigene_034051), ARF9 (Unigene_098607), and ETTIN/ARF3 (ETT/ARF3, Unigene_237143) ( Figure 7D). It was thus hypothesized that auxin may be the key factor affecting the morphology of ray florets. s 2023, 14, x FOR PEER REVIEW 9 of the developmental characteristics of different ray floret types, three gene sets are cons ered to be very important, including (ⅰ) genes differentially expressed at stages 5-6, 7 and 9-10; (ⅱ) genes differentially expressed at stages 7-8 and 9-10; and (ⅲ) genes diff entially expressed at stages 9-10 only. The three gene sets contained 478 transcription f tors, among which were 25 AP2/ERF family genes, 15 MADS-box family genes, 7 NACs ARFs, 5 AUX/IAA family genes, and 4 TCPs possibly involved in the regulation of m phological differences ( Figure 7C). Furthermore, DEGs in the plant hormone signal tran duction pathways of the three sets were screened, and STRING 11.5 (https://cn.strin db.org/, 6 June 2022) and Cytoscape software (version: v3.5.1) were used for interacti network analysis. It was found that DEGs were mainly auxin-related genes, includi TRANSPORT INHIBITOR RESPONSE (TIR1, Unigene_033158), IAA14 (Unigene_00997 IAA3/SHORT HYPOCOTYL 2 (IAA3/SHY2, Unigene_085389), MONOPTEROS/AR (MP/ARF5, Unigene_174210), ARF19 (Unigene_093663), IAA18 (Unigene_034051), AR (Unigene_098607), and ETTIN/ARF3 (ETT/ARF3, Unigene_237143) ( Figure 7D). It w thus hypothesized that auxin may be the key factor affecting the morphology of ray f rets. Based on the above transcriptome analysis of CVW and CVT, some candidate DE were obtained in this study. The expression levels and annotations of these genes in d ferent sequencing samples are shown in Table 1. Auxin related genes ARF and AUX/IA ethylene-related genes AP2/ERF, MADS-box family genes regulating flower develo Based on the above transcriptome analysis of CVW and CVT, some candidate DEGs were obtained in this study. The expression levels and annotations of these genes in different sequencing samples are shown in Table 1. Auxin related genes ARF and AUX/IAA, ethylene-related genes AP2/ERF, MADS-box family genes regulating flower development, TCP family genes regulating flower symmetry, and NAC family genes closely related to organ boundaries may be involved in the regulatory network of ray floret CTMD difference formation.

Validation of Key DEGs Related to Morphological Difference in Ray Florets
Based on the above results of the transcriptome analysis, it is hypothesized that the genes involved in plant hormone signal transduction, especially auxin-related genes, play an important role in regulating the development of different ray floret types. To determine the expression patterns of candidate genes in CVW and CVT, qRT-PCR analysis was performed in capitula at different developmental stages.
ARFs, as the core transcription factor connecting auxin signals and downstream genes, play a crucial role in the regulation of plant growth and development. Analysis of the expression pattern of CvARFs revealed that CvARF1, CvARF6, CvARF8, and CvARF9 were all expressed at significantly higher levels in CVT at different developmental stages than CVW. And transcript levels of CvARF2 and CvARF5 were higher in CVT than CVW at stages 7-8, when petal primordia of ray florets were flourishing. However, the transcription level of CvARF3 was higher in CVW, and this trend was also consistent with the results of a previous transcriptome ( Figure 8A). CvARF5 and CvARF6, which belong to class II of the ARF family, usually activate downstream gene transcription in response to auxin, and CvARF3, which belongs to class I, usually functions as a transcriptional repressor. Their differential expression may be the cause of different ray floret types formation.
pressed in CVT at stages 5-6, which was the ray floret primordia formation period, and these genes were generally expressed at higher levels in CVT than in CVW ( Figure 8B).
Transcription factors closely related to petal development were also screened from the transcriptome, including PISTILLATA (CvPI), APETALA3 (CvAP3), DEFH21(Cvdefh21), SEPALLATA 3 (CvSEP3), and CvCDM77 from the MADS-box family; AINTEGUMENTA (CvANT); AINTEGUMENTA-LIKE 6 (CvAIL6); and CvWRI1 from the AP2/ERF family; and NAC-activated by AP3/PI (CvNAP) from the NAC family. CvPI was expressed at the highest level in CVT at stages 9-10, while CvAP3 was at the highest level in CVT at stages 7-8. The expression of CvSEP3 and CvCDM77 in CVT at stages 7-8 and 9-10 was higher than in CVW. And CvWRI1 was expressed at high levels in CVT at all developmental stages compared with CVW ( Figure 8C).

Auxin and Auxin-Related Genes Are Involved in Regulating the Formation of Different Ray Floret Types
In this study, transcriptome sequencing and gene expression analysis revealed that a large number of genes related to auxin were differentially expressed in ray florets with different CTMDs. Additionally, we discovered that CvSAUR, CvGH3, and CvGH3.5 were differentially expressed among the flat, spoon, and tubular of the C. vestitum line CVZ [30]. All these results demonstrated that auxin and auxin-related genes were involved in regulating different ray floret types. Auxin plays a crucial role in petal formation and is essential for plant growth and development [42,43]. Auxin activity is necessary for the beginning of petal primordia [44,45]. The number and shape of single floral petals can be dramatically impacted by mutations in the auxin synthesis, polar transport, and response genes [46,47]. The mutation of YUCs related to auxin synthesis can lead to serious defects in flower development [48]. And disruption of auxin polar transport in pin1 and pinoid (pid) mutants can lead to anomalies in the number and positioning of organs or partial failure of floral organ initiation [49,50]. In addition, ARFs are the core transcription factors of auxin signal transduction and also affect the development process of petals. In Mimulus lewisii, the up-regulation of MlARF3 and MlARF4 resulted in unfused corollas [51]. Meanwhile, the study on the Asteraceae plant Matricaria recutita found that auxin affected the initiation of ray florets [24]. And it was found that CvARF3, CvARF5, CvARF6, and genes involved in auxin synthesis and transport were differently expressed between flat and tubular ray florets in this study. Therefore, we speculated that auxin may play an important role in ray floret morphogenesis, and the auxin-related DEGs screened in this Moreover, the expression patterns of genes related to auxin synthesis, transport, and response were analyzed in this study. CvYUCCA10 was expressed at the highest level in CVT at stages 7-8 and decreased at stages 9-10, similar to GRETCHEN HAGEN 3.5 (CvGH3.5), CvIAA14, and PINFORMED 1c (CvPIN1c). CvGH3, SMALL AUXIN UPREGU-LATED RNA 71 (CvSAUR71), CvAUX1, and LAX PANICLE (CvLAX) were highly expressed in CVT at stages 5-6, which was the ray floret primordia formation period, and these genes were generally expressed at higher levels in CVT than in CVW ( Figure 8B).
Transcription factors closely related to petal development were also screened from the transcriptome, including PISTILLATA (CvPI), APETALA3 (CvAP3), DEFH21 (Cvdefh21), SEPALLATA 3 (CvSEP3), and CvCDM77 from the MADS-box family; AINTEGUMENTA (CvANT); AINTEGUMENTA-LIKE 6 (CvAIL6); and CvWRI1 from the AP2/ERF family; and NAC-activated by AP3/PI (CvNAP) from the NAC family. CvPI was expressed at the highest level in CVT at stages 9-10, while CvAP3 was at the highest level in CVT at stages 7-8. The expression of CvSEP3 and CvCDM77 in CVT at stages 7-8 and 9-10 was higher than in CVW. And CvWRI1 was expressed at high levels in CVT at all developmental stages compared with CVW ( Figure 8C).

Auxin and Auxin-Related Genes Are Involved in Regulating the Formation of Different Ray Floret Types
In this study, transcriptome sequencing and gene expression analysis revealed that a large number of genes related to auxin were differentially expressed in ray florets with different CTMDs. Additionally, we discovered that CvSAUR, CvGH3, and CvGH3.5 were differentially expressed among the flat, spoon, and tubular of the C. vestitum line CVZ [30]. All these results demonstrated that auxin and auxin-related genes were involved in regulating different ray floret types. Auxin plays a crucial role in petal formation and is essential for plant growth and development [42,43]. Auxin activity is necessary for the beginning of petal primordia [44,45]. The number and shape of single floral petals can be dramatically impacted by mutations in the auxin synthesis, polar transport, and response genes [46,47]. The mutation of YUCs related to auxin synthesis can lead to serious defects in flower development [48]. And disruption of auxin polar transport in pin1 and pinoid (pid) mutants can lead to anomalies in the number and positioning of organs or partial failure of floral organ initiation [49,50]. In addition, ARFs are the core transcription factors of auxin signal transduction and also affect the development process of petals. In Mimulus lewisii, the up-regulation of MlARF3 and MlARF4 resulted in unfused corollas [51]. Meanwhile, the study on the Asteraceae plant Matricaria recutita found that auxin affected the initiation of ray florets [24]. And it was found that CvARF3, CvARF5, CvARF6, and genes involved in auxin synthesis and transport were differently expressed between flat and tubular ray florets in this study. Therefore, we speculated that auxin may play an important role in ray floret morphogenesis, and the auxin-related DEGs screened in this study, such as CvARF3, CvARF5, and CvARF6, may be the key genes regulating morphological differences of ray florets.

Multiple Plant Hormones Are Involved in Regulating the Morphology of Ray Florets
In addition to auxin, several plant hormones are usually involved in regulating petal development. It has been reported that the AUXIN-REGULATED GENE INVOLVED IN ORGAN GROWTH (ARGOS) acts upstream of ANT [52,53], which is an AP2/ERF transcription factor member of a subfamily that is associated with regulation of petal cell proliferation [54,55]. In addition, AIL6 is activated by MP/ARF5 [56], which has a redundant function with ANT and regulates petal development [57]. In this study, the qRT-PCR analysis revealed that CvARF5 was expressed at a higher level in CVT than in CVW, but CvANT and CvAIL6 did not show any significant differences between the two types. While CvWRI1, which belongs to the same gene subfamily as CvANT and CvAIL6, has a higher transcript level in CVT than in CVW, implying that it may be involved in regulating the formation of different ray floret types.
Studies on Asteraceae plants have also found that several plant hormones, including ethylene, CTK, GA, ABA, BR, and JA, also affect the ray floret types [58][59][60][61]. GhWIP2 in Gerbera acts as a transcriptional repressor, suppressing cell expansion and reducing the length of ray florets by modulating crosstalk between GA, ABA, and auxin [60]. In addition, CmJAZ1, a repressor of the JA signaling pathway, and CmBES1, the BR-related transcription factor, were involved in regulating the morphology of ray florets [25,26]. It is evident that numerous plant hormone signals and related genes are involved in the regulation of ray floret morphogenesis. The auxin-related genes, such as CvARF3, CvARF5, and CvARF6, and the AP2/ERF family gene CvWRI1 identified in this study, may mediate auxin and ethylene signals to jointly affect the morphological differences in ray florets.

MADS-Box Genes Affect the Morphogenesis of Ray Florets
MADS-box genes are widely involved in the regulatory network of floral organ development, in which class A, B, and E genes synergistically regulate petal development [62,63]. In this study, the expression patterns of B-class genes CvAP3, CvPI, and Cvdefh21 showed significant tissue specificity, with expression levels higher in CVT than in CVW, while transcript levels of E-class genes CvSEP3 and CvCDM77 were also significantly higher in CVT than in CVW. MADS-box genes control capitulum development in Asteraceae plants [64][65][66]. SEPALLATA-like genes GRCDs have been shown to contribute to meristem determinacy as well as flower type differentiation [67,68], and they also have functional redundancy in the regulation of organ identity [69]. Suppression of GGLO1, GDEF1, and GDEF2 expression resulted in the disappearance or degeneration of gerbera trans florets and outer whorls of ray florets [70]. The roles of some MADS-box genes in Chrysanthemum have been analyzed in previous studies [71,72]. Transcriptome and expression analysis revealed that some MADS-box genes influence the development of the ray floret and disc floret [73,74], and MADS-box genes involved in networks interact with CYC2 genes involved in networks to synergistically regulate floret differentiation [75,76]. The above studies showed that MADS-box genes were found to be extensively involved in the regulatory network of capitulum development. The B-class and E-class genes were differentially expressed in flat and tubular types in this study and could play important roles in regulating the formation of different ray floret types.

Conclusions
Based on floral development stage-specific transcriptome analysis combined with gene expression analysis, it was found that auxin-related genes, such as CvARF3, CvARF5, and CvARF6, and some transcription factors involved in floral development, such as AP2/ERF, MADS-box, CYC2, NAC, and bHLH, were differentially expressed in CVW and CVT. The results suggest that auxin plays a key role in ray floret development. Furthermore, CvARFs and homeotic genes involved in petal development work together to regulate ray floret differentiation. Future work will focus on the function and regulatory relationships of these genes. In conclusion, this study raises the prospects of auxin in regulating the morphological differences of ray florets and lays the foundation for further study of molecular mechanisms.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes14030766/s1. Figure S1: Schematic diagram of the capitulum development process of Chrysanthemum vestitum (Pu et al., 2020). The diagram was drawn according to the observation results of capitula obtained with paraffin sections and a scanning electron microscope. The ray floret primordia are initiating and developing at stages 5-6, the petal primordia of ray florets are developing at stages 7-8, and the morphological difference between a flat ray floret and a tubular ray floret has formed at stages 9-10. Figure S2: The length distribution of unigenes. The horizontal coordinate indicates the different length intervals of unigenes; the vertical coordinate indicates the number of unigenes in a certain length interval. Figure S3: Expression patterns of eigengenes in different modules. Expression patterns of 12 modules are shown in this figure, with lavenderblush2 and tan being shown in Figures 4A and 5A, respectively. Table S1: Primer sequences used in qRT-PCR experiments. Table S2: Summary statistics of clean reads in the transcriptome. Table S3

Data Availability Statement:
The datasets supporting the results presented in this study are included in this article and its additional files.