Comparative genome-wide analysis of WRKY transcription factors in two Asian legume crops: Adzuki bean and Mung bean

The seminal participation of WRKY transcription factors in plant development, metabolism and in the governance of defense mechanism implicated their gaining importance for genomic and functional studies. The recent release of draft genome sequences of two legume crops, Adzuki bean (Vigna angularis) and Mung bean (Vigna radiata) has paved the way for characterization of WRKY gene family in these crops. We found 84 WRKY genes in Adzuki bean (VaWRKY) and 85 WRKY genes in Mung bean (VrWRKY). Based on the phylogenetic analysis, VaWRKY genes were classified into three groups with 15 members in Group I, 56 members in Group II, and 13 members in Group III, which was comparable to VrWRKY distribution in Mung bean, 16, 56 and 13 members in Group I, II and III, respectively. The few tandem and segmental duplication events suggested that recent duplication plays no prominent role in the expansion VaWRKY and VrWRKY genes. The illustration of gene-structure and their encoded protein-domains further revealed the nature of WRKY proteins. Moreover, the identification of abiotic or biotic stress-responsive cis-regulatory elements in the promoter regions of some WRKY genes provides fundamental insights for their further implementation in stress-tolerance and genetic improvement of agronomic traits.

, respectively, has paved the way for the genomic and functional studies of these crops. Plants employ diverse approaches to acclimatize to environmental cues. Tuning the expression of genes deploying transcriptional regulators in respond to disparate environmental and physiological stimuli is one of the ways employed by eukaryotic organisms 16,17 . Plants dedicate approximately 7% of its genome to encode vital transcriptional regulators called transcription factors (TF), which facilitates sequence-specific binding to candidate genes through their conserved DNA-binding domains 18 . Adzuki bean genome and Mung bean genome allocates 2269 genes and 1850 genes encoding transcription factors, respectively 14,15 . Despite of their significance, none of the transcription factor family in these legumes crops has been explored yet.
WRKY transcription factor family is one such family known to endow their contribution in diverse biological phenomena in plants including defense responses 19 . Initially WRKY TFs were thought to be involved particularly, in plant-pathogen interaction, but recent functional studies implicate their gaining importance in abiotic stress responses as well 20 . WRKY TFs can act as positive as well as negative regulators of stress responses. For instance, one of the soy bean (Glycine max) WRKY genes, GmWRKY21 impart tolerance against cold stress in Arabidopsis thaliana, while GmWRKY54 confers salt and drought tolerance 21 . In rice (Oryza sativa), OsWRKY72 confers tolerance towards high salt and drought via ABA signaling 22 . OsWRKY8 provides osmotic stress tolerance 23 . OsWRKY11 induced by heat provides heat and drought tolerance 24 . AtWRKY26 is activated by a heat-induced ethylene-dependent response imparts heat resistance 25 . On the other hand, some WRKY proteins confer sensitivity in plants towards stress. Like, AtWRKY2 induced by NaCl and mannitol and AtWRKY40 induced by ABA have been proposed to act as a negative regulator of ABA-induced seed dormancy 26,27 . Similarly, overexpression of AtWRKY18 and AtWRKY60 make the plant more susceptible towards salt and osmotic stresses 28 . Moreover, a crosstalk between biotic and abiotic stress response components is reported to exist. Apart from its previously mentioned role, AtWRKY18 is induced by ABA as well as salicylic acid and pathogens 26,29 . OsWRKY89 induced by salinity, ABA and UV-B and also confers disease resistance 26,30 . AtWRKY53 which is a negative regulator drought tolerance, involved in senescence, also plays role in basal defense in plants [31][32][33] . AtWRKY70 is involved in disease resistance, basal defense as well as in senescence 34,35 . OsWRKY45 from rice enhances disease resistance, but also provides drought and cold tolerance in Arabidopsis 36,37 .
The WRKY family proteins contains one or two highly conserved domains of around 60 amino acids, comprising of a hallmark heptapeptide WRKYGQK at N-terminal, and a signature C-terminal zinc-finger motif 38,39 . They bind to the W box and a sugar-responsive cis-element SURE in the promoter of the target gene 40 . All known members of WRKY family can be classified into three groups viz. group I, II, and III, based on the number of occurrence of WRKY domains and the types of zinc finger motif. In the group I proteins, two WRKY domains can be found, whereas the group II and the group III possess only a single domain 38,39 . Group II is further divided into several subgroups based on the phylogenetic clades 19 .
In this study, we conducted a genome-wide survey to identify and designate nomenclature to WRKY TF family in Adzuki bean (VaWRKY) and Mung bean (VrWRKY) by implementing bioinformatics approach on the publicly available database of sequenced genome. We identified a comprehensive and non-redundant set of 84 VaWRKY genes and 85 VrWRKY genes encoding the WRKY transcription factor family in Adzuki bean and Mung bean, respectively and manually curated them. Subsequently, gene classification, exon-intron organization, chromosome distribution, gene duplication events, phylogenetic relationships, and conserved motifs were also studied, which lay a firm foundation for further comparative genomics studies. We also carried out promoter analysis to investigate stress-responsive cis-regulatory elements in 17 VaWRKY genes and 18 VrWRKY genes orthologous to strongly reported stress-responsive WRKY proteins of Arabidopsis, Rice and Soybean. Our results may provide a subset of potential candidates to be explored for stress-tolerance and genetic improvement of agronomic traits in Adzuki bean and Mung bean, as well as other related crops.

