Characterization of CsWRKY29 and CsWRKY37 transcription factors and their functional roles in cold tolerance of tea plant

Plants are vulnerable to environmental constraints and pathogen attacks throughout their life span^[1,2]. Plant growth and yield are pointedly influenced by these primary biotic and abiotic stressors. Being sessile organisms, the viability of plants profoundly relies on their ability to evolve complex genetic mechanisms through the activation of various physiological mechanisms and alterations in response to environmental stress via robust transcriptional control and accurate signal transduction pathways^[3]. One of the key factors that can mediate the transcriptional regulation processes in plants under stress conditions are transcription factors (TFs). TFs bind upstream of target genes using particular domains, thereby activating or suppressing the expression of target genes that govern various signals^[4]. When it comes to plant-specific transcription factor gene families, WRKY is among the most well-known TF family, comprised of a class of essential TFs involved in physiological change and response to multiple stresses^[5]. Proteins of this family generally consist of one or two widely known conserved WRKY domains including a total of 60 amino acids. The N-terminus generally consists of a highly conserved heptapeptide WRKYGQK domain, and a unique zinc finger-like motif C2H2 or C2HC located at the C-terminus. Gene expression can be controlled by binding the WRKY domain to an essential cis-acting region (TTGACC/T) known as the W-box, located in the target gene promoter. In addition to the W-box, research showed that some WRKY TFs could bind to other elements. For instance, the HvWRKY46 from barley act on the SURE sugar-responsive elements (TAAAGATTACTAATAGGAA)^[6], and OSWRKY13 from rice associates with PRE4 pathogen-responsive element (TGCGCTT)^[7]; such binding of WRKY TFs indicates the diversity in their functional mechanism resulting in the control of transcriptional activity that is synergistic with other factors.

WRKY proteins are knowingly implicated in an array of physiological activities, including seed germination, lateral root development, leaf senescence, flowering time, fruit ripening, and other crucial metabolic pathways^[5]. Moreover, recent studies have revealed the crucial roles of WRKY family members in multiple abiotic stress tolerance; for instance, wheat transcription factor TaWRKY46 was shown to alleviate drought tolerance via the ABA-dependent pathway^[8], while overexpression of TaWRKY2 in seedlings showed drought tolerance in transgenic wheat^[9]. Overexpression of tea plant CsWRKY26 improved drought tolerance of transgenic Arabidopsis^[10]. In another study, CsWRKY13 from Camellia sinensis was shown to be a negative regulator of lignin synthesis^[11]. In addition, the WRKY genes were also widely found to be involved in biotic stresses, such as pathogen defense. For example, the wooden strawberry transcription factor FvWRKY50 plays a vital role in fungal resistance^[12]. Similarly, SlWRKY46 isolated from tomato plants has been shown to negatively regulate defense against B. cinerea in transgenic plants^[13]. Apart from this, WRKY members have also been proven to be potentional candidates for improving cold tolerance in essential plant species. For instance, plants overexpressing PoWRKY1 in Arabidopsis under cold stress showed enhanced seed germination and growth of the root tissues as compared to wild-type plants^[14]. The overexpression of VbWRKY32 improved cold resistance of Verbena bonariensis^[15]. Therefore, it is worthy to explore the distinctive roles of novel WRKYs and potential functional roles in different economically important plant species, such as tea plants.

Tea is one of the best-known non-alcoholic beverages, utilised by a large number of people globally. The medicinal value (containing polyphenols, catechins, and caffeine), rich taste, and aroma make tea a popular non-alcoholic beverage^[16]. From an economic perspective, tea is of paramount importance; among more than 50 tea-producing countries, tea export in China alone reached 2.79 billion tons during the year 2019 (www.fao.org). However, the tea plant prefers a humid climate and grows in mountainous areas, making tea prone to frost and cold stress, which negatively impacts the yield and quality of harvested leaves in almost all types of tea cultivars^[17]. Hence, exploration of biotic and abiotic stress resistance genes for cultivar development with improved yield in economically important crops like tea has always been the main interest of plant breeders. More recently, the genomes of two key tea plant varieties, namely C. sinensis var. sinensis and C. sinensis var. assamica, have been sequenced and made available to the public^[18−20]. Having access to such high-quality genome sequences is a great way to find new genes in tea plants that produce resistant tea varieties to both biotic and abiotic stress. Nonetheless, the comprehensive information on the WRKY gene family from a genome-wide perspective and their potential functional roles in the tea plant have been overlooked. As a result, the discovery and characterization of new WRKY genes across the whole genome will aid in our understanding of how tea plants respond to cold stress and how stress resistance mechanisms work. Although few studies report the WRKY transcription factors in tea plants^{[10,11,21,22]}, deep investigation of the regulatory mechanism of WRKY members in tea plant is still needed, especially the functional roles of WRKYs against cold stress.

We herein performed the genome-wide identification of genes encoding WRKY transcription factors in tea plant. The proteins of WRKY transcription factors in C. sinensis var. sinensis species were analyzed for their phylogenetic relationships, as well as gene structure, putative motif, evolution trajectory, and expression analysis. The study was further extended to the silencing and generation of transgenic Arabidopsis lines expressing CsWRKY29 and CsWRKY37 genes with an aim to investigate their functional roles against cold stress. The current study provides knowledge about the evolutionary and biological significance of the tea WRKY transcription factor family, which will be beneficial in future resistant breeding in tea plants.

A total of 90 putative WRKY genes were identified in the tea plant genome by iTAK^[23]. Among them, 86 genes were selected for subsequent analysis after confirm the presence of WRKY domains (Supplemental Table S1). These 86 WRKY genes were named hereafter CsWRKY1 to CsWRKY86. All of the CsWRKY proteins comprised of at least one typical WRKY domain with lengths ranging from 87 to 1,080 amino acid residuals, with an average of 370 amino acids (Supplemental Table S1). Out of the 86 identified CsWRKY proteins, 21 contained the presence of two WRKY domains, whereas the rest of the proteins presented only one WRKY domain. The heptapeptide domain (WRKYGQK), which is highly conserved among WRKYs was found in 81 WRKY proteins. Additionally, the two variant heptapeptide domains were detected in the remaining five proteins, including WRKYGRK in three proteins (CsWRKY77, CsWRKY78, CsWRKY79), and WRKYGKK in CsWRKY48 and CsWRKY55. Out of 86 WRKY proteins, 79 have a zinc-finger motif of the C2H2 type, and the remaining seven proteins (CsWRKY80, CsWRKY81, CsWRKY82, CsWRKY83, CsWRKY84, CsWRKY85, and CsWRKY86) hold the C2HC type zinc-finger motif (Supplemental Table S1).

