The complete genome and proteome sequences of Arabidopsis were downloaded from the Arabidopsis Information Resource ¹ . The T. sinensis data reported in this study are available under Accession No. CNP0000958 in the CNGB Nucleotide Sequence Archive ² . The hidden Markov model (HMM) file of the WRKY domain (Accession Number PF03106) was downloaded from the Pfam database ³ (Mistry et al., 2021), and HMMER3.0 was used for the identification of WRKY genes with an E-value setting of 1e-5. In addition, SMARAT ⁴ and CCD ⁵ were used to confirm all the potential TsWRKY genes (Chen et al., 2018; Lu et al., 2020; Letunic et al., 2021). The molecular weights (Mw), instability index (II), aliphatic index (AI), the Grand Average of Hydropathicity (GRAVY), and isoelectric points (pI) of the identified WRKY proteins were predicated on the Expasy website ⁶ (Gasteiger et al., 2005). The subcellular locations were predicted using ProtComp – Version 9 from Softberry website ⁷ .

Multiple sequence alignments were created using ClustalW using default settings (Thompson et al., 2003), and then the WRKY proteins conserved domain sequences were modified in GeneDoc software (Nicholas, 1997). The distribution pattern of intron was analyzed by the Gene Structure Display Server (GSDS) (Hu et al., 2015) ⁸ . Conserved motif analysis of the identified T. sinensis WRKY proteins was carried out on the MEME online program (Bailey et al., 2015) ⁹ . The optimized parameters of MEME are as follows: the maximum number of motifs is 20, the motif width is between 8 and 50 aa, and the rest of the parameters are default. The promoters, which were extracted from 2000 bp upstream of the CDS region of TsWRKYs, were used for Cis-regulatory elements (CREs) prediction analysis by PlantCARE online software (Lescot et al., 2002) ¹⁰ .

The Circos (Krzywinski et al., 2009) was used to map the TsWRKY gene information based on T. sinensis genomic data. Tandem and segmental duplication of T. sinensis WRKY genes, as well as the synteny relationship between T. sinensis and six plant species genomes, were evaluated using MCScanX (Wang et al., 2012) and TBtools (Chen et al., 2020). The KaKs Calculator 2.0 (Wang et al., 2010) was used to estimate the non-synonymous (Ka) and synonymous (Ks) substitution of each duplicated WRKY gene. The sequence of WRKY proteins from Arabidopsis, tomato, citrus, maple, pineapple, and rice was downloaded from the NCBI ¹¹ .

The conserved domains from the predicted WRKY proteins sequences were confirmed using multiple sequence alignments. The amino acid sequences of WRKY proteins in A. thaliana and T. sinensis were alignment by ClustalW. The phylogenetic trees were constructed using the Neighbor-Joining (NJ) method in MEGA 7.0 (Kumar et al., 2016) with the following pa-rameters: p-distance, pairwise deletion, and 1000 bootstrap replications. The neighbor-joining tree construction method of T. sinensis and other green line species refers to the research of Rinerson (Rinerson et al., 2015). The green line species including: Micromonas pusilla, Ostreococcus tauri, Ostreococcus lucimarinus, Dunaliella salina, Chlamydomonas reinhardtii, Gonium pectorale, Volvox carteri, Physcomitrella patens, Selaginella moellendorffii, Brachypodium distachyon, Oryza sativa, Glycine max, Arabidopsis thaliana, and T. sinensis.

The T. sinensis var. ‘Heiyouchun’ used in this study is universally recognized as the best variety because of its nutritional value, good taste, and unique aroma. ‘Heiyouchun’ was grown in the field at the Forestry Nursery of the Taihe County, Fuyang City, Anhui Province, China (118°48’8’’ N and 32°3’52’’ E). Young and healthy leaves with at least six branches and 5–10 cm in length were collected in four different sampling periods from March 30 to April 20, 2021. All collected samples were immediately frozen in liquid nitrogen and stored at −80°C.

