Genome-Wide Identiﬁcation and Expression Analysis of WRKY Transcription Factors in Siraitia siamensis

: WRKY transcription factors, as the largest gene family in higher plants, play an important role in various biological processes including growth and development, regulation of secondary metabolites, and stress response. In this study, we performed genome-wide identiﬁcation and analysis of WRKY transcription factors in S. siamensis . A total of 59 SsWRKY genes were identiﬁed that were distributed on all 14 chromosomes, and these were classiﬁed into three major groups based on phylogenetic relationships. Each of these groups had similar conserved motifs and gene structures. We compared all the S. siamensis SsWRKY genes with WRKY genes identiﬁed from three diverse plant species, and the results implied that segmental duplication and tandem duplication play an important roles in the evolution processes of the WRKY gene family. Promoter region analysis revealed that SsWRKY genes included many cis -acting elements related to plant growth and development, phytohormone response, and both abiotic and biotic stress. Expression proﬁles originating from the transcriptome database showed expression patterns of these SsWRKY genes in four different tissues and revealed that most genes are expressed in plant roots. Fifteen SsWRKY genes with low-temperature response motifs were surveyed for their gene expression under cold stress, showing that most genes displayed continuous up-regulation during cold treatment. Our study provides a foundation for further study on the function and regulatory mechanism of the SsWRKY gene family.


Introduction
WRKY transcription factors are one of the largest families of transcriptional regulators in plants and are integral to signaling webs that modulate many important plant processes [1]. WRKY refers to a highly conserved protein domain, the WRKY domain, with nearly 60 amino acids [2]. The most significant feature of these proteins is that they all contain a WRKYGQK sequence in N-terminal and a C 2 H 2 or C 2 HC (C-X 4-5 -C-X 22-23 -H-X 1 -H or C-X 7 -C-X 23 -H-X 1 -H) zinc finger motif at the C-terminal [3,4]. Previous studies have shown that several variants of the WRKYGQK motif exist, including WRKYGEK and WRKYGKK. A typical WRKY motif has therefore been suggested, embracing the consensus sequence W(R/K)(K/R)Y [5]. According to the number of WRKY domains and the type of