To better understand the phylogenetic relationship among WRKY proteins in C. sinensis and A. thaliana, we performed phylogenetic analyses of the conserved WRKY domain sequences based on the Neighbor-Joining (NJ) method using Mega software^[24]. A total of 158 WRKY domains, including 72 AtWRKY and 86 CsWRKY domains, were used to construct the topology (Fig. 1). Based on the topology and structural feature of WRKY domains, all 86 CsWRKY proteins were clustered into three groups (group I, II, and III); group II was further categorized into five subgroups (Fig. 1 and Supplemental Table S2). In total, 21 CsWRKYs were assigned to group I, which contained two conserved WRKY domains (WRKYGQK) and one zinc-finger motif [C-C-H-H]. Group II comprised 58 CsWRKYs, which further consisted of five subgroups, II-a, II-b, II-c, II-d, and II-e, and each had six, 10, 18, nine, and 15 CsWRKY members, respectively. In group II, all members of CsWRKY proteins have a zinc-finger structure (C-X4/5-C-X22/23-HXH) and a highly conserved WRKY domain. Moreover, seven CsWRKY proteins from group III showed the presence of a single WRKY domain and an altered zinc-finger motif C-C-H-C (C-X7-C-X23-H-X-C) compared to groups I and II (Fig. 1 and Supplemental Table S2).

Figure 1. 

Phylogeny and classification of CsWRKY genes in the tea plants. The phylogenetic tree was generated from the conserved WRKY domains obtained from tea plants and A. thaliana using the neighbor-joining method. The circles and triangles represent tea plants and Arabidopsis, respectively. Different colors represent the classes and subclasses of each WRKY class. Green: group I; blue: II-a; purple: II-b; orange: II-c; dark blue: II-d; light green: II-e; red: III. The bootstrap values for some key nodes are presented beside the branches.

Motif composition, exon-intron organization, and cis-elements of CsWRKY genes

Using MEME, we were able to identify the diversity and similarity of the conserved motifs found in CsWRKY proteins, which helped us better comprehend their composition. The members share similar motifs, indicating that they perform comparable roles in general. In total, 10 motifs for each member were identified. Among the 10 identified motifs, motifs 1 and 3 both contained the heptapeptide stretch WRKYGQK (Supplemental Table S3) and were anticipated as the most conserved and fundamental characteristic of the WRKY family. The motifs within the subgroups were found quite similar; for example, subgroups II-a, and II-b both had motifs 1, 2, 4, and 7 except for three proteins (CsWRKY22, CsWRKY35, and CsWRKY36). Likewise, subgroup II-c had motifs 1, 2, and 5, while subgroup II-d and II-e had motifs 1, 2, and 4. It was also noted that some motifs were unique in some groups of CsWRKY protein, such as motifs 3 and 10 both were unique in group I, while motif 8 in group II-b and motif 9 in group II-d (Supplemental Fig. S1). They appeared as a pair in most cases, suggesting that they are functionally related in their corresponding groups.

In addition, we also analyzed the structural features of CsWRKY genes (Supplemental Fig. S2). Among the 86 identified CsWRKY genes, one had 7 exons and 6 introns, five had 6 exons and 5 introns, 14 had 5 exons and 4 introns, seven had 4 exons and 3 introns, nine had 2 exons, 34 had 3 exons and 2 introns. In contrast, the remaining genes had various structures, for instance single exon (CsWRKY76 , CsWRKY77 , CsWRKY78 , CsWRKY79 ) , 10 exons (CsWRKY07 ) and 12 exons (CsWRKY12 ) . Half of the CsWRKY genes (5 in 11) in group III had 3 exons, while the number of exons in group II ranged from 2 to 7 exons (CsWRKY22). Meanwhile, the number of exons in group I varied considerably, ranging from 4 to12.

Transcriptional regulation is accomplished by binding TFs to cis-regulatory elements in the 5' upstream regions of genes. Hence, the 1-kb upstream regulatory regions of CsWRKY genes were subjected to PLACE online webserver for the discovery of stress-related regulatory elements in the promoter regions. As a result, a wide range of stress-responsive elements was identified in the promoters of most CsWRKY genes linking to cold, defense, hormone, and so on (Supplemental Fig. S3). The light and hormone-responsive elements were the most abundant group of elements among the promoters of CsWRKY genes, followed by low temperature-responsive elements, suggesting their potential involvement in the response to low temperatures in tea plants.

Tandem duplication but not whole-genome duplication is likely to drive the amplification of CsWRKY genes in tea plants

We discovered that the tea plant genome encodes more WRKY proteins (86) than certain model plants such as Arabidopsis (72)^[25], cucumber (55)^[26], and grape (59)^[27] (Supplemental Fig. S4). Most (86.05%) of the CsWRKY genes are located in the 15 pseudo-chromosomes of the tea plant genome. Only 12 CsWRKYs genes reside on the unanchored scaffolds (Supplemental Fig. S5). The number of CsWRKY genes distributed on each chromosome ranges from 1 to 10, with an average of five genes for each chromosome. Among all the chromosomes, chromosome 6 contains the largest number of WRKY genes, including three members from subgroup I, six members of subgroup II and one member of subgroup III; followed by chromosome 1 which contains nine WRKY genes (two in subgroup I, five in subgroup II and two in subgroup III). Chromosome 9 contains the least WRKY genes, with only CsWRKY45 presented. It is worth noting that chromosome 15 contains three CsWRKY genes although it has the shortest chromosome length. In addition, we found that more than 30.23% (26/86) of CsWRKY genes tend to be clustered and are adjacent to each other on the chromosomes, suggesting their tandem duplication-driven nature of gene amplification (Fig. 2a).

