Characterization of 54 WRKY-encoding sequences in lotus.
Accumulating evidence suggested the important roles of WRKY genes in modulating the expression of genes related to plant’s stimuli response and development processes [45]. For instance, in soybean and Arabidopsis, studies suggested that GmWRKY21, 31, and 54, and AtWKRY25, and 33 confer differential tolerance to stressful environment respectively [46]. OsWRKY78 are reported to regulates stem elongation and seed development in rice [47]. However, WRKY gene family in lotus has so far been poorly described. Completion of the lotus genome by Ming et al.[48] makes it possible to carry out genome-wide study on WRKY gene family of lotus. WRKY proteins displayed highly conservation in thier WRKY domain [49] across different plant species. Thus, HMMER searches were performed to search local lotus protein data with HMM profile of conserved WRKY domains as queries [50]. A total of 54 lotus WRKYs (NnWRKYs) genes which contained complete WKRY domain were identified in lotus genome (Table 1). Using identified genes as queries, we further performed a homology search from amino acid sequences of lotus to identify new potentially WKRY encoding genes, which however resulted in no new WKRY genes. Domain analysis confirmed the presence of one or two WKRY domains in the predicted WKRY proteins. The properties of identified NnWRKY proteins were characterized and listed in Table 1, including protein length, molecular weight, isoelectric point, and motifs. The identified NnWRKYs contained a diverse number of amino acid from 106 (NNU_023323) to 897 (NNU_008419) in length with an average length of 407. The protein weight ranged from 12.21 (NNU_023323) to 98.71 (NNU_008419) kDa. The isoelectric point (PI) varied from 4.57 (NNU_021524) to 9.80 (NNU_020936), with 32 numbers have a pI < 7 and others have a pI > 7. All NnWRKY genes were found to scattered in 14 of 50 scaffolds. We found that most NnWRKY genes were located in the scaffold 1. This uneven distribution of NnWKRY genes on lotus scaffolds might be related to the scaffold length. The number of WKRY gene family were relative variable among different species, i.e. there were 102 WRKY genes in Oryza sativa L. ssp. Japonica [51], 116 in Gossypium raimondii [52], 37 in Physcomitrella Patens [49], which indicating that duplicated copies of WKRY genes are frequently retained after a whole-genome duplication event.
Table 1
WRKY transcript factor family in lotus.
NnWRKY (Genebank accession) | Deduced polypeptide | WRKY domain | Megascaffold |
Length (aa) | MW (kDa) | pI | Zinc- Finger Type | Domain Number |
NNU_005293 | 743 | 80.25 | 5.6 | C2H2 | 2 | 1 |
NNU_005576 | 461 | 50.53 | 9.1 | C2H2 | 2 | 3 |
NNU_009557 | 473 | 52.36 | 8.99 | C2H2 | 2 | 3 |
NNU_013081 | 565 | 61.64 | 6.74 | C2H2 | 2 | 3 |
NNU_019405 | 608 | 66.79 | 5.64 | C2H2 | 2 | 2 |
NNU_000906 | 538 | 58.50 | 7.41 | C2H2 | 2 | 4 |
NNU_008419 | 897 | 98.71 | 5.95 | C2H2 | 2 | 5 |
NNU_013102 | 629 | 69.51 | 7.66 | C2H2 | 2 | 3 |
NNU_014662 | 785 | 84.71 | 5.82 | C2H2 | 2 | 18 |
NNU_007725 | 588 | 64.27 | 6 | C2H2 | 1 | 3 |
NNU_012242 | 693 | 75.89 | 8.6 | C2H2 | 1 | 3 |
NNU_024662 | 496 | 55.04 | 5.42 | C2H2 | 1 | 10 |
NNU_005172 | 307 | 34.02 | 6.68 | C2H2 | 1 | 1 |
NNU_015799 | 493 | 55.50 | 6.42 | C2H2 | 1 | 1 |
NNU_002747 | 303 | 32.89 | 4.98 | C2H2 | 1 | 2 |
NNU_022759 | 296 | 32.11 | 6.06 | C2H2 | 1 | 9 |
NNU_019387 | 313 | 34.83 | 7.19 | C2H2 | 1 | 2 |
NNU_000764 | 576 | 62.37 | 6.22 | C2H2 | 1 | 4 |
NNU_022872 | 622 | 66.90 | 5.9 | C2H2 | 1 | 5 |
NNU_005007 | 613 | 66.19 | 6.35 | C2H2 | 1 | 1 |
NNU_024102 | 685 | 74.89 | 6.26 | C2H2 | 1 | 1 |
NNU_021907 | 497 | 54.44 | 7.31 | C2H2 | 1 | 33 |
NNU_025395 | 503 | 56.43 | 7.34 | C2H2 | 1 | 1 |
NNU_024663 | 526 | 57.63 | 5.86 | C2H2 | 1 | 10 |
NNU_000937 | 311 | 34.28 | 6.39 | C2H2 | 1 | 4 |
NNU_008364 | 308 | 33.93 | 6.78 | C2H2 | 1 | 5 |
NNU_001380 | 321 | 35.45 | 8.74 | C2H2 | 1 | 13 |
NNU_023785 | 327 | 36.04 | 7.35 | C2H2 | 1 | 8 |
NNU_009014 | 299 | 33.30 | 5.84 | C2H2 | 1 | 9 |
NNU_014445 | 310 | 34.51 | 6.21 | C2H2 | 1 | 2 |
NNU_010128 | 128 | 14.55 | 9.4 | C2H2 | 1 | 6 |
NNU_023323 | 106 | 12.21 | 9.22 | C2H2 | 1 | 1 |
NNU_003142 | 184 | 20.96 | 9.37 | C2H2 | 1 | 5 |
NNU_004050 | 183 | 20.82 | 9.49 | C2H2 | 1 | 4 |
NNU_010326 | 196 | 22.17 | 9.4 | C2H2 | 1 | 1 |
NNU_007649 | 303 | 33.67 | 6.21 | C2H2 | 1 | 3 |
NNU_011033 | 302 | 34.01 | 5.72 | C2H2 | 1 | 3 |
NNU_008704 | 216 | 24.62 | 8.19 | C2H2 | 1 | 2 |
NNU_025262 | 222 | 25.41 | 7.77 | C2H2 | 1 | 8 |
NNU_010097 | 223 | 25.88 | 9.19 | C2H2 | 1 | 6 |
NNU_014805 | 158 | 18.97 | 9.65 | C2H2 | 1 | 1 |
NNU_009041 | 282 | 32.07 | 4.94 | C2H2 | 1 | 9 |
NNU_014473 | 281 | 32.00 | 5.67 | C2H2 | 1 | 2 |
NNU_010321 | 365 | 40.11 | 5.07 | C2H2 | 1 | 1 |
NNU_021524 | 288 | 31.65 | 4.57 | C2H2 | 1 | 15 |
NNU_011577 | 358 | 39.41 | 9.51 | C2H2 | 1 | 1 |
NNU_020936 | 360 | 40.55 | 9.8 | C2H2 | 1 | 1 |
NNU_004768 | 326 | 36.50 | 6.14 | C2HC | 1 | 15 |
NNU_005350 | 329 | 37.06 | 5.46 | C2HC | 1 | 3 |
NNU_004769 | 355 | 38.99 | 5.37 | C2HC | 1 | 15 |
NNU_007278 | 686 | 78.85 | 8.72 | C2HC | 1 | 16 |
NNU_010416 | 362 | 40.55 | 6.24 | C2HC | 1 | 1 |
NNU_013507 | 363 | 40.37 | 5.27 | C2HC | 1 | 4 |
NNU_013509 | 362 | 40.56 | 5.51 | C2HC | 1 | 4 |
Phylogenetic analysis of lotus WKRY genes.
To study the evolutionary relationships between lotus WRKY genes, a phylogenetic tree of 54 NnWRKY genes was created based on multiple sequence alignment of nucleotide sequence data [53](Huelsenbeck et al., 2001). Using Modeltest 3.7 [32], we first test the best-fitting nucleotide substitution model for phylogenetic tree construction using Bayesian method, and GTR + G + I were selected. Based on the statistical support values and previous works [7, 10], 54 NnWKRY genes were classified into 5 separate monophyletic groups, termed WRKY I (14), WRKY IIa (10), WRKY IIb (6), WKRY IIc (17) and WRKY III (7) (Fig. 1). This classification was also consistent with the previous described classification based on the presence of WKRY and zinc finger domains [51, 52, 54]. For example, the members that contain two WRKY domains and feature the C2H2-type zinc finger motif (N: Cx4Cx22-23 HxH; C: Cx4Cx23HxH) were categorized into the group I. However, there were also some exceptions, i.e. some NnWKRY proteins in group I which carries one WRKY domain appeared to be more related to the group. The same phenomenon was also found in Arabidopsis [7], which might attribute to the loss of the WRKY domain, as domain gain or loss appear to have shaped the WKRY family [55]. These findings revealed the existence of five variants in NnWKRY genes recovery as five distinct clades in the Bayesian reconstruction.
Gene structure analysis of NnWRKYs.
As gene structures were an important part in assessing the evolution of NnWRKY gene family, we then examined the intron/exon distribution of NnWRKY genes [17, 19, 56]. In spite of the highly conserved WKRY domains, the length of these full sequences is highly variable and even more distinct in their gene sequences (Fig. 2). All NnWRKY genes contain at least one intron in their CDSs: 6 have one intron; 23 have two introns; 5 have three introns; 11 have four introns; 6 have five introns. Besides, NNU_013102, NNU_024102, and NNU_008419 were more complicated, which have 6, 8 and 10 introns respectively (Fig. 2). Introns were particularly abundant in members of WRKY IIa, which has an average intron number of 4.9, followed by it was WRKY I which have 4. In contrast, WRKY IIb, IIc and III only have an average number of introns approximate to 2 (Fig. 3). During the evolution of intron, intron loss was more frequently than intron gain following segmental duplication [57~59], which indicated that WRKY IIa might be the most ancestral that other NnWRKY genes might derive from it.
Conserved sequences in WRKY proteins.
A multiple sequence alignment of 54 full-length NnWKRY sequences showed that nucleotide identity among them varied from 23.53% (between NNU_016707 and NNU_022788) to 91.04% (Fig. 4). Their full length amino acid sequence was even more variable, ranged from 4.35 to 91.35%. However, despite their full gene or protein sequences were highly distinct, several residues in the WRKY domain were conserved in all NnWRKY genes, implied that corresponding genes encoding active enzymes. In contrast, amino acid sequences outside of the WKRY domain were the most divergence sequence.
We then calculated ratio (ω) of the synonymous substitution rate (dS) versus the non-synonymous substitution rate (dN) among different regions of WKRY protein, including full-length sequence, C-terminal WKRY domain, N-terminal WKRY domain and sequence outside WKRY domain. It was found that dN/dS ratios were always less than 1 for both WKRY domains, suggesting purifying selection on these regions which could maintain the function conservation of WRKY genes [60, 61]. In contrast, the ratio was much higher in regions outside WKRY domain with most dN/dS ratios higher than 1 (Fig. 5), implying the relaxed purifying or positive selection on this region [60, 61].
Further analysis of the sequence feature of NnWRKY proteins using software MEME [62] with default settings, predicted 20 conserved motifs. We listed the sequence of these motifs (Fig. S1) and visualized their distribution in Fig. 6. As shown in Fig. 6, motifs 1, 5, and 14 contained the characterized sequence (WRKYGQK) of WRKY domain, and at least one of them are presented in every sequence. Besides, one or more conservation motifs outside of the WKRY domain can be detected in a WKRY protein: Motifs 4 and 15 were found in most members of WRKY I; motif 9 was found in the WRKY IIa; motif 15 was identified in most members of WRKY IIc; motif 3 was unique in WRKY III. Generally, WRKY motif compositions were more similar between members in the same group, indicating functional similarity among them.
Analysis of functional divergence.
Functional divergence after gene duplication has been known as an important source of evolutionary novelty [63]. After gene duplication, one gene copy is free to change as it will not cause function deficit of the gene, and accumulating amino acid changes will contribute to functional divergence. To investigate whether amino acid substitutions in the NnWKRY sequences could have caused adaptive functional diversification, type I and type II functional divergence between groups of the WRKY gene family were estimated by posterior analysis [33]. Type I functional divergence represents altered functional constraints (e.g. different evolutionary rates) between duplicate genes, while type II indicate no altered functional constraints but radical change in amino acid (i.e., charge, hydrophobic) [33]. Our results showed that type I coefficients (θI) of functional divergence in WRKY I/WRKY IIb and WRKY IIa/WRKY IIb were significantly greater than zero (P < 0.05). Meanwhile, we did not find any evidence for type II functional divergence between groups (Table 2), indicated that type I functional divergence was the dominant pattern for the evolution of WRKY family in lotus. Further, we predicted amino acid residues which are potentially responsible for the functional divergence on the base of sites specific profiles. To greatly reduce false positive, a cut-off value, Q(k) ≥ 0.95, was used empirically as cutoffs for estimated key amino acid residues. Six amino acid residues (875, 878, 879, 881, 912 and 914) in all comparisons were identified as being most important for the functional divergence (Fig. 7). These amino acids were all identified in the WRKY domain and presumably have contributed to the functional divergence of lotus WKRY genes.
Table 2
Functional divergence between groups of the NnWRKY genes.
Type I | Type II |
| ΘI ± SEz | LRT | P | Qk≥0.80y | Qk≥0.95x | ΘII ± SE | P |
WRKY I/WRKY IIa | 0.354 ± 0.248 | 2.038 | 0.153 | 1 | 0 | -0.939 ± 1.240 | 0.449 |
WRKY I/WRKY III | 0.388 ± 0.379 | 1.051 | 0.305 | 0 | 0 | -1.103 ± 2.005 | 0.582 |
WRKY I/WRKY IIb | 0.55 ± 0.277 | 3.936 | 0.047 | 1 | 1 | -2.390 ± 3.233 | 0.46 |
WRKY I/WRKY IIc | 0.001 ± 0.022 | 0 | 0.964 | 0 | 0 | -1.028 ± 1.297 | 0.428 |
WRKY IIa/WRKY III | 0.762 ± 0.47 | 2.628 | 0.105 | 3 | 1 | -0.423 ± 1.357 | 0.755 |
WRKY IIa/WRKY IIb | 0.942 ± 0.463 | 4.148 | 0.042 | 11 | 5 | -1.613 ± 2.491 | 0.483 |
WRKY IIa/WRKY IIc | 0.289 ± 0.225 | 1.654 | 0.199 | 1 | 0 | -0.845 ± 0.908 | 0.352 |
WRKY III/WRKY IIb | 0.001 ± 0.022 | 0 | 0.964 | 0 | 0 | nan | nan |
WRKY III/WRKY IIc | 0.001 ± 0.022 | 0 | 0.964 | 0 | 0 | -0.297 ± 1.23 | 0.81 |
WRKY IIb/WRKY IIc | 0.685 ± 0.319 | 4.6 | 0.968 | 3 | 0 | -1.466 ± 2.351 | 0.533 |
zThe coefficient of functional divergence between two subgroups and its standard error. |
yThe number of critical amino acid residues with posterior prob-ability (Qk) > 0.80. |
xThe number of critical amino acid residues with posterior probability (Qk) > 0.95. |
Detection of positive selection and gene conversion events.
ML methods and two pairs of codon substitution models in PAML were chosen to test positive selection using likelihood ratio test (LRT) [38, 63, 64]. First, we compared one ratio model (M0) with the discrete model (M3) to evaluate the variations in dN /dS ratio among codon positions below each group. Overall, testing the selective pressure with the basic model M0 showed that the ML estimates of dN/dS values for all groups were close to zero, with a mean of 0.281, indicating that each group of NnWRKY genes have been under purifying selection (Table 3). Interestingly, the log likelihood differences between model M3 and M0 were statistically significant for all groups, revealing that the evolutionary forces acting on individual codon site are fluctuated and some genes might have been under positive selection. Then, the presence of codons that evolved under positive selection was further tested by comparing the M8 and M7 models. According to the criterion described previously[54] group IIa was believed to have undergone positive selection. We identified the positively selected sites under model M8 using the Bayesian method, and one positive selection sites (625G, corresponding to NNU_013102) were identified [65, 64].
Table 3
Detection of positive selection among codons of NnWKRY under site-specific model.
Group | Nz | dN/dS(ω) under M0y | 2△/M3 vs. M0x | 2△/M8 vs. M7x | M8 Estimatesw | Positively Selected Sitesv |
WRKY I | 14 | 0.356 | 477.1361** | 11.39814** | p1 = 0.05110 ω = 2.92059 p = 0.87912 q = 1.09104 | - |
WRKY IIa | 10 | 0.31614 | 1269.635** | 15.37348** | p1 = 0.03044,ω = 44.74768, p = 0.56713,q = 0.65863 | 625 G (NNU_013102) |
WRKY IIb | 6 | 0.16165 | 300.6672** | 0.00018 | p1 = 0.00001,ω = 2.03772, p = 0.71318,q = 2.67776 | - |
WRKY IIc | 17 | 0.27169 | 1358.758** | 177.983** | p1 = 0.00001,ω = 1.00000, p = 0.42827,q = 0.55990 | - |
WRKY III | 7 | 0.30279 | 228.188** | 2.541058 | p1 = 0.23920,ω = 1.00000, p = 0.77790,q = 1.87365 | -- |
zNumber of the sequences in the group. |
ydN/dS (ω) is the average ratio over sites under a codon model with one-ratio. |
xDouble asterisks represent P < 0.01 for chi-square (χ2) test, while single asterisk represents P < 0.05. |
wω is estimated under model M8; p1 is the inferred proportion of positively selected sites; p and q are the parameters of the beta distribution. |
vThe amino acid, which is listed in positively selected position, had a posterior probability > 0.95 under M8, df = 1. |
Gene conversion is an important population genetic process which also acted as an important force in evolution [66]. Possible gene conversion events during evolution of NnWRKY genes were analyzed using the program GeneConv [40]. Possible gene conversion event between ancestors of two sequences in the alignment and originated from outside of the alignment were all included in this analysis. The results demonstrated that WRKY I, IIa, IIc and III possessed at least one gene conversion event (Table 4). No gene conversion event was identified from the genes in WRKY IIb group. Most of the gene conversion events (47%) were identified within WRKY I group and they can be detected in 9 members of this group. All these results suggested that gene conversion events might have played significant role during the evolutionary history of lotus WRKY family, which is also found in WKRYs of other species as showed previously [10].
Gene expression analysis.
Previous work demonstrated that WRKY genes serve important functions in development and stress responses of plants [49]. To gain more information on the role of lotus WRKY gene family, in this study, available transcriptome data of lotus, related to different tissues, stresses, and developmental stages, were used to analyze expression patterns of WRKY genes [43, 44, 67, 68].
NnWKRY genes showed differential expression level in various lotus tissues (Fig. 8). Except NNU_004769 which was not detected in the petiole, no more tissue-specific WRKY genes were found. Among them, NNU_013509 showed the highest expression level while NNU_004769, NNU_025395, NNU_024662, and NNU_011033 showed the lowest transcript level (with average FPKM < 1) in most examined tissues. Genes from the same group frequently showed similar expression patterns within specific tissues. For example, NNU_005293, NNU_005576, and NNU_009557, which belonging to WRKY Ia, had relative lower expression level in old leaf and petiole (average FPKM < 10) and high expression level in root (average FPKM = 39.60) and rhizome (average FPKM = 72.79). These results indicated that NnWRKY genes might primarily have functions in sink tissues.
In response to heat, submergence and copper stresses, the total transcript abundance of NnWRKYs did not show apparent difference among different stresses, while Cu treatment induced significantly higher abundance of NnWRKY genes than other treatments (p < 0.05) (Fig. 9). In response to heat and submergence stresses, WRKY genes in lotus are tend to decrease, as evidenced by slightly decreased abundance of NnWRKY genes. In response to heat treatment, most NnWRKY genes were down-regulated. Some genes, i.e. NNU_019387, NNU_010416, and NNU_013507 were down-regulated by heat for over 2 folds. Only 14 of 54 NnWRKY genes were up-regulated by heat stress, with NNU_010326 being the highest upregulated gene (log2(heat/ck) = 1.72). Two NnWRKY genes (NNU_000937 and NNU_005007) were mostly induced (log2(heat/ck) > 2) by submergence treatment, while most NnWRKY genes (i.e. NNU_005172 and NNU_008364) were down-regulated. Transcript level of NNU_005172 was reduced to less of 14% of normal condition. It is interesting to note that the homologous gene of NNU_005172 in soybean, WKRY 20, plays important role in drought stress [69, 70]. Overexpression of soybean WKRY 20 in alfalfa and soybean could enhance drought tolerance of the plant, possibly by regulating ABA signaling [69, 70]. These results indicate that WKRY 20 might played as a switch to modulate plant response to different kinds of water stress via regulating ABA signaling. Some metabolic processes of plants use copper as an essential micronutrient, which acting as a constituent of enzymes and other proteins [71~75]. Excessive copper has been found to cause oxidation of biomolecules and damage cell membrane [76]. Under copper stress, the transcript abundance of most NnWRKY genes was increased. Numerous works reported the roles of WKRY gene family members in enhancing plant oxidative stress tolerance [77]. WRKY proteins domain have W-box binding site to bind W-box, which is the key regulatory element involving in stress response [1]. Therefore, these genes may be the important factors in the copper stress metabolic network of lotus.
Transcriptional expression of NnWRKY during lotus rhizome and flower bud development were also investigated (Fig. 10). We did not find significant difference among the overall abundance of NnWRKY genes in lotus during various development stages. During different development stage of rhizome, 24 NnWRKY genes showed obvious difference (fold change > 2 or specific expressed) from stage 1 to stage 3. Among them, 3 and 13 genes were down-regulated during the stage 2 and 3 in comparison with stage 1 respectively. One and 6 NnWRKY genes showed specific expression pattern, that 1 and 5 genes were not detected in stage 2 and 3 respectively. Only two genes in stage 3 were up-regulated for over 2 fold and 3 genes were expressed weakly in stages 2 and 3 (FPKM < 1). During development of lotus flower bud, 4 genes were obviously changed (log2fold > 2 or specifically expressed) in compared with stage when the flower bud was not started to form. The expression of NNU_008419 and NNU_024102 decreased for over 2 fold during flower bud development. NNU_013081 was not detected during flower develop stage. Only one gene (NNU_014445) was found to be significantly induced during flower bud development, which although expressed at very low abundance (FPKM < 1)
The gene expression patterns above pointed to potential co-regulated relationships between various stresses and development processes. Thus, we constructed the co-regulated network of NnWKRY to various stimuli using the transcriptome datasets described above (Fig. 11). The results showed that stresses up-regulated WKRY genes were commonly down-regulated by various development processes, i.e. NNU_005007, NNU_024102 and NNU_013509. The contrasting expression pattern under stresses and development processes indicated that these WKRY genes might play important roles in lotus energy metabolism and participate in stress response. GO analysis confirmed this conclusion, showed that these genes were mainly involved in plant stresses responses (Fig. S2). It is interesting to note that heat stresses downregulated a series of WRKY genes as development did. We speculated that heat stress at 37°C might have promoting effect on lotus development in short time. Some of these obviously changed genes in the network were then validated by semi-quantitative RT-PCR. As shown in Fig. S3, RT-PCR results of genes under different stresses were similar in magnitude to those obtained by deep sequencing, suggesting that the gene expression data generated from the RNA-Seq data applied in this study were reliable. The gene expression patterns of NnWKRY genes provides important clues for investigating gene function of lotus WRKY gene family.