Total RNA from samples was extracted using the RNA Extraction Kit 3.0 (Huayueyang Biotech, Beijing, China), with the RNase-free DNase I treatment to remove potential genomic DNA contamination. Qualified RNA was chosen as a template to produce the first-strand cDNA, as determined by gel electrophoresis and the A260/A280 ratio. Complementary cDNA was generated with SuperScript cDNA Synthesis Kit WX2050 (Huayueyang Biotech, Beijing, China). The specific TsWRKY gene primers were designed using Primer Premier 5, and TsActin gene served as the reference gene for normalization of the expression levels in different sampling periods. Supplementary Table S1 presents all the primer information. The qRT-PCR was performed with a 2×SYBR Green qPCR Mix (With ROX) (Sparkjade, Shandong, China), and amplification was performed using 96-well plates and CFX96 TouchTM RT-PCR system (Biorad, Los Angeles, CA, USA). Each reaction was performed in biological triplicates. The data from qRT-PCR amplification were analyzed using the 2^−ΔΔCt method. A calculation of Bonferroni’s multiple comparisons test was performed using SPSS statistical software version 25. A mean fold change greater than 2 and a p value less than 0.05 were considered significant differences between the two groups.

To thoroughly investigate the candidate WRKY genes in T. sinensis, 78 TsWRKY genes were finally identified, designated as TsWRKY1–TsWRKY78 based on the order of their HMM (Hidden Markov Model) search results, and were used for subsequent analysis. All extensive information TsWRKYs, including chromosomal location, subcellular localization prediction, protein length, molecular weight, GRAVY, instability index, and aliphatic index are provided on Table S2 . Among the 78 TsWRKY proteins, TsWRKY14 and TsWRKY39 proteins were determined to be the smallest and the largest proteins, with 116 and 1161 amino acids (aa), respectively. The proteins’ mo-lecular weights varied from 13.5 to 125.8 kDa, and their pI values were from 4.94 to 9.73. According to the expected subcellular localization results, 67 and 11 TsWRKY proteins were found in the nucleus and extracellular areas, respectively. All of the TsWRKY proteins have a GRAVY of less than 0, which means that they are all hydrophilic proteins, and more information is shown in Table S2 .

Multiple sequence alignments of the WRKY domains, which cover about 60 amino acids, were used to analyze the evolutionary relationships of TsWRKY proteins. The WRKY domains of seven distinct Arabidopsis WRKY proteins (AtWRKY1, 18, 6, 8, 7, 14, and 30) were randomly chosen as representatives for further evaluation. Figure 1 shows the substantially conserved WRKY domain sequences. The majority of the proteins in this family (76 out of 78) share the conserved WRKY domain WRKYGQK, while TsWRKY33 and TsWRKY62 differ by an amino acid.

The phylogenetic analysis using the WRKY genes of Arabidopsis as a reference revealed that these TsWRKY genes are more precisely classified into groups I, II (a–e), and III (Eulgem et al., 2000; Zhang and Wang, 2005). Among the 78 members of the TsWRKY family, group II contains the most TsWRKY proteins (54), followed by group I (11) and group III (11). Additionally, each category may be subdivided into numerous subcategories. There are 11 TsWRKY proteins in group I that contain two WRKY conserved domains that are classified as N-terminal WRKYs (IN) or C-terminal WRKYs (IC) depending on their locations on the protein. Group II of the TsWRKY proteins can be divided into five subgroups, including two subgroups IIa, six subgroups IIb, 27 subgroups IIc, 10 subgroups IId, and nine subgroups IIe ( Figure S1 , Table S2 ). The 11 TsWRKY proteins in group III contain a zinc finger motif of the form C-X7-C-X23-H-X-C, which is identical to the AtWRKY30 in Arabidopsis subgroup III ( Figure 1 ). We created a broader WRKY domain dataset for phylogenetic analysis to further understand the evolution of TsWRKYs familly. Figure 2 shows the neighbor-joining phylogenetic tree constructed by WRKY domains for 14 species. A scattered distribution of TsWRKYs was observed in groups I, II (a–e), and III. The TsWRKYs protein sequence branched away from algae, bryopsida, and pteridophyta but showed clustering with dicotyledons and monocotyledons within each subgroup. Interestingly, the TsWRKYs protein sequence branches tend to show close proximity between the two. For example, TsWRKY35-TsWRKY57 and TsWRKY25-TsWRKY51 pairs in group IIb, and TsWRKY13-TsWRKY15, TsWRKY21-TsWRKY67 and TsWRKY48-TsWRKY65 pairs in group IId ( Figure 2 ; Figure S2 ). Group III WRKY family members could be significantly subdivided into eight clades, but all TsWRKYs proteins were found on Clade 1, 4, 6, and 8.