Figure 2. 

Chromosome distribution and evolution of CsWRKY genes of tea plants. (a) Chromosome distribution of WRKY genes of tea plants. Only the 74 CsWRKY genes that reside on the tea chromosomes are shown. The outer circles represent the 15 chromosomes of the tea plant genome, while the inner circles exhibit the 725 genomic collinear blocks (green links) within the tea plant genome. (b) Distribution of the synonymous mutation rate (Ks) of 8,208 gene pairs located in 725 collinear blocks. The peak of Ks distribution represents the recent whole genome duplication (WGD) event that occurred in the tea plant genome. (c) Ks distribution of CsWRKY gene pairs. The red bar indicates the total number of CsWRKY gene pairs duplicated at the time of the WGD event.

Previous studies have shown that the tea plant genome went through a recent whole genome duplication (WGD) event about 30 million years ago (MYA), which led to a significant increase in the number of genes involved in tea quality and resistance to stress^[18,20]. We therefore further ask whether the amplification of CsWRKY genes is affected by the WGD event. To test this, we first identified 725 genomic collinear blocks in the tea plant genome and found 47 CsWRKY genes located in the blocks (Fig. 2a). We then further calculated the synonymous mutation rate (Ks) values of these gene pairs to estimate their duplication time. Results show that among all the 86 WRKY genes, only four of them have a Ks distribution peaked around 0.36 and were considered to be duplicated at the time of the recent WGD event of tea plants, suggesting fewer contributions of tea WGD event on the amplification of CsWRKY genes (Fig. 2b−c). Besides, we also found that the ratios of non-synonymous substitution rate (Ka) to Ks of all WRKY duplication pairs within genomic collinear blocks were less than 1, indicating that the tea WRKY genes are likely to experience strongly purify selection during their evolution (Table 1).

Expression analysis of CsWRKY genes in eight representative tissues of tea plants

We examined the expression pattern of 86 CsWRKY genes in eight tissues of tea plant, including root (RT), stem (ST), bud (AB), young leaf (YL), mature leaf (ML), old leaf (OL), flower (FL), and fruit (FR) using RNA-seq data. The results found that 71 of all 86 CsWRKY genes (82.5%) were expressed in all eight types of tissues studied (Supplemental Fig. 6a). 96.5% (83/86) of the CsWRKY genes were expressed in RT and FR. While 94.2% (81/86) of CsWRKY genes were expressed in AB and FL. 89.5% (77/86) of CsWRKY genes were expressed in ML and ST. 87.2% (75/86) of CsWRKY genes were expressed in FL and OL. Meanwhile 31.40% (27/86), 20.90% (18/86), 17.40% (15/86), 15.10% (13/86), 18.60% (16/86),15.10% (13/86), 4.60% (4/86) and 17.40% (15/86) of CsWRKY genes showed high expression levels (log₂^FPKM > 5) in RT, ST, AB, YL, ML, OL, FL and FR respectively. Some CsWRKY genes showed distinct differential expression patterns; for instance, CsWRKY35 could not be found in the leaves, and the CsWRKY46 could only be detected in roots. Likewise, CsWRKY01 was highly expressed in RT, and the expression level of this WRKY family member was 28 times higher than that in FL.

The expression of CsWRKY genes under drought and MeJA treatment

As described in Supplemental Fig. S6b, 10 genes, including CsWRKY01/09/23/37/41/43/55/58/74/76 were up-regulated distinctly at 24 h and then down-regulated slightly at 48 h after drought stress. Among the 10 CsWRKY genes mentioned above, nine CsWRKY genes declined continuously, while the remaining one (CsWRKY41) was up-regulated again in 72 h. Three CsWRKY genes (CsWRKY02/20/40) were up-regulated after 24 h, however, keeping high expression for either 48 or 72 h.

Plant hormones such as JA (Jasmonic acid) have well-established roles in modulating plant signaling networks^[28]. JA and its derivatives methyl jasmonate (MeJA) are considered to be related to biotic stress. In this study, the transcriptome profiles of CsWRKY genes at 0, 12, 24, and 48 h after MeJA treatment were analyzed by Illumina RNA-Seq. Overall, MeJA treatments resulted in a wide variety of CsWRKY gene expression profiles (Supplemental Fig. S6c). At the initial MeJA treatment (12 h), 22.10% (19/86) CsWRKY genes were up-regulated, while 59.30% (51/86) CsWRKY genes were down-regulated. Up to 24 h, the expression level of most genes does not change significantly, and 44.2% (38/86) of CsWRKY genes were slightly down-regulated. The expression of CsWRKY genes in MeJA-treated 48 h tea leaves was also diverse. Five CsWRKY genes (CsWRKY05/23/26/37/43) were up-regulated as before, particularly CsWRKY43, whose expression increased five times more than that of CK. Four CsWRKY genes (CsWRKY19, -37, -58, -86) those all under-expressed in CK, MeJA-treatment for 12 h and 24 h, were also up-regulated at this time. It is noteworthy that the expression of CsWRKY86 turned in performanceto be two times higher than CK, and 46 times higher than MeJA-treatment for 12 h.

Expression patterns of CsWRKY genes under cold treatment

We analyzed the expression patterns of the 86 CsWRKY genes during cold acclimation (Supplemental Fig. S6d). At the beginning of cold treatment (4 °C, 0−6 h), 62.80% (54/86) of CsWRKY genes were up-regulated, 10 of which (CsWRKY01/09/10/19/23/26/29/33/43/53) were increased fivefold. 12.80% (11/86) of CsWRKY genes were down-regulated at the same time. As the processing time progressed, eight CsWRKY genes (CsWRKY04/06/15/18/29/38/49/71) were up-regulated steadily in 7 d, and particularly CsWRKY29 presented an expression 10 times higher than CK. On the contrary, 37% (32/86) of CsWRKY genes were down-regulated equably. In the meantime, two of the aforementioned 11 down-regulated CsWRKY genes (CsWRKY41 and CsWRKY 71) were up-regulated likewise. When the temperature plummeted below 0 °C (−4 °C) for 7 d, three CsWRKY genes (CsWRKY18/29/49) were up-regulated as before. However, for the moment, five CsWRKY genes (CsWRKY01/19/27/37/85) which were down-regulated by chilling were up-regulated by freezing. Moreover, when the temperature rose eventually, 25.6% (22/86) of CsWRKY genes were slightly up-regulated, and the expression levels of some genes were restored to those before low-temperature treatment. It was speculated that these members of the CsWRKY gene family play a different role in cold acclimation. The most captivating genes (CsWRKY01/18/19/23/29/37/38/49/59/67), which involved chilling and freezing, would be a suitable focus for further research in the future.

