Genome-wide analysis of the WRKY genes and their important roles during cold stress in white clover

Identification and classification of the TrWRKY gene family in white clover

The white clover genome resource information was released from the previous study, and all files were provided by Stig Uggerhøj Andersen from Aarhus University (Griffiths et al., 2019). DNA, CDS, and protein sequences were retrieved from the white clover genome. The sequences of Arabidopsis WRKY family proteins were collected from the TAIR database and used as BLAST (version 2.9.0+) query sequences to search the genome of white clover (Altschul et al., 1997), with an evaluation setting of 1E-05 and the coverages were set as 80%. The HMM file (version 3.3) (PF03106) was downloaded from the Pfam database (Punta et al., 2011), and the HMMER (evalue: 0.01) was used to identify and confirm WRKY DNA-binding domain, which was characterized as candidate WRKY proteins (Finn, Clements & Eddy, 2011). Annotated information for all candidate WRKY genes was retrieved from the white clover genome, including genome position, protein length, intron numbers, and these WRKY genes were classified into groups based on similar WRKY genes in Arabidopsis.

Phylogenetic analysis of the TrWRKY genes in white clover

Multiple sequence alignment analysis of Arabidopsis thaliana and white clover WRKY proteins using MUSCLE (version 5.1.0) with default parameters (Edgar, 2004). The phylogenetic tree of the WRKY gene family was constructed using the MEGA (version 11; Tamura, Stecher & Kumar, 2021) method with the following parameters: (1) Neighbor-joining (NJ); (2) Poisson correction; (3) genetic distance; (4) pair-wise deletion; (5) bootstrap: 1,000 replications. The TrWRKY genes were classified into different groups and subgroups based on the phylogenetic tree of AtWRKY and TrWRKY sequences.

Motif composition distribution analysis of TrWRKY genes in white clover

The conserved motifs of white clover WRKY protein sequences were identified using MEME (Multiple EM for motif Elicitation, Version 4.8.1) with the following parameters: (1) minimum and maximum width: 10 and 50, respectively; (2) the maximum number of motifs, 10; (3) number of appearances of a single pattern distributed in the sequences with model: 0 or 1 per sequence (-modzoops) (Bailey et al., 2009). All results were displayed with TBtools (version 1.098; Chen et al., 2020a).

Cis-acting elements analysis in the promoters of members of the TrWRKYs gene family

