Genome-Wide Analysis of Members of the WRKY Gene Family and Their Cold Stress Response in Prunus mume

Prunus mume, which is a rosaceous arbor with very high ornamental, edible and medical values, has a distribution that is mainly restricted by low temperature. WRKY transcription factor genes play crucial roles in the growth, development, and stress responses of plants. However, the WRKY gene family has not been characterised in P. mume. There were 58 PmWRKYs identified from genome of P. mume. They were anchored onto eight link groups and categorised into three broad groups. The gene structure and motif composition were reasonably conservative in each group. Investigation of gene duplication indicated that nine and seven PmWRKYs were arranged in tandem and segmental duplications, respectively. PmWRKYs were discriminately expressed in different tissues (i.e., roots, stems, leaves, flowers and fruits) in P. mume. The 17 cold-related candidate genes were selected based on RNA-seq data. Further, to investigate the function of PmWRKYs in low temperatures, the expression patterns under artificial cold treatments were analysed. The results showed that the expression levels of the 12 PmWRKYs genes significantly and 5 genes slightly changed in stems. In particular, the expression level of PmWRKY18 was up-regulated after ABA treatment. In addition, the spatiotemporal expression patterns of 17 PmWRKYs were analysed in winter. These results indicated that 17 PmWRKYs were potential transcription factors regulating cold resistance in P. mume.


Introduction
Low temperature is one of the most important abiotic factors affecting plant growth and development, and ultimately affecting the geographical distribution of plants. Chilling stress (0-15 • C) disrupts cellular homeostasis by altering the fatty acid composition of the membrane and accumulating the reactive oxygen species (ROS) in some organelles, which can deactivate proteins and interfere with normal physiological processes, especially for photosynthesis [1][2][3]. Freezing stress causes