Cloning and nuclear localization of CsWRKY29 and CsWRKY37

In our preliminary experiments, we initially chose four CsWRKY proteins for functional validation of their roles under cold stress, including CsWRKY18, CsWRKY49, CsWRKY29, and CsWRKY37 (Supplemental Table S4). These genes were selected based on their high relative expression under cold stress (Fig. 3a, c & d), or differential expression levels in plant organs, for instance, higher expression levels of CsWRKY29 and CsWRKY37 in flower and fruit, while CsWRKY18, CsWRKY49 in mature leaf and stem (Supplemental Fig. S6d and Fig. 3f & g). Nonetheless, only two of them, namely, CsWRKY29 and CsWRKY37, were further selected for gene slicing and overexpression analysis to check their functions on cold response. The selection of CsWRKY29 and CsWRKY37 was made based on their upregulation expression profiles under cold and freezing treatment conditions (Supplemental Fig. S6). Cloning of CsWRKY29 and CsWRKY37 genes found that both of them contained a conserved WRKYGQK domain and C2H2 type zinc-finger motif, indicating they are the members of group II of WRKY families (Fig. 3b and Supplemental Figs S7 & S8). The qRT-PCR experiments showed that CsWRKY29 was highly expressed in mature leaves, while CsWRKY37 was highly expressed in flowers and roots (Fig. 3f, g). Both the genes were significantly up-regulated during cold treatment at 4°C, especially in the later stage of treatment (72 h); the expression of CsWRKY29 and CsWRKY37 were considerably up-regulated by three and 20 times respectively (Fig. 3c, d), and are potentially the key candidate genes for tea plants in response to cold stress. The subcellular localization of these two CsWRKY29 and CsWRKY37 proteins was then investigated via transient expression in N. benthamiana cells. The results indicated that both CsWRKY29 and CsWRKY37 reside in the nucleus (Fig. 3e).

Figure 3. 

Cloning and subcellular localization of two cold-responsive WRKY genes, CsWRKY29 and CsWRKY37, in tea plants. (a) Expression patterns of the 11 differentially expressed genes of tea plant during cold acclimation. CK: non-acclimated (25 °C day, 20 °C night); CA 1-6 h: fully acclimated (10 °C for 6 h); CA 1-7 d: 10 °C day and 4 °C night for 7 d; CA 2-7 d: cold response (4 °C day, 0 °C night for 7 d); DA-7 d: recovering (25 °C day, 20 °C night for 7 d). (b) Multiple sequence alignment of CsWRKY29, CAWRKY37, and their homologous in Arabidopsis (AtWRKY6). The red and blue boxes represent the WRKY domain and C2H2 type zinc-finger motif, respectively. (c) Relative expression levels of CsWRKY29 and (d) CsWRKY37 genes under cold treatments (4 °C) were revealed by qRT-PCR experiments. (e) Subcellular localization of CsWRKY29 and CsWRKY37 proteins. Expression patterns of CsWRKY29 (f) and CsWRKY37 (g) among different tissues of tea plants.

Silencing of CsWRKY29 and CsWRKY37 in vivo reduces cold tolerance of tea plants

We immersed 1-year-old tea seedlings in 1 ml of 20 μM antisense oligonucleotides (AsODN) solution to suppress the CsWRKY29 and CsWRKY37 gene expression, whereas the control was soaked with sODN (Fig. 4a). The results showed that compared with the control group of tea leaves fed with sODN, the expression levels of CsWRKY29 and CsWRKY37 decreased significantly after feeding with AsODN at 12 h and 24 h, respectively, and the inhibition rate was close to 52.30% on average (Fig. 4b, c). The control and treatment seedlings were then exposed to cold stress (0 °C) for 1 h. We found higher malondialdehyde (MDA) accumulation in CsWRKY29/37-silenced tea plant seedlings (AsODN) compared to control (Fig. 4d, e). In addition, the net photosynthetic rate and maximum photochemical efficiency of photosystem II (Fv/Fm) were measured (Fig. 4f, g). Under cold stress treatment, the Fv/Fm values in CsWRKY29/37-silenced tea plants (AsODN) were considerably lower than those observed in control tea plants (Fig. 4h, i). This suggests that CsWRKY29 and CsWRKY37 both protect photosystem II and respond to cold stress in the tea plant.

Figure 4. 

Antisense oligodeoxynucleotides- based gene silencing of CsWRKY29 and CsWRKY37 in tender leaves of tea plants. (a) A schematic diagram of AsODN silencing of CsWRKY29 and CsWRKY37. (b) Relative expression levels of CsWRKY29 and (c) CsWRKY37 feeding with AsODN compared to control. (d) MDA contents in CsWRKY29 and (e) CsWRKY37-silenced tea plant seedlings compared to control. (f) Determination of the net photosynthetic rate and maximum photochemical efficiency of photosystem II (Fv/Fm) in CsWRKY29- and (g) CsWRKY37-silenced and control tea plants. (h), (i) Fv/Fm values in CsWRKY29/37-silenced and control tea plants.

Overexpression of CsWRKY29 and CsWRKY37 gene confers cold tolerance with improved survival rate and reduced MDA content in transgenic plants

