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 short-day 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.

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).

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.

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

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).³ The known enriched motifs were identified via simple enrichment analysis using the MEME software (Bailey and Grant, 2021).⁴ Tomtom motif comparison tool was used to quantify the similarity between motifs (Gupta et al., 2007).

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).⁵ 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.

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™ 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™ 2000/2000c Spectrophotometers.

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.

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).

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.

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 II a, 10 members in II b, 10 members in II c, 10 members in II d, and 19 members in II e. 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 (Yu et al., 2016, 2018).

For the determination of intron-exon structure similarity, Gene Structure Display Server (GSDS)⁶ 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 II a varied from 3 to 9, group II b from 2 to 6, group II c from 1 to 2, group II d from 1 to 4, and group II e 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 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 cross-promoter 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.

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.

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 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).

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.