The 1,000 bp genomic sequence upstream of the transcription start site of the WRKYs gene family members was obtained from the white clover genome and the cis-acting elements in the promoter region of TrWRKYs gene family members were predicted by the PlantCARE online tool (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

Chromosomal location, gene duplication and Ka/Ks analysis of WRKY genes in white clover

All white clover proteins were compared to each other using the software BLASTP (version 2.9.0+), and the gene duplications were identified and characterized based on BLAST results by software MCSanX (version Python) with the default parameter (Wang et al., 2012). Based on the positional information of WRKY genes in the white clover genome and the duplication between genes, the software CIRCOS (version 0.69-8) was used to display the distribution of the WRKY gene family in the white clover genome (Krzywinski et al., 2009). Based on the phylogenetic tree and gene duplications results, a molecular evolutionary analysis of the TrWRKY genes was performed by calculating the nonsynonymous (Ka) to synonymous (Ks) substitution ratio of the duplicated gene pairs in S. tuberosum using the KaKs_Calculator in TBtools (Chen et al., 2020a).

Gene regulation network analysis of white clover WRKY gene family

The information on the Arabidopsis gene regulatory network (GRN) was extracted from the AraNet database (V2) (Lee et al., 2014), including 22,894 Arabidopsis genes and 895,000 interactions (links). All white clover proteins were subjected to BLAST searched with Arabidopsis proteins, with an e-value cut-off 1e-05, and the highest scoring hits were confirmed as homologous genes of the white clover gene. Meanwhile, all Arabidopsis proteins were also BLAST searched against white clover proteins with the same set, and hits with the highest scores were identified as homologous genes for Arabidopsis genes; the two genes were identified as homologous pairs with two BLAST results. Then, the GRN of white clover was constructed using the GRN of Arabidopsis based on homologous pairs. The sub-networks that contained the TrWRKY gene of white clover were searched and evaluated; the results were visualized using Cytoscape software (version 3.9.1), as well as the genes were annotated using GO information based on white clover genome annotation information (Shannon et al., 2003). Then, GO enrichment analysis of the sub-networks was performed using topGO software (version 2.50.0; Alexa & Rahnenfuhrer, 2019), the threshold level was set to 0.05 to show the most significant terms, and terms with high enrichment were assigned as GRN functions as described in the software protocol.

Expression analysis of white clover TrWRKY genes in response to cold stress

The RNA-seq data have been reported, the data included eight time points in response to cold stress; our previous work has described in detail, which could be assessed with accession numbers: PRJNA781064 (Zhang et al., 2022). All of the RNA-seq reads were mapped to the transcript sequences of white clover genomes using the Salmon software (version 0.12.0), and the expression level of each gene (FPKM value) was estimated using the subroutine quant of software Salmon (Patro et al., 2017). These expressional data were transformed using the “log2” function and they were centered using the “scale” function of the R program (version 4.2.1; R Core Team, 2022); then, all expression data were clustered and plotted using the “heatmap.2” function of the ggplots package (version 3.1.3).

Plant growth and qRT-PCR analysis

Seeds of white clover cv. Haifa were purchased from Barenbrug China Ltd. Com. (Beijing, China). All seeds were germinated and transferred to a mixture of perlite and sand each with a volume of 3:1, as our previously described (Zhang et al., 2022). In briefly, all seeds were growing in the pots, about 10–15 plants per pot. The growth temperature was 24 °C in light and 18 °C in darkness per day and irrigated with half-strength Hoagland solution once every 2 days. After 4 weeks, they were randomly divided into four groups for cold stress treatment. We collected at 0 min (control), 30 min, 1 h, and 3 h (4 time points in total) at a setting of 4 °C. For each group, three samples were randomly chosen of five seedlings were pooled to form a biological replicate. All samples were frozen in liquid nitrogen and stored at −80 °C. Total RNA was extracted from white clover seedlings at different time points of cold stress at 4 °C using the Total Plant RNA Extraction Kit (Tiangen, Beijing, China), divided into four time points (0 min (control), 30 min, 1 h, 3 h), and used Prime Script RT kit (Toyobo, Shanghai, China) was reverse transcribed into cDNA as a template for quantitative reverse transcription PCR. Primers were designed using Primer3 based on the nucleotide sequences of WRKY family genes (Table S1) (Untergasser et al., 2012). qRT-PCR was performed using a Light Cycler^® 96 system (Roche, Rotkreuz, Switzerland) and SYBR Premix Ex TaqTMII (Toyobo, Shanghai, China). Three replicates of each experiment were performed. The PCR conditions were set as follows. 95 °C for 2 min, 40 cycles, 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 1 min. Fold change values were calculated using expression abundance, which is based on the 2^−ΔΔCT method.

Identification of WRKY genes in white clover

A total of 145 TrWRKY genes were successfully identified from the white clover genome by multiple sequence alignment based on the amino acid sequence of the white clover WRKY gene family. These genes were proven to contain WRKY domains according to domain analysis and they were named TrWRKY1 to TrWRKY145 according to their chromosome locus and structure. All genomic information of these TrWRKY genes, including names, gene locus, chromosomal locations, group, introns and protein length (aa) were retrieved and summarized in Table 1. For these 145 TrWRKY genes, the largest protein was TrWRKY145, comprising 1,063 amino acids (aa), while the smallest one was TrWRKY65 (137aa). There are 16 chromosomes in the white clover genome; the TrWRKY genes were not evenly distributed across all chromosomes. The intron distribution was 1—11, of which the TrWRKY121 had the most introns (11 introns), while fifteen members contained only one intron; the results showed TrWRKY genes highly diverged.

Classification and phylogenetic analysis in white clover

An unrooted phylogenetic tree with 145 TrWRKY genes using neighbor joining methods (Fig. 1) was constructed to further explore the phylogenetic relationship of the WRKY transcription factor family in the white clover. This unrooted tree intuitively reflected the evolutionary status and grouping attribution of 145 members of the WRKY family. As shown in Fig. 1, the white clover WRKY proteins could be classified into three large groups (Group I-II) on the basis of the classifications of WRKYs in Arabidopsis. Specifically, the largest number of WRKY members in Group II was 84, while Group I and Group III had 37 and 24 members, respectively. In addition, Group II was further classified into five subgroups (SubGroup IIa-IIe). The most numerous subgroups were subgroup IIc, with 22 members. Next are IIb and IIe, with 21 and 20 members, respectively. Subgroup IId has 13 members. The subgroup with the lowest number of members in Group II is IIa, with only eight members. At the same time, the unrooted phylogenetic tree showed that the distribution of WRKY genes in white clover and Arabidopsis was highly consistent.

Motif composition distribution analysis of TrWRKY proteins in white clover

To further our understanding of the molecular structure and function of the TrWRKY gene family in white clover, we have analyzed the conserved motifs of WRKY gene family members and found the same subgroup had similar motif composition. Ten individual motifs were identified by the MEME tool, revealing the distinct regions of TrWRKYs (Fig. 2, and Figs. S1–S4). As Fig. 2 shows, the TrWRKY genes in Group I have contained nine motifs in total, most of them (except TrWRKY001, 003, 005, 015, 035) contained motif 3 and motif 6, while motif 3 has “WRKY” residues (see Fig. S1), which confirmed WRKY domain in these TrWRKY genes. Motif 1 is present in TrWRKY genes without motif 3, which also contained “WRKY” residues (see Fig. S1). The results showed WRKY domain has diverged in white clover. In addition, there are double WRKY domains identified in Group I members, even three copies of WRKY domains, which is also supported by domain search results. Group II and III have shown similar results, each TrWRKY gene contained motif 1 or motif 3, even two motifs, which consisted of BLAST and domain search results (Figs. S2 and S3). Meanwhile, the results of conservation motif composition also supported the results of sequence similarity and phylogenetic tree analysis, demonstrating clear structural motif differences between the three group. For example, most members of Group I contain motif 3, while members of Group III contain motif 1 (see Fig. S4), similar appearance was also discovered in Group II, each subgroup has diff motif composition patterns, see Fig. S2. In each subgroup, the proteins harbor a similar number and type of motif, which suggested the functional similarities of these TrWRKYs.

Cis-acting elements analysis of the TrWRKYs promoter

Promoter cis-elements influence the initiation of gene transcription. We performed a bioinformatics analysis to identify possible cis-elements in the promoter sequences of TrWRKYs. PlantCARE was used to identify putative cis-acting elements in the 1,000 bp upstream sequence of each TrWRKY gene promoter. A total of 15 stress response elements, consisting of TC-rich repeats (the cis-regulatory element for defense along with stress response), ACE (cis-regulatory element that engages in light response), LTR (cis-regulatory element that plays a role in low-temperature response), TCA-element (cis-regulatory element with a role in salicylic acid response), SARE (cis-regulatory element with a role in salicylic acid response), ABRE (cis-regulatory element associated with the abscisic acid response), AuxRR-core (cis-regulatory element with a role in auxin response), G-box (cis-regulatory element with a role in light response), CGTCA-motif (cis-regulatory element with a role in the MeJA-response), TGACG-motif (cis-regulatory element associated with MeJA-response), P-box (gibberellin-responsive element), GARE-motif (gibberellin-responsive element), WUN-motif (wound-responsive element), MBS (MYB binding site associated with drought-inducibility), and MRE (MYB binding site associated with light response), were identified (Fig. S5). All TrWRKYs had at least one stress response-linked cis-regulatory element. The cis-regulatory elements for hormone modulation consisting of CGTCA motifs, ABREs, AuxRR cores, P-boxes, TCA elements and TGA elements were also uncovered in numerous TrWRKY promoter regions. Overall, 64 TrWRKYs (44%) had more than one ABRE motif, which indicated the prospective abscisic acid response under stress conditions. Approximately 66 TrWRKYs (46%) had one or more CGTCA motifs that demonstrated the MeJA response potential, and the TCA element, TGACG motif, P-box, and AuxRR core were found in 30, 60, 12 and 8 TrWRKYs, respectively (Fig. S5). 68 G-box, 23 LTR, 29 MBS, and 25 TC-rich repeats were also found in TrWRKY promoter regions, which illustrated that these genes might play a role in cold, drought inducibility and defense responses.

Chromosome localization, gene duplication and Ka/Ks analysis of TrWRKY genes in white clover

To determine the evolution and expansion of WRKY genes, we used the MCScanX and Circos softwares to construct the distribution of WRKY genes across chromosomes. All TrWRKY genes were distributed across 16 chromosomes, but they were not uniformly located on these chromosomes, see Fig. 3. For example, the chromosomes TrChr1O, TrChr1P, TrChr8O and TrChr8P harbor more TrWRKY genes than other chromosomes, such as TrChr2O and TrChr2P. Based on BLAST results, there were 124 gene duplication events identified by MCScanX software, including 118 segment duplications (SD) and six tandem duplications (TD). The results have suggested chromosome doubling helps to bring about WRKY expansion in white clover, and distributions of TrWRKY were similar between some doubling chromosomes, for example, chromosomes TrChr6O and TrChr6P, chromosome TrChr5O and TrChr5P. However, there are some divergences between doubling chromosomes, such as chromosome TrChr3O and TrChr3P, chromosome TrChr7O and TrChr7P. The results suggested TrWRKY genes have undergone deep diversion in sub-genome evolution, some chromosomes have expanded TrWRKY members by gene duplication events, while some duplications have been purged, which caused some TrWRKY genes hots with numerous members clustering.

To estimate the divergence time of white clover WRKYs, synonymous (Ks) and nonsynonymous (Ka) substitutions between gene duplication pairs were calculated using the KaKs_Calculator in TBtool. When Ka/Ks <1, the genes experience purifying selection, which means the selection process could neutralize mutation to maintain the stability of the protein; in contrast, when Ka/Ks >1, the genes experience positive selection, which means great mutation happens in genes and eventually leads to a change in coded proteins. Our identified WRKY gene pairs had Ka/Ks values ranging from 0.06 to 0.89, proving that all of these genes experienced a purification selection process in white clover (Table S2).

Genetic regulation network analysis of white clover WRKY genes

Gene regulation networks (GRN) are increasingly used to explore the system-level functions of genes, we have reconstructed GRNs of TrWRKY and their interacting genes based public interaction database. The GRNs consisted of 349 genes and 463 interactions, as Fig. 4 shown. From GRNs, we have found most of TrWRKY have interacted with dozens of function genes, consisting of TrWRKY function on the transcription regulation process. For example, TrWRKY084 interacted with 37 genes, TrWRKY131 with 32 genes, and TrWRKY100 with 30 genes, the results indicated these TrWRKY genes played important roles in white clover lifespan. Gene Ontology (GO) annotation of these interacting genes was retrieved, and GO enrichment analysis was performed using the topGO package on the R platform. The results showed they were mainly distributed in the nucleus, see Fig. 5, which supported TrWRKY genes also functioned in the nucleus. In addition, molecular functions of these function genes were highly focused on transcription regulator activity, while they were mainly participating in the regulation of the transcription process, the results have confirmed the function of TrWRKY genes. It was notable that these genes were also enriched on terms “response to wounding” and “phosphorylay signal transduction system”, which are plant popular descriptions in response to abiotic stress, these results suggested TrWRKY genes probable function in response to abiotic stress.

Expression analysis of TrWRKY genes in response to the cold stress

In order to investigate TrWRKY genes function in response to abiotic stress, we have adopted our previous RNA-seq data under cold stress to assess their expression profiles. The white clover was treated with cold stress at 4 °C, and RNA-seq was analyzed at eight time points, including 0 h, 30 min, 1 h, 3 h, 6 h, 12 h, 24 h, and 72 h. All expressing TrWRKY genes (with FPKM value larger than 1) were collected, and their expressional value (FPKM) was grouped with a violin plot, and results showed TrWRKY genes were increased in response to cold stress, see Fig. S6. Especially, TrWRKY genes were rapidly activated at 30 min, and keeping a high expressional level in the following stages, the results suggested TrWRKY gene played critical roles in the early stage under cold stress. Expression profiles of these TrWRKY genes were clustered and displayed using the heatmap function, the results showed most of the TrWRKY genes were highly expressed at 30 min in response to cold stress, see Fig. 6. The finding consisted of violin plot analysis, which confirmed their rapid response to cold stress. Among these TrWRKY members, some were expressly high-expressing at 30 min, for example, TrWRKY017, 040, 049, 064, 065, 079, 085, 100, 102, 113, 138. Combining their hub-function in previous GRNs, such as TrWRKY100 and TrWRKY 113, they would be assigned rapid and critical regulation function in response to cold stress.

qRT-PCR validation of TrWRKY genes expression in response to cold stress

To validate the rapid response to cold stress, we have performed qRT-PCR analysis of seven TrWRKY genes with four time points, including 0 h, 30 min (0.5 h), 1 h, and 3 h. The qRT-RCR analysis results have confirmed their rapid and highly expressed in response to cold stress, see Fig. 7. All TrWRKY genes have shown a dramatic increase followed by a mild decrease in expression, the qRT-PCR and RNA-seq results are concordant, showing a similar pattern of TrWRKY genes expression in response to cold stress. These findings demonstrate that the TrWRKY genes are actively involved in the early response to cold stress and play critical roles in regulating gene expression in white clover.