Transgenic Arabidopsis plants were subjected to freezing stress to check the functions of CsWRKY29 and CsWRKY37 (Fig. 5a−c). Results revealed that the average survival rates of transgenic plants overexpressing CsWRKY29 and CsWRKY37 were > 76% greater than the wild type plants (48%) after the cold treatment (Fig. 5d−f). This suggested that transgenic plants overexpressing CsWRKY29 and CsWRKY37 were more tolerant to freezing stress. We then conducted a physiological test to see how CsWRKY29/37 overexpression contributes to cold tolerance. It was revealed that the MDA contents in the wild type plants spiked after 2 h of chilling stress. Nonetheless, there was significantly low MDA accumulation in the CsWRKY29 and CsWRKY37 overexpressing plants after the chilling stress than in the WT plants (Fig. 5g, h). Taken together, we showed that overexpression of both CsWRKY29 and CsWRKY37 from tea plant confers cold tolerance with improved survival rate and reduced MDA content in transgenic plants.

Figure 5. 

Overexpression of the CsWRKY29 and CsWRKY37 in Arabidopsis plants. (a) Screening and PCR validation of CsWRKY29/37 transgenic Arabidopsis plants. OE1-6: overexpression line 1−6; PC: positive control; WT: wild-type. (b, c) Expression levels of CsWRKY29/37 genes in six transgenic Arabidopsis plants. (d) The top panel shows the phenotypes of wild-type and transgenic plants growing normally for 2 weeks, while the bottom panel illustrates the phenotypes of wild-type and transgenic plants after cold treatments at −6 °C for 2 h. (e, f) Survival rates of wild-type and transgenic Arabidopsis plants (p < 0.05) after cold treatments. (g, h) MDA contents of wild-type and transgenic Arabidopsis plants (p < 0.05) after cold treatments.

Plants in their natural environment usually encounter a wide range of climatic extremes and challenging stress conditions. The evolutionary mechanism of physiological, biochemical, and developmental factors provides significant but not multiple insights to cope with such acute environmental drifts^[29]. In contrast, the molecular regulation of genes-encoding important transcription factor families in the plant may reveal the closest desirable output of many simultaneous stress responses^[30]. One of the important classes of transcription factors, the WRKY TF family, have been characterized for their significance in plant growth, phytohormone signaling, and the development of susceptibility against biotic and abiotic stress conditions^[5]. For instance, 11 OsWRKY genes were demonstrated to have a diverse response in coping with salt, drought, cold, or heat stresses in rice plants^[31]. Similarly, eight members of wheat WRKY encoding genes were preferentially expressed during their encounter with cold, heat, salt, and drought stress treatments^[32].

Tea plants are usually grown in mountainous areas, and the growth and productivity of this perineal evergreen woody crop are heavily hindered by extreme weather events, particularly cold stress^[33−35]. Although a handful of studies have shown the critical functions of the WRKY family in the model plant A. thaliana^[36], as well as in different plant species, including wheat^[32], maize^[37], banana^[38], and soybean^[39]. Yet, there is a dearth of information about WRKY genes and the putative roles of novel WRKYs under extreme weather conditions in the second most consumed beverage, the tea plant. Advancements in genome sequencing and the availability of high-quality genomes and well-annotated databases have allowed us to uncover the WRKY transcription factor family in the tea plant^[18,19,40]. Here, we have identified 86 putative WRKY transcription factor genes by screening the whole genome sequence of the tea plant. In addition, the functional roles of the two novel members of the WRKY family, CsWRKY29 and CsWRKY37, have proven to improve cold tolerance in transgenic Arabidopsis. This study provides the evolutionary information and potential functional roles of the WRKY genes identified in tea plant.

Although few studies have reported the WRKY members of the tea plant, however, these studies differ from our study in such a way that these studies are either limited to the transcript level or provide rare information about the important functional roles of the WRKY members, particularly under cold stress^{[10,11,21,22,41]}. For instance, Wu et al. used transcriptome level analysis to identify WRKY proteins in tea plants with the only prediction of their functional roles^[21]. Similarly, Wang et al. performed a systematic analysis of the WRKY gene family under stress conditions with no solid experimental validation of the functional roles under cold stress^[22]. Other notable differences between these studies and our findings are the total number of WRKYs identified, for instance, 56 and 50 genes respectively in earlier studies compared to 86 WRKYs identified in our study with additional information about whole genome duplication events in tea plant.

A zinc-finger motif-like structure and the conserved WRKY domain were used to classify all 86 potential WRKYs that encode proteins into three distinct categories. (I, II, III,) and subgroups (II-a, II-b, II-c, II-d, II-e). Further dissection of WRKY encoding proteins revealed that most members contained the most common WRKY domains, including the highly conserved WRKY heptapeptide domain (WRKYGQK); only a few members (5/86) contained the two variants heptapeptide domain, i.e., WRKYGRK, WRKYGKK. Importantly, 75/86 members of CsWRKYs showed the presence of the zinc-finger motif of the C2H2 type, and the remaining 11 proteins contained the C2HC type zinc-finger motif. The loss of zinc-finger motif C2HC motif from these candidate WRKY genes might involve functional diversity. Nonetheless, it is not yet clear whether the presence or absence of zinc-finger motif variants could modify the function and expression of WRKY genes^[42].