Figure 3B depicts the particular condition of the T. sinensis WRKY gene structures. The number of introns, with the exception of TsWRKY39 and TsWRKY59, ranges from 2 to 6, with an average of 3.59. The TsWRKY genes structure are composed of three exons and two introns in more than 60% (47 of the 78) of them. TsWRKY39, in particular, has the most exons and introns of all TsWRKYs, with 12 exons and 11 introns. Gene structures of genes in the same group, like IId and IIe, tend to be consistent in general.

The domain prediction results ( Figure 3C ) were validated by utilizing the MEME web server. By sequencing the TsWRKY proteins, 20 distinct motifs were discovered, comprising 8-50 amino acids. The majority of TsWRKYs contained motifs 1, 2, and 5, which corresponded to the DBD domain, while others contained motifs unique to each class. For instance, motif 8 is unique to group IId, while motif 10 is found only in groups IIa and IIb. Most importantly, each class had a distinct motif organization, and two genes that were tightly clustered on the tree usually exhibited identical motif patterns. Despite their heterogeneity in size and sequence, the projected WRKY domains and other conserved domains were cross-confirmed by the two combined approaches, implying that the group classifications are reliable.

CREs are genomic sequence motifs located in the 5’ upstream region of genes that bind to motif-specific proteins and function as regulatory switches for downstream genes (Korkuć et al., 2013). As shown in Figure 4 , the upstream 2000 bp regulatory regions of all TsWRKYs were extracted, several CREs were predicted using PlantCARE, and the 20 most common were visualized using TBtools software.

Our analysis revealed that T. sinensis contained many promoters’ core regulatory elements (CAAT-box, TATA-box), light responsive elements (Box 4, G-box, GT1-motif, AE-box, and TCT-motif), and W box elements. We observed a lot of abiotic stress responsive elements as well, such as wound-responsive elements (WUN-motif), drought-inducibility elements (MBS), dehydration, low-temp, salt stress responsive elements (DRE), low-temperature responsive elements (LTR and WRE3), and defense and stress responsive elements (TC-rich repeats). These are the hormone responsive elements: abscisic acid responsive elements (ABREs), methyl jasmonate (MeJA) responsive elements (CGTCA-motif and TGACG-motif), ethylene-responsive elements (EREs), auxin-responsive element (TGA-element), gibberellin-responsive elements (GARE-motif, P-box and TATC-box), and salicylic acid responsive element (TCA-element). Other CREs were also predicted, such as anaerobic responsive elements (AREs) and circadian control elements (circadian).

All TsWRKYs had at least one stress response-related CREs in this investigation. A total of 65 TsWRKY genes (83.3%) had one or more ABREs, which could be a sign that they have an ABA response when they are stressed. Additionally, more than 70% of TsWRKY genes have the EREs and AREs that have been speculated as having important promoter roles (Olive et al., 1991; Ohme-Takagi et al., 2000). TsWRKY27 and TsWRKY72 contained 16 out of 20 promoters in their promoter regions that surpass other TsWRKY genes. We also focused on CREs involved in wound, hypothermia, and drought responses, such as WUN-motif, LTR, MBS, and TC-rich repeats. The WRKY protein can be used efficiently in conjunction with W-box regions to activate or inhibit downstream target gene transcription (Jiang et al., 2017). It can form protein complexes with other active components, which improves transcription binding activity (Chi et al., 2013). Moreover, 62 TsWRKYs possessed one or even more W-boxes, implying that these WRKY genes are regulated by autoregulation or crossregulation (Rushton et al., 2010).

The 78 TsWRKY genes were dispersed randomly throughout the 28 T. sinensis chromosomes ( Figure 5 ). The bulk of the TsWRKYs were found at or at the ends of chromosomes. Of all TsWRKYs, 11 were identified on Chr24, scattered in several clusters, which is the largest number. On the contrary, there is only a single TsWRKY gene on Chr2, Chr4, Chr9, Chr10, Chr17, Chr20, Chr21, and Chr22. Tandem and segmental duplications both contribute to the generation of gene families throughout evolution (Cannon et al., 2004). Hence, we investigated the occurrences of TsWRKY genes duplication. On cross-referencing with Holub’s published research study, 20 TsWRKY genes (25.7%) were found to be tandem duplicated. Tandem duplication event is a chromosomal region within 200 kb, including multiple (two or more) members of a gene family (Holub, 2001). There were 10 distinct pairs of tandemly duplicated genes on Chr 1, 6, 11, 12, 13, 15, 16, 23, and 24. In addition to the tandem duplication events, 83 segmental duplication events involving 72 TsWRKY genes were discovered using the BLASTP and MCScanX approaches.