Identification and classification of VaWRKY and VrWRKY protein family based on the WRKY domain.
In this study, 91 proteins in Adzuki bean were identified possessing WRKY DNA-binding domain, whereas a total of 85 WRKY proteins were found in Mung bean (Supplementary Table I). Out of these, 78 VaWRKY and 68 VrWRKY proteins contain full-length WRKY domains of around 60 amino acids, while the remaining members possess partial WRKY domains. In case of VaWRKY proteins, six members have truncated N-terminal WRKY domain, while seven members have truncated C-terminal WRKY domain. Whereas, in case of VrWRKY proteins, six members lack intact N-terminal WRKY domain, whereas eleven members lack full length C-terminal WRKY domain (Table I). But as these proteins significantly possess typical conserved sequences, they were retained for their further analysis considering their putative functional roles. In case of Adzuki bean, seven proteins were isoforms, but no isoforms for WRKY proteins were observed in Mung bean. In total, 84 non-redundant VaWRKY genes and 85 VrWRKY genes were identified encoding the WRKY family of Adzuki bean and Mung bean, respectively. Although the presence partial proteins suggested the occurrence of pseudogenes encoding truncated WRKY protein. In total, 84 non-redundant putative VaWRKY genes and 85 putative VrWRKY genes were identified encoding the WRKY family of Adzuki bean and Mung bean, respectively. They were designated as VaWRKY1-84 and VrWRKY1-85, representing the WRKY members of Adzuki bean and Mung bean, respectively in Table I. The length of the proteins along with their theoretical pI (isoelectric point) and molecular weight were indicated in Supplementary Table I Generally, the WRKY proteins are classified in three groups and 5 sub-groups, depending on the number and nature of the WRKY domain, they possess. The members possessing two WRKY domains at N-terminal and C-terminal (Group-I NTWD and Group-I CTWD), were grouped into Group I, while the members having  only one WRKY domain belonged to Group II and Group III. In order to further classify the WRKY proteins, we  analyzed multiple sequence alignment of domain sequences, and constructed a phylogenetic of VaWRKY and  VrWRKY proteins, with WRKY sequences of a model legume crop, L. japonica, as reference, representing each group and subgroup. Moreover, a phylogenetic tree of WRKY domains, excluding the members with truncated domains was also constructed for further supporting the classification ( Supplementary Fig. S2). The phylogenies of the domains also indicate the consequence of the evolutionary selection and suggest their origins. As per our findings, the VaWRKY protein family have 17 members in Group I; five members in Group IIa, fourteen members in Group IIb, 24 members in Group IIc, seven members in Group IId and eleven members in Group IIe; and thirteen members in Group III ( Table 1). The distribution of genes encoding these 91 proteins, considering seven isoforms were as follows: fifteen members in Group I; five members in Group Ia, fourteen members in Group IIb, twenty members in Group IIc, seven members in Group IId and ten members in Group IIe; and thirteen members in Group III (Fig. 1A). While in Mung bean, where no isoforms were observed, we found sixteen genes in Group I; 1 gene in Group Ia, nineteen genes in Group IIb, twenty genes in Group IIc, seven genes in Group IId and nine genes in Group IIe; and thirteen genes encoding VrWRKY proteins in Group III. The distribution of WRKY genes in Mung bean was comparable to that of Adzuki bean (Fig. 1B). The distribution of Group I, Group IIc, Group IId, Group IIe and Group III was similar in Adzuki bean and Mung bean (17.9% and 18.8%, 23.8% and 23.5%, 8.3% and 8.2%, 11.9% and 10.6% and 15.5% and 15.3%, respectively). Although the numbers of Group IIa and Group IIb members differs between Adzuki bean and Mung bean, the combined composition of Group IIa and Group IIb was similar in both the crops, which is approximately 23%.
Gene exon-intron structure organization. The gene structural similarity and diversity plays salient role in the evolution of gene family. To study this, we generated exon-intron map of the VaWRKY and VrWRKY genes, using the software Gene Structure Display Server. The detailed representation of the coding region, introns and upstream or downstream regions of genes, was provided in Fig. 2. Introns, which are integral elements of eukaryotic genomes, actively participate in the genomic recombination leading to gene rearrangements and evolution. The group I genes have 2 to 8 introns. Most of the Group IIa genes and Group III possess 3 and 2 introns, respectively. Variable number of introns were found in Group IIb genes, 2 to 5 introns in VaWRKY genes, and 2 to 6 introns in VrWRKY genes. In Group IIc VaWRKY members, 0 to 3 introns were found with VaWRKY36 (Vang08g01570) possessing no intron, exceptionally. Whereas, 1 to 5 introns were found in case of Group IIc VrWRKY genes. In Group IId, the VaWRKY genes possess 2 to 5 number of introns, but the VrWRKY genes have 2 to 3 introns, showing less variation. Also, in Group IIe members, 1 to 5 and 2 to 4 introns observed in VaWRKY and VrWRKY, respectively. This variable distribution of introns in Group IIb, IIc and IIe members including Group IId members of Adzuki bean, indicate that both exon loss and gain has occurred during their evolution, which may explain why closely related WRKY genes falling in the same group can be diverse in function [31][32][33][34][35] .
Although, most of the WRKY genes were rich in introns, but the region encoding N-terminal WRKY domain possess only one intron, with the exception case of that region encoding N-terminal domain in Group I members lack introns ( Supplementary Fig. S2). Two types of introns, namely R (Arg)-type and V (Val)-type were found in most of the VaWRKY and VrWRKY genes (indicated in Fig. 2), discussed later in Section 3.4.
Chromosome location and gene duplication. Out of 84, only 83 of VaWRKY genes could be mapped on the chromosome. The precise location of one gene Vang04g16950 (VaWRKY71) could not be determined. As represented in Fig. 3A, most of the genes (about 60%) are located in chromosome 1 to 4, followed by chromosome 7, 9, 10 and 5. Only a few genes are located on chromosome 6, 8 and 11. Fifteen genes (two Group I, twelve Group II and one Group III) were mapped on chromosome 1; twelve genes (one Group I, nine Group II and two Group III) were mapped on chromosome 2; ten genes (three Group I, six Group II and one Group III) were mapped on Chromosome 3; twelve genes (four Group I, seven Group II and one Group III) were mapped on chromosome 4; five genes (two Group I, two Group II and one Group III) mapped on chromosome 5; one gene (Group II) was mapped on chromosome 6; eight genes (one Group I, three Group II and four Group III) were mapped on chromosome 7; four genes (one Group I and three Group II) were mapped on chromosome 8; seven genes (one Group I, four Group II and two Group III) were mapped on chromosome 9; seven genes (six Group II and one Group III) were mapped on chromosome 10; and two genes (Group II) were mapped on chromosome 11. The locus of 25 VrWRKY genes could be documented only on scaffolds, but not on chromosome, due to lack of information. As represented in Fig. 3B, most of the VrWRKY genes are located on chromosome 4 to 7. Four genes (Group II) were mapped on chromosome 1; one gene (Group II) was mapped on chromosome 1; four genes (two Group I and two Group II) were mapped on chromosome 3; six genes (one Group I, four Group II and one Group III) were mapped on chromosome 4; nine genes (two Group I, five Group II and two Group III) were mapped on chromosome 5; eleven genes (three Group I and eight Group II) were mapped on chromosome 6; eleven genes (two Group I, eight Group II and one Group III) were mapped on chromosome 7; three genes (one Group I and two Group II) were mapped on chromosome 8; five genes (three Group II and two Group III) were mapped on chromosome 9; three genes (one Group I and two Group II) were mapped on chromosome 10; and three genes (two Group II and one Group III) were mapped on chromosome 11. The detailed description of chromosome location and gene-length is mentioned in Supplementary Table S3.
It is well established that multiple members of WRKY gene family that form a large and extensively complex regulative network to control complicated physiological processes were expanded as a result of the long evolutionary history 40 . The individual genes undergo a series of genomic recombination and amplification during the process of evolution. The major player in this event are recent gene duplications, resulting in many paralogous pairs in different species 41,42 . Tandem duplication is one way of gene duplication, where a set of two or more genes located in the same chromosome within the range of 100-kb distance, separated by zero or few spacer genes 43  a block of genes is duplicated in a different chromosome. Furthermore, a gene cluster is defined as a chromosome region with two or more genes located within 200 kb sequence.
To unravel the mechanism of VaWRKY gene evolution, we explored gene duplication events. We found that, eleven VaWRKY genes formed five gene clusters. Chromosome 1, 3 and 7 contains 1 gene cluster. Whereas, two gene clusters were located on chromosome 2. Out of these five gene cluster, three were tandem duplications; the gene pairs Vang0228s00170 (VaWRKY78) and Vang0228s00230 (VaWRKY72) located on chromosome 2 with three spacer genes, Vang01g15570 (VaWRKY17) and Vang01g15560 (VaWRKY19) located on chromosome 3 with zero spacer gene, and Vang0459s00010 (VaWRKY82), Vang0459s00030 (VaWRKY79) and Vang0459s00020 (VaWRKY80) located on chromosome 7 with one and zero spacer genes (Fig. 3A). Further analysis of intron lengths and gene structure of the tandem duplicated gene revealed that, recombination might have taken place in addition to duplication event, which explains the variation in the encoded proteins. The gene pairs Vang05g09440 (VaWRKY18) and Vang05g09490 (VaWRKY20) located on chromosome 2 were also observed to be duplicated on chromosome 3 as gene pair Vang01g15570 (VaWRKY17) and Vang01g15560 (VaWRKY19). Similarly, 12 VrWRKY genes formed six gene clusters, with one gene cluster located on chromosome 1, 3, 4 and 5 each and two gene clusters located on scaffolds. Four of them were tandem duplications; Vradi04g07100 (VrWRKY65) and Vradi04g07130 (VrWRKY50) located on chromosome 4 with two spacer genes; Vradi05g05160 (VrWRKY78) and Vradi05g05170 (VrWRKY85) located on chromosome 5 with zero spacer genes; Vradi0214s00140 (VrWRKY80) and Vradi0214s00230 (VrWRKY79), and Vradi0338s00040 (VrWRKY83) and Vradi0338s00060 (VrWRKY84) located on scaffolds with zero spacer gene. In case of both Adzuki bean and Mung bean, the tandem and segmental gene duplication events are not that significant, suggesting that these phenomena do not play much significant role in the evolution of VaWRKY and VrWRKY genes.
Multiple sequence alignment and phylogenetic tree analysis.  Fig. 4. The R-type introns, which were phase-2 introns spliced exactly after the R position, before the zinc-finger motif, similar to the splicing position observed in Arabidopsis. However, the V-type introns were phase-0 introns, located just before the V position which is the sixth amino acid after the second cysteine residue in the C 2 H 2 zinc finger motif.
Based on the phylogenetic analysis of conserved WRKY domains, the VaWRKY and VrWRKY members were subdivided into 8 clades (Supplementary Fig. S2). In case of Group I members, the two WRKY domains, designated as Group I-NTWD and Group I-CTWD, clustered into two separate clades due to the divergence in their sequences. The clade of Group IIa was close to that of Group IIb, whereas the Group IId and IIe member were clustered, representing their origination from common gene ancestor or evolution under similar selection pressure. The Group III members are more similar to the Group I-NTWD members as compared to any other subgroup, suggesting that may have shared a common ancestor before their divergence into Group I and Group III, or the Group III members have been originated from Group I genes due to mutation in their zinc-finger domain after losing the C-terminal WRKY domain. Similarly, Group IIc members are more closely related to the Group I-CTWD members. This suggests that the Group IIc WRKY proteins might have been originated from the Group I proteins after the loss of the N-terminal domain. Moreover, two Group IIc members, Vang08g01570.1 (VaWRKY36) and Vang0051s00140.1 (VaWRKY35), and one member namely Vang0005s00450.1 (VaWRKY42) have been clustered with Group I-NTWD and Group I-CTWD members respectively, suggesting a common origin of their domains. Also in case of Mung bean, one Group IIc member namely Vradi04g07740.1 (VrWRKY42) clustered with Group I-CTWD member. Both of the two WRKY domains of Group I member Vradi05g10960.1 (VrWRKY3) clustered with the Group IIa member Vradi06g13520.1 (VrWRKY17), suggesting that VrWRKY3 has been originated due the duplication of VrWKY17 domain Supplementary Fig. S2. The phylogenetic relationships of the whole WRKY proteins sequence are illustrated in Fig. 5, where some of the Group IIc members have been clustered in Group I and Group III. This may be due to more sequence similarity in the region outside the domain of respective members.
Motif analysis. Apart from the conserved residue of 60 amino acids, other motifs also lie in the rest of the protein sequence, which may perform unknown functional or structural roles. To investigate further the similarity and diversity of motif, MEME analysis was performed to investigate conserved stretches of amino acids ranging from 6 to 50 amino acids 44 . The distribution of the 20 conserved motifs discovered significantly varied among different WRKY proteins (Fig. 6). In case of VaWRKY, five motifs, motif 1, motif 2, motif 3, motif 4 and motif 6, together comprise the WRKY domain, out of which motif 1 and motif 3 comprise the conserved heptapeptide (Supplementary Table S4). Whereas, six motifs, motif 1, motif 2, motif 3, motif 10, motif 15 and motif 16 makes the VrWRKY domains, with the WRKYGQK heptapeptide located in motif1 and motif 3. Out of the 15 remaining non-redundant motifs in VaWRKY and 14 non-redundant motifs in VrWRKY, the function of the majority of     (Fig. 7), suggesting the WRKY genes possessing similar domains express differently depending on the functional diversity. We also investigated the expression of closely located VaWRKY genes forming gene clusters. Their differential expression suggests their involvement in non-redundant signaling pathways.