To further elucidate the evolutionary relationship of CsWRKY encoding genes, we conducted a comprehensive phylogenetic analysis of CsWRKYs and AtWRKYs using conserved WRKY domain sequences. Owing to the topology and structural conservation of crucial WRKY domains, a group of 158 WRKY domains, including 72 AtWRKY and 86 CsWRKY domains, were clustered into three primary groups: Group I, Group II, and Group III. Furthermore, the third group (Group II) was then classified into (subgroups II-a, II-b, II-c, II-d, and II-e). Noticeably, the highest number of CsWRKY domains were allocated to group II covering a total of 56 WRKY members, which is the largest group compared to others. In addition, the CsWRKY members in group II were also found to contain a well-known conserved WRKY domain along with a zinc-finger structure (C-X4/5-C-X22/23-HXH). Previous studies also reported the occurrence of a relatively higher proportion of WRKY members in group II in many other plant species, including Caragana intermedia (64%)^[43] and Manihot esculenta (66%)^[44]. Similarly, it was also suggested that every individual group containing CsWRKY encoding genes retained a similar gene structure; however, other group members demonstrated the presence of certain features, in line with previous findings in chickpeas^[44], and common beans^[45]. Moreover, the structural diversity and organization of Exon-intron in WRKY-encoding plays an integral role in gene family evolution, which further supports the accuracy of the phylogenetic analysis^[46]. We, therefore, extensively investigated the structural features of Exon-intron in CsWRKY genes. The findings of this study indicated that the number of introns-exons in CsWRKY genes was unevenly distributed in different members of the CsWRKY gene family. For instance, one CsWRKY had 7 exons and 6 introns, five had 6 exons and 5 introns, 14 had 5 exons and 4 introns, seven had 4 exons and 3 introns, nine had 2 exons, 34 had 3 exons and 2 introns. Notably, the distribution of introns in rice and rubber tree were also identified in the range of 1−8 and 1−7, respectively^[47,48]. On the other hand, the other members, including CsWRKY76, CsWRKY77, CsWRKY78, and CsWRKY79 contained a single exon, CsWRKY07 10 exons, and CsWRKY12 12 exons. Most of the CsWRKY genes in group III contained 3 exons, while the number of exons in group II ranged from 2 to 7 exons (CsWRKY22). Meanwhile, there was a large amount of variation in the number of exons present in group I, which ranged anywhere from 4 to 12. Overall, a similar exon-intron structure was found in members of CsWRKY detected in the same cluster during phylogenetic analysis.

The expansion of gene families in plant genomes is usually driven by gene duplication mechanisms including tandem duplication, segmental duplication, and genome redundancy^[18]. In this study, the proportion of WRKY members identified from the tea genome sequence was compared to the genomes of various representative plant species. It was discovered that the tea plant genome encodes more WRKY proteins (86) than certain model plants like Arabidopsis (72)^[25], cucumber (55)^[26], and grape (59)^[27]. Most of the CsWRKY genes (86.05%) were found in the 15 pseudo-chromosomes of the tea plant genome, whereas only 12 CsWRKYs genes were found on the unanchored scaffolds. Further investigation suggested that chromosome 6 in the tea plant genome with the most WRKY genes, which has three members from subgroup I, six members from subgroup II, and one member from subgroup III; it is followed by chromosome 1, which showed the presence of nine WRKY genes (two in subgroup I, five in subgroup II and two in subgroup III). Investigation of the tea plant genome suggested that it underwent a whole-genome duplication (WGD) 30 million years ago. This led to an increase in genes related to tea quality and stress resilience^[18,20]. Nevertheless, our results demonstrated that only a relatively smaller number of CsWRKY genes (four pairs) were duplicated during the WGD event of tea plants, suggesting less contributions of tea WGD event on the expansion of CsWRKY genes. On the other hand, the ratios of Ka to Ks of 22 duplication pairs of CsWRKY genes showed less than one ratio, indicating that the tea WRKY genes are likely to experience strongly purify selection during their evolution. Recent studies have revealed that tandem duplication is the most essential and primary source of the WRKY gene family in Brassica napus^[49], which is also in line with our findings that tandem duplication events are the primary evolution force for WRKY gene family expansion in the tea plants.

Furthermore, cloning of two CsWRKY29 and CsWRKY37 genes revealed that both contained a conserved WRKYGQK domain and C2H2 type zinc-finger motif, indicating they are the members of group II of WRKY families. In general, group II and III WRKY proteins have been shown to have a single conserved WRKY domain, which differs in the type of zinc-finger motif^[36]. Protein localization in cells is critical for defining their biological activities^[50]. Previous research has shown that WRKY transcription factors are generally found in the nucleus with transcriptional activity^[29,36]. The subcellular localization of CsWRKY29 and CsWRKY37 proteins, using fused GFP with both WRKY proteins, agroinfiltrated into N. benthamiana leaves revealed the presence of both CsWRKY29-GFP and CsWRKY37-GFP fusion proteins in the nucleus, suggesting that both CsWRKY29 and CsWRKY37 are likely to be involved in important functional roles as nucleoprotein in tea plant.

Since CsWRKY29 and CsWRKY37 genes were significantly up-regulated during cold treatment at 4 °C, especially in the later stage of treatment (72 h), both WRKY37 and WRKY29 genes were further investigated for their physiological functions in cold tolerance by suppressing their expression level in tea plant leaves using gene-specific AsODN suppression method^[51]. Given the fact that the tea plant still lacks a stable genetic transformation system, we therefore used ASODN technology, coupled with in vitro overexpression technology, to validate the functions of CsWRKY29 and CsWRKY37 in response to cold stress. As expected, the expression levels of CsWRKY29 and CsWRKY37 decreased significantly after feeding with AsODN, with higher MDA accumulation and considerably lower Fv/Fm values in CsWRKY29/37-silenced tea plant seedlings, suggesting that CsWRKY29 and CsWRKY37 both protect photosystem II and respond to cold stress in the tea plant. In an additional experiment, transgenic Arabidopsis seedlings overexpressing CsWRKY29 and CsWRKY37 were subjected to freezing stress. The role of CsWRKY29 and CsWRKY37 was apparent in terms of higher survival rates of transgenic plants compared to wild-type plants after exposure to cold treatment. Previous studies demonstrated that CsWRKY46, a WRKY transcription factor of the tea plant, enhanced the survival rate of transgenic overexpressing Arabidopsis plants under chilling stress. MDA accumulation in the CsWRKY29 and CsWRKY37 overexpressing plants were also found significantly low, while in the WT plants, it spiked after 2 h of chilling stress. Low levels of MDA content have been linked to the degree of reactive oxygen species damage to cell membranes^[52]. We have shown that both CsWRKY29 and CsWRKY37 transcription factors are vital in regulating the cold response, either from a phenotypic survival or physiology and biochemistry. Overall, the overexpression of CsWRKY29 and CsWRKY37 in Arabidopsis confirm their functional roles in alleviating cold stress in tea plants.

