Healthy and mature propagules of K. obovata were collected from Yueqing Bay, Yuhuan city, Zhejiang Province, China (28°13′N, 121°10′ E). The propagules were cultivated in a growth chamber (25°C temperature, 75% humidity, 14 h light/10 h darkness photoperiod) in plastic pots containing sand and watered with 1/2-strength Hoagland’s nutrient solution weekly. At the six-leaf stage, the seedlings were cultivated under chilling stress (4 °C) for 0, 1, 3, and 12 h. All the treatments included three seedlings. The leaves were collected, frozen immediately in liquid nitrogen and then stored at −80 °C.

The complete genome sequence files of K. obovata were downloaded from the Genome Sequence Archive (GSA) webpage (https://bigd.big.ac.cn/gsa/browse/CRA002395). Information on the WRKY domain hidden Markov model (HMM) profile numbered PF03106 was retrieved from the Pfam protein family database (http://pfam.sanger.ac.uk/), and the A. thaliana WRKY (AtWRKY) sequence information was obtained from The Arabidopsis Information Resource (TAIR) ( www.arabidopsis.org ). The candidate WRKY protein sequences were identified by comprehensive research using HMMER (E-value cutoff <1E-5) and BLAST analyses (72 AtWRKYs were used as queries) of the K. obovata whole-genome and protein databases. These KoWRKY sequences were identified by checking the complete WRKY conserved domain with SMART (http://smart.embl-heidelberg.de/) and InterPro (http://www.ebi.ac.uk/interpro/) and were confirmed manually by checking all sequences containing WRKYGQK motifs and a typical C₂H₂ (CX_4-5Cx_22-23HXH) or C₂HC (CX₇CX₂₃HXC) zinc-finger domain. The sequences of the confirmed KoWRKYs were input into the ExPASy online tool website (http://web.expasy.org/protparam/) to calculate the physicochemical properties of the proteins, including the molecular weight (MW) and isoelectric point (pI). The subcellular localization of the KoWRKYs was determined by ProtComp 9.0 (http://linux1.softberry.com/berry.phtml?topic=protcomppl&group=programs&subgroup=proloc).

The chromosomal positions of all the KoWRKY genes were determined from the genome annotation file. The physical positions of the KoWRKY genes on the chromosome were mapped using MapChart (Voorrips, 2002). The exon/intron structure of the KoWRKY genes was determined by comparing their predicted coding sequence (CDS) with genomic sequences using the gene structure display server (http://gsds.cbi.pku.edu.cn/). The upstream 2 kb sequences from the transcription initiation site of all the KoWRKY genes were extracted, and the cis-elements in these regions were identified using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

The criteria used for identifying duplicate genes were as follows: the similarity of the two aligned sequences was greater than 75%, and the length of the shorter aligned sequence covered more than 75% of that of the longer sequence. Two or more adjacent duplicates on the same chromosome within 100 kb were defined as tandem duplications, while duplicates across different chromosomes or at a distance greater than 100 kb on the same chromosome were considered segmental duplications (Kan et al., 2021). The synonymous rate (Ks) and nonsynonymous rate (Ka) of the identified KoWRKY gene pairs were calculated using KaKs Calculator version 2 (Wang Y et al., 2012). The approximate date [million years ago (Mya])] of each duplication event was estimated using the mean Ks values with the formula T = Ks/2λ, in which the mean synonymous substitution rate (λ) was 6.1 × 10^–9 (Lynch and Conery, 2000; Zhu et al., 2014). OrthoFinder with default parameters was employed to identify homologous genes between two different genomes (K. obovata and A. thaliana), and a synteny graph was constructed with TBtools according to the identified homologous gene pairs (Chen et al., 2020a).

Multiple EM for Motif Elicitation (MEME) was used to analyze the KoWRKYs and to identify 15 possible conserved motifs (http://meme.nbcr.net/meme/intro.html). The parameters were as follows: the repetitive time was “any”, the maximum motif number was 15, and the motif width was between 5 and 50 residues. The MEME results were subsequently displayed with TBtools software (Chen et al., 2020a). Multiple sequence alignment of the KoWRKY domains was performed using MAFFT version 7 (Katoh and Standley, 2013). Based on the alignment of the WRKY domains of the KoWRKYs and AtWRKYs, a phylogenetic tree was constructed with FastTree version 2 via the generalized time-reversible (GTR) model in conjunction with the Shimodaira–Hasegawa (SH) test (Price et al., 2010). All the identified KoWRKYs were divided into different groups according to the classification of AtWRKY sequences.

The 2-kb DNA sequences upstream of the ATG start codon of all genes assembled from the K. obovata genome were used to identify WRKY TF-binding sites. The potential KoWRKY target genes that were predicted with PlantRegMap were used for further pathway enrichment analysis with the KEGG database (https://www.kegg.jp/kegg/pathway.html) (Tian et al., 2019). Hypergeometric Fisher’s exact test (p < 0.01) and the Benjamini test [false discovery rate (FDR])<0.05] were performed to detect statistically significantly enriched KEGG pathways. The R package ggplot2 was used to visualize the top 11 significantly enriched KEGG pathways.

Transcriptomic data of K. obovata were obtained from the National Center for Biotechnology Information (NCBI) publicly accessible database (Accession number: PRJNA678025). The expression levels of the KoWRKY genes were analyzed in various tissues, including root, stem, leaf, flower, pistil, stamen, sepal, and fruit tissues. The expression data were transformed into log₂ [transcripts per million (TPM)+1] values for differential expression analysis (Zeng et al., 2020). The resulting gene expression profiles were visualized as a heatmap via TBtools software (Chen et al., 2020a).

The expression profiles of the KoWRKY genes that responded to chilling stress were obtained from the NCBI transcriptomic database (Accession number: PRJNA678025). The expression of the KoWRKY genes was visualized as a heatmap using TBtools software (Chen et al., 2020a). For the transcriptome analysis of KoWRKYs in response to chilling stress, thresholds of p < 0.05 and |log₂ (fold-change) ≥ 1| were used to define differentially expressed genes (DEGs).

Total RNA was extracted using RNASimple Total RNA Kit (Tiangen, Beijing, China) and cDNA was obtained using TIANScript cDNA kit (Tiangen, Beijing, China) according to the manufacturers’ instructions. Following total RNA extraction and cDNA synthesis, quantitative analysis was performed via real-time PCR in conjunction with SuperReal PreMix Plus (SYBR Green, Tiangen, China) on a CFX96 Touch Real-Time PCR Detection System (Bio–Rad, United States) according to the manufacturers’ instructions. The PCR program included an initial denaturation step at 95°C for 30 s, followed by 40 cycles of 95°C for 10 s, 58°C for 10 s, and 72°C for 30 s. The relative transcript levels of the candidate genes were calculated according to the 2^−ΔΔCT method (Livak and Schmittgen, 2001). All the data were generated from averages of three independent replicates, and statistical significance was determined by one-way ANOVA followed by Duncan’s multiple range tests. All the primers used are listed in Supplementary Table S1.

Differentially expressed KoWRKY gene was used as analysis object. Pearson correlation coefficient (PCC) between KoWRKY gene and non-KoWRKY gene was calculated to determine the coexpressed gene pair. Only gene pair with |PCC| value greater than 0.95 (p < 0.05) was used to regard as the potential regulated gene. The network was visualized using Cytoscape software (version 3.6.1).

To identify the WRKY-encoding genes present in the genome of K. obovata, the hidden Markov model (HMM) profile PF03106 from the Pfam database was used as a query for an HMM search against the genome, and a local BLASTP search was performed for which all 72 A. thaliana WRKY protein sequences were used in the query. In total, 64 potentially WRKY-encoding genes were identified and annotated. The 64 genes were named as KoWRKY1-KoWRKY64 according to the order of their gene ID, and the location on the chromosome, length of CDS, and the MW, pI and subcellular localization of their encoded proteins were shown in Table 1. Each chromosome contained KoWRKY members, of which chromosome 12 had the most, with 9 members; chromosome 18 had the fewest members, with only one. The CDS length ranged from 351 (KoWRKY23) to 2,181 bp (KoWRKY35), with an average length of 1,203 bp. The MW ranged from 12.99 kDa (KoWRKY23) to 78.60 kDa (KoWRKY35), with an average of 43.92 kDa. The theoretical pIs varied from 4.46 (KoWRKY33) to 10.26 (KoWRKY41), and the predicted subcellular localization results indicated that the 64 KoWRKYs were located in the nucleus (Table 1).

To determine the distribution of the KoWRKY genes across the genome, all 64 identified KoWRKY mRNAs/open reading frames (ORFs) were mapped onto their corresponding chromosome via BLAST searches against the released K. obovata genome sequence. As shown in Figure 1, the KoWRKY genes were unevenly distributed across the 18 chromosomes, and the numbers on each chromosome were not related to chromosome length. Chromosome 12 had 9 WRKY genes (the majority), including one member from each of groups I, IIe and III of the KoWRKY gene family and 2 members from groups IIb, IIc, and IId, followed by 5 members on chromosomes 1, 5, 6, and 8; however, chromosome 18 contained only one WRKY gene (KoWRKY 64, in subgroup IIa).

To investigate the structural diversity of the KoWRKY genes, the localization of intron–exon interactions was analyzed. It was found that the number of introns in the KoWRKY genes ranged from 1 (KoWRKY23 and KoWRKY28) to 6 (KoWRKY9) in K. obovata. A total of 34 (53.13%) KoWRKY genes had 2 introns, followed by 8 (12.5%), 11 (17.19%), and 8 (12.5%) genes that contained 3, 4, and 5 introns, respectively (Figure 2). Most KoWRKY genes in group I had 3 to 5 introns, except KoWRKY18 and KoWRKY9, which had 2 and 6 introns, respectively. The number of introns in the KoWRKY genes in group II widely varied, ranging from 1 to 5. However, all 7 KoWRKY genes in group III had 2 introns (Figure 2).

To analyze the functional diversification of the WRKY members, the sequences of the 2-kb upstream promoter regions of the KoWRKY genes were retrieved and analyzed. A number of cis-elements, including 9 elements related to plant development and 17 motifs related to stress responses, were analyzed, and the 26 elements are represented in Supplementary Table S2. The cis-elements related to plant growth and development include light-responsive elements (Box4, G-box, Sp1, ACE, and GT1 motifs), meristem-specific activation elements (CCGTCC-boxes), meristem expression-specific elements (CAT-boxes), endosperm expression-specific elements (GCN4 motif), and circadian control-related elements (circadians). Cis-elements related to stress responses include eight hormone-responsive elements (EREs, ABREs, CGTCA motifs, TGACG motifs, GARE motifs, TATC-boxes, P-boxes, TCA elements, and SAREs), low-temperature–responsive elements (LTRs and DRE cores), wound-responsive elements (WUN motifs), auxin response elements (TGAs and AuxRR cores), anaerobic induction elements (AREs), As shown in Supplementary Table S2, each KoWRKY gene contained more than one cis-element in its promoter region, and most KoWRKY genes contained Box4s, G-boxes, EREs, ABREs, TGACG and ARE motifs. Our analysis suggested that KoWRKY genes play an important role during development and stress responses.

To elucidate the mechanism underlying KoWRKY gene family expansion in K. obovata, BLASTP and MCScanX were employed to identify gene duplication patterns. The results showed that among the 64 KoWRKY genes, 49 segmentally duplicated genes formed 34 pairs (Supplementary Table S3); however, no tandem duplications were detected for KoWRKY genes, implying that segmental duplication was the major driving force for the expansion of KoWRKY genes.

To examine the selective pressure of KoWRKY genes, the Ka/Ks ratios of the duplicated gene pairs were calculated. Of the 35 paralogous KoWRKY gene pairs, the Ka values ranged from 0.05 to 0.54, the Ks values varied from 0.27 to 3.02, and all duplicated KoWRKYs had a Ka/Ks < 1, ranging from 0.12 to 0.45 (Supplementary Table S3). Estimated by a universal substitution rate of 6.1 × 10^–9 mutations per site per year, duplications of KoWRKYs may have occurred at two time points, approximately 22.3–40.7 MYA and 194.6–247.8 MYA (Supplementary Table S3). Previous studies suggested that a Ka/Ks < 1, a Ka/Ks = 1, and a Ka/Ks > 1 indicate purifying selection, neutral evolution, and positive selection, respectively (Nan and Gao, 2019). In this study, none of the Ka/Ks ratios of repeated KoWRKY gene pairs in K. obovata were greater than 1, which indicates that they underwent purifying selection.

To further explore the synteny relationships of KoWRKY genes with A. thaliana, Orthofinder was used to find the orthologous genes and TBtools was employed to construct the comparative synteny map (Figure 3). It showed that orthologous relationships between 50 KoWRKY genes and 32 AtWRKY genes were detected, and 59 orthologous WRKY gene pairs were identified based on these genes and the syntenic loci in K. obovata and A. thaliana chromosomes, suggesting that WRKY genes in both species had a similar origin and evolutionary process.

To gain insight into the functional regions of the KoWRKYs, the MEME program was used to identify the conserved motifs among the 64 KoWRKYs. A total of 15 conserved motifs were identified, namely, motifs 1–15 (Figure 4). These 15 conserved motifs are indicated with colored boxes according to their scale, with sizes ranging from 15 to 49 amino acid residues (Supplementary Figure S1). Among these individual motifs, motifs 1 and 3 were found to encode the conserved WRKY domain, and motifs 2 and 4 were found to encode the conserved zinc-finger structure. In addition to these conserved WRKY and zinc-finger motifs, other conserved motifs (motifs 5–15) were also predicted to be present within the KoWRKYs. Each KoWRKY protein had at least two conserved motifs, with the maximum number being seven, which was the case for several KoWRKYs. The distributions of the conserved motifs varied among the different KoWRKY groups. For example, group I KoWRKYs had seven motifs (motifs 1, 2, 3, 4, 6, 10, and 12), each containing one motif 1 and motif 3 (except KoWRKY18) and one zinc finger (motifs 2 and 4). Motifs 5, 6, 9, 11, 13, and 15 were present only in members of groups IIb, I, IIb, III, IIe and IIc, respectively. Motif 8 was present in members of group III and subgroup IIe, and motif 14 was present in subgroup IIa and IIb members. On the whole, the analysis of KoWRKY motifs showed that the KoWRKY members of every group or subgroup had similar motif compositions that corresponded to the clustering results generated by the phylogenetic tree analysis.

On the basis of the multiple alignment results for K. obovata WRKY protein sequences, the 64 WRKYs could be categorized into three groups (I, II and III) (Supplementary Figure S2). Most members of group I contained two conserved domains, WRKYGQK and a CX₄CX_22-23HXH zinc-finger motif in the N-terminal or a CX₄CX₂₃HXH motif in the C-terminal region. However, in KoWRKY18, the WRKYGQK sequence was replaced by WRKYGEK in the C-terminal region, and KoWRKY18 lacked a C-terminal WRKY domain. The 39 KoWRKYs in group II had a conserved heptapeptide WRKYGQK sequence (only KoWRKY23 contained a WRKYGEK sequence) and a CX₄–₅CX₂₃HX₁H zinc-finger domain, while 7 KoWRKYs in group III had a conserved WRKYGQK sequence and a CX₇CX₂₃HX₁C zinc-finger domain.

To investigate the evolutionary relationships between KoWRKYs and WRKYs from A. thaliana, an unrooted maximum-likelihood phylogenetic tree based on multiple alignments of the predicted amino acid sequences of the WRKY domains from K. obovata and A. thaliana was constructed. According to the tree (Figure 5), the KoWRKYs could be classified into three primary groups (groups I, II and III), with 18 KoWRKYs in group I, 39 KoWRKYs in group II, and 7 KoWRKYs in group III. Moreover, the KoWRKYs in group II could be further classified into five subgroups (groups IIa, IIb, IIc, IId and IIe) containing 2, 6, 14, 10, and 7 members, respectively. Notably, among all these groups and subgroups, most KoWRKYs were the members of subgroup IIc, the case of which is similar to AtWRKYs.

In order to lay the foundation for the functional analysis of KoWRKY proteins, their potential target genes were subjected to KEGG enrichment analysis. KEGG pathway enrichment analysis showed that the potential target genes of KoWRKY proteins were significantly (p < 0.05) involved in 11 pathways (Figure 6). The three pathways of metabolism, lipid metabolism, and glycosyltransferases had the most enriched genes, which were 637, 113, and 86, respectively. In particular, 73 target genes were involved in environmental adaptation pathways (Supplementary Table S4), which may play an important role in how KoWRKY genes regulate the adaptation process of K. obovata to the environment.

The expression profiles of all 64 KoWRKY genes were investigated using a standard transcriptome analysis procedure based on publicly available transcriptomic data of different tissues of K. obovata, including root, stem, leaf, flower, pistil, stamen, sepal, and fruit tissues. Among the 64 KoWRKYs, 50 were expressed in samples above (TPM > 0). Some KoWRKY genes showed preferential expression across all the tissues tested. KoWRKY genes, especially KoWRKY33, KoWRKY39, KoWRKY45, KoWRKY50, KoWRKY52, KoWRKY56, KoWRKY57, KoWRKY59 and KoWRKY63, were most highly expressed in the roots. Several members, such as KoWRKY2, KoWRKY6, KoWRKY9, KoWRKY21, KoWRKY39, KoWRKY50, KoWRKY54, KoWRKY60, and KoWRKY62 presented higher expression levels in the fruits than in the other tissues. Among all the tissues, the fewest KoWRKY members were expressed the most in the sepals. In addition, KoWRKY19 and KoWRKY37 were barely expressed in any of the tissues tested (Figure 7A).

To investigate the potential functions of KoWRKYs in response to chilling stress, the expression of KoWRKYs after chilling treatment based on publicly available transcriptomic data was analyzed. The results showed that the expression levels of most KoWRKY genes, especially KoWRKY16, KoWRKY28, KoWRKY32, KoWRKY43, KoWRKY45, KoWRKY55, KoWRKY61, KoWRKY63, and KoWRKY64, were upregulated after chilling treatment, the expression of which increased more than two-fold. However, the expression levels of KoWRKY37, KoWRKY42 and KoWRKY53 were not different from those of the control after the first, second and fourth cold treatments (Figure 7B).

To confirm the candidate KoWRKY genes that are important for cold tolerance, 9 DEGs were selected, and their expression levels were quantified via qRT–PCR. As shown in Figure 8, the expression of all the selected KoWRKY genes was upregulated under 4°C. Under low-temperature stress, the expression of KoWRKY16 increased sharply but then decreased, peaking after 1 h, the level of which was significantly higher than that of the control. The expression of KoWRKY43, KoWRKY63 and KoWRKY64 peaked at 3 h after 4°C treatment; however, the expression of KoWRKY28, KoWRKY32, KoWRKY45, KoWRKY55 and KoWRKY61 peaked at 12 h. The largest increase in the expression level (approximately 90-fold) was detected for KoWRKY28 after 12 h of chilling treatment.

Genome-wide gene expression profiling of KoWRKY16, KoWRKY28, KoWRKY32, KoWRKY43, KoWRKY45, KoWRKY55, KoWRKY61, KoWRKY63, and KoWRKY64 and mRNAs from leaves of plants subjected to chilling stress was conducted to identify genes coexpressed with KoWRKYs. The potential target mRNAs were predicted via Pearson correlation test for the 9 above mentioned KoWRKY genes whose expression increased more than two-fold after chilling treatment. The results showed that 263 significantly expressed mRNAs were correlated (PCC >0.95, p < 0.001) with 9 KoWRKYs, and for all 315 coexpressed pairs, each KoWRKY and mRNA was positively correlated (Figure 9, Supplementary Table S5). Taken together, the results indicated that the KoWRKYs might positively regulate the response of these putative genes in chilling stress.