Characterization of WRKY Gene Family in Whole-Genome and Exploration of Flowering Improvement Genes in Chrysanthemum lavandulifolium

Chrysanthemum is a well-known ornamental plant with numerous uses. WRKY is a large family of transcription factors known for a variety of functions ranging from stress resistance to plant growth and development. Due to the limited research on the WRKY family in chrysanthemums, we examined them for the first time in Chrysanthemum lavandulifolium. A total of 138 ClWRKY genes were identified, which were classified into three groups. Group III in C. lavandulifolium contains 53 members, which is larger than group III of Arabidopsis. The number of introns varied from one to nine in the ClWRKY gene family. The “WRKYGQK” motif is conserved in 118 members, while other members showed slight variations. AuR and GRE responsive cis-acting elements were located in the promoter region of WRKY members, which are important for plant development and flowering induction. In addition, the W box was present in most genes; the recognition site for the WRKY gene may play a role in autoregulation and cross-regulation. The expression of the most variable 19 genes in terms of different parameters was observed at different stages. Among them, 10 genes were selected due to the presence of CpG islands, while nine genes were selected based on their close association with important Arabidopsis genes related to floral traits. ClWRKY36 and ClWRKY45 exhibit differential expression at flowering stages in the capitulum, while methylation is detected in three genes, including ClWRKY31, ClWRKY100, and ClWRKY129. Our results provide a basis for further exploration of WRKY members to find their functions in plant growth and development, especially in flowering traits.


