Genome-Wide Identification and Expression Analysis of the WRKY Gene Family in Cassava

The WRKY family, a large family of transcription factors (TFs) found in higher plants, plays central roles in many aspects of physiological processes and adaption to environment. However, little information is available regarding the WRKY family in cassava (Manihot esculenta). In the present study, 85 WRKY genes were identified from the cassava genome and classified into three groups according to conserved WRKY domains and zinc-finger structure. Conserved motif analysis showed that all of the identified MeWRKYs had the conserved WRKY domain. Gene structure analysis suggested that the number of introns in MeWRKY genes varied from 1 to 5, with the majority of MeWRKY genes containing three exons. Expression profiles of MeWRKY genes in different tissues and in response to drought stress were analyzed using the RNA-seq technique. The results showed that 72 MeWRKY genes had differential expression in their transcript abundance and 78 MeWRKY genes were differentially expressed in response to drought stresses in different accessions, indicating their contribution to plant developmental processes and drought stress resistance in cassava. Finally, the expression of 9 WRKY genes was analyzed by qRT-PCR under osmotic, salt, ABA, H2O2, and cold treatments, indicating that MeWRKYs may be involved in different signaling pathways. Taken together, this systematic analysis identifies some tissue-specific and abiotic stress-responsive candidate MeWRKY genes for further functional assays in planta, and provides a solid foundation for understanding of abiotic stress responses and signal transduction mediated by WRKYs in cassava.


INTRODUCTION
The WRKY family is a large family of transcription factors (TFs) found in higher plants (Rushton et al., 2010). SPF1, the first reported WRKY transcription factors, plays crucial roles in the regulation of gene expression (Ishiguro and Nakamura, 1994). WRKY TFs contain one or two WRKY domains which have a highly conserved WRKYGQK motif at the N-terminus and a zincfinger structure at the C-terminus (Llorca et al., 2014). Based on the variation in WRKY domain and the pattern of the zinc-finger motif, WRKY proteins can be divided into three major groups (1, 2, and 3) with several subgroups (Eulgem et al., 2000). The group 1 typically contains two WRKY domains including a C 2 H 2 motif, while group 2 and group 3 are characterized by a single WRKY domain. Group 2 also contains a C 2 H 2 zinc-finger motif and can be further divided into five subgroups (2a-2e) based on the phylogeny of the WRKY domains, whereas group three contains a zinc-finger-like motif ending with C 2 -H-C (Eulgem et al., 2000).
In recent years, accumulated evidence has confirmed that a large number of WRKY genes are induced by abiotic stresses and play important roles in the regulation of plant tolerance to abiotic stress (Seki et al., 2002;Rushton et al., 2010;Li et al., 2011;Scarpeci et al., 2013). In Arabidopsis, AtWRKY30 was induced by abiotic stress including treatments with methyl viologen (MV), H 2 O 2 , arsenic, drought, NaCl, and mannitol, and overexpression of AtWRKY30 increased plants tolerance to MV and salinity stresses (Scarpeci et al., 2013). WRKY46, another WRKY gene from Arabidopsis, was significantly induced by drought, salt, and H 2 O 2 , and wrky46 mutant was less tolerant to osmotic and salt stress than WT (Ding et al., 2014). WRKY25 and WRKY26 were induced under heat stress and were confirmed to play positive roles thermotolerance in Arabidopsis . Additionally, overexpression of WRKY25 or WRKY33 increased plant tolerance to salt stress and sensitivity to ABA (Jiang and Deyholos, 2009). Likewise, 41 out of 103 rice WRKY genes showed significant differences in their transcript abundance under abiotic stress (cold, drought and salinity; Ramamoorthy et al., 2008). Some rice WRKYs have been shown to be positive regulators of abiotic stresses, such as OsWRKY5 (Berri et al., 2009), OsWRKY7 (Ramamoorthy et al., 2008), OsWRKY11 (Wu et al., 2009), OsWRKY43 (Berri et al., 2009), OsWRKY45 (Qiu and Yu, 2009), and OsWRKY47 (Raineri et al., 2015). For example, overexpression of OsWRKY45 in Arabidopsis was found to increase plant tolerance to salt and drought, and to decrease sensitivity to ABA (Qiu and Yu, 2009). Overexpression of OsWRKY47 increased plant tolerance to drought and yield compared to WT (Raineri et al., 2015). This evidence demonstrated that the WRKY gene family may contain important regulatory factors involved in plant response to abiotic stress.
To date, genome-wide analysis has identified a large number of WRKY family members in several species with 74 WRKY genes in Arabidopsis (Arabidopsis thaliana; Ulker and Somssich, 2004), 103 in rice (Oryza sativa cv. Nipponbare; Ramamoorthy et al., 2008), 45 in barley (Hordeum vulgare; Mangelsen et al., 2008), 55 in cucumber (Cucumis sativus; Ling et al., 2011), 119 in maize (Zea mays; Wei et al., 2012), 182 in soybean (Glycine max; Bencke-Malato et al., 2014), and 109 in cotton (Gossypium aridum; Fan et al., 2015). However, there is currently no evidence regarding the WRKY family in the important tropical plant cassava. Cassava (Manihot esculenta Crantz) is the third most important crop after rice and maize in Africa, Asia, and Latin America, where it is an important food security crop (Oliveira et al., 2014). Cassava, a major staple crop, has the starchy roots that provide dietary carbohydrate for 800 million people across the tropical and sub-tropical world (International Cassava Genetic Map Consortium, 2014). Due to its high starch production and limited input, cassava is also a major producer of industrial starch and bioethanol (Zidenga et al., 2012;Perera et al., 2014). Cassava is particularly tolerant to drought and lowfertility soils when facing environmental stresses (International Cassava Genetic Map Consortium, 2014;Zeng et al., 2014). However, the mechanisms by which cassava responds to abiotic stress are poorly understood. Thus, understanding of the molecular mechanisms underlying the tolerance of cassava to abiotic stress may provide effective methods for genetic improvement of stress tolerance for cassava and other crops. The high-quality sequencing of cassava wild ancestor and cultivated varieties reported in our previous study have provided an excellent opportunity for genome-wide analysis of cassava genes (Wang et al., 2014a). Based on the significance of WRKYs involved in plant growth, development and adaption to the environment and on the lack of any genome-wide systematic analysis of cassava WRKY genes, the WRKY family was selected for a systematic analysis in cassava. In this study, 85 WRKY genes from the cassava genome were identified and detailed studies of their phylogeny, conserved motifs, gene structure, expression profiles in various tissues, and in response to drought, osmotic, salinity, cold, oxidative stresses, and signaling of ABA were performed. The current results may provide a novel insight into the future work on the function of WRKYs and abiotic stress responses in cassava.