In conclusion, we have performed a genome-wide investigation of tea WRKY transcription factor genes and examined their evolution and expression patterns in different organs and under abiotic stresses. In particular, the functional roles of two CsWRKY genes (CsWRKY29 and CsWRKY37) in cold tolerance were fully validated using both in vivo gene slicing and overexpression in transgenic Arabidopsis. The results reported in this study provides a strong foundation for studying the function and regulatory mechanisms of the WRKY TFs family in tea plants, which would help future hardiness-breeding programs.

To evaluate the functional roles of WRKY members, C. sinensis var. sinensis 'Shuchazao', a commercial tea cultivar of the region, whose genome was recently made accessible to the public was utilized as the plant material^[19,20]. Tea plants were grown in an open field at the Tea Garden in Anhui Agricultural University (Anhui, China). Different tea plant organs including, buds, flower, fruit, leaf (1^st, 2^nd, and 3^rd), stem, and root were collected and immediately frozen in liquid nitrogen and stored at −80 °C. Arabidopsis seeds were grown in a growth chamber under a controlled environment maintaining at 16 h day (25 ± 3 °C) and 8 h night (20 ± 3 °C) with 70% relative humidity. Two-year-old tea plant seedlings developed from stem cuttings were used to perform tissue-specific gene expression analysis for low-temperature stress treatments. While some of the plants were put into an incubator at a temperature of 4 °C, the other plants served as a control group. The tender leaves of the control plants, as well as the treated plants, were harvested at different times (0, 12, 48, and 72 h), promptly frozen in liquid nitrogen, and then stored at −80 °C for the purpose of RNA extraction. The Nicotiana benthamiana plants that were used in this study were grown in pots inside a growth chamber with a nighttime temperature of 20 °C and a daytime temperature of 25 °C, and the photoperiod was set at 16 h light and 8 h dark.

Identification, characterization, and evolutionary analysis of tea WRKY genes