INTRODUCTION
Chrysanthemum lavandulifolium is one of the original species of Chrysanthemum × morifolium and uses as a model plant. It is a major ornamental plant, originated in China, an important cut flower with diversified petal color and shape having a variety of flower structures, which makes it a highly valuable ornamental plant in the floral industry worldwide (Nguyen and Lim, 2019). Chrysanthemum is among the four most important ornamental plants in the world used for beautification and aesthetic value (Shinoyama et al., 2012). The flower of chrysanthemum is also known for its health benefits, and it contains important antioxidants and phenolic compounds (Lin and Harnly, 2010;Yue et al., 2018).
Due to a simple diploid (2n = 2x = 18) genetic background, C. lavandulifolium dominates over C. morifolium (Huang et al., 2012;Dong-ru et al., 2019). The wild ancestors having desired qualities against biotic and abiotic cues can be used for the genetic improvement of cultivated varieties. Identification and incorporation of such genes into commercial varieties reduces pesticide use, improves plant health, and induces photoperiodcontrolling flowering mechanisms (Wu et al., 2010;Fu et al., 2013;Yang et al., 2017).
Genetic architectural study is one of the major research directions in any plant. By narrowing the path further, transcription activation and silencing is the major field of concern. Many transcription factors (TFs) have been explored to switch on and off against certain environmental conditions. WRKY TF is one of the major families first identified in plants in 1994 by the name of DNA-binding protein in sweet potato, which is labeled as SWEET POTATO FACTOR1 (SPF1) (Ishiguro and Nakamura, 1994). WRKY TFs play important roles in resistance against biotic (fungal or bacterial pathogen) and abiotic stresses, growth, and development. Several WRKY TFs have been proven to provide resistance against biotic stresses exerted by fungal or bacterial pathogens by influencing other associated genes (Turck et al., 2004;Duan et al., 2007;Zhou et al., 2008Zhou et al., , 2011Li et al., 2009Li et al., , 2011Chen et al., 2010;Scarpeci et al., 2013;Phukan et al., 2016;Yokotani et al., 2018;Warmerdam et al., 2020). The name WRKY has been given to this family for their distinguishing characteristic, the WRKY domain, which is composed of highly conserved 60 amino acids and a zinc finger motif. WRKY TFs are classified into three groups based on the number of (WRKYGQK) domain, and zinc-finger motif. Group I has two domains, and groups II and III have one domain. Groups I and II have C2H2 type zinc finger motifs while group III has a C 2 HC zinc finger motif. Members of group II were further classified into subgroups (Rushton et al., 2010).
WRKY TFs are extensively studied in Arabidopsis and most other plants. WRKY7 and OsWRKY11 expression resulted in blooming delay in Arabidopsis and rice (Cai et al., 2014;Chen et al., 2019), WRKY13 has an adverse effect on flowering time, and WRKY12, GsWRKY20 and CsWRKY50 promote flowering, possibly through the influence of gibberellic acid under short-day conditions in Arabidopsis (Luo et al., 2013;Kumar et al., 2016;Li et al., 2016), GhWRKY22 regulating pollen development (Wang et al., 2019;Ramos et al., 2021). While different WRKY family TFs play an important role in chrysanthemum, they showed resistance to salt stress (Liu et al., , 2014Li et al., 2015b;Jaffar et al., 2016;Liang et al., 2017;Wang K. et al., 2017;He et al., 2018;Zhang et al., 2020), and CmWRKY1 and CmWRKY15 transgenic lines in C. morifolium interact with SA and ABA and altered plant development (Fan et al., 2015(Fan et al., , 2016Bi et al., 2021). CmWRKY48 showed resistance against aphid infestation (Li et al., 2015a). They also play a regulatory role in the flowering stages. The complex network of VQ20, WRKY2, and WRKY34 interaction influences the plant gametogenesis and plays a crucial part in the development of pollen and pollen tube growth, modulating flowering (Lei et al., 2017;Ma et al., 2020).
Flowering time is a complex character regulated by a network of different external and internal signals. Different pathways have been identified, including photoperiod pathway, vernalization pathway, gibberellin-regulated pathways, and autonomous floral initiation pathway (Soares et al., 2020). Therefore, flowering induction studies are more crucial in ornamental plants, especially in chrysanthemum, where all the ornamental characteristics are associated with flowering, i.e., flower color, shape, and size. Chrysanthemum is a short-day plant and needs a specific photoperiod to induce flowering, which makes it relatively expensive as well as limits the annual production of chrysanthemum flowers. We explored the WRKY gene family in C. lavandulifolium for the first time, prior, the WRKY gene family has only been studied in the C. morifolium . The size of the WRKY family in C. morifolium comprised 15 members, which were named from CmWRKY1 to CmWRKY15 (GenBank: KC615355-KC615369). As WRKY is a large family, influencing floral characteristics, we need to explore them deeply. Our present study was a part of such investigation in which we tried to explore the WRKY gene family in C. lavandulifolium comprised structural gene analysis, protein motifs, phylogenetic relatedness of C. lavandulifolium and Arabidopsis thaliana, cisacting elements, recognition of CpG rich regions, the methylation status of the genes, and expression of 19 important genes in different tissues at different stages to check their role in flowering-related traits. This investigation will lead us to identify variation in transcription regions in the WRKY family for future chrysanthemum research.

Plant Material and Treatment
The plant materials of C. lavandulifolium were obtained from Beijing Forestry University. The plants' G1 lines were grown under controlled environmental conditions in the growth chamber. All the plants used in this research came from the G1 line. The plants were propagated by cuttings and grown in pots. The plants were kept under long-day conditions (14 h L/10 h D) light/dark at a 25 • C temperature. The photoperiod conditions were changed to short-day conditions (12 h L/12 h D) when the plants attained maturity (more than 14 leaves). The temperature was kept constant. The plants were watered every 7 days. Stem, root, leaf, buds, and capitulum samples for quantitative real-time polymerase chain reaction (qRT-PCR) were taken at different growth stages, including young, mature, and flowering. The plants were considered at a young stage when they attained eight leaves. The plants were called into the adult stage when the number of leaves was 14, and they were transferred to shortday conditions. The flower develops in short conditions, and samples were taken at different stages of flower development (Supplementary Figure 1). The little buds (1-3 mm in diameter) appeared after 7 days in short-day conditions and were taken for RNA extraction. The median buds (5-8 mm in diameter) were taken after a few days before flower opening. The first flower samples were taken when the flowers had just opened.

Characterization of WRKY Gene Family in Chrysanthemum lavandulifolium
The C. lavandulifolium genome has been sequenced (Wen et al., 2022). All the protein and genomic sequences of ClWRKY members were obtained from PRJNA681093, National Center for Biotechnology Information (NCBI). The genomic information of C. lavandulifolium was completed by Henan University and Beijing Forestry University, and related papers have been published (Wen et al., 2022). The redundant sequences, which lack the WRKY domain, were removed after the sequence analysis using manual inspection in the Mega X software. A total of 138 definite ClWRKY sequences were taken for further study. The ClWRKY members were named from ClWRKY1 to ClWRKY138, according to the source provider (Supplementary Table 3).

Classification of ClWRKY Proteins
Multiple protein sequence alignments of C. lavandulifolium and A. thaliana WRKY were made using ClustalX and MEGA-X version 10.2.6 (Molecular Evolutionary Genetics Analysis) (Kumar et al., 2018). The phylogenetic tree based on amino acid (aa) sequences of ClWRKY and AtWRKY conserved WRKY domains was constructed following the neighbor-joining method with 1,000 bootstraps to determine the close association of both gene families.

Gene Structure and Protein Motif Analysis
The gene structure (intron-exon) of ClWRKY members, fulllength gene, and coding (CDS) sequences were analyzed in the online database, Gene Structure Display Server (GSDS) 1 as described (Hu et al., 2015). Analysis of the conserved protein motifs was performed using the NCBI CDD batch 2 and the TBtools software .

Cis-Acting Elements in the Promoter Region of ClWRKY Members
The promoter sequences of all 138 ClWRKY members were submitted to the PlantCARE software to analyze their promoter region for cis-acting elements (Lescot et al., 2002). For this purpose, a 2,000 bp sequence upstream of the gene was taken to identify cis-acting elements. The MEME suite software was used to identify de novo motifs and known enriched motifs (Bailey et al., 2015). The de novo motif analysis was performed to identify significant patterns (Bailey and Elkan, 1994). 3 The known enriched motifs were identified via simple enrichment analysis using the MEME software (Bailey and Grant, 2021). 4 Tomtom motif comparison tool was used to quantify the similarity between motifs (Gupta et al., 2007).

CpG Island Analysis of the Promoter Region of ClWRKY Gene Family
For CpG analysis, promoter regions from the WRKY gene family were sequenced. For the exploration of CpG islands, the "Methyl Primer Express version 1.0" software was used (Methyl Primer Express Software v1.0 Quick Reference Card). 5 All the sequences were submitted to the software one by one to identify CpG islands or rich regions. The minimum length of the island was kept at 200 bp, while the maximum length was 2,000 bp. The criterion of "C + G/Total bases" was 50%. The CpG observed/CpG expected ratio was 0.6.

Genomic DNA Extraction and Bisulfite Treatment
The total DNA was extracted from different tissues of C. lavandulifolium. The tissues included root, stem, leaf, and capitulum at three different growth stages, i.e., young stage, adult stage, and flowering stage. Total genomic DNA was extracted using the "Plant Genomic DNA Kit" (TOLOBIO), according to the manufacturer's instructions. DNA was quantified using Thermo Scientific NanoDrop TM 2000/2000c Spectrophotometers, and the quality of DNA was checked by 1.5% agarose gel electrophoresis.
Total genomic DNA from each sample was treated with sodium bisulfites using the "DNA Bisulfite Conversion Kit" (TIANGEN) according to the manufacturer's instructions. The conversion efficiency and quality were measured using Thermo Scientific NanoDrop TM 2000/2000c Spectrophotometers.

Methylation-Specific Polymerase Chain Reaction
Ten genes were selected on the basis of CpG analysis. Primers for all ten genes were designed using the "Methyl primer express version 1.0" software. For this purpose, two types of primers (methylated and unmethylated) were designed to check the possible methylated sites. The primer length ranged from 18 to 22 bp. The PCR amplicon length was between 100 and 175 bp. Methylation-Specific Polymerase Chain Reaction (MSP) was conducted to analyze the CpG island among the ten selected genes after treatment with bisulfite for methylation status. The MSP conditions were as follows: pre-denaturation for 3 min at 95 • C, followed by 34 cycles of the 30s at 95 • C and 30 s at 55 • C, 1 min at 72 • C, and a final extension for 5 min at 72 • C. PCR products were resolved onto 1.5% agarose gel and observed under UV transilluminator.

Quantitative Real-Time Polymerase Chain Reaction
Nineteen (Zhou et al., 2011) variable genes on the basis of CpG island and phylogenetic analysis were selected for gene expression patterns. Primers for all the 19 ClWRKY genes were designed using the "Primer Premier 5" software. RNA was extracted from different plant tissues, i.e., root, stem, and leaves at young, adult, and flowering stages, while the capitulum samples, which comprise little buds, median buds, and first flower (right at the time of flower opening), were taken at flowering stages. The RNA was extracted using the RNAprep Pure Plant Plus Kit (spin column) according to the manufacturer's instructions. The cDNA was synthesized using the "StarScript II First-strand cDNA Synthesis Mix with gDNA Remover." The expression of all the selected genes was observed in qRT-PCR. The qRT-PCR was carried out on light cycle "Quantageneq225 Fluorescent Quantitative PCR, " using "2 × RealStar Green Fast Mixture." The conditions of PCR were as follows: denaturation at 95 • C for 2 min, followed by 40 cycles of denaturation at 95 • C for 15 s, annealing at 60 • C for 20 s, and extension at 65 • C for 5 s. The actin gene of the chrysanthemum was used as an internal control. The data from qRT-PCR for relative expression were analyzed using the 2 − Ct method proposed by Livak and Schmittgen (2001).

Identification and Characterization of WRKY Gene Family in Chrysanthemum lavandulifolium
WRKY genes possess a conserved "WRKYGQK" domain. The "WRKYGQK" heptapeptide is the most defining characteristic of WRKY TFs. The C. lavandulifolium WRKY sequences were obtained from PRJNA681093, NCBI. A total of 153 WRKY gene members have been selected; among them, 138 were exhibited for sequence analysis, and 15 were considered defective. The WRKY members are named from ClWRKY1 to ClWRKY138 according to their location on the chromosome. Protein motif analysis showed homology in most of the conserved regions. Among 138 ClWRKY members, 85.51% of members have homology in the "WRKYGQK" core domain, and 19 members have two conserved WRKY domains. However, 14.5% of WRKY members showed variation in the conserved heptapeptide domain and had substitutions in the core domain. Out of 14.5% variants, 13 members (ClWRKY37,38,39,49,50,51,52,53,54,74,76,77,and 98) have "WRKYGKK, " substituting glutamine as a lysine aa. ClWRKY12 replaced histidine instead of lysine, two genes, ClWRKY20 and ClWRKY21, shared sixth and seventh as glutamine and asparagine substitutions. High diversity has been found in three genes, i.e., (ClWRKY27) "WKKYGEQK, " (ClWRKY97) "WKKYGEKK, " and (ClWRKY131) "WRKNGQN" WRKY domains, respectively (Figure 1).
For the determination of important properties such as aa length, molecular weight, and isoelectric point, we thoroughly analyzed the WRKY gene. The average length of WRKY proteins was 313 aa, while 74 aa was the minimum and 697 aa was the maximum length. The WRKY members with the least and highest lengths were ClWRKY58 and ClWRKY36, respectively. The molecular weight ranged from 0.02 kDa (ClWRKY55) to 76.26 kDa (ClWRKY36). The pIs (isoelectric point) ranged from 4.73 (ClWRKY75) to 10.63 (ClWRKY97). Among 138 ClWRKY members, 82 were acidic (having less than 7 pI), while 56 were basic (having more than 7 pI). More details of ClWRKY features are available in Supplementary Table 1.

Phylogenetic Analysis and Classification of ClWRKY Proteins
A phylogenetic tree was constructed between C. lavandulifolium and A. thaliana to investigate the relatedness of different WRKY genes and to understand the classification of different groups among WRKY members. Arabidopsis is a model plant for genetic study and has been explored for the function of WRKY genes. The MEGA-X software was used for agglomerative construction following the neighbor-joining method to find the correlation between WRKY members of C. lavandulifolium and A. thaliana. The WRKY members are divided into three groups or seven subgroups. Group I has 27 members, and group II has 58 members with 9 members in IIa, 10 members in IIb, 10 members in IIc, 10 members in IId, and 19 members in IIe. Group III has 53 members and was found to be the largest group in this classification (Figure 2). The majority of the ClWRKY members clustered in the corresponding groups. Some of the WRKY members (ClWRKY20, 21, 55, 56, and 58) had one WRKY domain, but still, clustered in group I; by the rules, group I members should have two domains. These WRKY members may evolve from two WRKY domain genes and lose one domain in the process of evolution. Previously, such exceptional cases of WRKY members with one WRKY domain clustering in group I have been reported (Bi et al., 2016).
Some of the ClWRKY members clustered with important Arabidopsis WRKY genes that helped us identify the genes that are associated with flowering. ClWRKY36 is clustered with Arabidopsis AtWRKY2 and AtWRKY34 (Figure 2), both of which play an important role in pollen development through interaction with the VQ20 protein (Lei et al., 2017). Similarly, ClWRKY45 clustered with AtWRKY71 (Figure 2), which play a key role in the regulation of flowering, directly and indirectly, by influencing other flowering-associated genes .

Gene Structure of WRKY Members and Protein Motif Analysis
For the determination of intron-exon structure similarity, Gene Structure Display Server (GSDS) 6 database has been used. Variation among ClWRKY members was observed in terms of a number of introns. The minimal number of introns in a gene was one, which occurred in seven members including ClWRKY27, ClWRKY58, ClWRKY59, ClWRKY76, ClWRKY86, ClWRKY92, and ClWRKY125 while ClWRKY90 had 9 introns, which is the maximal in the whole family. The overall average number of introns in the ClWRKY family was three (Figure 3). The number of introns in group I ranged from 1 to 5, group IIa varied from 3 to 9, group IIb from 2 to 6, group IIc from 1 to 2, group IId from 1 to 4, and group IIe from 1 to 4. Group III varied from 1 to 6 (Supplementary Table 2). The minimal number of exons was two, which appeared in ClWRKY27, ClWRKY58, ClWRKY59, ClWRKY76, ClWRKY86, ClWRKY92, and ClWRKY125, while the maximal number of exons was 10, which occurred in ClWRKY87 and ClWRKY90 (Supplementary Table 2). The overall average number of exons in the ClWRKY family was four (Figure 3). The exon numbers varied within the groups; group I ranged from 2 to 6, group IIa from 4 to 10, group IIb from 3 to 7, group IIc from 2 to 3, group IId from 2 to 5, group IIe from 2 to 4, and group III from 2 to 7, respectively. The average number of exons within the groups is as follows: group I has 4, group IIa has 6, group IIb has 5, group IIc has 3, group IId has 3, group IIe has 3, and group III has 3, respectively (Supplementary Table 2). Protein motif analyses were performed on the ClWRKY gene family. The first two motifs represented by green and yellow colors are WRKY protein motifs (Figure 4). At least one or two WRKY protein motifs have been identified among all the gene sequences. This suggests that WRKY is a conserved domain.

CpG Island Analysis of the Promoter Region of ClWRKY Gene Family
CpG island is the region of DNA that is rich in CpG sites and is often located near the promoters of genes. Ten members were found to have CpG islands; among them, seven members contained CpG in their promoter region, two in the crosspromoter and gene region, and one member was found to have it in the gene body. All the ten members had one CpG island. ClWRKY115 had the shortest island with 438 bp, while ClWRKY83 had the longest island with a length of 1,065 bp ( Table 1). The genes with the CG-rich content are the obvious means of DNA methylation and a prominent source of epigenetic improvement.

Identification of DNA Methylation
Total genomic DNA of C. lavandulifolium was extracted from the root, stem, leaf, little buds, median buds, and first flowers (flowers at just opening time). These tissues were taken at different growth stages of the plant, which include the young stage, adult stage, and flowering stage. The methylation and unmethylation primers were designed for the ten different genes using the "Methyl Primer Express version 1.0" software. These genes were selected due to the presence of CpG islands. The MSP was carried out for the following ten genes: ClWRKY9, ClWRKY31, ClWRKY33, ClWRKY41, ClWRKY83, ClWRKY97, ClWRKY100, ClWRKY115, ClWRKY129, and ClWRKY130. The methylation status detected by MSP was analyzed through gel electrophoresis. The MSP results showed the presence of methylation in ClWRKY31, ClWRKY100, and ClWRKY129 at different stages in different tissues. These genes showed significant methylation status in all tissue samples at all growth stages of the plant (Figures 5A-C). However, the methylation status detection through bisulfite sequencing can give us more clear results compared with the MSP method.

Cis-Acting Elements in the Promoter Region of ClWRKY
The cis-acting elements in the promoter region are really important for gene expression (Rombauts et al., 1999). The 2,000 bp sequence upstream of the genes was submitted to the PlantCARE software and searched for cis-acting elements. A range of cis-acting elements was identified in the promoter region of the ClWRKY family, which were highly diverse.
A total of 110 cis-acting elements were found in the promoter region of ClWRKY members. Most of these elements have a diverse role in the plant life cycle. According to our objectives, we screened out the elements that can play a role in plant growth and development and particularly in flowering induction. Therefore, we selected 15 cis-acting elements on the basis of their possible role in flowering-related traits and the overall development of FIGURE 2 | A phylogenetic tree of 138 ClWRKY proteins of Chrysanthemum lavandulifolium (blue) and 74 AtWRKY proteins of Arabidopsis thaliana (red) on the basis of amino acid sequences. The tree was constructed in MEGA X. The domains were clustered into three groups, namely, group I, group II, and group III. The group II was divided into five subgroups (a-e). the plant. The detailed features of the 15 cis-acting elements are provided in Supplementary Table 4.
Most of the genes have important cis-acting elements, such as auxin responsiveness (AuR) and gibberellin responsive elements (GRE), associated with flowering-related traits (Figure 6). W-box was most widely available among the ClWRKY family. ClWRKY45 and ClWRKY36, the genes that resulted in high expression at the capitulum stage, had some important elements involved in endosperm expression, elements involved in auxin responsiveness, and ethylene-responsive elements (AuR, EE, and ERE), which are supposed to play a role in flower development. Methyl jasmonate (MeJAR), which plays important role in flowering induction, was present in most of the ClWRKY members, especially in group III members, which indicates that WRKY members may play a role in the overall development of a plant.
The de novo motif analysis in MEME identified 15 motifs (Supplementary Figure 2). Tomtom tool comparison showed homology with some important motifs, such as motifs similar to MEME12 and MEME14, which are involved in seed germination, pollen tube growth, and WRKY binding sites (Supplementary Table 5). The enriched motifs identified in MEME are mostly different from the motifs found in the PlantCare software. The simple enrichment analysis (SEA) in MEME identified 95 known enriched motifs. Among them, 49 motifs encoded ethylene, which is essential for plant growth and development. Other important motifs included are WRKY binding sites, TEOSINTE BRANCHED, cycloidea and PCF (TCP) binding sites, and so on (Supplementary Table 6).

qRT-PCR Analyses
Different CIWRKY genes having special characters were selected on the basis of previous analysis. Total RNA was extracted from different tissues (root, stem, leaves, and flowers) of C. lavandulifolium at different stages (young, adult, and flowering). The expression pattern of 19 WRKY genes was observed in all tissues at different stages of the plant's life span (Figure 7). Among them, 10 genes were selected due to the presence of CG-rich regions, and nine genes were selected as they showed close association with important Arabidopsis WRKY genes that are associated with floral traits. Most of the genes exhibited higher expression at the flowering stage in the capitulum tissues. The highest expression in the capitulum tissues at the flowering stage was exhibited by ClWRKY36, ClWRKY45, ClWRKY57, ClWRKY59, and ClWRKY82 (Figure 7). The expression of these genes was higher in the flowering stage, especially in the first flowering (FF) stage and median bud stage, which suggests their role in flowering regulation. These genes were selected for qRT-PCR as they clustered with flowering-related AtWRKY genes (Figure 2). On the basis of previous analysis, ClWRKY36 clustered with CIWRKY71, so it might have a basic role in the control of flowering time and plant development, as WRKY71 homolog was found to promote flowering in woodland strawberries (Lei et al., 2020). ClWRKY31, ClWRKY33, ClWRKY41, ClWRKY97, ClWRKY100, ClWRKY129, and ClWRKY130, which were selected on the basis of having CpG rich regions (Supplementary Table 1), showed higher expression in capitulum tissues (Figure 7). In addition, ClWRKY31, ClWRKY100, and ClWRKY129 were detected with DNA methylation (Figure 5), which is linked with flowering regulation and is one of the basic sources of epigenetics. Interestingly, ClWRKY83 remains stable across all the tissues and stages as it is constantly expressed in all samples of C. lavandulifolium. It might play a role in adaptability. Moreover, ClWRKY83 showed close association with important AtWRKY7 of Arabidopsis. The other WRKY member's expressions showed less or non-significant changes at different stages in all tissues.
We also investigated CILFY (leafy), as it not only promotes flowering time but also influences other genes that regulate flowering-related traits (Yang et al., 2016). The expression of this gene was high in the median bud and adult leaf stages, but low in stem and root (Figure 7). Besides, ClLFY may promote flowering time in C. lavandulifolium by making association with ClWRKY36 and ClWRKY45.

DISCUSSION
WRKY TFs are one of the largest transcription families, playing an important role in plant development and resistance against different biotic and abiotic stresses. WRKY TFs are spread broadly across higher plants . WRKY family investigation is limited to C. morifolium with 15 genes . Being a model plant for chrysanthemums (Dongru et al., 2019), the characterization of the WRKY gene family in C. lavandulifolium is of great importance. For the first time, we identified 138 WRKY members in C. lavandulifolium. Out of 138 members, 118 were properly conserved in the "WRKYGQK" domain. Substitution of glutamine occurs instead of lysine in thirteen WRKY heptapeptide domains. This WRKYGKK sequence is reported to be a major variant in many studies Song et al., 2016a,b). ClWRKY12 (WRKYGHK) "Q" replaced by "H" ClWYKY20 and 21 (WRKYGQN) "K" replaced by "N, " ClWRKY27 (WKKYGEQK) "K" replaced by "Q, " ClWRKY97 (WKKYGEK) "Q" replaced by "E, " and ClWRKY131 (WRKNGQK) "Y" replaced by "N" as their respective domains were also reported in carrot (Nan and Gao, 2019). The variation in the heptapeptide domain of WRKY genes has been largely demonstrated (Yamasaki et al., 2005;Brand et al., 2013;Rinerson et al., 2015;Yang et al., 2020); however, variations in aa sequence were different from our findings. The variation in the WRKYGQK conserved domain mainly occurs from Q to K aa. The same kind of substitution was observed in the WRKY gene family of Camelina sativa . The 138 ClWRKY genes were analyzed for CG-rich regions using the "Methyl Primer Express version 1.0" software. The ClWRKY gene family is composed of 138 genes, which is bigger than the AtWRKY gene family (Rushton et al., 2010). The ClWRKY family was classified into three groups. The classification of ClWRKY members was according to the classification of the WRKY family in Arabidopsis. Most of the ClWRKY members clustered in the corresponding groups. However, some of the ClWRKY members with one WRKY domain clustered in group I, which included, ClWRKY58,95,60,56,20,21,and ClWRKY55. Among them, ClWRKY20, 21, and 55 possessed a unique pattern "WNSIVV" which was conserved in these three members. Clustering of one domain WRKY of factors in group I is consistent with the study of the WRKY TF family in Linum usitatissimum (Yuan et al., 2021). The same results are reported in Cicer arietinum (Kumar et al., 2016).
Some WRKY members with one conserved domain and clustered into group I have been reported in peanut (Arachis hypogaea L.) (Zhao et al., 2020).
The intron-exon structure seems to be diverse, starting from 1 to 9 introns. Diversification in intron-exon structure may play an important role in the evolution of genes (Xu et al., 2012). In this study, we have identified different intron constitutions within gene groups. Group I has four introns on average, group IIa has 5, group IIb has 4, and all other groups have an average of three introns. This suggests that the latter groups originated from group I and earlier subgroups of group II. The same trend of intron presence in different groups of WRKY members in eggplant was reported (Yang et al., 2020). A total of six motifs were detected; among them, two motifs were considered to be WRKY domains. The conserved domain "WRKYGQK" was almost included in all of the members; however, some of the WRKY members showed a slight variation as "WRKYGKK" domain as suggested by Liu et al. (2021) and Qu et al. (2021).
The cis-acting elements are involved in plant development and flowering induction. There are several elements, such as auxin, gibberellins, and jasmonic acid, which affect floweringrelated traits. Jasmonates, cytokinins, and auxin control second bud growth in chrysanthemum . Auxin controls plant growth and development (Singla et al., 2006). The petal is the most important factor influencing the chrysanthemum's ornamental value. Auxin signaling is associated with petal growth in chrysanthemums (Wang J. et al., 2017). AuR is found in most ClWRKY members and may play a role in plant development and flowering induction. GRE was abundant in ClWRKY, an important hormone associated with flowering induction in chrysanthemum (Dong et al., 2017). CmBBX24 regulates flowering in chrysanthemum in conjunction with gibberellin synthesis (Zhu et al., 2020). In addition, the key recognition site for the WRKY gene, the W-box containing the TGAC sequence, was observed in most ClWRKY members. The W-box contributes under abiotic stress conditions (Dhatterwal et al., 2019) by influencing the early senescence of rice flag leaves . The presence of flower-associated cis-acting elements in most genes indicates the potential of ClWRKY to play a role in flowering control in C. lavandulifolium.
DNA methylation has an impact on growth and development, with CG content, or CpG, being the most abundant source of DNA methylation (Dong-ru et al., 2019). Various methods such as methylation-sensitive amplification polymorphism and MSP have been applied to methylated genes (Herman et al., 1996;Coronel et al., 2018). In this study, we have identified ten genes with strong CpG content in the promoter region (Table 1), which may cause epigenetic changes in C. lavandulifolium (Chip and Dresselhaus, 2019). Three genes having methylated sites, i.e., ClWRKY31, ClWRKY100, and ClWRKY129, were identified on the basis of MSP results. These genes can regulate gene expression and may play a role in flowering-related traits, as they are expressed significantly in capitulum tissues at the flowering stage. CmMET1 can affect methylation and regulate candidate genes associated with flowering in C. morifolium (Kang et al., 2021). DNA methylation changes in the promoter and gene body regions can influence gene expression and function of the protein FIGURE 6 | Cis-acting element analysis in the promoter region of the ClWRKY family. The PlantCARE software was used to search the cis-acting elements. Abbreviation of important cis-acting elements are as follows: EE, endosperm expression; ARE, auxin-responsive elements; GR, gibberellin responsiveness; LR, light responsiveness; REEAI, regulatory element essential for anaerobic induction; AuR, elements involved in auxin responsiveness; SSR, seed-specific regulation; MeJAR, methyl jasmonate responsiveness; ERE, ethylene-responsive elements; GRE, gibberellin responsive elements; ESNE, endosperm specific negative expression; LR2, elements involved in light response; PAP2L, plant-AP2 such as AuR2, part of auxin responsiveness; and W-box, recognition site of WRKY TFs.
Frontiers in Plant Science | www.frontiersin.org FIGURE 7 | Expression analysis of the selected ClWRKY genes in different tissues (capitulum, leaf, stem, and root) at different stages (young, adult, and flowering). Expression of 19 genes: 10 genes were selected due to the presence of CpG island, nine genes were selected based on their close association with important Arabidopsis genes found in the phylogenetic tree built between C. lavandulifolium and A. thaliana. Abbreviations of different growth stages and tissues are as follows: YSR, young stage root; YSS, young stage stem; YSL, young stage leaf; ASR, adult stage root; ASS, adult stage stem; ASL, adult stage leaf; FSR, flowering stage root; FSS, flowering stage stem; FSL, flowering stage leaf; FSLB, Flowering stage little buds; FSMB, Flowering stage buds; FSF, Flowering stage first flower. (Miguel and Marum, 2011), the morphology of floral characters, photoperiod-related sensitivity, and fruit ripening (Sun et al., 2014;Xiao et al., 2020). We studied only four stages to detect methylation, and it will be better to explore the methylation status in more stages in the future.
This study focuses on the expression of ClWRKY members to find their role in plant development, especially their influence over flower induction. The expression analysis revealed variations of different genes in different organs. ClWRKY45 showed significant expression at the flowering stage in the capitulum tissue, which may be related to some morphological characters and floral organ development. In the same fashion, ClWRKY45 clustered with Arabidopsis AtWRKY71, a prominent gene that controls flowering characteristics , providing strong evidence for our results. Moreover, ClWRKY36 and ClWRKY57 are significantly expressed in capitulum at the flowering stages and have close homology with Arabidopsis AtWRKY2 and AtWRKY34, indicating the role of CIWRKY36 and ClWRKY57 in flowering, as Arabidopsis homologs are the key regulators of pollen development (Lei et al., 2017). ClWRKY59 is expressed significantly in the first flower tissue and can positively regulate flowering, as ClWRKY59 homolog, AtWRKY75, promoted flowering in Arabidopsis . ClWRKY82 exhibited higher expression in the median bud stage, which indicates its potential to play a role in flowering regulation. AtWRKY7, a ClWRKY82 homolog, regulates plant growth (Chen et al., 2019). The qRT-PCR results indicate that ClWRKY83 is a stable gene as it expresses across all the tissues at different stages of plant life and may play a role in adaptability. The higher expression of three genes in capitulum tissues, including ClWRKY31, ClWRKY100, and ClWRKY129, which are methylated as well, can play an important role in floweringrelated traits. These results suggest that ClWRKY33, ClWRKY36, and ClWRKY45, along with ClWRKY31, ClWRKY100, and ClWRKY129, can play a crucial role in the genetic and epigenetic improvement of the chrysanthemum. ClWRKY45, ClWRKY36, and ClWRK129, in particular, have the potential to alter flowering-related traits. The overall investigation of the CIWRKY domain is providing theoretical insight into C. lavandulifolium transcriptional activation and inactivation mechanisms for the improvement of plant development.

CONCLUSION
This study aimed to identify the WRKY gene family in C. lavandulifolium and explore the potential of WRKY genes for genetic improvement. A total of 138 WRKY genes were analyzed in C. lavandulifolium. All the data results are unique to explore the underlying mechanism of the WRKY TF in C. lavandulifolium plant improvement. The gene structure analysis suggested the vast variation among WRKY members. The WRKY domain remains conserved across all the members of this family. Ten WRKY genes possessed CG-rich content in the promoter region where three genes exhibit methylation status.
This research can provide a basis to study the role of WRKY genes in the improvement of chrysanthemum plants, especially in flowering traits.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
WZ designed the research project, coordinated, and supervised the study. MA performed all the experiments and wrote the manuscript. MA and KD did most of the experimental work and analyzed the data. WuY, AP, and WaY performed specific experiments during the research. All authors read and approved the final manuscript.