Plant Materials and Treatments
W14 (M. esculenta ssp.flabellifolia) is an ancestor of the wild cassava subspecies with a strong tolerance to drought stress (Wang et al., 2014a). South China 124 (SC124) is a widely planted cassava cultivar in China (Zeng et al., 2014). Argentina 7 (Arg7) adapts to a geographical high-latitude region of Argentina (Zhao et al., 2014). All plants were grown in a glass house of the Chinese Academy of Tropical Agricultural Sciences (Haikou, China). Stem segments with three nodes were cut from 8 months old cassava plants and inclined into pots with a mixture of soil and vermiculite (1:1) where they were regularly watered (Hu et al., 2015). The plants were grown from April to July 2013 during which time the temperature in the glass house ranged from 20 to 35 • C. The transcripts of cassava WRKY genes in different tissues, including stems (90 days after planting), leaves (90 days after planting), and middle storage roots (150 days after planting) were examined with wild subspecies (W14) and cultivated variety (Arg7) under normal growth conditions. Ninety-days-old leaves and roots were sampled from Arg7, SC124 and W14 to study the transcriptional response of cassava WRKY genes under 12 days drought stress. After 60 days of normal cultivation, the Arg7 seedlings similar in growth vigor were used in the following treatments. For abiotic stress and signal molecule treatments, Arg7 seedlings were subjected to 200 mM mannitol for 14 days, 300 mM NaCl for 14 days, 100 µM abscisic acid (ABA) for 24 h, 3.27 M (10%) H 2 O 2 for 24 h and low temperature (4 • C) for 48 h, respectively. According to Scarpeci et al. (2013) and Ding et al. (2014), 20 mM H 2 O 2 can induce oxidative stress in Arabidopsis. In this study, high concentration of H 2 O 2 (3.27 M) was used to strongly induce oxidative stress due to the woody feature of cassava.