Prediction of putative stress-responsive WRKY in Adzuki bean and Mung bean and their promoter analysis. The interaction between transcription factor and the stress-inducible cis-acting reg-
ulatory elements present in the promoter modulates the expression of gene regulatory networks response to the respective physical, environmental and biological stress. As these elements are highly conserved among orthologous or paralogous and co-regulatory genes, we investigated stress-responsive homologs of AtWRKY, OsWRKY and GmWRKY proteins reported to be involved in various abiotic and biotic stresses. Subsequently, the −1.5 kb promoter region of those putative stress-responsive candidates were analyzed using two different promoter analysis tool, PlantCare and PLACE, to identify stress-responsive cis-elements. Based on the phylogenetic data ( Supplementary Fig. S7), the VaWRKY and VrWRKY members which clustered closest to the reported stress-responsive AtWRKY, OsWRKY and GmWRKY proteins, used as references, were selected as putative homologous stress-responsive candidates ( Supplementary Fig. S7 and Table S8). A total of 17 VaWRKY and 18 VrWRKY proteins were chosen, which were also represented as phylogenetic tree along with their respective homolog members in Arabidopsis, Rice or Soybean, forming six different clades, depending on their respective groups and subgroups (Table 2 and Fig. 8). The stress-responsive elements recognized in the promoter of these genes are listed Supplementary Table S9. Exceptionally, Vradi07g15410 (VrWRKY18), a putative homolog of OsWRKY72 ( Supplementary Fig. S7), clustered in Clade III rather than Clade IV which accommodates rest of We found cis-elements involved in various stresses, for instance, 'DPBFCOREDCDC3' induced by abscisic acid, ' ABRELATERD1' involved in abscisic acid signaling and drought stress, 'ERELEE4' responsive to ethylene, 'PREATPRODH' involved in osmotic stress, 'DRECRTCOREAT' induced by high salt, cold and drought etc. Interestingly, in both the Adzuki bean and Mung bean, we found stress-inducible elements similar to their respective homologs in Arabidopsis and Rice (Supplementary Table S8). The location of these elements are depicted in Fig. 9. The detailed interpretation of these stress-responsive elements and their functions are mentioned in the Supplementary Tables S10 and S11.
The clade I members (as indicated in Fig. 8), predominantly consist of cis-elements responsive to osmotic stress, drought, and senescence. The clade II members possess elements mainly induced by ABA, drought, and pathogens. Few members also possess cold and senescence responsive elements. The clade III members possess elements generally induced by ABA, salt, salicylic acid including pathogens. The clade IV and V members predominantly possess elements responsive to drought, ABA and biotic stress, which were also found in stress-responsive OsWRKY members. Whereas clade VI containing the Group III members chiefly contained elements involved in the heat, cold, drought, osmotic stress, senescence and biotic stress, similar to as found in AtWRKY and OsWRKY.