Identification and Characterization of SsWRKY in S. siamensis
To identify potential WRKY genes, the S. siamensis genome database (https://dataview. ncbi.nlm.nih.gov/object/SRR22947134?reviewer=qmfs7mc075pohiv011jqjceh94, accessed on 2 January 2023) was screened using HMM3.0, and a total of 61 candidates were identified. Of these, two candidate SsWRKY with incomplete WRKY domains were not considered to be SsWRKY. These were classified as SsWRKY-like genes, and were removed from our further analysis. The remaining 59 SsWRKY genes were named SsWRKY1-SsWRKY59 according to their position on chromosome of S. siamensis. The detailed characterization of 59 SsWRKY genes is listed in Table 1. The protein length ranged from 128 (SsWRKY14) to 745 (SsWRKY39) amino acids, and the corresponding predicted molecular weights varied from 14.59 kDa to 81.91 kDa, with an average of 40.46 kDa. The theoretical pI ranged from 4.83 (SsWRKY51) to 9.92 (SsWRKY23). Subcellular localization prediction assigned all 59 proteins to the nucleus. The great majority of proteins were considered unstable, except SsWRKY3, SsWRKY10, SsWRKY13, SsWRKY22, SsWRKY23, and SsWRKY50 proteins.

Classification and Phylogenetic Relationships of SsWRKYs
Combining 70 AtWRKY protein sequences with 59 SsWRKY protein sequences, MEGA7.0 software was employed to construct the phylogenetic tree and classify the sequences into three groups based on the features of the WRKY gene family in A. thaliana ( Figure 1). Eleven SsWRKY proteins were classified into Group I, forty SsWRKY proteins were classified into Group II and were further divided into five subgroups, with 4, 5, 16, 7, and 8 members in subgroups IIa-IIe, respectively. Group III consisted of the remaining eight SsWRKY proteins.
Multiple sequence alignments were used to identify the structural characteristics of SsWRKY proteins. As is shown in Figure 2, the majority of the 59 SsWRKY proteins (48/59, 81%) had a single conserved WRKY domain, and the remaining 11 proteins (19%) contained two WRKY domains. The highly conserved WRKY motif represented by 58 WRKYGQK and only one WRKYGKK (SsWRKY10), was identified in 59 SsWRKY proteins, similar to results of studies from other species including Isatis indigotica [37], O. pumila [5], and tomato [38]. All Group I members contained two complete WRKY domains, including a WRKYGQK sequence and a C 2 H 2 -like zinc-finger motif. All 40 Group II proteins contained one WRKY domain and a C 2 H 2 motif (C-X 4-5 -C-X 22-23 -H-X 1 -H). Eight members of Group III had a WRKYGQK sequence and the C 2 HC zinc finger motif (C-X 7 -C-X 23 -H-X 1 -H). Multiple sequence alignments were used to identify the structural characteristics of SsWRKY proteins. As is shown in Figure 2, the majority of the 59 SsWRKY proteins (48/59, 81%) had a single conserved WRKY domain, and the remaining 11 proteins (19%) contained two WRKY domains. The highly conserved WRKY motif represented by 58 WRKYGQK and only one WRKYGKK (SsWRKY10), was identified in 59 SsWRKY proteins, similar to results of studies from other species including Isatis indigotica [37], O. pumila [5], and tomato [38]. All Group I members contained two complete WRKY domains, including a WRKYGQK sequence and a C2H2-like zinc-finger motif. All 40 Group II proteins contained one WRKY domain and a C2H2 motif (C-X4-5-C-X22-23-H-X1-H). Eight members of Group III had a WRKYGQK sequence and the C2HC zinc finger motif (C-X7-C-X23-H-X1-H). thaliana WRKY proteins were aligned by ClustW, and the phylogenetic tree was constructed using the neighbor-joining method and MEGA 7.0 with 1000 bootstrap replicates. The WRKY proteins from different groups are indicated with different colors, proteins from A. thaliana and S. siamensis are denoted by red circles and blue squares.

Gene Structure and Conserved Motif Analysis
To better illustrate the similarity and diversity of motif composition among SsWRKYs, the MEME program was employed to analyze their conserved motifs. A total of ten conserved motifs were identified in the WRKY proteins of S. siamensis, ranging from 15 to 38 amino acids (Figure 3d), as shown in Figure 3 and Table S1. Motifs 1 and 2 were extensively distributed in all SsWRKY proteins, corresponding to the conserved WRKY domain. Motifs 3, 5, and 10 were unique to Group I. Motif 6 was found mainly in Groups IIa and IIb; Motifs 7 and 8 were mainly distributed in Group IIe, Group IId, and Group III; Motif 9 was found only in Groups IIa and IIb (Figure 3b). In general, almost all the SsWRKYs that belonged to the same subgroups displayed similar motif composition, implying that these SsWRKY proteins have analogous functions.  To help clarify the evolution of SsWRKY family genes and their genetic structural diversity, the exon-intron organization of the coding sequences of SsWRKY genes was examined. As is shown in Figure 3c, the number of the exons in SsWRKYs varied from two to seven exons (7 with two exons, 29 with three exons, 8 with four exons, 10 with five exons, 4 with six exons), with SsWRKY21 containing the highest number of predicted exons (n = 7). These results show that the genetic structural diversity among SsWRKY genes may be the result of exon-loss and -gain events during the evolution of the SsWRKY gene family. SsWRKY genes from the same group showed similar exon/intron structures, as in Group IIe. In addition, SsWRKY genes with similar structures clustered together, as in Group III and Group IIe which contained three exons and two introns (or three introns) except for SsWRKY24, SsWRKY11, and SsWRKY36.

Gene Structure and Conserved Motif Analysis
To better illustrate the similarity and diversity of motif composition among SsWRKYs, the MEME program was employed to analyze their conserved motifs. A total of ten conserved motifs were identified in the WRKY proteins of S. siamensis, ranging from 15 to 38 amino acids (Figure 3d), as shown in Figure 3 and Table S1. Motifs 1 and 2 were extensively distributed in all SsWRKY proteins, corresponding to the conserved WRKY domain. Motifs 3, 5, and 10 were unique to Group I. Motif 6 was found mainly in Groups IIa and IIb; Motifs 7 and 8 were mainly distributed in Group IIe, Group IId, and Group III; Motif 9 was found only in Groups IIa and IIb (Figure 3b). In general, almost all the SsWRKYs that belonged to the same subgroups displayed similar motif composition, implying that these SsWRKY proteins have analogous functions.  To help clarify the evolution of SsWRKY family genes and their genetic structural diversity, the exon-intron organization of the coding sequences of SsWRKY genes was examined. As is shown in Figure 3c, the number of the exons in SsWRKYs varied from two to seven exons (7 with two exons, 29 with three exons, 8 with four exons, 10 with five exons, 4 with six exons), with SsWRKY21 containing the highest number of predicted exons (n = 7). These results show that the genetic structural diversity among SsWRKY genes may be the result of exon-loss and -gain events during the evolution of the SsWRKY gene family. SsWRKY genes from the same group showed similar exon/intron structures, as in Group IIe. In addition, SsWRKY genes with similar structures clustered together, as in Group III and Group IIe which contained three exons and two introns (or three introns) except for SsWRKY24, SsWRKY11, and SsWRKY36.

Chromosomal Location and Collinearity Analysis of SsWRKY Genes
As displayed in Figure 4, all 59 SsWRKY genes were mapped to 14 chromosomes of S. siamensis. Remarkably, although WRKY genes were not distributed evenly on every chromosome, they were present on each chromosome. There was no apparent relationship between the number of SsWRKY genes and chromosome length ( Figure 4). SsWRKY genes were most abundant on Chr 3, with eight genes, and least abundant on Chr 12, 13, and 14, each with two genes. Collinear blocks within the S. siamensis genome were investigated to identify the relationship between SsWRKY genes and gene duplication events ( Figure 5a). Twenty-four gene pairs were detected in the S. siamensis genome and were located on different chromosomes, suggesting that segmental duplications in these locations potentially contributed to expansion of the SsWRKY family. In addition, only one tandem duplication gene pair was identified, SsWRKY8 and SsWRKY9 located on Chr 2

Chromosomal Location and Collinearity Analysis of SsWRKY Genes
As displayed in Figure 4, all 59 SsWRKY genes were mapped to 14 chromosomes of S. siamensis. Remarkably, although WRKY genes were not distributed evenly on every chromosome, they were present on each chromosome. There was no apparent relationship between the number of SsWRKY genes and chromosome length ( Figure 4). SsWRKY genes were most abundant on Chr 3, with eight genes, and least abundant on Chr 12, 13, and 14, each with two genes. Collinear blocks within the S. siamensis genome were investigated to identify the relationship between SsWRKY genes and gene duplication events (Figure 5a). Twenty-four gene pairs were detected in the S. siamensis genome and were located on different chromosomes, suggesting that segmental duplications in these locations potentially contributed to expansion of the SsWRKY family. In addition, only one tandem duplication gene pair was identified, SsWRKY8 and SsWRKY9 located on Chr 2 ( Figure 4). Some reports have indicated that tandem gene duplication is an important reason for gene clustering [39]. To enquire into the phylogenetic mechanisms of the SsWRKY gene family, three collinear maps of S. siamensis were constructed using three representative species, one monocotyledon (O. sativa) and two dicotyledons (A. thaliana and C. sativus) (Figure 5b-e). A total of 56 SsWRKY genes showed a collinear relationship with C. sativus (n = 56), followed in terms of frequency by A. thaliana (n = 50), and O. sativa (n = 20) (Table S2). There were 85, 72, and 29 orthologous pairs between the other three species, respectively (C. sativus, A. thaliana, and O. sativa). Syntenic pairs were identified with three other species (with 13 SsWRKY genes, SsWRKY7, 8,11,12,16,26,36,38,40,45,47,49), implying that these gene pairs may have been present before ancestral divergence and played a key role in the evolution of the SsWRKY gene family. sativus (n = 56), followed in terms of frequency by A. thaliana (n = 50), and O. sativa (n = 20) (Table S2). There were 85, 72, and 29 orthologous pairs between the other three species, respectively (C. sativus, A. thaliana, and O. sativa). Syntenic pairs were identified with three other species (with 13 SsWRKY genes, SsWRKY7, 8,11,12,16,26,36,38,40,45,47,49), implying that these gene pairs may have been present before ancestral divergence and played a key role in the evolution of the SsWRKY gene family.

Analysis of Cis-Acting Elements of SsWRKY Gene Family
To further our understanding of the potential function and regulation of WRKY transcription factors in S. siamensis, nearly 2.0 kb of sequence regions upstream of the translation initiation site was extracted to identify the cis-acting elements, using PlantCARE software. Almost all SsWRKY genes had at least one cis-acting element in the promoter region. These cis-acting elements were further classified into three categories, plant growth and development, phytohormone-responsive, and abiotic and biotic stress ( Figure 6). The first category, plant growth and development, included light-responsive and developmentrelated elements. Various light-responsive elements were presented in the promoter region of SsWRKY genes, containing Box4, G-box, GT1-motif, GATA-motif, Sp1, etc. The Box 4 element (203) was the most abundant cis-acting element found in SsWRKY promoters. There were also many elements related to development, including CAT-box, MSA-like, and O 2site. The second category contained phytohormone-responsive elements, composed of the TGACG-motif, TCA-element, ABRE, GARE-motif, etc. For example, the promoter region of SsWRKY1 had ABRE, AuxRR-core, CGTCA-motif, GARE-motif, and TGACG-motif, indicating that these cis-acting elements are related with ABA, auxin, methyl jasmonate (MeJA), gibberellin, and salicylic acid responsiveness. It is speculated that SsWRKY1 may therefore be involved in multiple hormone responses. The third category included cisacting elements related to abiotic and biotic stress, including antioxidant response element (ARE), low temperature response (LTR) element, drought inducibility elements (MBS), wound responsive elements (Wun-motif), and TC-rich repeats. In S. siamensis, 139 ARE elements were present in 53 promoters (90%) of SsWRKY genes. In contrast, Wun-motif was the least abundant element in S. siamensis. LTR elements were found in 15 promoters of SsWRKY. Based on the strong induction in response to cold-stress treatment, SsWRKY44 and SsWRKY56 had four copies of the LTR element. Remarkably, the SsWRKY37 gene did not contain any elements related to abiotic or biotic stress. In general, these results demonstrate that the SsWRKY genes might be associated with growth and development, hormones, and stress responses. The analysis of the putative cis-acting elements of the SsWRKY gene family might improve understanding of stress response, especially low temperature stress in S. siamensis.

Analysis of Cis-Acting Elements of SsWRKY Gene Family
To further our understanding of the potential function and regulation of WRKY transcription factors in S. siamensis, nearly 2.0 kb of sequence regions upstream of the translation initiation site was extracted to identify the cis-acting elements, using PlantCARE software. Almost all SsWRKY genes had at least one cis-acting element in the promoter region. These cis-acting elements were further classified into three categories, plant growth and development, phytohormone-responsive, and abiotic and biotic stress ( Figure 6). The first category, plant growth and development, included light-responsive and developmentrelated elements. Various light-responsive elements were presented in the promoter region of SsWRKY genes, containing Box4, G-box, GT1-motif, GATA-motif, Sp1, etc. The Box 4 element (203) was the most abundant cis-acting element found in SsWRKY promoters. There were also many elements related to development, including CAT-box, MSAlike, and O2-site. The second category contained phytohormone-responsive elements, composed of the TGACG-motif, TCA-element, ABRE, GARE-motif, etc. For example, the promoter region of SsWRKY1 had ABRE, AuxRR-core, CGTCA-motif, GARE-motif, and TGACG-motif, indicating that these cis-acting elements are related with ABA, auxin, methyl jasmonate (MeJA), gibberellin, and salicylic acid responsiveness. It is speculated that SsWRKY1 may therefore be involved in multiple hormone responses. The third category included cis-acting elements related to abiotic and biotic stress, including antioxidant response element (ARE), low temperature response (LTR) element, drought inducibility elements (MBS), wound responsive elements (Wun-motif), and TC-rich repeats. In S. siamensis, 139 ARE elements were present in 53 promoters (90%) of SsWRKY genes. In contrast, Wun-motif was the least abundant element in S. siamensis. LTR elements were found in 15 promoters of SsWRKY. Based on the strong induction in response to cold-stress treatment, SsWRKY44 and SsWRKY56 had four copies of the LTR element. Remarkably, the SsWRKY37 gene did not contain any elements related to abiotic or biotic stress. In general, these results demonstrate that the SsWRKY genes might be associated with growth and development, hormones, and stress responses. The analysis of the putative cis-acting elements of the SsWRKY gene family might improve understanding of stress response, especially low temperature stress in S. siamensis.

Protein Interaction Network and Gene Co-Expression of SsWRKYs
String software was used to construct a protein-protein interaction network of 18 SsWRKY proteins based on the A. thaliana association model. In the interaction network diagram (Figure 7a and Table S2), connected genes might have close functional relationships. The thicker the connecting line, the stronger is the predicted interaction between the two proteins. Among the 18 SsWRKY proteins, there were three proteins from Group I, eleven proteins belonging to Group II, and four proteins in Group III. In addition, SsWRKY53, SsWRKY54, SsWRKY11, SsWRKY12, SsWRKY38, SsWRKY47, SsWRKY3, and SsWRKY39 were shown to have stronger predicted interaction networks with other proteins. These proteins might play a key role in the whole process of WRKY regulation. The results of SsWRKYs co-expression showed that most genes were co-expressed, with the red color representing higher co-expression levels. As shown in Figure 7b, the cluster members SsWRKY4, 5, 27, 28, 41, 42, 44, 45, 46, 47 had high co-expression levels. Among them, the co-expression levels between SsWRKY44-47 and SsWRKY4-5, SsWRKY27-28, SsWRKY41-42 were significantly higher than those among other members.

Expression Profiles of the SsWRKY Genes in Different Tissues
The expression levels of the 59 SsWRKY genes were investigated in four different tissues of S. siamensis (roots, stems, tubers, and leaves), using the S. siamensis transcriptome database (https://dataview.ncbi.nlm.nih.gov/object/PRJNA917212?reviewer=qmfs7mc0 75pohiv011jqjceh94, accessed on 2 January 2023). All 59 SsWRKY genes were found to have significant expression in at least one tissue, and changes in expression were found between the different tissues surveyed, except for SsWRKY12 (Figure 8). The results showed that SsWRKY genes had the highest expression levels in the roots; 37% (22/59), 31% (18/59), 29% (17/59), and 2% (1/59) of the SsWRKY genes displayed their highest expression levels in roots, tubers, stems, or leaves, respectively. These differences in the expression profiles of the WRKY genes indicate potentially different roles in these tissues. Plants 2023, 12, x FOR PEER REVIEW 12 of 24 Figure 6. Cis-acting elements analysis in the promoter region of SsWRKY genes. The number of cisacting elements are shown by the digits. The greater the quantity, the darker the color. Light blue, yellow, and green represent three categories of cis-acting elements related to plant growth and development, phytohormone response, and abiotic and biotic stress.

Protein Interaction Network and Gene Co-expression of SsWRKYs
String software was used to construct a protein-protein interaction network of 18 SsWRKY proteins based on the A. thaliana association model. In the interaction network diagram (Figure 7a and Table S2), connected genes might have close functional relationships. The thicker the connecting line, the stronger is the predicted interaction between the two proteins. Among the 18 SsWRKY proteins, there were three proteins from Group I, eleven proteins belonging to Group II, and four proteins in Group III. In addition, SsWRKY53, SsWRKY54, SsWRKY11, SsWRKY12, SsWRKY38, SsWRKY47, SsWRKY3, and SsWRKY39 were shown to have stronger predicted interaction networks with other proteins. These proteins might play a key role in the whole process of WRKY regulation. The results of SsWRKYs co-expression showed that most genes were co-expressed, with the red color representing higher co-expression levels. As shown in Figure 7b, the cluster members SsWRKY4, 5,27,28,41,42,44,45,46,47 had high co-expression levels. Among them, the co-expression levels between SsWRKY44-47 and SsWRKY4-5, SsWRKY27-28, SsWRKY41-42 were significantly higher than those among other members.

Expression Profiles of the SsWRKY Genes in Different Tissues
The expression levels of the 59 SsWRKY genes were investigated in four different tissues of S. siamensis (roots, stems, tubers, and leaves), using the S. siamensis transcriptome database (https://dataview.ncbi.nlm.nih.gov/object/PRJNA917212?re-viewer=qmfs7mc075pohiv011jqjceh94 , accessed on 2 January 2023). All 59 SsWRKY genes The color variation from blue to red represents low to high expression. Our findings will provide references for the function and regulation mechanism of the WRKY gene family in other plants.

Expression Analysis of SsWRKY Genes under Cold Stress
In order to examine the expression patterns of SsWRKY genes potentially associated with responses to low temperatures, 15 SsWRKY genes containing LTR cis-acting elements in their promoter were selected and surveyed for their expression levels during different stages of induced low-temperature stress (4 • C) (0, 6, 12, 24, 36 and 48 h). The expression of the selected SsWRKY genes could be classified into three groups according to their expression patterns, as shown in Figure 9. The highest expression levels in the majority of selected SsWRKY genes (SsWRKY12, SsWRKY13, SsWRKY16, SsWRKY20, SsWRKY29, SsWRKY30, SsWRKY47, SsWRKY49, SsWRKY57) were found after exposure to low temperature for 48 h. SsWRKY44, SsWRKY51, and SsWRKY56 showed their highest levels at 6 h, and then gradually decreased, returning to basal expression levels after 48 h cold treatment, whereas SsWRKY41, SsWRKY42, SsWRKY53 were highly expressed at 24 h, falling again to initial expression levels after 48 h cold treatment. Figure 8. Expression patterns of SsWRKY genes in four different tissues of S. siamensis. RPKM value was assessed to investigate the expression levels of SsWRKY genes in stems, roots, tubers, and leaves. The color variation from blue to red represents low to high expression. Our findings will provide references for the function and regulation mechanism of the WRKY gene family in other plants.

Expression Analysis of SsWRKY Genes under Cold Stress
In order to examine the expression patterns of SsWRKY genes potentially associated with responses to low temperatures, 15 SsWRKY genes containing LTR cis-acting elements in their promoter were selected and surveyed for their expression levels during different stages of induced low-temperature stress (4 °C) (0, 6, 12, 24, 36 and 48 h). The expression of the selected SsWRKY genes could be classified into three groups according to their expression patterns, as shown in Figure 9. The highest expression levels in the majority of selected SsWRKY genes (SsWRKY12, SsWRKY13, SsWRKY16, SsWRKY20, SsWRKY29,  SsWRKY30, SsWRKY47, SsWRKY49, SsWRKY57) were found after exposure to low temperature for 48 h. SsWRKY44, SsWRKY51, and SsWRKY56 showed their highest levels at 6 h, and then gradually decreased, returning to basal expression levels after 48 h cold treatment, whereas SsWRKY41, SsWRKY42, SsWRKY53 were highly expressed at 24 h, falling again to initial expression levels after 48 h cold treatment.

Discussion
S. siamensis is important for its economic and medicinal value [31]. However, little research has been conducted at the molecular level due to a lack of genome and transcriptome data. Studies in many plant species have shown that the WRKY gene family is the Figure 9. qRT-PCR analysis of SsWRKY genes under cold stress. The Y-axes represent the relative expression and X-axes represent the duration of the cold treatment. Data represents the mean ± SD of three technical repetitions. * p < 0.05, ** p < 0.01.

Discussion
S. siamensis is important for its economic and medicinal value [31]. However, little research has been conducted at the molecular level due to a lack of genome and transcriptome data. Studies in many plant species have shown that the WRKY gene family is the largest transcription factor family, and members often participate in important plant processes, especially stress response [40]. Based on whole genome sequences, the structure and function of the WRKY gene family has been widely investigated in several species including A. thaliana [41], O. sativa [42], C. sativus [26], and Cucumis melon [27]. However, there have been no equivalent studies in S. siamensis. Genome sequencing of S. siamensis has been recently completed (unpublished), providing a means for performing genome-wide analysis. In this study, we provide identification and analysis of the SsWRKY gene family for the first time.
Based on the number of the conserved WRKY domains and the feature of the zinc finger motif, the WRKY gene family was classified into three groups. Following phylogenetic analysis with AtWRKY genes, the 59 SsWRKYs were divided into three main groups, 11 SsWRKY proteins belonging to Group I, 40 SsWRKY proteins classified into Group II, accounting for the largest proportion (68%), indicating that this group may have experienced gene duplication events during its evolution, and the remaining eight SsWRKY proteins belonged to Group III. Similar results have been reported from Zea mays [43], C. sativus [26], Coffea canephora [44]. Furthermore, variation of the conserved WRKY domain (WRKYGKK) was observed in the SsWRKY10 protein sequence from the subgroup IIc. A similar situation has also been described in other plants such as O. pumila (OpWRKY16, OpWRKY27, OpWRKY28, OpWRKY30, and OpWRKY47), and C. sativus (CsWRKY10 and CsWRKY47) [5,26]. Previous research has shown that the WRKYGQK domain can bind to W-box cis-elements in the promoter region to interact with target genes [22], and changes of the WRKYGQK domain may influence their binding specificity. Therefore, the function and DNA-binding specificities of the SsWRKY10 protein should be further analyzed. Based on the analysis of conserved motifs, WRKY members in the same group or subgroup were found to have similar motif composition. Almost all SsWRKY contained motif 1 and motif 2, which may be a core element retained during evolution. The conserved motifs of IId and IIe subgroups were similar, which may suggest their genetic relationship through evolution.
Gene duplication events, including tandem, segmental, and whole-genome duplications, are important for acquiring new genes and enabling expansion on gene families in organisms [45]. Among 69 SsWRKY genes, only SsWRKY8 and SsWRKY9 displayed tandem duplication (Group IIa). Among the SsWRKY genes, the number of segmental duplications exceeded the number of tandem duplications. Twenty-four gene pairs were detected in the S. siamensis genome and were located on different chromosomes, suggesting that segmental duplications in these locations potentially contributed to expansion of the SsWRKY family. Moreover, three collinear maps of S. siamensis were constructed with three representative species, and it was shown that 56 SsWRKY genes (nearly 95%) were orthologous between CsWRKY genes (Table S2), implying that segmental duplication of WRKY genes might have occurred in diploid plants.
A gene's expression profile can often be related to its function. In this study, 59 SsWRKYs were analyzed in four different tissues of S. siamensis; root, stem, tuber, and leaf. Among SsWRKYs, 22 SsWRKY genes were expressed in the roots (Figure 7). These results concur with research reported for other species, such as A. thaliana [46], A. truncatum [19], and C. sativus [26]. It may be that roots, as an important organ for absorbing water and nutrients, respond first to stress when a plant is subjected to drought, high salt, or other stresses. For example, AtWRKY33 is a negative regulator that mediates Pi-deficiencyinduced remodeling of root architecture [47]; AtWRKY6 and AtWRKY23 were shown to regulate root growth and development [48]. SsWRKY12, SsWRKY21, and SsWRKY52 are homologous genes of AtWRKY6, AtWRKY23, and AtWRKY33. These genes may play key roles in root growth and development, and may perform a similar function in S. siamensis as they do in A. thaliana. SsWRKY46 is more highly expressed in stems than in other tissues, indicating that some SsWRKY genes may be involved in specific organ responses including stress. These results imply that SsWRKY genes showing different patterns of expression in different tissues are involved in numerous biological metabolism processes in S. siamensis.
The possible regulatory function of SsWRKY protein can be better understood using interaction network analysis of WRKY proteins from model plants. Five SsWRKYs (Ss-WRKY11, SsWRKY12, SsWRKY38, SsWRKY53, SsWRKY54) with high sequence similarity with AtWRKY62, AtWRKY33, AtWRKY30, AtWRKY40, and AtWRKY46 were identified as potential central nodes of the interaction network. Studies have reported that AtWRKY40 (homolog of SsWRKY53) can modulate the expression of stress-responsive nuclear genes encoding mitochondrial and chloroplast proteins [49]. Also, AtWRKY40 can alter resistance to pathogens. Notably, SsWRKY12 is also a node in the interaction network diagram. AtWRKY33 (homolog of SsWRKY12) has been proposed to interact with AtWRKY30 and AtWRKY15 (homolog of SsWRKY38, 22). Related studies showed that overexpression of AtWRKY33 in Arabidopsis could improve plants' salt-tolerance [50]. Furthermore, Ss-WRKY47 was predicted to interact with SsWRKY11, SsWRKY54, AtWRKY70 (homolog of SsWRKY47), and AtWRKY46 (homolog of SsWRKY54), and to play a key role in the regulation of brassinolide in plant growth, development, and drought tolerance [16]. In addition, a strong interaction between SsWRKY3 and SsWRKY39 was predicted. AtWRKY1 (homolog of SsWRKY3) negatively regulates plants' defense responses to Pst DC3000 through SA signaling pathways [51]. It is known that homologous proteins with similar sequences from different plants may have similar functions. In general, based on previous studies and prediction results of protein-function interaction, it can be inferred that the corresponding SsWRKY gene in S. siamensis may also have similar functions, and that closely related SsWRKY genes may have close functional relationships.
The WRKY gene family is especially associated with responses to biotic and abiotic stress [52]. Previous studies have shown that WRKY genes function in response to cold stress in many plants. For example, AtWRKY34 negatively effects the CBF-mediated cold response pathway [53]. In Vitis vinifera, cold stress induced the most rapid upregulation of VvWRKY genes of all the tested abiotic stress treatments [54]. In Brassica campestris, BcWRKY46 is induced by cold stress and ABA, improving the plant's tolerance of low temperature by activating related genes in the ABA signal pathway [55]. VbWRKY32 positively regulates the expression level of genes in response to cold, which increases antioxidant activity and maintains membrane stability, improving survival ability under cold stress [56]. In this study, 15 SsWRKY genes that contained one LTR cis-acting element in the promoter region were surveyed for their response to cold-temperature stress in leaf tissue (Figure 8). The results showed that all had altered expression throughout the experiment, and that nine SsWRKY genes were maximally expressed under cold stress after 48 h. Similar expression patterns were also found in V. vinifera [57]. In addition, several SsWRKY genes were highly expressed at 24 h, returning to the initial expression level after 48 h cold treatment, indicating that these SsWRKY transcription factors may function variously at different periods of the stress response. Although WRKYs have been observed to function in many plants in response to low temperature, the mechanism of how WRKYs respond to cold signals and regulate the expression of downstream genes remains largely unknown. Further research is required to demonstrate the function of these genes in relation to low temperatures and their involvement in cold signal pathways.

Plant Materials and Stress Treatment
S. siamensis plants used in this research were collected from Guangxi Academy of Agricultural Sciences (108.24 E;22.84 N) (Nanning, China). Tissue culture seedlings of S. siamensis were placed in a chamber with a mean temperature of 25.0 ± 1.0 • C, relative humidity of 60% ± 10%, and a day/light cycle of 12/12 h. For the cold treatment, S. siamensis seedlings were placed in low-temperature refrigerator at 4 • C and samples were gathered at 6, 12, 24, 36, and 48 h with 0 h as control. The leaves of six tissue culture seedlings were collected and mixed to provide one sample. The samples were snap frozen in liquid nitrogen and then stored at −80 • C freezer to extract total RNA.

Multiple Protein Sequence Alignment and Phylogenic Tree Construction
Based on the conserved domain of the SsWRKY protein, the family was divided into different groups using DNAMAN7.0 software. The ClusterW tool in MEGA7.0 software was applied to perform multiple protein sequence alignment analysis of WRKYs from S. siamensis and A. thaliana, using default parameter values and the neighbor-joining method to construct a phylogenetic evolution tree, with 1000 bootstrap iterations [61]. Subsequently, the graphic was transformed using Illustrator 2020 (v24.0.0.330). Based on the reported classification of WRKYs from A. thaliana, the SsWRKYs were classified into relevant groups.

Gene Structure and WRKY Protein Conserved Motif Analysis
The DNA and CDS sequences of the SsWRKY gene family were screened from the whole genome sequence and gene annotation file of S. siamensis. Then, the exon-intron structures of SsWRKYs were analyzed and visualized using TBtools software. The MEME v5.1.1 online program (https://meme-suite.org/meme/, accessed on 23 October 2022) was used to predict the conserved motifs of WRKY protein, with the maximum number of motifs set to 10 [62]. The diagrammatic sketches of gene structures and conserved motifs were re-edited using Photoshop 2022 software.

Chromosomal Location, Gene Duplication and Collinearity Analysis
All SsWRKY genes were mapped to the 14 chromosomes of S. siamensis according to their positions on the S. siamensis genome sequence, which was displayed using Circos (http://circos.ca/, accessed on 23 October 2022). Blastp was used to analyze the gene duplication events with default parameters and MCSanX to detect collinearity relationship between S. siamensis and three other species (Oryza sativa, A. thaliana, C. sativus) [63,64].

Promoter Cis-Acting Element Analysis of SsWRKY Gene Family
TBtools software (GTF/GFF3 sequences extractor) was applied to extract the 2000 bp sequence upstream of the CDS sequences of S. siamensis, and the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/, accessed on 23 October 2022) was used for prediction of cis-acting elements in the promoter regions [65].

Protein Interaction Network and Gene Co-Expression of SsWRKYs
For the protein interaction network, String (http://string-db.org/, accessed on 23 October 2022) was employed to analyze the SsWRKY protein interaction with reference to AtWRKY protein, with the parameter threshold set to 0.15 [66]. A. thaliana was selected as the model plant to analyze the co-expression of SsWRKY. The network was visualized using Cytoscape version 3.7.0 [67].

Expression Analysis of SsWRKY in Different Tissues
In order to analyze the expression patterns of SsWRKY genes in different tissues, the RNA-seq data of SsWRKY genes were gathered from the transcriptome database (unpublished), which includes gene expression levels in roots, tubers, leaves, and stems. The values for fragments per kilobase of transcript per million fragments mapped (FPKM) were utilized to calculate the genome-wide expression of the SsWRKY genes, which was then visualized on a heatmap using TBtools software, with the expression levels shown by different colors from blue to red. 4.9. Total RNA Extraction, cDNA Synthesis, and qRT-PCR Validation Samples consisting of 50-100 mg of tissue in liquid nitrogen were fully milled to a fine powder, and total RNA was extracted using the Trizol method as previously reported [68]. The purity and concentration of RNA was measured using a Micro Drop system (BIO-DL, Shanghai, China), and the RNA integrity was checked on 1% agarose gels. Total RNA (1000 ng) was used for reverse transcription using an Evo M-MLV RT mix kit for qPCR, with gDNA Clean (Accurate Biology, Changsha, China) to remove genomic DNA contamination, at 20 µL volume according to the manufacturer's instructions. The cDNA samples were diluted 5-fold with RNase-free water to use as templates for qRT-PCR analysis. Specific primer pairs were designed using Beacon Designer 7.0 software and were restricted to primer sequences of 18-25 bp, amplicon length 80-150 bp, melting temperature (Tm) 55-60 • C, and GC content of 45-55%; the primers used for this study are shown in Table S3.
The reaction was carried out in 96-well plates with a DiceTM Real Time System III thermal cycler (Takara, Japan) and a 20 µL reaction system: 2 × SYBR Green Pro Taq HS premix 10 µL (Vazyme, China), 2 µL cDNA, 0.4 µL each of forward primer and reverse primer, 7.2 µL RNase free water. The amplification reaction procedure was: 95 • C, 30 s; 95 • C, 5 s; 60 • C, 30 s, for 40 cycles. Each qRT-PCR analysis was performed with three technical replicates and three biological replications, the relative expression level was calculated using the 2 −∆∆CT method. The relative expression under cold-stress treatment of SsWRKYs containing LTR-motifs was visualized using GraphPad Prism 9 software.

Conclusions
We identified 59 WRKY genes distributed on 14 chromosomes in S. siamensis based on a genome-wide analysis. These 59 SsWRKY genes were divided into three major groups based on their phylogenetic relationships. The comparison of WRKY genes in S. siamensis with three diverse species helped to visualize the expansion and contraction of the SsWRKY gene family through its evolution. Expression profiles based on the RNA-seq data from four different tissues showed that these genes have complex expression profiles and may therefore have diverse but important roles in roots, stems, tubers, and leaves. Furthermore, a protein-protein interaction network was predicted for 18 of 59 SsWRKYs. Finally, 15 SsWRKY genes with LTR cis-acting elements revealed altered and complex expression in response to cold stress. Our study describes a functional framework for WRKY genes, which can provide the basis for further exploration of the function and regulatory mechanism of WRKY genes in S. siamensis.