Genome-wide analysis of the WRKY gene family in drumstick (Moringa oleifera Lam.)

Sequence database searches

The complete genome and proteome sequences and General Feature Format (GFF) file for Arabidopsis were downloaded from TAIR (http://www.arabidopsis.org). The annotated drumstick genome sequences were provided by Yunnan Agricultural University. A WRKY-domain Hidden Markov Model (HMM) Profile, which was downloaded from Pfam (http://pfam.xfam.org/), was used as a query with which to search all of the annotated proteins in the drumstick genome with an E value cut-off of 1E−5. The candidates selected using HMMER were examined to determine whether they had typical features of WRKY proteins by employing the Pfam database. Finally, the CD-HIT program and the Pfam database were used to eliminate duplicate and incomplete sequences. Non-overlapping WRKY protein sequences were used for further analysis.

Multiple sequence alignment and phylogenetic analyses

The conserved WRKY domains of MoWRKY genes obtained using manual inspection in the Pfam program were aligned using ClustalX 1.83 software. Phylogenetic analysis including seven representative domains from Arabidopsis was carried out to obtain better classifications of the different clades by applying the Neighbor-Joining method with 1,000 bootstrap replicates using MEGA 6 software.

Gene structure and motif composition analysis

Analysis of the exon-intron organization of Mo WRKYs was performed by comparing the coding sequences of Mo WRKYs with their corresponding genomic sequences using GSDS software (Hu et al., 2015). Conserved motifs in each WRKY protein were investigated using the Multiple Expectation Maximization for Motif Elucidation (MEME) online program: http://meme-suite.org/. The following parameters were employed in analysis: maximum number of motifs 20; minimum motif width 6; maximum motif width 50.

Promoter cis-acting elements analysis of MoWRKYs

The promoter sequences, 1.5kb upstream of the translation start site, of the MoWRKY genes were obtained from drumstick genome. PlantCARE (Lescot et al., 2002) was used to analyse the MoWRKY gene promoters and identify their cis-acting elements.

Tests for selective pressure

The multiple sequence alignment of drumstick MoWRKY proteins was carried out using ClustalW with default parameters. Then the sequences were trimmed to reduce gap penality. DNAMAN was used to search for nucleotide sequences encoding additional WRKY proteins, with the aligned MoWRKY protein sequences as guides. The synonymous (K_S) and nonsynonymous (K_a) substitution rates were calculated with the YN00 program in PAML4.9 with default parameters (Yang, 2007).

Expression analysis

To investigate the patterns of expression of MoWRKY genes under normal and abiotic stress conditions, seedlings of drumstick were cultivated in potting soil at 25 °C under 14: 10 h light: dark conditions in a growth chamber for 20 days before treatment. For salt and oxidative stress treatments, seedlings were sprayed for 12 h with, respectively, 150 mM NaCl and H₂O₂ solution. Cold and heat stress were applied by transferring plants to a climate chamber at, respectively, 4 °C and 42 °C for 12 h. Drought stress was induced by withholding water for 2 weeks. Each treatment consisted of three replicates. After stress treatments, total RNA was isolated from leaf, stem, stem tip and root tissues of each seedling using a Total RNA Kit (OMEGA, Guangzhou, China). Total RNA was reverse transcribed into cDNA using a PrimeScript RT Master Mix (Perfect) Real Time Kit (Takara, Dalian, China). Gene specific primers were designed using Primer 5.0 and the RPL gene was used as a reference (Deng et al., 2016). Expression of all MoWRKY genes was examined by RT-PCR and products from each sample were analyzed using a 1% agarose gel. Among all MoWRKY genes, nine genes belonging to different subgroups were selected for analysis of gene expression levels using qRT-PCR according to the method described in Wei et al. (2016) and RPL was amplified as a reference gene (Deng et al., 2016). Relative expression levels were evaluated using the 2^−ΔΔCT. Three technical replicates were conducted for test and reference genes of each sample to obtain precise and reproducible results. Statistical analysis was carried out using SPSS 19.0 software (SPSS Inc., Chicago, IL, USA), Duncan’s multiple range test was used to detect differences among means. A p-value <0.05 was considered significant.

Identification of WRKY family members in the drumstick genome

To identify all the WRKY genes in the drumstick genome, we employed the HMM profile of the WRKY domain (PF03106) as a query to search against the drumstick genome database using HMMER 3.0 and BLAST. A total of 54 nonredundant genes (Table 1) were identified as WRKY genes and a unique name was assigned to each drumstick WRKY gene, consisting of two italic letters denoting the source organism and sequential numbers: MoWRKY1 to MoWRKY54. All the putative 54 WRKY genes were further analyzed to confirm the presence of the WRKY domain and all of them were annotated with gene ontology (GO) terms (File S1). Fifty-three MoWRKY genes containing complete WRKY domains were identified; only one gene (MoWRKY50) lacked a complete domain. The highly conserved domain WRKYGQK was present in 52 of the MoWRKY proteins, whereas the remaining one (MoWRKY24) contained a WRKYGKK domain. The lengths of the MoWRKY proteins ranged from 106 (MoWRKY24) to 834 (MoWRKY3) amino acids; the average length was 391 amino acids.

Phylogenetic relationship and classification of MoWRKY genes

The most prominent structural feature of WRKY genes is a conserved WRKY domain; there is also a zinc-finger motif. Among the 54 MoWRKY proteins identified, nine MoWRKY proteins contained two WRKY domains; since one MoWRKY protein did not have a complete WRKY domain, a total of 62 WRKY domains were found in this study. In each protein that contained two WRKY domains, we designated these domains by the WRKY name plus N or C for the N-terminal or C-terminal domain respectively. In order to examine phylogenetic relationships and classify all 62 MoWRKY domains, a phylogenetic tree based on conserved WRKY domains was constructed. Representative WRKY domains from Arabidopsis were used in our analysis, and the candidate domains were obtained from Diao et al. (2016) and Li et al. (2016). Figure 1 shows a multiple sequence alignment of the 62 WRKY domains. Three major groups were identified, as previously described in poplar (He et al., 2012), pepper (Diao et al., 2016) and carrot (Li et al., 2016). Additionally, several subgroups were apparent on the basis of the phylogenetic analysis.

Group I contained 10 WRKY proteins, of which all contain two WRKY domains except for MoWRKY10. This member might have lost the N-terminal WRKY domain during evolution, since its single WRKY domain showed high similarity to MoWRKY1C, which is located in the C-terminal WRKY domain clade, suggesting a common origin for these two domains. Group II had the largest numbers of WRKY proteins and was divided into five major subgroups: IIa, IIb, IIc, IId and IIe. Subgroup IIa (three members) and IIb (eight members) were two subgroups in the same branch, while subgroup IId (five members) and IIe (seven members) were derived from one clade. Subgroup IIc, with 14 members, was more similar to group I than to any other subgroups according to the phylogenetic analysis. Furthermore, six WRKY domains belonged to group III, which is widely considered to be the most advanced in terms of evolution and the most relevant to adaptability (Dou et al., 2016; Kalde et al., 2003; Huang et al., 2016). Comparing the two phylogenetic trees, constructed for MoWRKY domains and genes, similar groups and subgroups were identified, though the classifications of a few members were different (Fig. 1 and Fig. S1), indicating that the conserved WRKY domain is an important unit in WRKY proteins.

Structure analysis of MoWRKY genes

Intron/exon organization and numbers of introns are typical imprints of evolution within some gene families. In this study, we analysed the structure of MoWRKY genes to gain further insight into evolutionary events that had shaped them and found that all MoWRKY genes contain introns (Fig. 2A). The number of introns varies among genes, with the minimum, one intron, identified in five MoWRKY s (MoWRKY50, MoWRKY44, MoWRKY51, MoWRKY19 and MoWRKY11) of subgroup IIc and the maximum, 10 introns, being present in MoWRKY22. Gene structure analysis revealed that genes with similar structures always clustered in the same class. For example, six members of group III all contained three exons and two introns. Similarly, five exons and four introns were present in MoWRKY2, MoWRKY3, MoWRKY4, MoWRKY5 and MoWRKY6, which belonged to group I. However, the other five MoWRKY s in group I exhibited different gene structures.

Motif composition analysis of MoWRKY proteins

The conserved motifs of WRKY proteins in drumstick were investigated using the MEME online software suite (http://meme-suite.org/) to better understand the similarity and diversity of motif compositions. Twenty distinct motifs were identified and a schematic overview of these motifs is provided in Fig. 3. For MoWRKY proteins, motif 1 was broadly distributed in all MoWRKY proteins, which was corresponded to WRKY domain. Motif 3 was only detected in the type I group. Motifs 5, 6, 7, 8, 9, 16, 17, 18 and 20 were only detected in the type II group; among them, motifs 5 and 6 were only detected in subgroup IIa and IIb, motifs 8 and 9 were only detected in subgroup IIb, motif 17 was only detected in subgroup IIc, and motif 16 was only detected in subgroup IId. Motifs 12 and 15 were only detected in the type III group. Generally, proteins with similar motif compositions were clustered in the same class indicating that members of the same class may have similar functions.

Rapid expansion of group III WRKY genes in land plants

Group III WRKY genes have only been characterized in flowering plants, and a large number of duplications and diversifications in this group appear to have resulted from different selection challenges (Dou et al., 2016; Huang et al., 2016; Kalde et al., 2003). To explore the evolutionary relationships of group III WRKY genes across drumstick and other land plant species, we performed a multiple sequence alignment among the 81 group III WRKY proteins from drumstick and another seven species. A phylogenetic tree was constructed from the results of the alignment using the neighbor-joining method (Fig. 4). The marked difference in group III WRKY gene size among different species suggests that group III WRKY gene expansion occurred after the divergence of monocotyledons and dicotyledons. MoWRKY clearly shared more sequence similarity with VvWRKY and PaWRKY than with other WRKYs.

The cis-acting elements analysis of MoWRKYs

For further understand the possible functions of MoWRKY genes, the cis-acting elements in all MoWRKY genes promoters were analyzed using PlantCARE software based the drumstick genome data. Various types of cis-acting elements were found and all MoWRKY genes contained several cis-acting elements in their promoter regions. The 10 most common elements were summarized in Table 2. These elements included three hormone responsive elements (ABRE, CGTCA motif and TGACG motif), an essential element for the anaerobic inductio (ARE), a drought stress responsive element (MBS), a heat stress responsive element (HSE) and four light responsive elements (Sp1, Box 4, G box and GT1 motif).

Divergence in selective pressure between subgroups

The ratio (ω) of the non-synonymous substitution rate (Ka) to the synonymous substitution rate (K_S) provides a sensitive measure of selective pressure acting on a protein-coding gene. Homologous genes with ω ratios of 1, <1, or >1 are usually assumed to be evolving under neutral evolution, purifying selection, or positive selection, respectively. To test for deviations in the substitution rates of MoWRKY genes, we calculated ω values across all pairwise comparisons within the 54 WRKY genes using the YN00 program in the PAML software package. The frequency distribution of ω values is shown in Fig. 5A. The results suggested that the WRKY gene family evolved mainly under strong purifying selection. However, there are several ω values greater than 1, such as those for the comparison between MoWRKY8 and MoWRKY9 and that between MoWRKY8 and MoWRKY10, indicating that positive selection acted on these genes. Only 0.5% of the ω values approximated to 1, indicating that no selective pressure acted on these genes.

To test whether the rate of evolution among the subgroups of WRKY genes was identical, we calculated ω values across all pairwise comparisons within each of the subgroups; the results are shown in Fig. 5B. The average ω values of each subgroup were different. In order (highest first) they were: IIc, III, I, IIe, IId, IIb and IIa, indicating that different subgroups had evolved at different rates and that IIc had evolved the fastest.

Expression patterns of WRKY genes in drumstick under normal growth conditions and abiotic stress conditions

To investigate the responses of MoWRKY genes to stresses, we examined the expression patterns of all 54 full-length MoWRKY s under normal growth conditions and under five abiotic stresses (heat, cold, drought, salt and oxidative) in different tissues (leaves, roots, stems, stem apex) using RT-PCR. As shown in Fig. 2B, among the 54 MoWRKY genes, 13 genes were expressed in all tissues under all growth conditions. In contrast, six genes, including MoWRKY24, the only gene with a variant WRKY domain (WRKYGKK), were not expressed in any tissue or in response to any of the treatments applied in this study. Thus, these six WRKY genes are expressed at undetectably low levels, or they are only induced in response to treatments and/or in tissues not examined in our study, or they are pseudogenes. The other 35 WRKY genes were expressed selectively in a specific tissue and/or in response to a specific treatment. Six of these genes were not expressed in any tissue under normal growth conditions but were expressed under stress conditions, suggesting that they play specific roles during stress conditions. At the same time, some genes, such as MoWRKY46, were only expressed in specific tissues under normal growth conditions but were expressed in all tissues under certain stress conditions, indicating that these genes may also play specific roles under stress conditions.

Nine MoWRKY genes from different subgroups were selected and their expression profiles were analyzed in root tissue under normal growth conditions and five abiotic stresses using qRT-PCR. As shown in Fig. 6, these selected MoWRKY genes were sensitive to abiotic stresses. All nine exhibited a high level of transcript accumulation under cold stress, especially MoWRKY30 (GO: 0006950), followed by MoWRKY54 (GO: 0006950 and GO: 0080134). Interestingly, the genes that were most strongly up-regulated under cold treatment were always up-regulated in response to heat and salt treatments. In drought stress, MoWRKY22, which had the most introns and MoWRKY3, which was the longest in MoWRKY gene family were found to be slightly upregulated, whereas weak expression were found for the other seven genes. The expression levels of almost all the nine MoWRKY s were decreased under oxidative stress. MoWRKY49, MoWRKY53 and MoWRKY54, which all belonged to group III, have similar gene structures and the same motifs. But the expression levels of the three genes under abiotic stresses were slightly different. They were evidently upregulated in cold and salt to different degrees; MoWRKY53 (GO: 0006950 and GO: 0080134) and MoWRKY54 (GO: 0006950 and GO: 0080134) were also responsive to heat. Overall, the expression patterns of MoWRKY s under various conditions suggest that different MoWRKY genes may be involved in different signaling and stress responses, and that an individual MoWRKY gene can also participate in multiple signaling and stress process.