To identify WRKY genes from the tea plants, the whole genome sequences, coding sequences, protein sequences, as well as functional annotations of C. sinensis var. sinensis (cultivar: 'Shuchazao') were extracted from the tea plant information archive (TPIA; http://tpia.teaplant.org/)^[40]. Furthermore, A. thaliana WRKY transcription factors protein sequences were retrieved from the Arabidopsis genome database as queries versus the tea plant genome. The iTAK software was implied to identify the sequences that contain the WRKY domain from the tea plant protein database^[23]. Moreover, the online tool NCBI CD-search^[53] and EMBL SMART^[54] were used to confirm all the potential CsWRKY genes, at this stage, the genes that were incorrectly predicted were discarded, and ultimately the non-redundant and only confirmed genes were obtained as putative WRKYs in tea plant.

The full-length WRKY protein sequences of the tea plant were multiple-aligned with A. thaliana WRKY transcription factors using MUSCL with default parameters^[55]. The alignment was further trimmed by trimAl with the heuristic method^[56]. To further elucidate the phylogenetic relationship of CsWRKYs and aid their classification, the phylogenetic tree of the tea plant and Arabidopsis WRKY gene was established based on the conserved WRKY domain sequences following the method of neighbor-joining with the help of MEGA software^[24]. The stability of internal nodes was assessed with a bootstrap value of 1,000 interactions. The phylogenetic tree was viewed using an online program Evolview^[57]. The identified genes were finally named 'CsWRKY' with continuous numbers and were divided into different groups and subgroups based on the phylogenetic tree constructed.

The structural information of the CsWRKY genes was extracted from the tea plant gene annotation file, and the diagram of the exon-intron structure was generated using the web-based program Gene Structure Display Server (GSDS) to compare the respective full-length sequence’s structure^[58]. In addition, all non-redundant candidate tea WRKY proteins were also assessed to confirm whether they contained conserved motifs using the online program Multiple Expectation maximization for Motif Elicitation, with the maximum number of motifs set to 10. In contrast, the other parameters were set to default values^[59].

Expression analysis of CsWRKY genes among different tissues and abiotic stresses

To examine the effects of three abiotic stress applications on the expression of the CsWRKY genes, as well as the spatial and temporal patterns of expression in eight different tissues during development, the raw sequences from transcriptome data sets of tea plant root (RT), stem (ST), bud (AB), young leaf (YL), mature leaf (ML), old leaf (OL), flower (FL) and fruit (FR) (SRP056466)^[20], as well as under drought (PRJEB11522)^[60], cold (PRJNA387105)^[61], and methyl jasmonate (MeJA) treatments (PRJNA288922)^[62] were downloaded from the NCBI database. The raw sequencing datasets were pre-analyzed for exemption of the low quality bases (phred score < 20) and short length reads (length < 25 bp) by SolexaQA package^[63]. The remaining high-quality data obtained after pretreatment was then aligned against the tea plant reference genome using TopHat to calculate gene expression levels (FPKM value)^[64]. The Edger software tool in R was used to compare the expression patterns of CsWRKYs during cold stress^[65]. The differentially expressed genes were then screened. The number of up-regulated and down-regulated genes across the two treatment periods was counted using the Edger method with fold change ≥ 2 and adjusted p-value ≤ 0.01. The expression patterns of CsWRKY genes among different tissues and abiotic stress were visualized using pheatmap method in the R program.

Isolation and bioinformatics analysis of CsWRKY29 and CsWRKY37

The sequences of two WRKY genes, namely CsWRKY29, and CsWRKY37 were downloaded from TPIA database^[40]. The total RNA was obtained from the leaf tissues of the tea plant and then the first-strand cDNA templates were prepared using a PrimeScript RT reagent Kit (Takara Biotech, China). Afterwards, gene specific primers were designed for cloning using Premier3 software^[66], and then synthesized by General biological system (Anhui) Co., Ltd (China). The PCR product was purified from agarose gels with AxyPrep DNA Gel Extraction Kit (Axygen, USA) and ligated into pGEM-T Easy vector (Promega, China) after confirming the gene sequence by sequencing (General biological system (Anhui) Co., Ltd, China). The amino acid length, molecular weight, GC content, and isoelectric point of cloned genes were predicted using DNASTAR (www.dnastar.com). The multiple sequence alignments among CsWRKY29, CsWRKY37, and AtWRKY6 proteins were performed using MUSCLE^[55]. The open reading frame of CsWRKY29 and CsWRKY37 genes were predicted using NCBI Open Reading Frame Finder (www.ncbi.nlm.nih.gov/orffinder/).

Quantitative RT-PCR assays

Total RNA was isolated from the different tea plant organs including bud, flower, fruit, 1^st, 2^nd, 3^rd leaf, stem, and root. Quantitative real-time PCR (qRT-PCR) reaction was carried out using 2 μL of cDNA and 0.4 μM of each primer in a 20 μL reaction volume using SYBR Premix Ex Taq™ II (Takara, Japan). The qRT-PCR reaction program was performed under the following conditions: 95 °C for 5 min, 40 cycles at 95 °C for 10 s, and 60 °C for 30 s. The Actin (NCBI accession: HQ420251) was used as an internal control. The specific primers of Actin and selected CsWRKY genes used for qRT-PCR analysis were presented in Supplemental Table S5. Three replicates of each PCR run were performed. According to the threshold cycle (Ct), the relative changes in gene expression were quantified using the previously reported 2^−ΔΔCᴛ method^[67].

Protein subcellular localization studies

The primers of plant vectors were designed based on the Gateway method (Supplemental Table S6). The full-length sequences of CsWRKY29 and CsWRKY37 were ligated into pK7WGF2.0 vector without terminator codon, containing the green fluorescent protein (GFP) reporter gene, to generate 35s::CsWRKY-GFP fusion constructs. Then, the recombinant plasmid (35s::CsWRKY-GFP) and the control (35s::GFP) were transformed into Agrobacterium tumefaciens EHA105 using the electroporation method. The cells of A. tumefaciens were cultured in LB medium containing spectinomycin (100 μg/mL) and rifampicin (20 μg/mL) at 28 °C for 16−20 h until the optical density at ƛ = 600 nm (OD600) was 1.0. This medium added with bacteria was then centrifuged at 4000 rpm for 10 min. The sediment was resuspended with suspension buffer containing 10 mM MgCl₂, 10 mM MES, and 100 μM acetosyringone (AS), with OD600 adjusted to 0.6. The suspension was injected into well-developed N. benthamiana leaves. After infiltration, the tobacco plants were cultured in darkness for 8−12 h at 24 °C and then put in light for 2−3 d. The GFP signals were observed by a confocal laser scanning microscope (OLYMPUS, Japan).

Gene silencing

Candidate antisense oligonucleotides (AsODN) were selected using SOLIGO software with CsWRKY29 and CsWRKY37 as input sequence^[68] (Supplemental Table S7). Synthesized AsODNs were purchased from General Biosystems Company. To silence CsWRKY29 and CsWRKY37 in the tea leaves, isolated plants with one bud and two leaves were soaked in 1 ml of 20μM AsODN CsWRKY29, and CsWRKY37 solution, while the sense oligonucleotides (sODN) were used as a control^[51]. After 6 h, 12 h and 24 h of incubation, the leaves were harvested and kept at −80 °C before analysis. According to the best silencing time point, the silencing treatment was followed by a cold treatment of 0 °C for 1 h, and the darkness was restored for 30 min at a room temperature of 25 °C. After the recovery period, both net photosynthetic rate and maximum photochemical efficiency of photosystem II (Fv/Fm) were measured. We also determined the content of malondialdehyde (MDA) according to the kit (SolarBio, Beijing, China), with no less than three repeats of low-temperature treatment.

Freezing tolerance assays

Following the previously reported protocol^[69], the freezing assay was performed with minor modifications. Briefly, Arabidopsis seeds were grown on Petri dishes containing half-strength MS medium at 22 °C (day 16 h/night 8 h) with 70% humidity for 2 weeks. For freezing stress, Petri dishes containing plants were opened and exposed to a freezing chamber (RuMED4001), which was set to 4 °C for 2 to 3 h to evaporate the water from the surface of Petri dishes. The freezing chamber was then set at −1 °C and programmed to cool at 1 °C after every hour. After being exposed to −6 °C for 2 h, the Petri dishes containing plants were removed from the chamber. After the freezing procedure, the plants were incubated at 4 °C in the dark for 12 h before being transported to 22 °C in the light. After 4 d, the seedlings' survival rate was assessed visually^[70].

Functional analysis of CsWRKY29 and CsWRKY37 under cold tolerance

The ORFs of the CsWRKY29 and CsWRKY37 genes were cloned into the Pb2GW7 vector under the control of the CaMV 35S promoter. The constructs were then transformed into the Agrobacterium strain GV3101 using the freeze-thaw method, followed by transforming the Arabidopsis plants via the floral dip method. Arabidopsis ecotype Columbia-0 (col-0) genetic background was used as the wild type in the study. Six lines of transformed plants were selected based on their resistance to Basta. Genomic DNA was extracted and verified through PCR amplification using 35S and gene-specific primers (Supplemental Table S6). We further quantitatively analyzed the transcriptional levels in six transgenic lines overexpressing CsWRKY29 and CsWRKY37, the internal control was Atactin2 (F: 5'- TTCGTTTTGCGTTTTAGTCCC -3', R: 5'- GGGAACAAAAGGAATAAAGAGGC -3'). CsWRKY29 was more strongly expressed in the overexpressing line OE-1 and OE-5, while CsWRKY37 was more highly expressed in the OE-2 and OE-3. Therefore, we selected these overexpressing lines for further analysis. Transgenic lines were then exposed to −6 °C for 2 h and then grown under 25 °C conditions for 4 d, the survival rate was recorded. Malondialdehyde content was determined using a similar method used in a previous study^[71] using kit (SolarBio, Beijing, China).

Data analysis

All experiments were performed in triplicate, the numerical data generated from three experimental repeats were subsequently subjected to statistical analyses using Microsoft-EXCEL 2010 and R software. The statistical significance between groups was further analyzed using student’s t-test.

The coding sequences, protein sequences and expression data of all 86 WRKY genes are available in tea plant information archive (TPIA: http://tpia.teaplant.org).