Genome-wide analysis of WRKY transcription factors in wheat (Triticum aestivum L.) and differential expression under water deficit condition

Database search and identification of WRKYs

The protein sequences of 20 plants from nine different major taxonomic lineages were retrieved from several public databases. All of the amino acid sequences were obtained from the following sources: the eudicots Arabidopsis thaliana (At) and Populus trichocarpa, the monocots Brachypodium distachyon, Oryza sativa, Sorghum bicolor, T. aestivum (Ta), and Zea mays, basal magnoliophyta Amborella trichopoda, the bryophyte Physcomitrella patens, the lycophyte Selaginella moellendorffii, the chlorophytes Ostreococcus lucimarinus, and the rhodophytes Cyanidioschyzon merolae from the Ensembl Plants database (http://archive.plants.ensembl.org/info/website/ftp/index.html); the eudicots Cucumis sativus and chlorophytes Coccomyxa subellipsoidea C-169, Micromonas pusilla CCMP1545, and Volvox carteri and glaucophyte Cyanophora paradoxa from the JGI database (PhytozomeV9, http://genome.jgi.doe.gov/pages/dynamicOrganismDownload.jsf?organism=Phytozome#); the chlorophyte Ostreococcus tauri, gymnosperm Picea sitchensis, and rhodophyte Galdieria sulphuraria from NCBI (http://www.ncbi.nlm.nih.gov/protein/). The evolutionary relationship of these 20 species were obtained from NCBI (https://www.ncbi.nlm.nih.gov/Taxonomy/CommonTree/wwwcmt.cgi), and visually displayed by phylogenetic tree using FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

To identify the WRKY TFs in various species, the HMM profile of the WD (PF03106) downloaded from the Pfam database (http://pfam.xfam.org) (Finn et al., 2016) was applied as a query to search against the local protein database using HMMsearch program (HMMER3.0 software: http://hmmer.janelia.org/) (Finn, Clements & Eddy, 2011) with an E value cutoff of 1.0. The sequences obtained were then submitted to the Pfam database to detect the presence of WDs. The protein sequences containing complete or partial WDs, which may be pseudogenes, incomplete assemblies, sequencing errors, or mispredictions (Rinerson et al., 2015) were both considered as putative WRKYs. The physical and chemical properties including number of amino acids (NA), molecular weight (MW), theoretical pI, grand average of hydropathicity (GRAVY), aliphatic index (AI), and instability index (II) of putative TaWRKY proteins were calculated using the online ExPASy-ProtParam tool (http://web.expasy.org/protparam/).

Phylogenetic analysis

MEGA7.0 program was employed to construct the unrooted phylogenetic tree of identified WRKY protein domains in T. aestivum L. and A. thaliana L. using the maximum likelihood method (Kumar, Stecher & Tamura, 2016). The parameters of the constructed trees were: test of phylogeny: bootstrap (1,000 replicates), gaps/missing data treatment: partial deletion, model/method LG model, rates among sites: gamma distributed with invariant sites (G). Only bootstrap values greater than 60 could be displayed on the tree.

Chromosomal location of TaWRKY genes

To map the locations of WRKY gene transcripts in T. aestivum L., MapInspect software (http://www.softsea.com/download/MapInspect.html) was employed to visualize the chromosomal distribution of deduced TaWRKY genes according to their initial position and length of chromosome. The chromosomal location information of TaWRKYs was obtained from Ensembl Plants database (http://archive.plants.ensembl.org/Triticum_aestivum/Info/Index).

To detect the gene duplication, the CDS sequences of WRKY genes in wheat were blasted against each other (E value <1e⁻¹⁰, identity > 90%) (Song et al., 2014). Tandem duplicated TaWRKY genes were defined as two or more adjacent homologous genes located on a single chromosome, while homologous genes between different chromosomes were defined as segmental duplicated genes (Bi et al., 2016).

Characterization of gene structure, conserved motif, and putative cis-acting elements

The exon–intron structures of TaWRKY genes were obtained by mapping the CDS to DNA sequences using the Gene Structure Display Server2.0 (http://gsds.cbi.pku.edu.cn/) (Hu et al., 2015). CDS and genomic sequences in T. aestivum L. were retrieved from Ensembl Plants database (ftp://ftp.ensemblgenomes.org/pub/release-31/plants/fasta/triticum_aestivum/).

To discover motifs in TaWRKY protein sequences, the online tool Multiple Expectation Maximization for Motif Elication (MEME) 4.11.2 (http://meme-suite.org/) was utilized to identify the conserved motifs in full-length TaWRKYs (Bailey et al., 2009). The optimized parameters were as follows: distribution of motifs, 0 or 1 occurrence per sequence; maximum number of motifs, 10; minimum sites, 6; maximum width 60.

The 1.5-kb upstream of the transcription start site (−1) of all identified TaWRKY transcripts was extracted as promoter to predict cis-acting elements using the PlantCARE online (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002). Then statistics derived from hits of various cis-acting elements in all TaWRKY transcripts were constructed and displayed by diagram.

Plant materials, water deficit condition, and qRT-PCR

A hexaploid winter wheat (T. aestivum L.) cv. Zhengyin1 (St1472/506) was taken in our experiment, carried out from October 2015 to June 2016 in a greenhouse. Seeds were sown in each plastic pot filled with 7 kg of soil, earth-cumuli-orthic Anthrosols collected in northwest China. An equivalent of 0.447 g (urea)/kg⁻¹ (soil) and 0.2 g (K₂HPO₃)/kg⁻¹(soil) were mixed in soil with a net water content of 29.2% at the largest field water capacity. Water control was carried out from anthesis (April 17, 2016). Normal water supply and artificial soil desiccation were implemented with 70% to 75% (control group) and 45% to 50% (moderate water-deficit stress) of the largest field capacity, respectively. Spikes and flag leaves of wheat collected at 0, 1, 3, 5, 10, 15, and 25 days after anthesis (DAA), were immediately frozen in liquid nitrogen and then stored at −80 °C for subsequent analysis. Spikes were separated as glume, lemma, grain, palea, and rachis, and only glume and lemma were used in the experiment.

Total RNA was extracted from wheat tissues using the Trizol reagent (Tiangen, Biotech, Beijing, China) following the manufacturer’s instruction, and then digested with RNase-free DNase I. The quantity and concentration of RNA was evaluated by UV spectrophotometry. The first-strand cDNA was generated using PrimeScriptTM RT Reagent Kit (TaKaRa, Dalian, Liaoning Sheng, China), and the synthesized cDNA products were diluted 1:9 with nuclease-free water to use in qRT-PCR. Primer Primer 5.0 and AllelelID 6.0 (http://www.premierbiosoft.com/index.html) were used to design gene-specific primers (Table S1). Wheat Tublin was used as the reference gene.

The qRT-PCR was carried out using SYBR GreenSYBR^® Premix Ex Taq™(TaKaRa, Kusatsu, Shiga, Japan) according to the manufacturer’s instructions with Bio Rad CFX96TM real-time PCR detection system (BioRad, Hercules, CA, USA). Reaction parameters for thermal cycling were: 95 °C for 30 s, followed by 39 cycles of 95 °C for 5 s and 60 °C for 30 s, and finally a melting curve (65 °C to 95 °C, at increments of 0.5 °C) generated to check the amplification. The gene expression levels were calculated with the 2^−ΔΔCT method (Livak & Schmittgen, 2001), and three biological replicates were used.

Identification of WRKYs in wheat and comparative analysis

To comprehensively analyze and identify WRKY TFs in plants, 20 plants representing the nine major evolutionary lineages were chosen for analysis. After searching by HMMER and detecting WDs by Pfam database, a total of 1113 WRKY TFs were obtained (Table S2). The evolutionary relationships of various species and the number of WRKY TFs are shown in Fig. 1. Most terrestrial plants, including Monocots, Eudicots, and Bryophytes, contained 86 to 171 WRKY proteins, while Picea sitchenis, belonging to the Gymnosperms, carried only eight WRKY TFs, which could be due to incomplete sequencing. However, six or less WRKY proteins were found in aquatic algae, of which no WRKY TFs were found in Rhodophytes and Glaucophytes. In general, the number of WRKY TFs in many higher plants was more than that in lower plants, which suggested that the WRKY TFs may play an important role in the process of plant evolution. The number of WRKY proteins increased as plants evolved, possibly because of genome duplication.

In wheat, a total of 174 WRKY proteins were searched using the HMM search program. Subsequently, all obtained sequences were verified by Pfam database, which resulted in the identification of 171 WRKY TFs (Table S3). Based on their chromosome locations, we named 115 TaWRKYs as TaWRKY1 to TaWRKY115, and another 56 sequences were called TaWRKY116 to TaWRKY171 as they were anchored in the scaffolds. Among the 171 TaWRKYs, manual inspection showed that some were partial sequences, in which WDs or zinc-finger structures were incomplete or nonexistent.

The parameters used to describe the TaWRKY proteins were shown in Table S3. Molecular weight, theoretical pI, and aliphatic index could not be computed in sequences containing several consecutive undefined amino acids. The lengths of TaWRKY proteins ranged from 44 (TaWRKY121) to 1,482 residues (TaWRKY78), whereas the PI ranged from 4.96 (TaWRKY164) to 10.73 (TaWRKY40). This suggested that different TaWRKYs might operate in various microenvironments (Wang et al., 2014). The values of grand average of hydropathicity were all negative, which indicated that TaWRKY proteins were all hydrophilic. Almost all TaWRKYs were defined as unstable proteins, and only 30 TaWRKYs with instability index less than 40 were considered to be stable proteins.

Classification and phylogenetic analysis of TaWRKYs

To categorize and investigate the evolutionary relationship of the TaWRKY proteins in detail, we constructed an unrooted maximum-likelihood phylogenetic tree with 262 putative WDs in Arabidopsis and wheat (Fig. 2). Basing on the classification of AtWRKYs and primary amino acid structure feature of WRKY (Eulgem et al., 2000), we classified TaWRKYs into three major groups (Groups I, II, and III). The 30 TaWRKYs possessing two WDs and C ₂H₂-type zinc finger motifs (C–X_3–4–C–X_22–23–H–X₁–H) were classified into group I. Group II comprised 95 sequences, and each protein contained a single WD and C₂H₂-type zinc finger structure (C–X_4–5–C–X₂₃–H–X ₁–H). We further divided Group II into five subgroups, including IIa, IIb, IIc, IId, and IIe with 11, 7, 50, 17, and 10 members, respectively. Finally, 45 TaWRKYs with a single WD were assigned to Group III because of their C₂HC zinc-finger structure (C–X_6–7–C–X_23–28–H–X ₁–C).

As shown in Table S3, besides the highly conserved WRKYGQK motifs, we found three variants in TaWRKYs, namely WRKYGKK (10), WRKYGEK (11), and WSKYGQK (1), which were distributed in subgroup IIc, III, and TaWRKY157, respectively. In addition, two zinc-finger form variants, C–X₆–P–X₂₃–H–X–C and C–X6–F–X23–H–X–C were identified in TaWRKY80 and TaWRKY166, respectively. TaWRKY157, the most unique among all putative TaWRKYs, contained two WDs but with C₂HC-type zinc finger structure (C–X₇–C–X_23––H–X₁–C). The “Group I Hypothesis” sees all WRKY genes evolving from Group I C-terminal WDs (Rinerson et al., 2015). Therefore, TaWRKY157 could preliminarily be taken as an intermediate member of Groups I to III, although it was classified into Group III in the phylogenetic tree.

In this study, Group II was found to be the largest group of WRKY TFs in wheat. The members of Group II accounted for approximately 55.6% of all putative TaWRKYs, which was consistent with Musa balbisiana (Goel et al., 2016), pepper (Diao et al., 2016), and soybean (Song et al., 2016a). Subgroups IIa and IIb were separated from one clade, and IId and IIe were clustered to a branch, which is similar to previous studies in wheat (Okay, Derelli & Unver, 2014; Zhu et al., 2013).

Chromosomal location of TaWRKY genes

Among the 171 TaWRKY genes, 115 were mapped onto the 21 wheat chromosomes, and the other 56 were anchored in the scaffolds (TaWRKY116–171) (Table S3, Fig. 3). More TaWRKY genes were relatively distributed in Chromosomes 3B (18, 15.7%), 5B (11, 9.57%), 2A (9, 7.8%), and 5D (9, 7.8%). In contrast, chromosomes 6B, 7A, and 7B contained only one TaWRKY gene (0.870%). In general, most identified TaWRKY genes were observed in distal regions of chromosomes and only a few were observed in proximal regions. This phenomenon suggested that the TaWRKY genes were mapped on the all chromosomes with a significantly non-random and uneven distribution. The TaWRKY genes density in each chromosome ranged from 0.004/Mb (7B) to 0.056/Mb (5D) (Fig. S1).

Duplication events of WRKY genes have been found universally in a number of plants, such as peanut (Song et al., 2016b), white pear (Huang et al., 2015), and Brassica napus (He et al., 2016). In this study, we identified 79 TaWRKY gene duplication pairs which corresponded to 85 genes (Table S4, Fig. 3). This phenomenon indicated that some of the TaWRKY genes have more than one duplicated gene, which could be due to the multiple rounds of whole genome duplication in wheat. As shown in Fig. 3, two WRKY tandem duplication clusters (TaWRKY59/TaWRKY60, TaWRKY113/TaWRKY114/TaWRKY115) were identified on chromosomes 3B and 7D, respectively. In addition, 80 genes were found to have undergone segmental duplication, which were paralogs of WRKY genes on different chromosomes (Bi et al., 2016).

Gene structure analysis of WRKY genes in wheat

The exon–intron distribution was analyzed to further detect structural features of TaWRKY genes. Figure S2 showed that the number of introns in TaWRKY family genes varied from 0 to 5, while 0 to 8 in rice (Xie et al., 2005) and 0 to 22 in Musa acuminate (Goel et al., 2016), respectively. This phenomenon suggested that WRKYs in wheat show lower gene structure diversity. A total of 72 (42.11%) TaWRKY genes with two introns accounted for the largest proportion, followed by 44 (25.73%), 23 (13.45%), 15 (8.77%), 14 (8.19%), and 3 (1.75%) genes, possessing 1, 3, 4, 0, and 5 introns, respectively. The distribution pattern of introns and exons was group specific, which was similar to cassava (Wei et al., 2016) and carrot (Li et al., 2016), and TaWRKY gens belonging to the same subfamily shared a similar exon–intron structure. For example, TaWRKYs in Group III contained 0–5 introns, while approximately 91.11% (41/45) possessed 1–2 introns.

Two types of introns (V-type and R-type) were located in the WD characterized based on their splice site (Bi et al., 2016; Wang et al., 2014; Xie et al., 2005). V-type introns (phase 0) have a splice site before the V (valine) residue in C₂H₂ zinc finger structure, and R-type introns (phase 2) on the R (arginine) residue of the WD (Bi et al., 2016; Wang et al., 2014; Xie et al., 2005). In our study, all of the TaWRKY genes (17) in Groups IIa and IIb only contained V-type introns except TaWRKY165, which had no intron. However, R-type introns were mostly observed in all the other groups (Groups I, IIc, IId, IIe, and III) (Fig. S2). This phenomenon indicated that the intron phases were significantly conserved within the same group but remarkably different between groups (Chen, 2014). These results provided additional evidence to support the phylogenetic groupings and TaWRKYs classification.

Motif composition analysis of TaWRKYs

The conserved motifs of WRKY proteins in wheat were analyzed to explore the similarity and diversity of motif compositions. A total of 10 distinct motifs, named motifs 1–10, were detected using the MEME online program (Fig. 4). Among these 10 motifs, motifs 1 and 4 contained a WRKYGQK sequence, which is a basic feature of TaWRKYs. At least one of them contained almost all deduced TaWRKYs, except several incomplete proteins, such as TaWRKY11, 43, and 162. Motif 1 was observed almost in all groups, whereas motif 4 dispersed in Group I mostly.

As displayed schematically in Fig. 4, TaWRKYs within the same group or subgroup shared similar motif compositions. For instance, motifs 6 and 10 were unique to Group I, whereas motif 7 is specific to Group III. The motif unique to a particular group is likely to be involved in specific biological process in plants. Therefore, each family or subfamily of WRKY genes might be responsible for the specific biological process (Goel et al., 2016; Lippok et al., 2003). Furthermore, members of subgroups IIa and IId showed almost identical motif distribution patterns, indicating functional similarity among them. Interestingly, these two subgroups were also clustered to a branch in the phylogenetic tree. Likewise, the same phenomenon was also observed in subgroup IId and IIe. These results further validated the categorization of TaWRKYs and phylogenetic relationships.

Variety of cis-acting elements in promoter regions of wheat WRKY  genes

cis-acting elements in the promoter are crucial to gene expression, which is an essential part of its function (Dehais, 1999; Lescot et al., 2002). The 1.5 kb upstream promoter regions of all TaWRKYs were used to predict cis- acting elements using the online database PlantCARE. Here, various cis- acting elements were found in 142 out of 171 TaWRKY genes, while the remaining WRKYs could not be detected because of short sequence in their upstream regions (Table S5 , Fig. 5). Many cis- acting elements were related to response of hormones and biotic stresses, including MeJA, abscisic acid (ABA), SA, gibberellins (GA), auxin, zein, and fungus. MeJA-responsive elements with the largest portion were found in the promoter regions of 109 TaWRKY genes. Additionally, some elements involved in various abiotic stresses, such as light, wound, cold, heat, anaerobic induction, and drought, were identified in a large number of TaWRKY genes. A total of 44 light-responsive elements were almost distributed in all of the TaWRKYs. Some elements also observed in genes may regulate expression of different tissues (seed, root, shoot, leaf, phloem/xylem, endosperm, and meristem) in wheat development. Interestingly, a total of 76 TaWRKYs contained W-box (TTGACC), which regulates gene expression by binding WRKY, indicating these genes may auto-regulated by itself or cross-regulated with others (Chi et al., 2013; Jiang et al., 2014). MBSI, a MYB binding site involved in flavonoid biosynthetic genes regulation, only existed in TaWRKY87 and TaWRKY142, which suggested that these two genes may regulate flavonoid metabolism. Two unique genes were found, TaWRKY58 and TaWRKY94, which might respond to cold and water-deficit stresses for containing a cold and dehydration responsive element, C-repeat/DRE. Another special MYB binding site MBS that participated in drought response, were identified in 103 genes, indicating that most TaWRKYs seem to be involved in drought stress response (Table S5). Notably, all members analyzed contained more than one cis-element. Our analysis and previous studies both suggested that TaWRKY genes are involved in transcriptional regulation of plant growth and stress responses (Bakshi & Oelmüller, 2014; Ding et al., 2015; Raineri et al., 2016; Rushton et al., 2010).

Expression profiles of TaWRKY genes under water-deficit condition

With the exception of two shorter sequences (TaWRKY122 and 169), 12 out of 171 transcripts were selected as candidate drought responsive genes according to their orthologous WRKYs in Arabidopsis, which are involved in water deprivation, using the Biomart (http://plants.ensembl.org/index.html) (Table S6). The AtWRKYs responding to water-deficit stress were obtained based on function annotation in the TAIR database (http://www.arabidopsis.org/index.jsp). To validate these candidate drought-response genes, we determined their expression pattern in flag leaves, glumes, and lemmas using qRT-PCR. In our study, expression of genes could be detected at the transcript level in almost all selected tissues during the grain-filling period except TaWRKY8 (Fig. 6).

As shown in Fig. 6, we found that TaWRKY genes in glumes and lemmas share a more similar expression pattern compared with that in flag leaves. A relatively large group of genes, including TaWRKY1, 20, 31, 112, 123, 142, and 149 were significantly up-regulated in flag leaves at 0, 3, or 5 DAA, which suggested that these genes were highly induced at the early grain-filling stage (0–8 DAA). Among them, some genes, like TaWRKY123 and 142, were slightly up-regulated initially and then were restrained followed by an increase in the last point under water-deficit condition in glumes. In addition, peaks in the expression of several members (TaWRKY1, 20, 123, and 142) were mostly found at 5DAA in lemma, which lags behind other two tissues. Furthermore, water-deficit stress induced the most rapid up-regulation of some genes, and showed differences in three tissues. For example, both TaWRKY31 and TaWRKY149 were induced quickly in flag leaves after the onset of the water-deficit stress (0 DAA), approximately increasing up 6.67- and 8.22-fold, respectively. However, later induction was observed in glumes and lemmas. The immediate transcription response observed upon water-deficit stress appeared to be related to a more rapid perception of the drought (Eulgem et al., 2000; Rushton et al., 2010). Genes in another interesting cluster, composed of TaWRKY90, 97, 120, and 133, were down-regulated or slightly changed in glumes and lemmas during the early grain-filling stage, and induced during the middle (9–15 DAA) or late grain-filling stage (16–25 DAA) under water-deficit stress. However, two genes (TaWRKY120, 133) in flag leaves were strongly induced 3 days after water-deficit stress. This phenomenon indicated that TaWRKY120 and TaWRKY 133 genes were predominantly expressed in flag leaves at the early-filling stage. Our data suggested that the tissue-specific expression of TaWRKYs existed in wheat, and it appeared to be consistent with their role in tissues.