Discussion
It is well established that classification of genes is not only essential but also a pre-requisite for the functional analysis of a gene family. In general, the gene families of transcription factors which bind DNA in a sequence-specific manner, contain highly conserved characteristic DNA binding domains or motifs, which are crucial for their biological functions. Furthermore, domain is considered as functional as well as an evolutionary unit of the protein, whose coding sequence can be duplicated and undergo recombination. It makes genome-wide analysis of DNA-binding domains of transcription factor pre-eminent. WRKY TFs are one of the largest families of  It can be elucidated from Fig. 1, that no correlation exists between the number of WRKY genes and the size of various crop genome. Although, the number of WRKY genes is proportional to the genome size in case of the three legume crops, viz. V. angularis, V. radiata and L. japonica. Furthermore, there are differences in the distribution of WRKY genes in respective groups or subgroup (except subgroup IIx), especially in Group III, among various species (Fig. 1). However, the percentage gene distribution among groups/subgroups, shows that the distribution of WRKY genes in closely related Vigna species, i.e. Adzuki bean and Mung bean is highly comparable (Fig. 1B). In most of the dicotyledonous legume crops, the Group I or Group IIc hold the most substantial number of members. In case of Adzuki bean and Mung bean, the Group IIc is the largest. However, in rice (monocot), the Group III possessing 36 members represents the largest group, indicating that evolution is more active in Group III and may be the members have more function in monocots. Such variation in the distribution of WRKY genes among dicots and monocots suggests that Group III members have been evolved independently after the dissection of monocots and dicots.
The average gene length of VaWRKY and VrWRKY was found to be 2.8 kb and 3 kb, respectively. Any intron conserved in the gene is considered ancient intron. Most of the WRKY domains contains two types of introns with conserved splicing positions, known as R-type intron and V-type intron. In our study, we found the V-type introns in Group IIa and Group IIb members and the R-type introns in rest of the group members ( Supplementary  Fig. S2). The above finding indicates that WRKY gene family encoding WRKY domain, was formed due to duplication of ancient genes carrying an intron, followed by divergence instead of formation of similar genes as a result  of convergence events. The N-terminal WRKY domain is surprisingly intron-less. Although some other group members also suffered a loss in introns, perhaps during evolution. Because of their substantial contribution to various physiological processes, it is likely that the WRKY family in angiosperms has expanded dramatically during evolution. Recent gene duplication events have been reported to be more prevalent in the expansion of WRKY genes in many crops like Arabidopsis, Rice 50 , Populus trichocarpa 54 , etc. However, in some cases, for instance, in Lotus Japonica 45 , recent duplications seemed to play no significant role in WRKY gene expansion. In our case also, we found few tandem and segmental gene duplications in Adzuki bean, but not that significant as compared to Arabidopsis and rice. The recent gene duplications still need to be study in Mung bean.
The motif analysis by MEME revealed interesting facts regarding the gene evolution. The conserved motif 4 of Group I proteins, occurring just before the WRKYGQK residue containing motif 1 of the C-terminal WRKY domain, were also found in Group IIc proteins (Fig. 6). It indicates that the Group IIc genes have been originated from the loss of the N-terminal WRKY domain of Group I genes, which is also evident by the phylogeny of WRKY domain regions (Supplementary Fig. S2). Moreover, the typical zinc-finger type Group IIc similar to that of Group I CTWD, and the clustering of the Group I CTWD and Group IIc domains in the same clade, further confirm the evolution pattern. The phylogenetic closeness and the conserved V-type intron of Group Ia and Group IIb indicate their evolution from a common origin. The Group IId and Group IIe also seemed to share a common ancestor. Similarly, the phylogenetic relationships between Group III domains and Group I-NTWD suggest their common origination (Section 3.4). Thus, our data support the theory of evolution of WRKY genes that the Group I is the oldest group, and Group II and Group II have been evolved from Group 59 .
The cis-acting regulatory elements present in the promoter regions are important molecular switches involved in the transcriptional regulation of a gene via controlling an extensive network of gene involved in various biological phenomenon including stress responses and developmental processes. Furthermore, it is evident that defined cis-elements can successfully contribute to the genome-wide screening of ABA and abiotic stress-responsive genes 60 . The identification of prominent cis-regulatory elements in the promoter region of Adzuki bean and Mung bean genes suggest the putative involvement of the respective genes in environmental stresses like drought, salinity, heat, osmotic stress, senescence, ABA signaling as well as pathogen resistance. Interestingly, these elements were also found in the WRKY genes of Arabidopsis and Rice, reported to be involved in various abiotic and biotic stress. In Arabidopsis and Rice, fusion genes containing a C-terminal WRKY motif and a NBS-LRR (nucleotide binding site-leucine-rich repeat) motif in the R gene were identified 50 . The R gene mainly confers resistance Figure 9. Cis-regulatory stress-responsive elements identified in the 1.5 kb upstream promoter region of (A) Stress-responsive VaWRKY candidates (B) Stress-responsive VrWRKY candidates. The elements commonly identified by both PlantCare and PLACE, and those involved in major stresses were chosen for pictorial representation. The sequence and position of these elements and some additional elements, not mentioned here, are described in Supplementary Tables S10 and S11. The first scale represents the location of cis-elements mentioned in the promoter analysis data. The second scale denotes the position of those elements from the 'start codon' as zero reference point. The candidates possessing isoforms are indicated by an asterisk (*).