Constructing syntenic graphs between T. sinensis and several typical species allows us to investigate the evolutionary clues for the T. sinensis WRKY gene family ( Figure 6 ). The representative species consist of four dicots, Citrus clementina, Acer yangbiense, A. thaliana, and Lycopersicon esculentum, and two monocots, O. sativa and Ananas comosus. Syntenic links were found between 77 TsWRKY gene members and those in citrus (73), maple (70), tomato (70), Arabidopsis (63), pineapple (51), and rice (34). There were 136, 129, 113, 101, 73, and 45 orthologous pairings between the six species (tomato, citrus, maple, Arabidopsis, pineapple, and rice), respectively. In general, the TsWRKYs comprised more syntenic gene pairs in dicots than in monocots. C. clementina and A. yangbiense, well-known members of Sapindales, show greater synteny with T. sinensis, which belonged to Sapindales. Notably, in the interactive Venn map of WRKY genes across species ( Figure 7A ), 30 TsWRKY genes shared syntenic WRKY gene pairings with all six species, implying that these orthologous pairs existed prior to the ancestral split. Certain TsWRKY genes were shown to relate to 3, 4, or 5 collinear gene pairs (between T. sinensis and maple/citrus/tomato WRKY genes), indicating the possibility that these TsWRKY genes have significant roles in the evolution of the WRKY gene family. Syntenic gene pairings between T. sinensis and other species may be important for elucidating WRKY gene evolution. The Ka/Ks (non-synonymous substitution/synonymous substitution) ratios of the WRKY orthologous gene pairs of six species were computed to assess the evolutionary constraints operating on the T. sinensis WRKY gene family. Figure 7B shows that almost all TsWRKY orthologous gene pairs had Ka/Ks < 1. As a result, we speculated that the T. sinensis WRKY gene family may have been subjected to significant purifying selection forces throughout evolution (Hurst, 2002).

Six main expression patterns of TsWRKY genes were observed ( Figure 8 ). TsWRKY8 and 74 were upregulated at first and then downregulated, while TsWRKY12 and TsWRKY65 showed the opposite trend. Suppressed expression patterns were seen in 18 TsWRKY genes (7, 16, 17, 21, 26, 44, 46, 51, 55, 58, 59, 60, 61, 67, 69, 76, 77, and 78), while upregulated expressions were observed in 22 TsWRKY genes (5, 6, 18, 19, 22, 24, 32, 33, 35, 38, 39, 40, 41, 43, 45, 53, 57, 62, 63, 64, 66, and 68). Furthermore, the expression of nine TsWRKY genes (2, 13, 14, 15, 23, 28, 31, 34, and 56) were first upregulated, then downregulated, and then upregulated to a high level. The other TsWRKY genes (1, 3, 4, 9, 10, 11, 20, 25, 27, 29, 30, 36, 37, 42, 47, 48, 49, 50, 52, 54, 70, 71, 72, 73, 75, and 78) were not significantly altered (fold change ≥ 2) at different stages.

We investigated the expression patterns of the genes involved in terpenoid biosynthesis in order to unravel the regulatory mechanism of terpenoid accumulation patterns in various developmental stages of T. sinensis. All of the key genes, except for TsFPPS, showed significant changes (fold change ≥ 2) in expression during development ( Figure 9B ). After April 6, the expressions of TsAACT, TsHMGS, TsHMGR, TsDXS, and TsDXR all changed distinctly. Notably, TsIDI, TsDXS, and TsDXR are more than five-fold changes in expression over four periods.

