Characterization of the WRKY gene family reveals its contribution to the adaptability of almond (Prunus dulcis)

Screening of WRKY gene family members

The GDR database (https://www.rosaceae.org/) was accessed to download the genome data for Prunus dulcis cv Texas (Alioto et al., 2020). The WRKY domain hidden Markov model (PF03106) dataset was downloaded from the Pfam database (Sonnhammer et al., 1998). The whole-genome protein data of Prunus dulcis were searched and compared with regard to the WRKY domain by using the Hmmer tool, and the E-value was set to ≤1e⁻⁵ (Finn, Clements & Eddy, 2011). The protein molecular weight, isoelectric point, and hydrophobicity of the protein were assessed. ExPASy (http://au.expasy.org/tools) and Wolfe PSORT II (https://www.genscript.com/wolf-psort.html?src=leftbar) were used for protein subcellular localization analysis (Gasteiger et al., 2003; Paul, Keun & Takeshi, 2007).

Multiple sequence alignment and phylogenetic tree construction

Cluster Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) ClusterW multiple sequence alignment was performed on PdWRKY protein sequences (Sievers & Higgins, 2018). The Arabidopsis thaliana AtWRKY protein sequences were downloaded from the UniProt protein database (https://www.uniprot.org/). Multiple alignments of PdWRKY and AtWRKY proteins were performed using the MEGA X tool selection MUSCLE method (Kumar et al., 2018; Liu et al., 2018). The aligned sequences were analysed by the neighbour-joining (NJ) method to build a phylogenetic tree based on the following parameters: Poisson model, pairwise deletion, and 1,000 bootstrap replications.

Analysis of the motifs and gene structures of PdWRKY family members

The MEME (https://meme-suite.org/meme/tools/meme) tool was used to perform conserved motif analysis on the PdWRKY protein sequences. The number of search motifs was set to 10, and the remaining parameters were set to the default values (Bailey, Mikael & Buske, 2009). The numbers of exons/introns in the PdWRKY gene sequences across the whole genome were extracted for gene structures analysis. Tbtools was used to combine the conserved protein motifs and gene structures with the PdWRKYs phylogenetic tree (Chen, Chen & Zhang, 2020).

Chromosome location, collinearity and Ka/Ks analysis of the PdWRKY genes

According to the chromosome location information for the PdWRKY genes, TBtools software was used to draw the location distribution map (Chen, Chen & Zhang, 2020). MCScanX was used to analyse the collinearity of the PdWRKY genes (Wang, Tang & Debarry, 2012), and the Circos tool was used to produce a chromosome collinearity map (Martin et al., 2009). In addition, we further selected PdWRKY and WRKY gene family members from Arabidopsis thaliana and Rosaceae plants (Malus domestica, Prunus persica, Prunus avium, Prunus armeniaca and Rosa chinensis) to study the collinear gene distribution relationship, and we used the Ka/Ks Calculator tool to calculate the Ka/Ks ratios of the collinear genes (Wang, Zhang & Zhang, 2010).

Analysis of cis-elements in PdWRKY genes

The 2,000-bp upstream sequences of the identified PdWRKY1 genes were extracted using TBtools. The extracted sequences were then submitted to the PlantCARE website (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to predict the cis-elements in the promoter regions (Lescot, 2002). A diagram of the cis-elements within the PdWRKY1 genes was constructed using TBtools and modified via Adobe Illustrator CC 2019 (Chen, Chen & Zhang, 2020), and statistical analyses were performed (Liu et al., 2018).

Expression pattern analysis

Transcriptome data from four almond tissues (leaves, buds, flowers, and pistils) were downloaded from the NCBI SRA database to obtain the corresponding FPKM values, and the expression patterns of PdWRKY genes in different tissues were analysed.

Material handling and qRT–PCR analysis

To analyse the transcription level of each PdWRKY genes after different abiotic stress treatments, 1-year-old ‘Zhipi’ almond seedlings were collected from the almond resource nursery in Shache County, Xinjiang. These seedlings were placed in Murashige and Skoog (MS) liquid medium containing 300 mmol/L NaCl and 20% polyethylene glycol (PEG; mass fraction) and placed in an artificial climate chamber at 4 °C to simulate low-temperature stress (Alam et al., 2019). For each treatment, the leaves were collected at 0, 6, 12 and 24 h.

To verify the expression pattern of 1-year-old almond seedlings, we selected 10 representative PdWRKY genes for qRT–PCR analysis. Three independent biological repeats containing three independent plants were used for qRT–PCR detection. The qRT–PCR primers of the selected PdWRKY genes were designed by Primer Premier 5 (Table S1). Total RNA was extracted from frozen embryonic samples by using the RNAprep Pure Plant Kit (TIANGEN, Beijing, China). The reaction system (25 μL) consisted of the following: 10 μL of SYBR Green PCR Master Mix, 0.4 μL of forward primer (10 μmol/L), 0.4 μL of reverse primer (10 μmol/L), 5 μL of diluted cDNA (50 ng/μL), and 25 μL of RNase-free water. The reaction conditions were as follows: predenaturation at 95 °C for 3 min; 40 cycles of 95 °C for 10 s, 55 °C for 20 s, and 72 °C for 20 s; and finally, 75 °C for 5 s. Plate reads were then performed to detect fluorescence signals (40 cycles). Melting curve analysis was conducted at a temperature range of 65–95 °C with a temperature increment of 0.5 °C/5 s. The obtained cycle threshold (CT) values were quantitatively analysed by the 2^−ΔΔCT method (Livak & Schmittgen, 2001) (Table S2).

Construction of the PdWRKY proteins interaction network

The PdWRKY protein sequences were uploaded to the STRING database (https://string-db.org/) for comparison, and the relationships between important proteins were predicted based on Arabidopsis thaliana protein interactions. Gene Ontology (GO) annotation of PdWRKY genes was performed using TBtools.

Identification and annotation of target genes of PdWRKY genes

To obtain the possible downstream target genes regulated by WRKY, TBtools was used to extract the 2,000-bp gene promoter sequences for almond genes. The conserved motif of the WRKY DNA-binding site (MA1298.1) was obtained from the JASPA_CORE database, which contains eukaryotic TF binding profiles (http://jaspar.genereg.net) (Aziz et al., 2018). Then, the Motif FIMO (https://meme-suite.org/meme/) tool was used to detect the motifs in the almond promoters that bound to WRKY. The final candidate target genes were determined based on the screening criteria of P < 1 × e⁻⁶ (Bailey, Mikael & Buske, 2009). Target genes of candidate WRKYs were functionally annotated using GO and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Analysis and visualization were performed using the Omicshare online platform (https://www.omicshare.com/tools/). Finally, protein domain annotation of candidate target genes was performed with the InterPro database (https://www.ebi.ac.uk/interpro/).

Analysis of the characteristics of PdWRKY family members

A total of 62 WRKY genes were screened from the whole genome of almond and renamed PdWRKY1–PdWRKY56. Six pairs of homologous genes were renamed PdWRKY7a/7b, PdWRKY11a/11b, PdWRKY28a/28b, PdWRKY311a/31b, PdWRKY38a/38b and PdWRKY48a/48b. The protein lengths of the PdWRKY family members ranged from 170 aa (PdWRKY40) to 748 aa (PdWRKY42); the molecular weight ranged from 19.55 kDa (PdWRKY40) to 80.85 kDa (PdWRKY42); and the isoelectric point ranged from 4.99 (PdWRKY38a) to 10.03 (PdWRKY13). The PdWRKY proteins GRAVY (grand average of hydropathicity) values ranged from −0.532 (PdWRKY13) to −1.096 (PdWRKY25), indicating hydrophilicity. The predicted results of subcellular localization indicated that all examined proteins were located in the nucleus (Table S3).

Multiple sequence alignments results and phylogenetic relationships of PdWRKY family members

Almond PdWRKY family members are divided into three groups, I, II and III, according to the domains found in AtWRKY family members (Fig. 1). Among them, 12 PdWRKYs were distributed in Group I, seven PdWRKYs were distributed in Group III, and 43 PdWRKYs were distributed in Group II. In addition, Group II was further divided into five subgroups, IIa, IIb, IIc, IId and IIe, with four, nine, 14, nine and seven PdWRKYs, respectively. In addition, we further constructed a maximum likelihood (ML) phylogenetic tree. Compared to the neighbour-joining (NJ) phylogenetic tree, the number of PdWRKYs in group I was reduced to seven and the number of PdWRKYs in group II was increased to 48. Also, members of group I and subgroup IIc did not cluster on the same branch. Notably, the five genes in group I were PdWRKY42/43/53/48a/48b and they were divided into subgroup IIc. Therefore, the clustering results of the neighbour-joining phylogenetic tree were better. The WRKYGQK heptapeptide is considered a hallmark of the WRKY protein. The multiple alignment method was used to analyse the conservation of domains across PdWRKYs. The PdWRKYs all contained WRKYGQK and zinc finger domains (Fig. S1). The members of Group I all had two WRKYGQK domains, while those of Group II and Group III had only one WRKYGQK domain, and all examined proteins had conserved C-C and H-H/C domains. Additionally, PdWRKY7 and PdWRKY37 had a single amino acid substitution: K changed to Q.

Motifs and gene structures of PdWRKY family members

PdWRKY proteins have different motif types (Fig. 2B). Motifs 1 and 2 are typical WRKY domains (Fig. 2D). All PdWRKY proteins contain Motif 1, except for PdWRKY8. Moreover, the motif types were conserved within the same group, indicating that members in the same group have similar functions. Members in Group I all have two copies of Motif 1, which is consistent with the multiple sequence alignments results.

To explore the gene structures of PdWRKYs, we constructed an exon/intron structures map (Fig. 2C). The results showed that PdWRKYs have at least one intron, and four, 15, five, 30 and eight members have six, five, four, three and two exons, respectively. The genes in the same group have similar structures; for example, all genes in Group III have three exons and two introns. In addition, we compared the distribution of WRKY sequences and exons for each PdWRKY and found that most PdWRKYs have WRKY domains covering or spanning one exon, Members in Group I have two WRKY domains, often spanning two exons.

Gene mapping, collinearity and Ka/Ks of PdWRKY family members

The gene chromosome mapping results showed that 62 PdWRKY members were unevenly distributed on eight chromosomes (Fig. 3). Chromosomes 1–8 contained 14 (22.58%), 8 (12.90%), 11 (17.4%), 7 (11.29%), 6 (9.68%), 13 (20.97%), 1 (1.61%) and 2 (3.2%) members. These results indicate that PdWRKY genes can mainly be found on chromosomes 1–6.

The collinearity results for genes within the genome showed that the PdWRKYs had 24 pairs of segmental duplications (Fig. 4) and nine pairs of tandem replications (Fig. 3; Table S4). The segmental duplications were mainly distributed on six chromosomes: Chr1 (six), Chr2 (four), Chr3 (nine), Chr4 (five), Chr5 (four) and Chr6 (five). Chr8 has only one PdWRKY gene, PdWRKY55, Chr7 has no PdWRKY gene, and most of the segmental duplications of PdWRKY genes occur in Group II. There were two or more pairs of tandem repeats in Groups I, II and III, which were mainly distributed on chromosomes 1 and 6 (Fig. 3).

To explore the evolutionary relationship of the PdWRKY gene family members, we constructed a comparative collinear map of six species to depict the collinearity of WRKY gene members in almond, Arabidopsis thaliana and five representative Rosaceae species (Fig. 5; Table S5). The results showed that Prunus dulcis had 38, 82, 55, 54, 46 and 51 collinear gene pairs with Arabidopsis thaliana, Malus domestica, Prunus persica, Prunus avium, Prunus armeniaca and Rosa chinensis, respectively. Compared with Arabidopsis thaliana, Prunus dulcis and Rosaceae have more collinear WRKY genes.

This study further calculated the Ka/Ks ratios between the PdWRKY duplications and collinear WRKY genes in almond and six other species to explore the PdWRKY gene family and WRKY gene evolution among species under selection. The results showed Ka/Ks ratios <1 for all duplicated PdWRKY genes and Ka/Ks ratios <1 for collinear WRKY genes among all species. Thus, the PdWRKY gene family members and WRKY genes in different species have experienced strong purifying selection during the evolution process (Fig. 6; Table S6).

Upstream cis-regulatory elements

The promoter sequences 2,000-bp upstream of the PdWRKY genes were extracted for cisacting element analysis, and the numbers of target elements were counted (Fig. S2; Tables S7 and S8). The results showed that in addition to a large number of basic light, TATA- and CAAT-box elements, additional annotated elements were related to plant hormones, stress resistance, defence and stress reactions, drought induction, anaerobic induction, and growth and development, including the auxin, gibberellin, abscisic acid, and MeJA pathways, indicating that PdWRKY genes are widely involved in regulating a variety of life processes such as almond growth, development, metabolism and stress resistance.

Expression pattern analysis of the PdWRKY family members

The transcriptome data for almond leaves, buds, flowers and pistils were selected for expression analysis across different tissues. The results showed that five PdWRKY genes were not expressed in these four tissues (Fig. 7; Table S9). These genes may be pseudogenes or genes with special temporal and spatial expression, or they may be expressed in other tissues. The other 57 PdWRKY genes were expressed in at least one tissue, and most of the PdWRKY genes were only expressed at low levels. This result may be related to the interaction of WRKY transcription factors with other genes or proteins during plant growth and development. There were 13, 20 and 22 PdWRKY genes expressed in leaves, buds and pistils, respectively, with at least one member of each group expressed in each tissue. In particular, the FPKM value of PdWRKY53 in buds was as high as 153.85. This finding indicates that this gene may play an important role in the development of almond buds. In addition, flowers expressed only seven PdWRKY genes as follows: Group I (four), Group IIc (two) and Group IId (one). Notably, 57 PdWRKYs were expressed in different tissues, indicating that PdWRKY genes have significant tissue-specific expression characteristics.

Expression analysis of PdWRKY family members

Ten PdWRKY genes were selected for qRT–PCR analysis to determine their expression levels in almond leaves at four time points, 0 (CK), 6, 12 and 24 h, under low-temperature, drought and salt stress (Fig. 8). After 6 h of low-temperature stress, five genes, PdWRKY1/7a/10/16/39, were highly expressed, especially PdWRKY39. At 12 h, only the expression of PdWRKY10 and PdWRKY25 was increased, while that of the other genes decreased, and at 24 h, four genes, PdWRKY1/16/25/54, were expressed at higher levels. When drought stress was applied for 6 h, only three genes, PdWRKY1/19/54, were highly expressed; at 12 h, four genes PdWRKY1/19/39/54, were highly expressed; and at 24 h, the expression of all examined genes was inhibited. After 6 h of salt stress treatment, only PdWRKY19 was highly expressed; at 12 h, three genes, PdWRKY1/10/19, were highly expressed; and at 24 h, the expression of all genes was inhibited. The results show that the members of the PdWRKY family have different temporal expression characteristics under different stresses.

PdWRKY protein interactions and GO annotations

The STRING database was used to predict potential interactions between PdWRKY proteins (Fig. 9). In the PdWRKY protein interaction network, 27 members connected to each other to form 13 nodes; some nodes have multiple PdWRKY members, and each node interacts with other nodes. Some proteins exhibited direct interactions, such as PdWRKY51 and PdWRKY44, while others exhibited more complex polygenic interactions, such as PdWRKY16/39, PdWRKY19/49, and PdWRKY43/48a/48b. Notably, multiple members, such as PdWRKY43/48a/48b and PdWRKY7a/7b/8/23, were predicted as central nodes and were connected to multiple genes. In addition, GO annotation and enrichment analyses of PdWRKY proteins were performed (Table S10). The results showed that nine binding-related protein functions were annotated in the biological process category, such as nucleic acid binding, DNA binding and organic cyclic compound binding. Forty-eight functions were enriched in the molecular function category and were mainly related to metabolism and gene expression, such as nitrogen compound metabolic process, biosynthetic process and gene expression.

Identification and annotation of PdWRKY target genes

To identify the possible downstream target genes regulated by the PdWRKYs and determine their functions, we searched the promoter sequence 2,000-bp upstream of the PdWRKYs using the conserved WRKY motif in the JASPAR database (Fig. S3). A total of 321 target genes were identified. Of these 321 genes, 128 had GO annotations, and 100 were mapped to the KEGG database (Table S11). GO analysis showed that multiple protein functions, such as metabolic process (GO:0008152), catalytic activity (GO:0003824) and catalytic activity (GO:0003824), were highly enriched in the target genes. The top 20 GO terms with the most target gene enrichment are shown in Fig. 10. In the cellular component category, most annotations were related to cell and membrane functions, including the terms cell (GO:0005623), intracellular (GO:0005622) and membrane (GO:0016020) (Fig. 10A). Binding (GO:0005488) function was the most enriched among the molecular function categories (Fig. 10B), while metabolic process (GO:0008152) was the most enriched among the biological process categories (Fig. 10C). Similarly, in the KEGG analysis, 86 target genes were annotated into metabolism-related pathways (Fig. 10D). The top 20 KEGG pathways for target gene enrichment included glycolysis/gluconeogenesis (ko00010), pentose and glucuronate interconversions (ko00040), fructose and mannose metabolism (ko00051), and alanine, aspartate and glutamate metabolism (ko00250). These results suggest that PdWRKY can affect multiple pathways by regulating target genes. In addition, more than 100 types of protein domains were found for the target genes, such as ribonucleases, UDP-glycosyltransferases and zinc fingers, indicating that WRKYs regulate a wide range of target genes to affect the growth and development of almond.