Exon-intron organization and chromosome location.
To illustrate the exon-intron organization of the genes, Gene Structure Display Server 2.0 (http://gsds.cbi.pku.edu.cn/) was used 65 . The tool required the gene sequences and corresponding coding sequences as input. The phylogenetic tree of the genes belonging to different groups, created by MEGA 6.0 66 , was also uploaded with the gene and coding sequences, to generate the gene structure. Chromosome location of the VaWRKY genes was determined using the Vigna Genome Server (http://viggs.dna. affrc.go.jp/) 67 . The physical location of the genes was documented manually.
VaWRKY and VrWRKY protein classification, phylogenetic reconstruction and conserved motif analysis. To classify the VaWRKY and VrWRKY proteins in their respective groups, combination of two approaches, multiple sequence alignment and phylogenetic analysis were employed. First, the domain sequences were aligned by Clustal Omega (http://www.ebi.ac.uk/Tools/ msa/clustalo/) using default settings 68 , and the conserved amino acid residues were displayed with GeneDoc software 69 . The domains were screened manually to classify them in respective groups and subgroups, based on their typical zinc-finger-motif and sequence similarity. Subsequently, the phylogenetic tree of VaWRKY and VrWRKY proteins, along with LjWRKY used as reference, were created using Clustal Omega to confirm the classification 68 . To further support the protein classification, phylogenetic analysis of the WRKY domain sequences was also performed. After classification the WRKY domains were analyzed with MEME suite 4.11.1 (http://meme-suite.org/tools/meme), for the identification of conserved motifs, with optimum search parameters as: minimum motif width = 6; maximum motif width = 50; maximum number of motif = 20; minimum sites per motif = 2; maximum sites per motif-600) 44 .
RNA sequence data analysis. The raw RNA sequence data of Adzuki bean and Mung bean (SRR1652394, SRR3406553, SRR1407784 and SRR1653637) were downloaded from NCBI Sequence Read Archive and their SRA files were saved in FASTQ format [46][47][48][49] . To acquire quality reads, we performed quality trimming and adapter removal of raw sequencing reads by Trimmomatic 0.32 with the following options: ILLUMINACLIP:adapters. Prediction of putative stress-responsive VaWRKY and VrWRKY and their promoter analysis. Using the reference of six AtWRKY, five OsWRKY and two GmWRKY proteins reported to be involved in various abiotic and biotic stresses [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37] , phylogenetic trees of VaWRKY and VrWRKY protein sequences were created with Clustal Omega tool. The proteins closely clustered to the reference WRKYs in the tree, were selected as stress-responsive homologs in Adzuki bean and Mung bean. The 1.5 kb upstream promoter region sequence from the transcription start site of the homologous VaWRKY and VrWRKY genes were retrieved using the 'Jbrowse' tool of Crop Genomics Lab database (http://plantgenomics.snu.ac.kr/) in order to perform promoter analysis. To investigate the stress-responsive cis-elements present in the 1.5 kb upstream region of the promoters, two different tools PlantCare (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) and PLACE (https://sogo.dna.affrc.go.jp/cgi-bin/sogo.cgi) were used 73,74 . The cis-elements identified were then mapped at their respective positions manually.