The WRKYs usually control the expression of terpenoid synthesis genes by activating or repressing their promoters, thereby regulating the accumulation of terpenoids. We constructed a co-expression network of TsWRKYs with terpenoid synthesis genes ( Table S3 ; Figure 10 ). The results showed that TsFPPS, TsIDI, TsMTPS, TsWRKY9, TsWRKY24, TsWRKY35, TsWRKY38, TsWRKY39, TsWRKY62, and TsWRKY64 are possible core members of the terpenoid synthesis co-expression network. TsFPPS, TsWRKY36, and TsWRKY75 were negatively correlated with other genes. By contrast, TsWRKY24, TsWRKY35, TsWRKY39, and TsWRKY40 are almost all positive correlated factors.

Variation in the number of gene family members is a key mechanism for shaping adaptive natural variation during the evolution of species (Guo, 2013). We discovered 78 proper TsWRKY genes in this investigation. The results of the neighbour-joining phylogenetic tree of 14 species constructed with the WRKY protein sequence’s conserved domain indicate that they can be divided into three major groups (I–III) and five subgroups (IIa–IIe). WRKY TFs diverge early in the green lineage and TsWRKYs are more closely related to dicotyledons and monocotyledons. The number of presumed TsWRKYs are comparable to the count of WRKY genes in Fagopyrum tataricum (78) (He et al., 2019), and it is somewhat lower than those in S. lycopersicum L. (81) (Huang et al., 2012), but much lower than those in O. sativa L. ssp. indica (102) (Ross et al., 2007) and Glycine max (174) (Yang et al., 2017). These findings corroborated previous research that suggested herbaceous plants tend to have a larger number of WRKY genes than woody plants (Wu et al., 2016).

The variety of gene structures reflects the historical evidence of gene family evolution and serves as the foundation for phylogenetic categorization (Xiao et al., 2017). Wheat and tea plants have 0–5 introns and 0–11 introns, respectively, whereas TsWRKYs have 2–12 introns, suggesting that TsWRKYs have abundant gene structural variation (Ning et al., 2017; Wang et al., 2019). The intron-exon distribution pattern is comparable across members of the same subfamily, which is the basis for functional similarity among members of the same evolutionary group (Li et al., 2020a). For example, the number of introns in almost all members of group III is 2. Furthermore, TsWRKYs from group I contain many more introns than other groups, which implies that it is more likely that other groups came from group I (Chen et al., 2017a). Twenty conserved motifs were discovered in 78 TsWRKY proteins ( Figure 3 ), with motifs 1 and 2 belonging to the WRKY conserved domains. Partiular motifs that only arose in one group, such as motif 10 in group IIa and group IIb proteins and motif 8 in subgroup IId, which has yet to be characterized for some roles, should be given additional attention.

According to a recent study, dicotyledons have experienced less evolutionary loss of the WRKY conserved domain than monocotyledons (Wei et al., 2012). This occurrence was validated in this investigation, with the majority of TsWRKYs having the conserved heptapeptid WRKYGQK motif, despite two TsWRKY genes having the variants WKKYGQK (TsWRKY33 and TsWRKY62, Group WKKY) ( Figure 1 ). The WRKY proteins demonstrate a propensity for binding to W-box elements, and WRKYGQK motif changes may affect DNA-binding interactions with downstream genes, as previously discovered. As a result, additional exploration of the functional and binding properties of these two WRKY proteins is required (Chen et al., 2017a).