Identification of P. mume WRKY Genes
To identify the WRKY family in P. mume, BLASTP and BLASTX were performed in the P. mume genomic database with the WRKY sequences of Arabidopsis thaliana that were downloaded from Plant Transcription Factor Database (http://planttfdb.cbi.pku.edu.cn/index.php). The hidden Markov model (HMM) searches were conducted in P. mume protein sequences using the WRKY domain HMM profile (PF03106). Candidate protein sequences without a WRKY domain were abandoned by screening of the SWISS-MODEL (http://swissmodel.expasy.org/). The sequences of genomic DNA, coding sequences (CDS) and protein sequences of PmWRKYs were obtained from National Center for Biotechnology Information (NCBI) (gene ID were shown in Table 1). ExPASy (http://web.expasy.org/compute_pi) and Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/) were used to separately predict the molecular weight (MW), isoelectric point (pI) and protein subcellular localization.

Genomic Localization and Gene Family Expansion Pattern Analysis
Based on the positional information in the P. mume genome database, genomic localization of PmWRKYs was drawn by MapDraw2.1 and Photoshop CS6.0. Research on gene family expansion pattern focused on tandem and segmental duplications. The tandemly arrayed genes were characterised as tandem duplications when they exhibited close phylogenetic relationships and were located at the same chromosomal location (within 100 kb) according to the criteria reported by Kong [24]. Segmental duplication regions were identified by the method of the Plant Genome Duplication Database [25]. First, the search for potential anchors was conducted by BLASTP (E < e −6 ). Subsequently, MCscan was applied to identify the homologous regions. Finally, syntenic blocks were evaluated using ColinearScan, and alignments with an E value < e −10 were deemed to be significant matches.

Gene Structure and Conserved Motif Analysis
The genomic sequences and CDS of PmWRKYs (Additional File 5) were obtained from the P. mume genome database. The exon-intron structures were identified with the Gene Structure Display Server2.0 (GSDS, http://gsds.cbi.pku.edu.cn/). It should be noted that since the database only contains the information of gene coding region, the untranslated region (5 UTR and 3 UTR) cannot be displayed by mapping the gene structure. PmWRKY protein sequences were submitted to the Multiple Em for the Motif Elucidation program (MEME) (http://meme.nbcr.net/meme/) to identify the conserved motif with the following parameters: any number of repetitions; the motifs number was set to 20; motif width was set to 6-200. SMART (http://smart.embl-heidelberg.de/) and Pfam were employed to annotate the identified motifs.

Heat Map Analysis by Transcriptome Data
The RNA-seq data of the five tissues (young roots, young stems, leaves, flowers, and immature fruits) (NCBI, Sequence Read Archive (SRA): SRP014885) and leaf buds collected before and after freezing stress in winter of P. mume cultivar 'Zhusha' (NCBI, SRA: SRP131731) have been used to draw the heatmaps of the PmWRKY genes. The 'before freezing' samples were collected when the lowest temperature of the day was higher than 5 • C, and the 'after freezing' samples were collected when the highest temperature of the day is lower than 0 • C in 2014. The time interval between the two samplings was one month. HemI (Heat map illustrator) was employed to draw heat maps with the default value [26].

Chilling and ABA Treatments
The annual branches were collected from the P. mume cultivar 'Yudie' for artificial treatments. Before the chilling treatment, the braches were exposed to 22 • C overnight by water culture and then transferred to 4 • C for 0, 0.5, 1, 2, 4, 8, 12, 24, 48, and 72 h in dark. For the ABA treatment, the branches were sprayed with 100 µM ABA for 0, 0.5, 2, 6, 12, 24, and 48 h. The stems were collected after the artificial treatments. The first time point (0 h) served as a control.
The stems and buds of the cultivar 'Yudie' were collected in winter from 5 November, 2017 to 4 March, 2018, which were planted in the open air. The daily land surface temperatures were recorded. Three replications of each sample were collected and all of the test samples were stored at −80 • C before total RNA isolation.

Gene Expression Analyses
Total RNA of each sample was extracted using the EASYspin Plus Plant RNA Extraction Kit (Aidlab, Beijing, China). The first-strand cDNA was obtained using a TIANScript RT Kit (KR107, Tiangen, Beijing, China). The specific primers were designed by Beacon Designer 8 based on cDNA sequences (Additional File 1). The expression levels of PmWRKYs, PmCBF1 (LOC103333423), PmCBF5 (LOC103337424), PmCBF6 (LOC103344251), PmLEA10 (LOC103340137) and PmLEA29 (LOC103321165) during the artificial low temperature and exogenous ABA treatments were examined using qRT-PCR with the following programme: 95 • C for 30 s, 40 cycles of 95 • C for 5 s and 60 • C for 30 s. Using the Actin gene of P. mume (ID: LOC103332029) as an internal control gene, the relative expression levels were calculated by the 2 −∆∆Ct method. Three independent experiments were performed with similar results.

Identification and Classification of WRKY Genes in P. mume
Two strategies were applied to identify PmWRKY genes: an HMM search and BLAST using well-characterised WRKY sequences from A. thaliana as queries. After SWISS-MODEL analysis, 58 sequences, which contained at least one WRKY domain, were obtained ( Table 1). The 58 genes were named PmWRKY01-PmWRKY58 according to the order of the gene ID and their corresponding protein sequences varied in length, MW as well as pI, as shown in Table 1. The length ranged between 162 and 884 AA, the MW ranged from 18.46 to 98.18 kDa and the pI varied from 4.92 to 9.68. According to the results of Plant-mPLoc, almost all of the PmWRKYs were predicted to be localised in the nucleus with high reliability except PmWRKY51, which was predicted to be targeted to chloroplasts besides the nucleus.
The relationships among the 58 PmWRKYs are shown in the phylogenetic tree ( Figure 1) produced by MEGA5.1; they were phylogenetically clustered into three main groups, which were similar to that of AtWRKYs [7]. Ten WRKYs, which contained WRKY domains and C 2 H 2 -type zinc-finger motifs in the N-terminal and the C-terminal were aggregated into group I. Group II consisted of 40 PmWRKYs, which were categorised into five subgroups using 12 diverse AtWRKY proteins, as references. Three members were clustered in subgroup II a, nine in II b, fourteen in II c, six in II d, and seven in II e. Furthermore, subgroup II c showed higher divergence than other subgroups, which was similar to those reported in the similar studies on other species [14,[27][28][29]. It is worth mentioning that PmWRKY51 was not clustered into group II although it contained a WRKY domain and a C 2 H 2 -type zinc-finger motif. Since the protein length and domain locations were considerably different, PmWRKY51 was classified into another group. In addition, there were also eight WRKY proteins, which contained a C 2 HC-type zinc-finger motif in group III.
Genes 2019, 10, 911 8 of 21 Figure 1. Tree of PmWRKYs. The unrooted phylogenetic tree of PmWRKY proteins was constructed using MEGA5.1 program by the neighbor-joining method with 1,000 bootstrap replicates. AtWRKYs were used as references to categorise PmWRKYs. The tree was divided into seven phylogenetic subgroups, designated as I, IIa-e, and III. The black solid points denote AtWRKYs, and the hollow points denote PmWRKYs.

Genomic Localization and Duplication of the PmWRKY Genes
Based on the genomic database, 57 PmWRKYs were distributed on all eight link groups of P. mume randomly and unevenly ( Figure 2), which leaves PmWRKY58 on the scaffolds. There were twelve PmWRKY genes on LG 1 and eleven genes on LG 2. In contrast, LG 8 contained only one gene.
LG 6, LG 7, LG 3, LG 4, and LG 5 contained 3, 5, 7, 9 and 9, respectively. It should be noted that group III genes were only located on LG 1, LG 5 and LG 7. Similar phenomenon occurred in P. persica [16], whereas the distribution of WRKY group III genes was even across all chromosomes in Brassica rapa ssp. Pekinensis [29]. The homologues of PmWRKYs in P. persica are shown in Table 1.
Gene duplication events fall into three categories: whole-genome duplication, tandem repeat and segmental duplication, which resulted in gene functional diversity, family evolution and plant adaptations [24]. Among the 58 PmWRKY genes, nine members were found in tandem repeats: PmWRKY08 to PmWRKY10, PmWRKY17 and PmWRKY18, PmWRKY28 and PmWRKY29, and PmWRKY41 and PmWRKY42 ( Figure 2). According to Holub's definition, 3 or more genes within 200 kb are considered a gene cluster [30]. The only gene cluster, which consisted of three group III genes (PmWRKY08, PmWRKY09, and PmWRKY10), was on LG 1. Twenty-eight genes were located on Figure 1. Tree of PmWRKYs. The unrooted phylogenetic tree of PmWRKY proteins was constructed using MEGA5.1 program by the neighbor-joining method with 1,000 bootstrap replicates. AtWRKYs were used as references to categorise PmWRKYs. The tree was divided into seven phylogenetic subgroups, designated as I, IIa-e, and III. The black solid points denote AtWRKYs, and the hollow points denote PmWRKYs.

Genomic Localization and Duplication of the PmWRKY Genes
Based on the genomic database, 57 PmWRKYs were distributed on all eight link groups of P. mume randomly and unevenly ( Figure 2), which leaves PmWRKY58 on the scaffolds. There were twelve PmWRKY genes on LG 1 and eleven genes on LG 2. In contrast, LG 8 contained only one gene. LG 6, LG 7, LG 3, LG 4, and LG 5 contained 3, 5, 7, 9 and 9, respectively. It should be noted that group III genes were only located on LG 1, LG 5 and LG 7. Similar phenomenon occurred in P. persica [16], whereas the distribution of WRKY group III genes was even across all chromosomes in Brassica rapa ssp. Pekinensis [29]. The homologues of PmWRKYs in P. persica are shown in Table 1.
Gene duplication events fall into three categories: whole-genome duplication, tandem repeat and segmental duplication, which resulted in gene functional diversity, family evolution and plant adaptations [24]. Among the 58 PmWRKY genes, nine members were found in tandem repeats: PmWRKY08 to PmWRKY10, PmWRKY17 and PmWRKY18, PmWRKY28 and PmWRKY29, and PmWRKY41 and PmWRKY42 (Figure 2). According to Holub's definition, 3 or more genes within 200 kb are considered a gene cluster [30]. The only gene cluster, which consisted of three group III genes (PmWRKY08, PmWRKY09, and PmWRKY10), was on LG 1. Twenty-eight genes were located on duplicated segments; among them, three genes (PmWRKY21, PmWRKY38 and PmWRKY56) were considered to be close relatives. PmWRKY21, PmWRKY38 and PmWRKY56 were located on LG 2, LG 4 and LG 7, respectively, and there was a triplet relationship among the three chromosomes [31]. Therefore, we surmised that these three PmWRKYs were generated along with a genome triplication process of Prunus. In addition, PmWRKY25 and PmWRKY34, as well as PmWRKY42 and PmWRKY53, were segmental duplicate pairs. The absence of the remaining 21 PmWRKYs' corresponding relatives may be due to the gene evolution of the P. mume WRKY gene superfamily. No duplication event was observed in group I. Therefore, 15.5% of the PmWRKY genes can be explained by tandem duplication, 12.1% can be accounted for by segmental duplication, whereas 72.4% were monogenes, which indicated that the formation of PmWRKY genes may not rely primarily on gene duplication.
Genes 2019, 10, 911 9 of 21 duplicated segments; among them, three genes (PmWRKY21, PmWRKY38 and PmWRKY56) were considered to be close relatives. PmWRKY21, PmWRKY38 and PmWRKY56 were located on LG 2, LG 4 and LG 7, respectively, and there was a triplet relationship among the three chromosomes [31]. Therefore, we surmised that these three PmWRKYs were generated along with a genome triplication process of Prunus. In addition, PmWRKY25 and PmWRKY34, as well as PmWRKY42 and PmWRKY53, were segmental duplicate pairs. The absence of the remaining 21 PmWRKYs' corresponding relatives may be due to the gene evolution of the P. mume WRKY gene superfamily. No duplication event was observed in group I. Therefore, 15.5% of the PmWRKY genes can be explained by tandem duplication, 12.1% can be accounted for by segmental duplication, whereas 72.4% were monogenes, which indicated that the formation of PmWRKY genes may not rely primarily on gene duplication. There was a WRKY gene (PmWRKY58) that could not be clearly located on the chromosomes, but could be identified on scaffolds.

Gene Structure and Conserved Motif Analysis of PmWRKYs
To make a thorough inquiry into the structural similarities and differences of PmWRKYs, we obtained the exon/intron structure diagrams and conserved motifs on the basis of the phylogenetic tree ( Figure 3A). First, according to the genome sequences and CDS of the PmWRKYs, we found that all of the genes contained at least one intron ( Figure 3B). Generally, the closest genes had similar structures, which only varied in the length of the intron and exon, whereas some genes exhibited different exon/intron arrangements. For instance, PmWRKY03 contained three exons, whereas its nearby paralogous gene PmWRKY40 had four exons and three introns even though their evolutionary relationships reached a 99% bootstrap value. Finally, all of the members of the subgroups II d and III were consistent with the number of exons. However, there was no significant consistency in the number of exons within other subgroups. The scale refers to a 5 Mb chromosomal distance. Genes in tandem repeats are underlined in black. Segmental duplicate genes are linked by blue lines. There was a WRKY gene (PmWRKY58) that could not be clearly located on the chromosomes, but could be identified on scaffolds.

Gene Structure and Conserved Motif Analysis of PmWRKYs
To make a thorough inquiry into the structural similarities and differences of PmWRKYs, we obtained the exon/intron structure diagrams and conserved motifs on the basis of the phylogenetic tree ( Figure 3A). First, according to the genome sequences and CDS of the PmWRKYs, we found that all of the genes contained at least one intron ( Figure 3B). Generally, the closest genes had similar structures, which only varied in the length of the intron and exon, whereas some genes exhibited different exon/intron arrangements. For instance, PmWRKY03 contained three exons, whereas its nearby paralogous gene PmWRKY40 had four exons and three introns even though their evolutionary relationships reached a 99% bootstrap value. Finally, all of the members of the subgroups II d and III were consistent with the number of exons. However, there was no significant consistency in the number of exons within other subgroups.
The MEME online tool predicted 20 individual motifs and revealed the specific regions of PmWRKYs ( Figure 3C and Additional File 2). An analysis of the 20 motifs revealed that the lengths of PmWRKY motifs ranged from 7 to 113 amino acid residues and the number of motifs varied from 1 to 11 in each PmWRKY protein. As shown in Figure 3C, motif composition was similar among the same subgroup, which suggested functional similarities of these PmWRKYs, whereas those of different subgroups had no common conserved motifs except for the C-terminal conserved motifs. For instance, motifs 3, 12, and 20 were conserved in group I. Motifs 9,11,16,19 were specific to subgroup II b. Motif 10 existed only in group II d. Motifs 1 and 2 were commonly shared by nearly all of the members of groups I and II, which were part of the WRKY domains. The results indicated that the conserved motifs mentioned above may bear special functions. According to the evolutionary analysis, II a and II b are two adjacent subgroups with near genetic distance, three unique motifs (motif 6, motif 7, and motif 14) nearly exist in all of the members of subgroups II a and II b, which supported the classification of PmWRKYs. In groups I and II, motif 1 and motif 8 both corresponded to the conserved heptapeptide domains, whereas motif 2 or motif 3 denoted the C 2 H 2 -type zinc-finger domain. However, in group I, motif 2 and motif 3 existed as a part of the N-terminal WRKY domain (NTWD) and the C-terminal WRKY domain (CTWD), respectively, which illustrated that the two WRKY domains belonging to family members of group I may be different in origin or function differentiation. It was also worth mentioning that there are three specific motifs (motif 4, motif 15, and motif 18) in PmWRKY08, PmWRKY09 and PmWRKY10 of group III. Although the function of the major motifs in the PmWRKYs was still indefinite; the PmWRKYs with the same conserved motifs may have similar functions. The MEME online tool predicted 20 individual motifs and revealed the specific regions of PmWRKYs ( Figure 3C and Additional File 2). An analysis of the 20 motifs revealed that the lengths of PmWRKY motifs ranged from 7 to 113 amino acid residues and the number of motifs varied from 1 to 11 in each PmWRKY protein. As shown in Figure 3C, motif composition was similar among the same subgroup, which suggested functional similarities of these PmWRKYs, whereas those of different subgroups had no common conserved motifs except for the C-terminal conserved motifs. For instance, motifs 3, 12, and 20 were conserved in group I. Motifs 9,11,16,19 were specific to subgroup II b. Motif 10 existed only in group II d. Motifs 1 and 2 were commonly shared by nearly all of the members of groups I and II, which were part of the WRKY domains. The results indicated that the conserved motifs mentioned above may bear special functions. According to the evolutionary analysis, II a and II b are two adjacent subgroups with near genetic distance, three unique motifs (motif 6, motif 7, and motif 14) nearly exist in all of the members of subgroups II a and II b, which supported the classification of PmWRKYs. In groups I and II, motif 1 and motif 8 both corresponded to the conserved heptapeptide domains, whereas motif 2 or motif 3 denoted the C2H2-type zinc-finger domain. However, in group I, motif 2 and motif 3 existed as a part of the N-terminal WRKY domain (NTWD) and the C-terminal WRKY domain (CTWD), respectively, which illustrated that the two WRKY domains belonging to family members of group I may be different in origin or function differentiation. It was also worth mentioning that there are three specific motifs (motif 4, motif 15, and motif 18) in PmWRKY08, PmWRKY09 and PmWRKY10 of group III. Although the function of the major motifs in the PmWRKYs was still indefinite; the PmWRKYs with the same conserved motifs may have similar functions.

Multiple Sequence Alignment of the WRKY Domains of PmWRKYs
The WRKY domain is an important functional and evolutionary unit of WRKY transcription factors and the conserved heptapeptide WRKYGQK near the N-terminal is regarded as the core

Multiple Sequence Alignment of the WRKY Domains of PmWRKYs
The WRKY domain is an important functional and evolutionary unit of WRKY transcription factors and the conserved heptapeptide WRKYGQK near the N-terminal is regarded as the core sequence of genes; their variation often leads to a decline or loss of DNA binding activity, which means that WRKY gene mutations may no longer have the original biological function [6]. Multiple sequence alignment analysis using WRKY domains of PmWRKYs found that mutations only occurred in PmWRKY04 and PmWRKY16 of subgroup II c. The mutation sites were both changed from Q to K, which formed the WRKYGKK sequence ( Figure 4). Interestingly, this type of mutation also occurred in the genome of Populus trichocarpa, O. sativa and A. thaliana [32][33][34]. We identified a CX 4 CX 22-23 HXH zinc finger motif in the N-terminal of subgroup I genes, a CX 4 CX 23 HXH motif in the C-terminal of subgroup I genes and II c genes, a CX 5 CX 23 HXH motif in II a, II b, II d and II e genes, and a CX 7 CX 23 HXC motif in subgroup III genes (Figure 4). From the multiple sequence alignment of the WRKY and zinc finger domains, we discovered that the homology of group I and II genes was higher than that of group III genes.
Through an analysis of the CDS and DNA sequences, we found that there were two types of introns in the WRKY domains. One is spliced exactly after the R position, similar to the splicing position found in Arabidopsis [7], and is designated as R-type intron. Another is spliced before the V position at the sixth AA after the second C residue of the C 2 H 2 -type zinc-finger motif [35]. We designate this type of intron as the V-type intron. Interestingly, the R-type intron is located before the zinc finger motif region in WRKY domains in subgroups I, II c, II d, II e and III, whereas the V-type intron is in subgroups II a and II b (Figure 4).

Expression Profiles of PmWRKY Genes in Different Tissues
According to the transcriptome data of five different organs (roots, stems, leaves, flowers and fruits) of P. mume, the heat map has been drawn by the RPKM values of PmWRKY genes ( Figure 5 and Additional File 3). As shown in the expression patterns of PmWRKYs in different subgroups in Figure 5A, the expression of WRKYs exhibited a wide range of diversity in five tissues. Most PmWRKYs were expressed in the roots at a high expression level. The expression of PmWRKY09 and PmWRKY10 (which all belonged to group III) could not be detected in the stems, leaves, flowers or fruits, and 12 genes (PmWRKY03, PmWRKY04, PmWRKY13, PmWRKY21, PmWRKY29, PmWRKY38, PmWRKY40, PmWRKY41, PmWRKY44, PmWRKY45, PmWRKY48 and PmWRKY56) lacked expression in one or two tissues. The expression levels of the other 44 genes were discovered in all of the detected tissues, which suggests that these PmWRKYs probably play vital roles in the developmental process of various tissues. We also hold the belief that the functions of duplicated genes diverged after the duplication event, which is supported by evidence that the expression patterns of several duplicated genes (PmWRKY08, PmWRKY09 and PmWRKY10, PmWRKY28 and PmWRKY29, PmWRKY41 and PmWRKY42, and PmWRKY21, PmWRKY38 and PmWRKY56) were quite different. In general, the group I genes had wider expression scopes than groups II and III genes. All of the genes in group I were expressed in all five tissues. However, a minority of groups II and III genes showed tissue-specific expression patterns. All of the subgroups II a and II b genes were hardly expressed in the stems, which meant that these genes might not participate in stem development. Furthermore, in the subgroups II d and II e, all of the members were expressed in five tissues, while PmWRKY48 showed a tissue-specific expression pattern, which demonstrated the functional change of genes in the same subgroups during evolution. The similarity of the expression patterns of subgroup II a, II b and subgroup II d supported the phylogenetic analysis.
As shown in the hierarchical clustering in Figure 5B, 58 PmWRKY genes were classified into six clusters according to their expression patterns. It was found that the genes with a closer phylogenetic relationship were more likely to be clustered into the same group in the heat map. The finest examples are three subgroup II c genes (PmWRKY16, PmWRKY35, and PmWRKY55) in cluster 6 and four subgroup II c genes (PmWRKY04, PmWRKY23, PmWRKY26, and PmWRKY32) in a subcluster of cluster 3, which show similar expression patterns. Among these clusters, clusters 1 and 6 contained the genes whose transcripts were detected in every tissue. However, the genes in other clusters presented tissue-specific expression patterns: cluster 2 genes abundant in the stems and flowers, cluster 3 genes in the roots, leaves and fruits, cluster 4 genes in the roots and leaves, and cluster 5 genes in the roots.

Expression Analysis of PmWRKY Genes under Chilling Treatment
To further examine the functioning of PmWRKY genes in cold tolerance, the expression patterns of PmWRKY genes before and after freezing stress in winter were analysed using transcriptome data. In Figure 6 and Additional File 4, after undergoing short freezing stress in P. mume 'Zhusha', 6 (PmWRKY18, PmWRKY23, PmWRKY32, PmWRKY37, PmWRKY44 and PmWRKY56) and 11 (PmWRKY03, PmWRKY04, PmWRKY06, PmWRKY08, PmWRKY13, PmWRKY14, PmWRKY27, PmWRKY28, PmWRKY42, PmWRKY52 and PmWRKY55) PmWRKYs were found up-or downregulated over two-fold, respectively. Then, these 17 PmWRKY genes, which may play a role in the cold tolerance of P.mume in winter, were selected for further study.

Expression Analysis of PmWRKY Genes under Chilling Treatment
To further examine the functioning of PmWRKY genes in cold tolerance, the expression patterns of PmWRKY genes before and after freezing stress in winter were analysed using transcriptome data. In Figure 6 and Additional File 4, after undergoing short freezing stress in P. mume 'Zhusha', 6 (PmWRKY18, PmWRKY23, PmWRKY32, PmWRKY37, PmWRKY44 and PmWRKY56) and 11 (PmWRKY03, PmWRKY04, PmWRKY06, PmWRKY08, PmWRKY13, PmWRKY14, PmWRKY27, PmWRKY28, PmWRKY42, PmWRKY52 and PmWRKY55) PmWRKYs were found up-or down-regulated over two-fold, respectively. Then, these 17 PmWRKY genes, which may play a role in the cold tolerance of P. mume in winter, were selected for further study.
To clarify the roles WRKYs played in cold resistance, the expression patterns of selected PmWRKY genes in different stages of artificial low temperature (4 • C) (0, 0.5, 1, 2, 4, 8, 12, 24 and 48 h) and exogenous ABA (100 µM) (0, 0.5, 2, 6, 12, 24 and 48 h) treatments were determined by qRT-PCR using cold resistant cultivar 'Yudie'. As shown in Figure 7A, the expression of PmCBF1, PmCBF5, PmCBF6 and dehydrin genes PmLEA10 and PmLEA29, which were reported to be involved in cold resistance previously [36], was induced significantly after chilling treatment. The expression levels of 17 genes were changed with different patterns during artificial chilling treatment. The expression levels of PmWRKY18, PmWRKY23 and PmWRKY 32 were increased with prolonging of the treatment time in the cold. The largest increase of the expression level (approximately 25-fold) was detected in PmWRKY18 after cold treatment for 24h. PmWRKY23 and PmWRKY32 also had the greatest up-regulation of more than 9.5-and 6-fold, respectively, at 48 h after being exposed to the cold condition. The three genes mentioned above were up-regulated continuously during low temperature treatment, whereas PmWRKY37 showed an irregular pattern. The expression of PmWRKY37 had fluctuations before and after it reached the maximum (nearly 2.7-fold) at 8 h. PmWRKY04 and PmWRKY06 had the greatest down-regulation of nearly 5-and 4-fold, respectively, and then slight up-regulation. The expression levels of six genes (PmWRKY03, PmWRKY08, PmWRKY14, PmWRKY42 and PmWRKY55) gradually decreased over time and ultimately changed over 2-fold, but the expression levels of the other genes (PmWRKY13, PmWRKY27, PmWRKY28, PmWRKY44, and PmWRKY56) only changed slightly. To clarify the roles WRKYs played in cold resistance, the expression patterns of selected PmWRKY genes in different stages of artificial low temperature (4 °C) (0, 0.5, 1, 2, 4, 8, 12, 24 and 48 h) and exogenous ABA (100 μM) (0, 0.5, 2, 6, 12, 24 and 48 h) treatments were determined by qRT-PCR using cold resistant cultivar 'Yudie'. As shown in Figure 7A, the expression of PmCBF1, PmCBF5, PmCBF6 and dehydrin genes PmLEA10 and PmLEA29, which were reported to be involved in cold resistance previously [36], was induced significantly after chilling treatment. The expression levels of PmWRKY56 were the same, showing a fluctuating increase over winter. In addition, the expression levels of PmWRKY04, PmWRKY06, PmWRKY08, PmWRKY13, and PmWRKY18 did not change significantly during winter. The results suggest that the transcription factor PmWRKYs may play an important role in overwintering survival of P. mume by different types of mechanisms. PmWRKY08 and PmWRKY28 may regulate negatively, while PmWRKY18, PmWRKY23 and PmWRKY44 may play a positively role in overwintering survival of P. mume.  Based on previous research, the transcriptional regulation of plant cold resistance can usually be categorised into either ABA-independent or ABA-dependent signal pathways; the latter also response to dehydration stress [37]. To explain how the PmWRKYs respond to ABA and whether the PmWRKYs are involved in the ABA-dependent cold signaling pathway or not, annual branches of P. mume were treated with exogenous ABA. Among these 17 genes, only the expression level of PmWRKY18 was continuously up-regulated, peaked (6-fold) at 6h and then gradually decreased during exogenous ABA treatment ( Figure 7B), whereas that of the others were lower than 2-fold and had no remarkable changes throughout, which suggested that PmWRKY18 may take part in cold adaptation in an ABA-dependent manner.

Expression Analysis of PmWRKY Genes in Winter
To verify the function of PmWRKY genes in cold tolerance, the expression patterns were analysed in winter using qRT-PCR. The daily land surface temperature was recorded from 12 Oct. to 4 Mar. As shown in Figure 8A, daily minimum and maximum temperatures below zero first appeared on 4 Nov. and 4 Jan., respectively, while the lowest temperature (−15.224 • C) was recorded on 26 January, after which the temperature began to increase with fluctuations. The stem and bud are the main aboveground organs of P. mume that overwinter. The expression of PmCBFs and PmLEAs was detected first ( Figure 8B). The expression levels of PmCBF1, PmLEA10 and PmLEA29 in stems increased during winter and peaked on 4 February then decreased throughout the remaining winter, while the expression levels of PmCBF5 and PmCBF6 decreased during winter. The expression patterns differed in buds and stems. The expression patterns of PmCBF1, PmCBF6, PmLEA10 and PmLEA29 in buds were similar; they increased as the temperature decreased, peaking on 3 December or 4 February The expression pattern of PmCBF5 in buds was similar in stems with both decreasing over winter. The expression patterns of PmWRKY03 and PmWRKY44 in stems were similar to those of PmCBF1, PmLEA10 and PmLEA29, which increased over winter, peaked on 4 February and 7 January then decreased. In contrast, the expression levels of PmWRKY08 and PmWRKY52 decreased over winter, were at their lowest levels on 7 January, then increased. The expression levels of PmWRKY32, PmWRKY37 and PmWRKY55 first decreased then increased during winter, peaking on 7 January or 4 February, then decreasing. The expression patterns of PmWRKY04, PmWRKY06, PmWRKY18, PmWRKY23, PmWRKY27 and PmWRKY28 in stems were similar. Their expression levels increased as temperature decreased with two peaks on 3 December and 4 February The expression levels of PmWRKY14 and PmWRKY42 increased gradually while PmWRKY13 decreased during winter. In addition, the expression of PmWRKY56 in stems did not change significantly in winter. The expression of PmWRKYs in buds also had different patterns. The expression patterns of PmWRKY23 and PmWRKY37 were similar to PmCBF1, PmCBF6, PmLEA10 and PmLEA29 in buds, increasing during winter and peaking on 3 December and 4 February The expression pattern of PmWRKY28 was similar to PmCBF5 in buds, with the highest levels on 5 November then decreasing during winter. In contrast, the expression level of PmWRKY14 increased during winter. The expression pattern of PmWRKY32, PmWRKY42, PmWRKY44 and PmWRKY55 were the same in buds. Their expression levels increased in winter time, and reached the highest peaks at 3 December or 7 January, and then decreased. On the contrary, the expression patterns of PmWRKY03 and PmWRKY27 decreased during winter and increased with increasing temperature. The expression patterns of PmWRKY52 and PmWRKY56 were the same, showing a fluctuating increase over winter. In addition, the expression levels of PmWRKY04, PmWRKY06, PmWRKY08, PmWRKY13, and PmWRKY18 did not change significantly during winter. The results suggest that the transcription factor PmWRKYs may play an important role in overwintering survival of P. mume by different types of mechanisms. PmWRKY08 and PmWRKY28 may regulate negatively, while PmWRKY18, PmWRKY23 and PmWRKY44 may play a positively role in overwintering survival of P. mume.

Discussion
With the development of molecular biology, especially sequencing technology, the genomic sequencing of more and more plants has been completed, which provides convenient conditions for research from the genomic level in a systematic and global way. From the P. mume genome databasebased analysis, 58 PmWRKY genes were identified (Table 1). Though 58 PpWRKY genes were

Discussion
With the development of molecular biology, especially sequencing technology, the genomic sequencing of more and more plants has been completed, which provides convenient conditions for research from the genomic level in a systematic and global way. From the P. mume genome database-based analysis, 58 PmWRKY genes were identified (Table 1). Though 58 PpWRKY genes were reported in P. persica [16], these genes in two species do not exactly correspond one by one, PmWRKY28 and PmWRKY29 are the two closest homologues of PpWRKY19. The homologue of PmWRKY51 was not described, and no homologous gene was found for PmWRKY16 in P. persica. Gene duplication, which consists of whole-genome duplication, tandem repeat and segmental duplication, plays a vital role in the amplification and evolution of the plant gene family. The study on the WRKY gene family of O. sativa concluded that the production of 80% OsWRKYs resulted from gene duplication events [35]. However, as Figure 2 shows, there were only 16 PmWRKYs (27.6%) involved in seven gene duplication events, which is much lower than the rate of duplication genes in O. sativa. We speculated that the PmWRKY gene family may mainly originate from two reasons: One is that the formation of PmWRKY genes may not depend on gene duplication. Another, which occurs in the evolution process, is the loss-function or redundant genes in the PmWRKY gene superfamily may be lost in the genome to reduce energy consumption. It was shown that the gene duplication of WRKYs in A. thaliana and O. sativa mainly occurred in group III [35]; in contrast, most of the PmWRKY gene duplication events were found in group II (13/16). In addition, genomic localization indicated that 58 PmWRKY genes were unevenly mapped to all eight chromosomes (Figure 2). This indicated that the WRKY gene families in different species originated from different evolution patterns.
The analysis of PmWRKY domain structure showed that there were two kinds of conserved introns in PmWRKY domains: R-type intron and V-type intron. The R-type introns only existed in the PmWRKY domains of subgroups I, II c, II d, II e and III, whereas the V-type introns were only detected in subgroups II a and II b and this was consistent with previous research on O. sativa [35]. In addition, another form of intron, which was located at the fourth AA residue (K) after the second C residue of the zinc-finger motif, was found in subgroups II a and II b of A. thaliana and C. sativus [7,38]. It is not clear whether or not the insertion of conserved introns in different loci affect the function of WRKYs. The heptapeptide sequence WRKYGQK is the conserved structure of WRKY genes. It can be found in multiple alignment analysis of WRKY domains in P. mume. Some mutations occurred in subgroup II c (Figure 4), which indicated that the selection pressure and evolution pattern of PmWRKY genes vary between subgroups. As is well known, theWRKY domain contains a zinc-finger motif C 2 H 2 or C 2 HC besides the conserved heptapeptide WRKYGQK [7]. Our study demonstrated that the zinc-finger motif C 2 HC only existed in the PmWRKY domains of group III whereas C 2 H 2 was owned by PmWRKY domains of other groups. We found that this result is slightly different from that of previous research on O. sativa [35], which announced that zinc-finger motif C 2 H 2 and C 2 HC simultaneously occurred in the N-termial of OsWRKYs of subgroup I.
RNA-seq data and qRT-PCR analysis provide a good clue to identify the important PmWRKY genes in cold resistance of P. mume. In our study, the expression of some PmWRKYs was changed under cold stress in P. mume cultivar 'Zhusha' based on the RNA-Seq ( Figure 6). There were 6 and 11 PmWRKYs that were found to be up-or down-regulated over two-fold, respectively, after being subjected to cold stress in winter. Since some factors were uncontrolled during winter time, artificial low temperature treatment was performed to verify the cold response of PmWRKYs. Based on the qRT-PCR results, the expression levels of the 12 selected PmWRKYs were significantly different and 5 slightly changed under the artificial chilling treatment ( Figure 7A). In addition, spatiotemporal expression patterns of the 17 PmWRKYs candidate genes varied during winter. The expression patterns of most PmWRKYs differed in stems and buds, which may result from the asynchronous icing of the organs ( Figure 8B). However some factors, for example developmental stage, air humidity and day length, which were uncontrollable in winter time, may also affected gene expression. The results suggest that these genes potentially took part in the cold resistance of P. mume. Bud dormancy is the important way plants in Prunus genus respond to the cold winter. It was reported that PpWRKY11 (Prupe.1G071400), PpWRKY18 (Prupe.1G393000), PpWRKY33 (Prupe.6G286000), PpWRKY46 (Prupe.2G185100), PpWRKY48 (Prupe.1G114800) and PpWRKY53 (Prupe.2G307400) may play a role in bud dormancy in P. persica (Table 1) [16]. PmWRKY18, PmWRKY55 and PmWRKY27, the homologous genes of PpWRKY11, PpWRKY33 and PpWRKY48 in P. mume, showed low temperature response in different degrees. These genes may enhance cold tolerance of P. mume by regulating bud dormancy in a low temperature induction pattern.
Low temperature signal transduction pathways are mainly categorized into either ABA-dependent or ABA-independent manners [37], both of which ultimately regulate the expression of functional genes, COR47, RD29A and KIN1, and enhance the chilling tolerance of plants. WRKY transcription factors have been reported to be key nodes in ABA signalling [39]. It is reported that AtWRKY18, AtWRKY40 and AtWRKY60 genes in A. thaliana can be induced by ABA and act as negative regulation factors in the ABA signaling pathway [40,41]. AtWRKY6, the homologous gene of PmWRKY37, functioned as a positive regulator of ABA signalling during the germination process [42]. In this work, we detected the response of 17 selected genes to exogenous ABA treatment. The results showed that only PmWRKY18 was significantly up-regulated after ABA treatment ( Figure 7B); whereas other genes had no significant changes compared with the control (data were not shown). We speculated that PmWRKY18 may take a part in an ABA-dependent cold signaling pathway.

Conclusions
Although cold tolerance related WRKY genes were studied in various species of plants, the cold signal response mechanism of WRKY genes and the down-stream target genes are mostly unclear. Research on the molecular biological function of PmWRKY genes in response to low temperature needs to be conducted. In this study, the characterisation and the expression analysis of PmWRKYs can create a foundation for further gene isolation and function analysis to clarify the role of WRKY genes in the cold resistance of P. mume.