Identification and Phylogenetic Analyses of the WRKY Gene Family in Cassava
The whole protein sequence of cassava was obtained from the cassava genome database (http://www.phytozome.net/cassava. php). Sequences of the AtWRKY and OsWRKY genes were downloaded from UniPort (http://www.uniprot.org/) and RGAP databases (http://rice.plantbiology.msu.edu/), respectively. To identify the cassava WRKY family genes, two different approaches were used as follows: firstly, the local Hidden Markov Model-based searches (HMMER: http://www.ebi.ac.uk/Tools/ hmmer/) built from known WRKYs to search the cassava genome database (Finn et al., 2011); secondly, BLAST analyses with all the Arabidopsis and rice WRKYs as queries were employed to check the predicted WRKYs in cassava database. With the help of CDD (http://www.ncbi.nlm.nih.gov/cdd/) and PFAM databases (http://pfam.sanger.ac.uk/), all the potential cassava WRKY genes identified from HMM and BLAST searchs were only accepted if they contained the WRKY domain, then using multiple sequence alignments to confirm the conserved domains of predicted WRKY sequences. Additionally, Clustal X 2.0 and MEGA 5.0 were used to constructed a bootstrap neighbor-joining (NJ) phylogenetic tree based on amino acid sequence of WRKY domains of cassava WRKY members and selected Arabidopsis WRKYs with 1000 bootstrap replicates (Larkin et al., 2007;Tamura et al., 2011). Furthermore, to better exhibit the characteristic of MeWRKY gene structure and conserved motifs, a NJ phylogenetic tree was created based on the full amino acids of cassava WRKYs.

Protein Properties and Sequence Analyses
The online ExPASy proteomics server (http://expasy.org/) was used to investigate the molecular weight (MW) and isoelectric points (pI) of presumed WRKY proteins. The conserved motifs in full-length WRKY proteins were identified using the MEME program (http://meme.nbcr.net/meme/cgi-bin/meme. cgi). Parameters employed in the analysis were: maximum number of motifs was 10 and the optimum width of motifs was set from 15 to 50 (Tao et al., 2014). Furthermore, all identified motifs were annotated according to InterProScan (http://www.ebi.ac.uk/Tools/pfa/iprscan/). The gene structures were identified by gene structure display server program (GSDS, http://gsds.cbi.pku.edu.cn/). Exon/intron organization was further checked by alignment of coding sequence and genomic DNA sequence of each WRKY gene.

Transcriptomics Analysis
Total RNA was extracted from stems, leaves and storage roots in Arg7 and W14 under normal growth conditions, and was also extracted from leaves and roots of Arg7, SC124 and W14 under normal conditions and 12 days drought treatment. Total RNA was isolated using plant RNeasy extraction kit (TIANGEN, China) following manufacturer's instructions and the concentration and purity were evaluated by NanoDrop 2000c (Thermo Scientific, USA). Three µg total RNA of each sample were used to construct the RNA pools according to the Illumina instructions, and subsequently sequenced by Illumina GAII following Illumina RNA-seq protocol. A total of 610.70 million 51-bp raw reads was generated from the 18 samples. Adapter sequences were removed from raw sequence reads using FASTXtoolkit (version 0.0.13, http://hannonlab.cshl.edu/fastx_toolkit/). Sequence quality was examined using FastQC (http://www. bioinformatics.babraham.ac.uk/projects/fastqc/) and low quality sequences (including reads with unknown base pairs "N") were removed, which produced 583.82 million clean reads. Clean reads were mapped to cassava reference genome (version 4.1) derived from the phytozome website (ftp://ftp.jgi-psf.org/pub/ compgen/phytozome/v9.0/Mesculenta/) using Tophat v.2.0.10 (http://tophat.cbcb.umd.edu/) (Trapnell et al., 2009), and 88.7% reads were aligned. The resulting alignment files are provided as input for Cufflinks to generate transcriptome assemblies (Trapnell et al., 2012). Gene expression levels were calculated as FPKM according to the length of the gene and reads count mapped to this gene: FPKM = total exon fragments/[mapped reads (millions) × exon length (kb)]. DEGseq was applied to identify differentially expressed genes with a random sampling model based on the read count for each gene (Wang et al., 2010).

Quantitative RT-PCR Analysis
Expression of MeWRKY genes in response to various abiotic stress (osmotic, salt, cold, and oxidative stress) and ABA signaling were examined by qRT-PCR analysis with Stratagene Mx3000P Real-Time PCR system (Stratagene, CA, USA) using SYBR R Premix Ex Taq ™ (TaKaRa, Japan) according to the manufacturer's instructions. Total RNA was extracted from leaves of control and treated samples. Two hundred ng Poly(A) + mRNA from each treatment was converted into cDNA using AMV Reverse Transcriptase (Promega, Madison, WI, USA) at 42 • C in a 20 µL reaction volume that subsequently served as the template for qRT-PCR. The amplification conditions used for all PCRs were implemented as follows: 10 min at 95 • C, and followed by 40 cycles of 10 s at 95 • C, 15 s at 50 • C, and 30 s at 72 • C. The relative expression of the target genes was determined using the 2 − Ct method (Livak and Schmittgen, 2001). The specific primers were designed according to the WRKY gene sequences by Primer 5.0 software (Table S1). Subsequently, reaction specificities for each primer pair was tested using qRT-PCR melting curve analysis, agarose gel electrophoresis, and sequencing PCR products. Amplification efficiencies of gene-specific primers ranged from 90 to 110%. β-tubulin gene (TUB) and elongation factors 1α gene (EF1) verified to be constitutive expression and suitable as internal controls were used as internal references for all the qRT-PCR analyses (Salcedo et al., 2014). Each treated sample contained a corresponding regularly-watered control and each sample was performed with three independent biological replications. Then, the treated and control plants at each time point were sampled to perform expression analysis. The relative expression levels of MeWRKY genes in each treated time point were compared with corresponding regularly-watered control (Wang et al., 2014b). Statistical difference were performed by Duncan's multiple range test (n = 3). Means denoted by the same letter do not significantly differ at P < 0.05.

Identification and Phylogenetic Analysis of Cassava WRKYs
To identify the WRKY family members in cassava, both BLAST and HMMER searches were performed to search the cassava genome with Arabidopsis and rice WRKY sequences as queries. After these searches, a total of 85 putative members of the WRKY family were detected in the complete cassava genome. Conserved domain analysis further confirmed that all the WRKYs contain single or double WRKY domains at the N-terminus, which are the basic characteristics of WRKY family. The 85 predicted WRKY proteins ranged from 149 (MeWRKY22) to 737 (MeWRKY64) amino acids (aa) in length with an average of 369.4 aa, the relative molecular mass varied from 17.19 kDa (MeWRKY22) to 79.76 kDa (MeWRKY64), and the pIs ranged from 4.91 (MeWRKY59) to 9.89 (MeWRKY1) with 38 numbers pI >7 and others pI <7 ( Table S2). cDNAs of all 85 MeWRKY genes have been submitted to GenBank and their accession numbers in GenBank are shown in Table S3.
To study the evolutionary relationships between cassava WRKY proteins and known WRKYs from Arabidopsis, an unrooted neighbor-joining phylogenetic tree was created based on multiple alignments of the predicted amino acid sequences of the WRKY domains from cassava and Arabidopsis. As shown in Figure 1, 85 MeWRKY proteins were classified into three major groups, among which group 2 was subdivided into five subgroups together with WRKYs from Arabidopsis. This was in accordance with the classification of WRKY family in Arabidopsis (Eulgem et al., 2000), cucumber (Ling et al., 2011), maize (Wei et al., 2012), and soybean (Bencke-Malato et al., 2014). Groups 1, 2, and 3 contained 17, 56, and 12 MeWRKY proteins, respectively. A total of 5, 14, 20, 8, and 9 proteins were assigned to subgroups 2a, 2b, 2c, 2d, and 2e, respectively. Generally, group 1 contained two WRKY domains, but there were a few MeWRKY proteins that contained only one WRKY domain, such as, MeWRKY9, -13, -35, and -82. The same phenomenon was also found in Arabidopsis (Eulgem et al., 2000) and maize (Wei et al., 2012). According to one previous report (Wei et al., 2012), the loss of WRKY domain seems to be more common in monocotyledons than in dicotyledons. It can be deduced that group 1 might contain the original genes of other groups and that MeWRKY9, -13, -35, and -82 emerged later during evolution.
Phylogenetic analysis also showed that there were some closely related orthologous WRKYs between cassava and Arabidopsis (MeWRKY42 and AtWRKY51; MeWRKY47 and AtWRKY48; MeWRKY48 and AtWRKY9; MeWRKY69 and AtWRKY27; MeWRKY35 and AtWRKY49; MeWRKY70/MeWRKY72 and AtWRKY12; MeWRKY37/MeWRKY57 and AtWRKY35; MeWRKY1/MeWRKY44 and AtWRKY7; MeWRKY45 and AtWRKY21; MeWRKY26 and AtWRKY3/AtWRKY4), suggesting that an ancestral set of WRKY genes existed prior to the divergence of cassava and Arabidopsis and that WRKYs from cassava generally have close relationship with the proteins from Arabidopsis. MeWRKY1 and MeWRKY44 showed a high degree of similarity with AtWRKY7, which was reported to negatively regulate plant defense against bacterial pathogens (Kim et al., 2006). MeWRKY69 shared considerable similarity with AtWRKY27 that is also involved in the regulation of plant defense against the bacterial pathogens by regulating the expression of nitrogen metabolism and nitric oxide (NO) generation genes (Mukhtar et al., 2008). AtWRKY51, which showed a high degree of similarity with MeWRKY42, was reported to mediate jasmonic acid (JA) signaling and partially alter resistance to virulent pathogens (Gao et al., 2011). These results suggested the possible functions of WRKY genes in cassava.

Conserved Motifs and Gene Structure of Cassava WRKYs
To further detect the structural features of cassava WRKYs, conserved motifs and intron/exon distribution were analyzed according to their phylogenetic relationships. A total of 10 conserved motifs in cassava WRKYs were found using MEME software and further annotated by InterPro Scan 5 (Figure 2; Figure S1). Results showed that three (1-3) of 10 motifs were annotated as WRKY DNA-binding, which is a basic characteristics of the WRKY family. All MeWRKYs contained at least one of them, indicating that the cassava WRKYs identified in this study had conserved features of the WRKY family. Notably, all the MeWRKYs contain at least two motifs, except for three members (MeWRKY8, -28, and -45) only containing motif 2 in cluster E. In cluster A, all the MeWRKYs, except for MeWRKY5, -22, -24, and -42, contained motifs 1, 2, 6, and 9. Interestingly, most of the MeWRKY members in cluster A specially showed motifs 8 and 9 in comparison to MeWRKYs in other clusters. In cluster B, all the MeWRKYs, except for MeWRKY82, contained motifs 1 and 7, and motif 10 was uniquely dispersed in four members (MeWRKY4, -9, -13, and -73). In cluster C, all members contained motifs 1, 2, and 4, except for MeWRKY35 which did not contain the motif 4. In cluster D, all members contained motifs 1 and 2, and five members (MeWRKY3, -39, -51, -80, and -15) also contained motif 4 in addition to motifs 1 and 2. In cluster E, all members contained motifs 1 and 2, except for the closely related MeWRKY8, -28, and -45, which only contained motif 1. Generally, WRKY members in the same cluster commonly shared similar motif compositions, indicating functional similarity among them.
FIGURE 1 | Phylogenetic analysis of WRKY proteins from cassava and Arabidopsis. The NJ tree was constructed with WRKY domains of WRKYs from cassava and Arabidopsis using ClustalX 2.0 and MEGA5 with 1000 bootstrap. Branches with less than 50% bootstrap support were collapsed. The WRKY proteins are grouped into three groups (1-3) and five subgroups (2a-2e). Group 1 proteins with the suffix "N" or "C" indicates the N-terminal WRKY domains or the C-terminal WRKY domains. "AtWRKYs" are the WRKY proteins from Arabidopsis. "MeWRYKs" indicate the WRKY proteins from cassava.
Exon-intron structural diversity, an important part in the evolution of gene families, provides additional evidence supporting phylogenetic groupings (Shiu and Bleecker, 2003;Wang et al., 2014c). Intron/exon distribution was analyzed to better understand phylogenetic relationship and classification of cassava WRKYs. As shown in Figure 3, the number of introns in MeWRKY genes varied from 1 (MeWRKY9, to 5 (-17, -21, -29, -32, -33, -34, and -63). However, in rice and rubber tree, the number of introns varied from 0 (OsWRKY10 and OsWRKY44) to 8 (OsWRKY41.D1 and OsWRKY41.D2) and 1 (HbWRKY22, -34, -35 and -36) to 7 (HbWRKY15), respectively (Xie et al., 2005;. These results indicated that WRKYs in cassava have less gene structure diversity than that in rice and rubber tree. Additionally, 42 out of 85 MeWRKY genes each had two introns. The same phenomenon was also observed in rice and rubber tree with 42 of 92 and 40 of 81 WRKY genes containing two introns each, respectively (Xie et al., 2005;. Cluster A  contained 2-5 introns; cluster B contained 1-5 introns; cluster C contained 1-3 introns; all cluster D MeWRKYs contained 2 introns, except for MeWRKY46 with only one intron; and cluster E MeWRKYs contained two introns. According to a previous report (Nuruzzaman et al., 2010), the rate of intron loss is faster than the rate of intron gain after segmental duplication in rice. Consequently, it can be concluded that clusters A and B might contain the original genes, from which those in other clusters were derived. Generally, MeWRKYs in the same cluster of the phylogenetic tree show similar exon-intron structures.

Expression Profiles of MeWRKY Genes in Different Tissues
To provide some clues on the roles of MeWRKY genes in cassava growth and development, the expression profiles of MeWRKY genes from different organs, including stems, leaves and storage roots were tested in a wild subspecies (W14) and cultivated variety (Arg7) using transcriptomic data. W14, a wild cassava subspecies, has a low rate of photosynthesis, tuber root yield, and starch content in root tubers, but strong tolerance to drought stress (Wang et al., 2014a). Arg7, a cultivated variety, can tolerate moderate drought stress (Zhao et al., 2014). Expression analysis of MeWRKY genes in these two accessions will provide insight into cassava development between wild subspecies and cultivated variety. Seventy-two of 85 MeWRKY genes were captured from the corresponding transcriptomic data (Figure 4; About 20.8% (15/72) of MeWRKY genes with high expression levels (value >5) in all three tested organs in Arg7 and W14, suggesting that MeWRKY genes may be involved in organ development. Transcriptomic data also showed that 56 MeWRKY genes had a constitutive expression pattern that expressed in all the tissues of the two accessions, suggesting that these genes might play a role in plant growth, development, and cellular homeostasis. The remaining 16 MeWRKY genes exhibited differential expression patterns, with specific to some particular tissues, such as MeWRKY16, MeWRKY20, and MeWRKY23. This phenomenon was also observed in rice (Ramamoorthy et al., 2008), cucumber (Ling et al., 2011), rubber tree  and grape (Wang et al., 2014c), indicating that the functions of the WRKYs are diverse in both monocotyledon and dicotyledon.
There were 33 MeWRKY genes that showed higher expression levels in leaf and stem tissues in Arg7 than that in W14. However, 25 MeWRKY genes had higher expression levels in storage roots in W14 than that in Arg7. Interestingly, MeWRKY8, -18, -34, -45, -54, -80, and -83 showed higher expression levels in Arg7 than in W14 in leaf and stem tissues, but opposite result was observed in storage roots. These MeWRKY genes have strong expression levels for special tissues in different accessions, indicating their key roles in tissue development or tissue functions.

Expression of MeWRKY Genes in Response to Drought in Different Accessions
Accumulated evidence has suggested that WRKY family genes play a significant role in plants' response to drought or osmotic stress (Ramamoorthy et al., 2008;Ren et al., 2010;Rushton et al., 2010;Ling et al., 2011;Tripathi et al., 2014). Thus, there is need to examine the expression patterns of WRKY genes in response to drought stress, which may provide important clues for further understanding the mechanisms of cassava involved in strong tolerance. For this reason, 3-month-old cassava seedlings (a wild subspecies W14 and two cultivated varieties Arg7 and SC124) were deprived of water for 12 days, and then the leaf and root tissues were collected to extract RNA for subsequent RNA-seq analysis. Heatmap representation of expression profiles of 78 MeWRKY genes under drought stress conditions were captured from the corresponding transcriptomic data (Figure 5; Table S5).
The transcriptomic data given above showed that there were significantly more WRKY genes upregulated by drought at the transcription level in roots than in leaves in W14, but there were fewer in roots than in leaves in Arg7 and SC124. There were also more WRKY genes significantly induced by drought (value >1) in roots than in leaves in W14 but fewer in roots than in leaves in SC124. W14 showed stronger tolerance to drought stress than SC124 and Arg7, two varieties commonly cultivated in China and Southeast Asia, respectively (Wang et al., 2014a). Cassava can form deep root systems (soil depth below 2 m), which is beneficial for penetrating into deeper soil layers and absorbing water stored in the soil (Okogbenin et al., 2013). Moreover, numerous studies have confirmed that the WRKY family genes play a positive role in the drought stress response in various species (Qiu and Yu, 2009;Ren et al., 2010;Jiang et al., 2012;Ding et al., 2014;Raineri et al., 2015). Therefore, these findings indicate that cassava WRKY genes might play an important role in water uptake from soil by roots, and hence maintaining strong tolerance to drought stress in W14 subspecies.

Temporal Expression Profiles of MeWRKY
Genes upon Exposure to Various Stress and Related Signaling WRKY genes have been reported to play pivotal role in the regulation of plant tolerance to various stress and related signaling transduction in various species (Rushton et al., 2010(Rushton et al., , 2012Tripathi et al., 2014;Banerjee and Roychoudhury, 2015). Hence, to investigate the roles of MeWRKY genes in response to various environmental stresses and related signaling, the expression profiles of MeWRKY genes under these treatments were analyzed. Nine MeWRKY genes (MeWRKY8, distributed in different subgroups and up-regulated by drought stress as indicated by RNA-seq data in different cassava accessions were selected for further examination of their transcriptional response to osmotic, salt, cold, ABA, and H 2 O 2 treatments. Under NaCl treatment, MeWRKY18 was induced after 2-6 h and 14 days treatment with significant up-regulation at 6 h. MeWRKY65 and MeWRKY83 were significantly induced at 6 h. MeWRKY28 was significantly induced at 14 days, while MeWRKY85 was visibly down-regulated at 14 days. MeWRKY37 showed down-regulation at all the treated timepoints. Other three WRKY genes (MeWRKY49, -62, and -76) did not display obvious trends during salt treatment (Figure 6). In Arabidopsis, some WRKY genes, including AtWRKY8  (Ramamoorthy et al., 2008;Yu et al., 2010;Tao et al., 2011). Accumulating evidence has suggested that some WRKY genes play a positive role of in response to salt stress, such as AtWRKY8 (Hu et al., 2013), AtWRKY25 (Jiang and Deyholos, 2009), AtWRKY30 (Scarpeci et al., 2013), and AtWRKY33 (Jiang and Deyholos, 2009). However, other WRKYs, including AtWRKY18 , OsWRKY45-2 (Tao et al., 2011), AtWRKY46 (Ding et al., 2014), and AtWRKY60  were found to act as negative regulators in salt stress response in Arabidopsis and rice. These studies indicated that MeWRKY genes may be involved in the salt stress response.
As shown in Figure 7, under mannitol treatment, MeWRKY18, -28, -37, -49, and -76 were induced during 2 h 14days treatment and showed significant induction at 6 h, 6 h, 14 days, 6 h, and 2 h, respectively. MeWRKY62 and MeWRKY83 expression were induced during 2-6 h and 14 days treatment with significant up-regulation at 6 h. MeWRKY65 showed induction at 2 h treatment. MeWRKY85 did not show obvious trends during mannitol treatment. Notably, MeWRKY18 showed up-regulation at all treated points and reached the highest expression level (value >40) at 6 h, indicating its possible roles in osmotic/drought stress responses. In Arabidopsis, some WRKY genes, including AtWRKY57  and AtWRKY63/ABO3 (Ren et al., 2010), have been reported to positively regulate drought stress tolerance. However, some WRKY genes, including AtWRKY18 , AtWRKY46 (Ding et al., 2014), AtWRKY53 (Sun and Yu, 2015), AtWRKY54 , AtWRKY60 , and AtWRKY70 , which showed significant induction during drought stress, have been reported to negatively regulate drought stress tolerance. MeWRKY18, showing high similarity with AtWRKY54, may represent a functional gene involved in drought tolerance in cassava. In rice, 23 WRKY genes have been reported to be induced under drought treatment (Ramamoorthy et al., 2008;Qiu and Yu, 2009;Wu et al., 2009;Shen et al., 2012;Raineri et al., 2015), among which OsWRKY11 (Wu et al., 2009), OsWRKY30 (Shen et al., 2012, OsWRKY45 Yu, 2009), andOsWRKY47 (Raineri et al., 2015) have been confirmed to function as positive factors in the regulation of plant tolerance to drought/osmotic stress. In cucumber, the expression of 4 WRKY genes (CsWRKY2, was found to be upregulated after drought treatment (Ling et al., 2011). In cotton (Gossypium hirsutum) roots, 15 out of 26 GhWRKs (GhWRKY9, and 7 out of 26 GhWRKs (GhWRKY12, were up-and down-regulated, respectively, under dehydration conditions . Together, these results indicate the important roles of these WRKY genes in response to osmotic/drought stress. Cold stress, a common environmental stress, affects plants growth and crop productivity, especially in tropical and subtropical origin (Wang et al., 2012). However, little is known about the mechanisms underlying the action of WRKYs in cold stress response. In Arabidopsis, WRKY34 was reported to be significantly induced by cold treatment and act as a negative regulator to cold response (Zou et al., 2010). In rice, 2 and 15 WRKY genes were up-and downregulated by cold treatment, respectively (Ramamoorthy et al., 2008;Yokotani et al., 2013). Among them, overexpression of OsWRKY76 increased tolerance to cold stress (Yokotani et al., 2013). Under cold treatment, MeWRKY18, -49, -65, and -83 showed up-regulation at all the treated timepoints, with significant up-regulation at 15, 48, 2, and 48 h, respectively. MeWRKY37, -62, and -85 showed significant up-regulation at 5, 48, and 48 h, respectively. However, MeWRKY28 and -76 expression was repressed during all the treated time points (Figure 8). The expression levels of MeWRKY49 and -83 increased as treatment time continued, suggesting their possible function in cold response. They could be used in further functional characterization. Cassava, an important tropical crop, is distributed in tropical areas all over the world. Cold stress significantly restricts plant growth, agricultural productivity, and the development of cassava. Research on WRKY-mediated cold response in cassava may benefit further functional characterization of WRKY genes and investigations of the mechanisms underlying the cold response in cassava.
H 2 O 2 , a well-known toxic molecule, plays a key role in several biotic and abiotic signaling pathways and its accumulation has been found to be induced by environmental and developmental stimuli (Costa et al., 2010). In Arabidopsis, several WRKYs, including AtWRKY6, -22, -28, -30, -46 (Scarpeci et al., 2008), and AtWRKY25 (Jiang and Deyholos, 2009), are rapidly and highly induced after oxidative stress treatment. Among them, AtWRKY28 (Babitha et al., 2013) and AtWRKY30 (Scarpeci et al., 2013) were found to positively regulate oxidative stress tolerance, whereas AtWRKY25 (Jiang and Deyholos, 2009) acts as a negative regulator of oxidative stress response. In other species, some evidence has suggested that WRKY genes play a positive role in response to oxidative stress; for example, silencing of SlDRW1, a WRKY gene from tomato plants (Solanum lycopersicum), increased the sensitivity of transgenic plants to H 2 O 2 with less chlorophyll content in leaf discs . Overexpression of ThWRKY4 in Arabidopsis, a WRKY gene from tamarisk (Tamarix hispida), enhanced tolerance to oxidative stress (Zheng et al., 2013). To determine whether cassava WRKY genes play a role in oxidative stress response, the expression of 9 MeWRKY genes in response to H 2 O 2 was examined. Results suggested that MeWRKY18 and MeWRKY37 showed significant up-regulation at 2 and 6 h treatments, respectively. MeWRKY76 was significantly induced at 24 h, while MeWRKY49, -83, and -85 were seriously down-regulated at 24 h. MeWRKY28 and MeWRKY62 were strongly repressed at all the treated timepoints. MeWRKY65 did not show obvious trends after H 2 O 2 treatment (Figure 9). These results suggest that cassava WRKYs are likely to be involved in oxidative stress response.
Overall, the patterns in the expression of MeWRKYs under various conditions suggest that different MeWRKY genes may be involve in different signaling and stress responses, and that a single MeWRKY gene also participates in multiple signaling and stress processes. Moreover, most of the cassava WRKY genes can be quickly and significantly induced by multiple stressors, ABA, and H 2 O 2 treatments, indicating that WRKY FIGURE 11 | Expression profiles of MeWRKY genes in leaves under various stresses and ABA treatments. Log2 based values from three replicates of qRT-PCR data were used to create the heatmap. The scale represents the relative signal intensity values. Relative expression values for each gene after various treatments are provided in Figures 6-10 and Table S6.
Frontiers in Plant Science | www.frontiersin.org genes may function on multiple transduction pathways in cassava (Figure 11; Table S6).

CONCLUSIONS
In this study, 85 WRKY genes from the cassava genome were identified and their basic classification and evolutionary characteristics were established. This information may provide abundant resources for functional characterization of WRKY genes. The differential expression patterns of MeWRKYs in tissues of the wild subspecies and cultivated varieties revealed that they play different roles in cassava development, and a large number of them exhibited tissue-specific expression, thus assisting in understanding the molecular basis for genetic improvement of cassava. In addition, transcriptomic analysis of different cassava accessions associated with drought stress indicated that the majority of MeWRKYs in the root of W14 subspecies were activated in response to drought, which may contribute to its strong tolerance to drought. Furthermore, analysis of the expression of MeWRKY genes after various treatments suggested that they have a comprehensive response to osmotic, salt, ABA, H 2 O 2 , and cold, implying that cassava WRKYs may represent convergence points of different signaling pathways. These data may facilitate further investigation of WRKY-mediated signaling transduction pathways. Taken together, this work would provide a solid foundation for future functional investigation of the WRKY family in cassava.

Availability of Supporting Data
The cassava WRKY genes identified in this study was submitted to GenBank and the accession number was shown in Table S3. The transcriptomic data was submitted to NCBI and the accession number was shown in Table S7.

AUTHOR CONTRIBUTIONS
WH, KL, HS, and ZX conceived the study. YW, ZX, WT, ZD, YY, WW performed the experiments and carried out the analysis. WH, YW, and HS designed the experiments and wrote the manuscript. All authors read and approved the final manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 00025 Figure S1 | Sequence logos for conserved motifs identified in MeWRKYs by MEME analysis.
Table S1 | Primers used in qRT-PCR analysis.