How did the number of members of the WRKY family expand from 1 in the unicellular green algae to 78 in the T. sinensis? The analysis of segmental and tandem duplications contributed to revealing the number and function of the TsWRKY gene family. Based on the chromosome distribution and synlinearity analysis results of the TsWRKYs, 83 segmental duplication events within 72 TsWRKY genes were observed, while tandem duplication events existed. We found that 20 out of the 78 genes (25.6%) in this family are tandem repeats in T. sinensis, suggesting that the abundance of tandem repeats may be a possible reason for the larger number of TsWRKYs. In addition, whole-genome duplication (WGD) events often lead to the growth of gene families, which is common in the evolution of angiosperms. Previous study revealed that T. sinensis has a large number of gene duplications and that WGD events happened approximately 7.8 and 71.5 million years ago (Mya) (Ji et al., 2021). We determined that, whereas some TsWRKY genes are the result of tandem duplication, segment duplication events are the driving force behind gene family evolution. The three basic evolutionary mechanisms are segmental duplication, tandem duplication, and transposition events like retroposition and replicative transposition. Individual gene duplication, chromosomal segment duplication, and even complete genome duplication supply the fresh materials required for gene generation (Yu et al., 2005). Gene duplication, which is associated with the generation of new gene functions, is one of the primary driving forces underlying genome evolution and also is essential to plant adaptive evolution (Moore and Purugganan, 2003; Kong et al., 2007). Tandem duplications produce highly diverse duplicates that have lineage-specific functions (Ezoe et al., 2020), and new research tea plants suggests that tandem duplication of genes plays an active part in flavor accumulation (Wang et al., 2021). Additionally, we speculated that the TsWRKY gene family may have been subjected to significant purifying selection forces throughout evolution because almost all TsWRKYs orthologous gene pairs had Ka/Ks<1. Divergence is a key feature of the evolution of paralogous homologous genes and DNA segments that make up fixed repeats. However, selection on copies of paralogous homologs would be relaxed in the case of full redundancy, which is when any number of functional copies of a gene give the same fitness. Therefore, this negative selection of TsWRKYs are associated with the post-fixation evolution of gene duplications.

The main volatile aromatic compounds of T. sinensis are terpenes (isopentene, monoterpenes, and sesquiterpenes), phenyl/phenylpropanes, and fatty acid derivatives. In T. sinensis, the terpenes have been demonstrated as the most important volatile compounds. The aroma of T. sinensis leaves is an important factor in determining its quality and an important criterion for measuring its economic value. During the germination and maturation process of T. sinensis tender buds, aromatic substances are gradually synthesized and accumulated. The results of real-time quantitative PCR revealed that 67.9 percent of TsWRKY had significantly varied transcript levels during the four stages of budding. Although little is known about the transcriptional regulatory network controlling terpenes synthesis, most of the identified regulators are WRKY TFs (Patra et al., 2013). It is reasonable to conclude that TsWRKYs expression patterns and the buildup of aromatic compounds are related. Therefore, we further analyzed the expression patterns of terpenoid synthase genes and detected that multiple genes of the MVA and MEP pathways are involved in regulating terpenoid accumulation. TsWRKYs showed a high correlation with the expression trends of 10 terpene synthesis genes. For example, TsFPPS, TsIDI, TsMTPS, TsWRKY9, TsWRKY24, and TsWRKY35 may coordinately regulate terpene synthesis.

WRKY genes are essential regulators of secondary metabolite production in plants (Luo et al., 2022), while their regulatory functions varied substantially due to distinct binding mechanisms (Li et al., 2020b). Evidence suggests that certain WRKYs, alone or in concert with other transcription factors, govern the biosynthesis of valuable natural products (Hsin et al., 2022). There have been several studies on the regulatory effects of WRKYs upon the activation or repression of genes involved in plant terpenoid production. Individual WRKY can be associated with a number of regulatory mechanisms, as SlWRKY73 transactivates the SlTPS3, SlTPS5, and SlTPS7 monoterpene synthase genes in tomato (S. lycopersicum) (Spyropoulou et al., 2014). Gossypol (sesquiterpene phytoalexins) in Gossypium arboretum (Xu et al., 2004), DP (diterpenoid phytoalexin) in rice (Akagi et al., 2014), Artemisinin (a type of sesquiterpene lactone) in Artemisia annua (Chen et al., 2017b), ginsenosides (a group of triterpene) in Panax quinquefolius (Sun et al., 2013), and tanshinone (one category of bioactive diterpenes) in Salvia miltiorrhiza are all regulated by WRKYs (Cao et al., 2018). From correlative analysis in sweet Osmanthus fragrans, it is speculated that the OfWRKY gene participates in aroma synthesis by regulating the synthesis of monoterpene volatiles, and that the expression of OfWRKYs are closely related to monoterpene synthesis (Ding et al., 2019). Heterologous expression of WRKY and MYC2 in Salvia sclarea causes coactivation of MEP-biosynthetic genes and accumulation of abietane diterpenes (Alfieri et al., 2018). The co-expression network between several key genes for terpene biosynthesis and TsWRKYs provides important insights into the terpene biosynthesis pathway in T. sinensis. It helps to further characterize the functions of candidate WRKY gene families in T. sinensis and provide new ideas for agronomic genetic improvement and quality variety breeding.