Characterization and expression analysis of the WRKY gene family in moso bamboo

The WRKY family of transcription factors (TFs) is one of the ten largest families of TFs in higher plants and has been implicated in multiple biological processes. Here, we identified 121 WRKY TFs in moso bamboo, including five novel members that were not annotated in the Phyllostachys edulis genomic database. Estimation of the divergence time of paralogous gene pairs revealed an important role of the recent whole-genome duplication in the expansion of the WRKY family. Expression analysis based on quantitative reverse-transcription polymerase chain reaction (qRT-PCR) data revealed that a large number of PheWRKY genes varied significantly under cold or drought stress treatments, which could be defined as abiotic stress-responsive genes. The overexpression of PheWRKY72-2 in Arabidopsis resulted in a decreased sensitivity to drought stress during early seedling growth. PheWRKY72-2 may enhance plant tolerance to stress by functioning as a positive regulator of stoma closure. Our study provides a theoretical foundation and some experimental evidence for further functional verification of the PheWRKY family of TFs.


ID WRKY group
. However, to date, no studies have addressed the WRKY gene family and its function in the organ growth and development of P. edulis. Here, for the first time, we performed a detailed analysis of the subgroup classification, gene structure and conserved motif composition of 121 WRKY TFs in the P. edulis genome. We analysed the expression patterns of PheWRKYs from floral expression profiles and fast-growing shoot expression profiles to determine whether the expression patterns of PheWRKYs were related to fast-growing shoots and flower development in P. edulis. We further studied changes in expression of the PheWRKY genes in response to abiotic stress (drought or cold) in P. edulis seedlings. Furthermore, a drought stress and cold stress-inducible gene, PheWRKY72-2, was characterized in transgenic Arabidopsis to comprehensively understand its role in the response to abiotic stress. Our study provided a theoretical foundation and some experimental evidence for further functional verification of PheWRKYs.

Results
Identification of the WRKY family of TFs. We identified 121 WRKY members in the P. edulis genomic database (alternative splicing variants were not considered; Table 1), including five novel members that were not annotated in the P. edulis genomic database (Table S1). Because all the proteins and CDS sequences obtained from moso bamboo genome database were annotated by their putative orthologous sequences of O. sativa, so we renamed each moso bamboo WRKY sequences based on the individual similarity with rice WRKY proteins by blast method, and the different proteins of moso bamboo which corresponded to the same sequences in rice were suffixed with '−1' , '−2' , '−3' , '−4' depending on levels of sequences similarity (Table S2) 29 . The polypeptides encoded by the PheWRKY genes ranged from 159 to 1,168 amino acids in length, with predicted pI values ranging from 4.7 to 10.1 and molecular weights ranging from 19.77 to 131.06 kD (Table S3).

Structure and classification of PheWRKYs.
According to the zinc-finger structure (Table 1), sequence alignment ( Figure S1), conserved motif ( Fig. 1 and Figure S2) and gene structure ( Figure S3), the PheWRKYs can be classified into three main groups. Twenty PheWRKY proteins that contained two WRKY domains were assigned to group 1. Fifty of these contained two intact WRKY domains. However, the other five members of group 1 (PheWRKY35-2, 85-1, 96-1, 96-2 and 96-3) had only one complete WRKY domain in their N-terminals, and their C-terminal domains lacked a zinc-finger. PheWRKY125, which should have been placed in group 3 because of its C 2 HC type zinc-finger, was assigned to group 1 because it contained two WRKY domains. Except for PheWRKY125, the zinc-finger structures of the PheWRKY proteins in group 1 were of C 2 H 2 type with a C-X 4 -C-X 22-23 -H-X 1 -H motif (Table 1). Sixty-three PheWRKY proteins had a single WRKY domain, but in 62 members, the motif was C-X 4-5 -C-X 23 -H-X 1 -H; one member had the motif C-X 5 -C-X 88 -H-X 1 -H (PheWRKY9-1). All 62 members were assigned to group 2 and were further classified into five subgroups based on their conserved motif, their gene structure, and phylogenetic analysis. The subgroups included the following: group 2a (5 members), group 2b (11 members), group 2c (28 members), group 2d (7 members), and group 2e (12 members). The zinc-finger structure of group 3 was C 2 HC (27 members), with the motif C-X 6-7 -C-X 23-30 -H-X 1 -C (Table 1). Eleven proteins contained only partial WRKY domains and did not fit in any group. Most group 2c members shared motif 2 with most group 1 protein sequences, which also primarily harboured motifs 7  appeared only in group 2a, whereas motifs 8, 9 and 10 were found exclusively in group 2b. Group 2d harboured motifs 5, 6 and 11 (Fig. 1b). Structural analysis of the genes revealed that consisted with previous report 30 , gene family members within the same group shared similar gene structures in terms of their intron numbers and intron phases ( Figure S3). The intron phases of most sequences in group 1 were 0 at the first two introns and 2 at the last two introns. Except for PheWRKY97-1, the intron phases of all other members from groups 2a and 2b were 0. For groups 2c, 2d, 2e and 3, most introns phases were 2. Group 2b had more introns (average: 4.375) than the other groups. In contrast, almost all members of groups 2a, 2c, 2d, and 2e had only two introns. Generally, the motif and gene structure analysis revealed that members from different species assigned to the same group always shared similar motifs, intron numbers, and phases ( Fig. 1b, Figure S3).
Although the WRKYGQK motif was highly conserved, we found several sequence variations in 15 P. edulis WRKY proteins, most of which belonged to groups 3 and 2c (Table S4). We identified WRKYGKK as the most common variant in six domains, whereas WRKYGEK was common to five domains. The other four variants-WRKYGKK, CRKYGQA, WRKYGQQ, and WRKFGQK-were found in a single WRKY domain in PheWRKY7-3, PheWRKY80-1 (N-terminal), PheWRKY80-2 (N-terminal), and PheWRKY61 (C-terminal), respectively.  To confirm the specific zinc-finger structure of PheWRKY9-1 (C-X 5 -C-X 88 -H-X 1 -H), we designed a pair of primers at the 5′-UTR and 3′-UTR of PheWRKY9-1. Different RNA samples isolated from different tissues, including one-year-old leaves, two-year-old leaves, flowers, one-year-old culms, shoots, seeds, roots, and seedlings under various abiotic stress, were used as clone templates. Although we did not obtain the CDS sequence of PheWRKY9-1, which is annotated in the moso bamboo CDS database, from any of the tissues, we did identify a gene with a normal zinc-finger structure (C-X 5 -C-X 23 -H-X 1 -H). Compared with the annotated PheWRKY9-1 sequence from the moso bamboo CDS database, this gene lacked two exons ( Fig. 2a and b), as revealed by sequencing studies. In addition, whether the specific zinc-finger exists or is merely a mis-annotation could not be determined. edulis, A. thaliana, O. sativa, and Brachypodium distachyon, we performed phylogenetic analyses of the WRKY domain sequences from all four species based on a neighbour-joining method using Mega 6 software. Because the N-terminal and C-terminal domains form distinct clusters, we designated the two domains as 1 N and 1 C, respectively (Fig. 3). The phylogenetic tree indicated a divergence between monocotyledons and dicotyledons. PheWRKY clearly shared more sequence similarity with OsWRKY and BdWRKY than with AtWRKY. We then used the bidirectional best hit (BBH) method, which is restricted to a 1:1 ratio of orthologues 31 , to arrange possible orthologues for PheWRKY from the three sequenced species mentioned above. Of the 121 sequences, we identified orthologues for 62 PheWRKY genes from at least one of the three plant species. The highest number of orthologues (52 members) was identified in O. sativa, whereas the lowest number (10 members) was detected in A. thaliana (Table S2) Expression patterns of PheWRKY genes in different seedling tissues. We detected the expression of 81 PheWRKY genes in different tissues of three-month-old seedlings by RT-PCR. These selected genes covered all subgroups and included both up-regulated genes and down-regulated genes in flower development and shoot   Figure S4).

Expression analysis of the PheWRKY genes' roles in flower development and shoot growth.
We divided the PheWRKY genes into two groups (Fig. 4b) based on their expression in shoot and one-year-old culm. The expression levels of genes from cluster 1 were more abundant in the culm than in the shoot. In contrast, the expression of cluster 2 was significantly elevated in the shoot growth stage than in the one-year-old culm (CK). The genes PheWRKY28, 29-3, and 62 were significantly up-regulated during the shoot growth stage, suggesting that they play roles in the rapid elongation of the bamboo shoot (Table S6).
The expression profiles from the floral digital gene expression (DGE) data were divided into three cluster groups. Cluster 1 had genes with high expressions at stage F4. Cluster 2 genes exhibited steady expression with little fluctuation from stage F1 to F4. The genes from cluster 3 were initially up-regulated, reached their peak expression at stage F2 or F3, and declined thereafter. Compared with CK (leaf), PheWRKY48-1 appeared to be preferentially expressed throughout the flowering process ( Fig. 4a and Table S6). To validate the reliability of RNA-Seq and DGE data, quantitative reverse-transcription polymerase chain reaction (qRT-PCR) assays were randomly performed on four selected genes ( Figure S5). As expected, in most cases, the expression trends of the selected genes corresponded to the presented data, indicating that the DGE and RNA-seq data were highly reliable ( Figure S5).
Expression profiles of PheWRKY genes under abiotic stress conditions. We detected the expression of 81 PheWRKY genes under drought and cold stress by qRT-PCR. Of these, 32 genes showed different degrees of up-regulation with at least one stress treatment, including 28 genes responding to cold treatment and 22 genes responding to drought treatment (Fig. 5). Interestingly, the most up-regulated genes in drought treatment (17 of the 22) were always up-regulated in response to cold treatment. The expression levels of PheWRKY4, 26-2, 48-2, 49, 62, 69-1, 72-3 and 96-1 increased gradually and peaked at 12 h in response to cold treatment, whereas those of PheWRKY11-2, 111, 34-2, 45-1, 114-1, 69-2, 70-2 and 76 rapidly accumulated at 1 h and then decreased to low levels. For drought treatment, PheWRKY45-1, 62, 69-1, 88-1, 96-1, 96-2   Our expression analysis revealed that 54 genes were significantly induced by at least one abiotic stress (cold stress or drought stress) or physiological process (flower development or shoot growth). Of these WRKY genes, 7, 4, 17 and 3 were regulated by only cold stress, drought stress, flower development, and shoot growth, respectively ( Figure S7).

Overexpression of PheWRKY72-2 in Arabidopsis. Subcellular location analysis indicated that green
fluorescent protein (GFP)-tagged PheWRKY72-2 was located in the nucleus, in accordance with its function as a TF (Fig. 6).
To further investigate the potential role of PheWRKY72-2, which was up-regulated under both drought and cold stress, PheWRKY72-2 was transformed into Arabidopsis (wild type [WT]). Following qRT-PCR analysis of PheWRKY72-2 in WT and transgenic lines confirmed that PheWRKY72-2 performed high expression level in transgenic lines but no expression could be detected in WT (Fig. 7). Under cold stress treatment, no significant differences could be found between the WT and transgenic plants (data not shown). For simulated drought stress, transgenic Arabidopsis seeds were surface sterilized and germinated on Murashige and Skoog (MS) agar medium containing 10% and 20% polyethylene glycol (PEG) (Fig. 7b and e). The overexpression of PheWRKY72-2 (T3 generation) resulted in significantly longer roots than those of WT plants on medium supplemented with 10% PEG (Fig. 7b and c). We investigated the root hairs further but found no visible difference between those of the WT and transgenic lines (Fig. 7d). However, when the PEG concentration reached 20%, the WT line began to exhibit etiolation symptoms, unlike the transgenic lines (Fig. 7e).
We then compared the stomatal apertures under drought treatment because water loss mainly depends on stomatal regulation. To this end, we measured the stomatal apertures of WT and transgenic plants under control and drought treatments (3 h). The length and width of the stomata were determined, and the length/width ratio was used to reflect the degree of stomatal closure. Under normal growth conditions, no significant differences were detected (Fig. 8b). Under drought stress, the PheWRKY72-2ox lines showed a higher length/width ratio of stomata than the WT plants ( Fig. 8a and b). These results suggest that stomatal closure in PheWRKY72-2ox plants is more sensitive to drought stress than that in WT plants. Under the drought treatment, OE-PheWRKY72-2 suppressed the expression of AtABI2, AtABI5 and AtABA2 and significantly increased the expression of AtABI4, AtAREB and AtLEA (Fig. 8c).

Discussion
Characterization of the P. edulis WRKY gene family. To facilitate the analysis of P. edulis WRKY gene evolution in group 1, we further classified P. edulis group 1 WRKY genes into two subgroups based on their zinc-finger motifs 32 . Except for PheWRKY61, whose zinc-finger motif of C 2 HC-C 2 HC was assigned to group 1b, all the other members of group 1 could be allocated to group 1a. PheWRKY61 probably evolved through an intra-molecular duplication event of a group 3 WRKY domain that had already evolved to have the C 2 HC type zinc-finger. PheWRKY35, 85-1, 96-1, 96-2 and 96-3 all contained two WRKY domains but lacked a complete zinc-finger in their C-terminal domains. The closing genetic relationship of the N-terminal WRKY domains among PheWRKY85-1, 96-1, 96-2 and 96-3 determined based on the phylogenetic tree suggests that these sequences underwent some similar evolutionary events. Domain acquisition and domain loss events appear to have shaped the WRKY family of proteins 33 . Thus, these five genes might have arisen when a two-domain WRKY gene lost the zinc-finger motif on its C-terminal domain during evolution. Small variations in the WRKYGQK sequence are described for some OsWRKY proteins (accounting for 11 OsWRKY domains) 32 . These variations are also very common (11 domains) in B. distachyon, but only three A. thaliana members had variant heptapeptide sequences (WRKYGKK) ( Table S4). In moso bamboo, we found that amino acid substitutions in the conserved heptapeptide signature (15 domains) were much more common than in the three model plants.
In previous reports, bamboo was considered to have a tetraploid origin. Approximately 7-12 mya, P. edulis experienced a long progression from tetraploidy to diploidy 4 . P. edulis carries two duplicates as that of O. sativa gene model sets 4 . For PheWRKY genes, a similar phenomenon was always observed. Although P. edulis suffered a large-scale gene loss event after the whole-genome duplication, 36 putative paralogous pairs of PheWRKY genes remained, far exceeding the numbers in A. thaliana, O. sativa and B. distachyon. Furthermore, most duplication events of WRKY paralogous pairs in this species occurred 6-15 mya, much later than those of the three model plants mentioned above. All of our results suggest that the recent whole-genome duplication, which was likely linked to polyploidy events, played an important role in the WRKY family expansion.  35 . Similarly, 55 VvWRKY genes differentially respond to at least one abiotic stress treatment 36 . We profiled the expression of 80 PheWRKY genes in plants subjected to cold and drought stress and detected 31 P. edulis WRKY genes with significantly higher expression in response to at least one abiotic stress using real-time PCR, indicating that, in P. edulis, the PheWRKY genes function in the abiotic stress response. The overexpression of PheWRKY72-2 enhanced plant tolerance to drought stress. In Arabidopsis, AtWRKY54 and AtWRKY70 enhance plant tolerance to osmotic stress by regulating the size of the stomatal aperture 37 . Various environmental stresses result in the rapid accumulation of ABA, leading to stomatal closure, which reduces water loss by transpiration 38 . Recent reports showed that AtWRKY46 participated in the feedforward inhibition of lateral root inhibition via regulation of ABA signaling and auxin homeostasis under osmotic/salt stress treatment 39,40 . Besides, WRKY46 is regulated by light and modulates starch metabolism and ROS levels to inhibit stomatal closing under osmotic stress treatment by another independent way. Unlike the overexpression line of AtWRKY46 whose stomatal closing is impaired under osmotic SCIenTIfIC RepoRts | 7: 6675 | DOI:10.1038/s41598-017-06701-2 stress treatment 40 , the transgenic lines of PheWRKY72-2 showed enhanced stomatal closure under drought stress by comparing with the WT line. ABA controls stomatal movement via a dual mechanism. Regulation can occur via the biochemical effects of ABA on guard cells 38 . Salt stress is often accompanied by drought stress, and both cause water deprivation through ABA-dependent and ABA-independent pathways 41 . Our results showed that OE-PheWRKY72-2 significantly promoted the expression of AtLEA, AtAREB and AtABI4 and repressed the expression of AtABI2, AtABI5 and AtABA2. Plants with up-regulated ABI1/2 expression exhibited decreased stress tolerance, whereas plants with down-regulated ABI1/2 expression exhibited increased stress tolerance 42 . Reactive oxygen species (ROS), such as superoxide, hydrogen peroxide and hydroxyl radicals, play a dual role in plants: They act as necessary signalling molecules but can cause damage to plant cells when overproduced under stress conditions 43 . The constitutive expression of AtAREB1 in Arabidopsis modulates ROS accumulation and endogenous ABA levels to improve drought tolerance 43 . In plants, the presence of late embryogenesis abundant (LEA) proteins has been associated with cellular tolerance to dehydration, which may be induced by freezing, saline conditions, or drying 44 . Thus, up-regulating AtLEA might improve stress tolerance in transgenic plants.

Diverse PheWRKYs are involved in various biological processes. WRKY
Compared with their expression in nonflowering moso bamboo leaves (CK), the significant up-regulation of PheWRKY19-1, 19-2, 35-2, 48-1 and 65-1 during flower development indicated their potential roles in morphogenesis and organ development. WRKY TFs associated with senescence [45][46][47] and stress responses 48 were significantly up-regulated in the panicles of P. edulis. Under normal circumstances, P. edulis rarely flowers, and the plants often flower concurrently and die collectively after flowering. However, collective death after flowering in a large area is usually associated with external environmental changes, such as climatic variation or the over-exploitation of bamboo 1 48 . The involvement of WRKY factors in floral induction may be induced by signalling hormones, such as salicylic acid 7, 50 , jasmonic acid 45,46 or gibberellic acid 10 . We found that PheWRKY 28,49, 62, 80-3 and 111 have significant expression levels in growing shoots but not in culms. OsWRKY78 has positive effects on stem elongation and cell length, as determined by transgenic studies 51 . The up-regulation of several WRKY genes in shoot growth and development indicate that many WRKY genes function in the rapid growth and elongation of P. edulis shoots. In summary, PheWRKY genes play varied roles in different biological processes in P. edulis.
The functional conservative and divergence of orthologous genes. We further analysed the correlations of orthologous pairs under various processes. In Arabidopsis, the expression pattern of AtWRKY53 was identified as a leaf senescence signal 44 . Constitutively overexpressing AtWRKY53 under the control of the CaMV35S promoter accelerated the senescence symptoms, whereas WRKY53 RNAi and knock-out lines exhibited delayed development and senescence 45 . The orthologues, PheWRKY15-1 (4E-039) and PheWRKY74-1 (2E-31), were always found in high abundance in leaves collected from flowering plants (flowering implies senescence and death in P. edulis) compared to their abundance in leaves collected from other growth stages ( Figure S8), indicating that there are conserved roles between WRKY homologues in leaf senescence. Moreover, we compared the expression trends of PheWRKY genes with those of their possible orthologues from A. thaliana, O. sativa 34 , and B. distachyon 46 under cold and drought stress. Although the sample times and experimental conditions were different, our data revealed that some P. edulis homologues in other model plants showed an entirely opposite expression trend. PheWRKY85-1 shares 87% sequence similarity with OsWRKY85; however, OsWRKY85 is up-regulated under drought stress, whereas PheWRKY85-1 is down-regulated. The function of these segregated duplicated genes (PheWRKY85-1 and OsWRKY85) may have been altered during evolution 34 .

Conclusions
In this study, a total of 121 WRKY genes were identified in P. edulis. The expression profiles derived from DGE, RNA-seq and qRT-PCR analyses indicated that PheWRKY involved in various abiotic stress. The overexpression of PheWRKY72-2 in Arabidopsis increased drought tolerance by functioning as a positive regulator of stomatal closure. Better understanding this gene family and its potential involvement in growth, development and stress responses will facilitate further research on the evolutionary history and biological functions of the PheWRKY gene family.

Materials and Methods
Plant materials, growth conditions, and treatments. Samples of P. edulis culms were harvested in Huoshan County (116° 10′; 31° 12′), Anhui province from spring to autumn of 2014. Flowering plant samples from different developmental stages of the flowering process (F1-floral bud formation, F2-inflorescence development, F3-anthesis, and F4-embryo formation) were collected in Guilin city (110° 17′; 25° 048′), Guangxi Zhuang Autonomous Region, from April to August in 2014. Plants were sampled from sites suitable for bamboo growth that were free of insect pests and artificial destruction. Samples were collected according to previously described methods 1,3 .
For the abiotic stress treatment, P. edulis seeds were germinated in culture dishes lined with soggy filter paper for 1 week and then transferred into four seed pots containing approximately 0.5 kg of vermiculite. The seedlings were grown in an artificial climate chamber with long-day conditions (16-hrs light, 8-hrs dark) at 26 °C in the light and 18 °C in the dark. Drought stress was created by irrigating the plants with media containing 18% (v/v) PEG. Cold stress was created by placing plants in a 4 °C lighted growth chamber; the other conditions were the same as described above. The leaves from three seedlings were harvested at 0, 0.5, 1, 3, 6, 12 and 24 hrs after abiotic stress treatment, rapidly frozen in liquid nitrogen and stored at −80 °C prior to use. RNA isolation, reverse transcription and gene expression analysis. Total RNA was extracted using Trizol (Invitrogen, USA). First-strand cDNA synthesis was conducted with approximately 1 μg of RNA using the reverse transcriptase AMV (Promega, Madison, Wisconsin, USA). qRT-PCR was performed using the fluorescent intercalating dye Light Cycler 480 SYBR Green I Master Mix (Roche, Mannheim, Germany) on a Light Cycler 480 (Roche, Rotreuz, Switzerland). Primers were designed using Primer 3 software and Oligo dT7 (Table S7). According to the manufacturer's instructions, the 20-μL reaction mixture contained 0.4 μL (10 μM) of each primer, 1.5 μL (30 ng) of cDNA and 10 μL of SYBR Green I Master Mix. TIP41 (tonoplast intrinsic protein 41) 52 and β-tubulin were selected as internal controls. All reactions-technical and biological-were performed in triplicate. The ΔCT and ΔΔCT values were calculated by the formulas ΔCT = CT target − CT reference and ΔΔCT = ΔCT treated sample −ΔCT untreated sample, respectively.
Database searches. The  To avoid missing potential members of the WRKY family, we created a local blast database using the CDS sequence that we obtained from the P. edulis CDS database and novel transcripts that we obtained from RNA-seq data. The novel transcripts were assembled in a GTF file (exported by Cufflink based on a series of transcriptome sequencing results), and the P. edulis genome data were analysed using Perl. Published AtWRKY sequences and OsWRKY were retrieved and used as queries in blastn searches against the CDS database, and sequences were selected as candidate genes if their E-value was ≤ −10. The Pfam (http://pfam.sanger.ac.uk/search) and Smart (http://smart.emblheidelberg. de/) databases were used to confirm each predicted WRKY sequence. The pI and molecular weight were estimated using the Compute pI/Mw tool from ExPASy (http://web.expasy.org/compute pi).
Multiple sequence alignment and phylogenetic tree construction. Multiple alignments of the amino acid sequences were performed using the ClustalX 2.1 program with default settings 53 . A phylogenetic tree based on sequence alignment was generated using MEGA 6.0 (http://www.megasoftware.net/) by the neighbour-joining method 54 . In addition, the BBH method 31 was used to arrange possible orthologues and paralogues.
Gene structure and conserved motif analysis. The gene structure based on full-length CDS alignments with relevant genomic sequences was investigated using the online service of the Gene Structure Display Server 55 . MEME 56 was used to identify motifs in the PheWRKY sequences.
Estimation of the duplication time in paralogous pairs. The K a and K s values of the paralogous genes were computed by the DNASP program. The synonymous substitution rate (K s ) was considered as a proxy for time to estimate the dates of the segmental duplication events. The formula T = K s /2λ was used to o calculate the approximate date of the duplication event and the λ in formula represented for clock-like rates of synonymous substitution. The estimated clock-like rate for P. edulis, O. sativa and B. distachyon were 6.5 × 10 −9 substitutions/ synonymous site/year and that that for A. thaliana was 1.5 × 10 −8 substitutions/synonymous site/year 57, 58 . DGE data, transcriptome data analysis, and co-expression network generation. The DGE data of four different developmental stages of flowering (floral bud formation, inflorescence growing, blooming and embryo formation) and the transcriptomic data of fast-growing shoots at two developmental stages (shoots mixed by six unearthed heights [10, 50, 100, 300, 600, and 900 cm] and mature culm) were sequenced in our previous study 1, 3 and used for PheWRKY gene expression analyses. For floral development, the DGE data of nonflowering, moso bamboo leaves (CK) were regarded as the reference gene database 1 . For shoot growth, the transcriptome data of one-year-old culms constituted the reference gene database 3 . The expression levels of PheWRKY genes at different growth stages were hierarchically clustered based on Euclidean distance with complete linkage found in Cluster 3.0. A gene was considered up-regulated if it had a false discovery rate (FDR) ≤0.05 and│log 2 fold change│ ≥ 2.
The microarray data and RNA-seq data were used for coexpression analysis via the WGCNA package. The co-expression correlation was then calculated using Pearson's correlation with R 2 ≥ 0.95. The network picture was created using Cytoscape.

Overexpression of PheWRKY72-2.
The complete open reading frame (ORF) of the PheWRKY72-2 (FP101056.1) gene cloned by Chen 59 in our previous work was inserted into the pGEM-T Easy plasmid vector using the following primers: forward, 5′-TGCTCTAGAATGGAAGCCTACCCTATGCT-3′; reverse, 5′-CCGGAATTCTCAGTGGAACCGGCCAGAC-3′. TIANGEN DNA Polymerase (TIANGEN biological company, China) was used to amplify the PheWRKY72-2 gene. The PCR parameters were 94 °C for 5 min, followed by 28 cycles of 94 °C for 30 s, 62 °C for 1 min and 72 °C for 1 min. The PCR products was digested with BamH I and Xba I and then cloned into pGEM-T Easy (Promega, USA) after gel extraction. The coding region of the PheWRKY72-2 cDNA was cloned into the pCAMBIA 2300 vector under the control of the CaMV 35 S promoter and CAMV terminator. The constructed plasmid was named pCAMBIA2300-PheWRKY72-2 and confirmed by the chain-termination method on an ABI 3100 automated sequencer (USA).
The construct vector was introduced into WT Arabidopsis plants (Columbia-0) by Agrobacterium-mediated transformation. Transgenic plants were selected on kanamycin, and the first generation of transgenic plants was examined in terms of their phenotypes. At least 10 independent transgenic plants exhibiting severe phenotypes were selected, subjected to phenotypic characterization and screened for the transgene. Transgenic Arabidopsis seeds (T3 generation) and WT seeds were cultured with 10% and 20% PEG solution on MS agar medium. One-week-old seedlings were used to observe the phenotype. Three-week-old transgenic plants and WT which growth in normal growth condition were treated with PEG solution for 3 h and further used for RT-PCR experiment and stomatal closure observation. Subcellular localization. The subcellular localization of PheWRKY72-2 was performed by transfecting GFP-tagged PheWRKY72-2 into Arabidopsis sheath protoplasts. The full-length cDNA of PheWRKY72-2 was fused in frame with the GFP cDNA and ligated between the CaMV 35 S promoter and the nopaline synthase terminator. The fluorescence signals in transfected protoplasts were examined by a confocal laser scanning microscope (Leica Microsystems).
Availability of supporting data. All sequencing data were deposited in the Short Read Archive at NCBI under accession number SRR961047, SRR1187864 and SRR1185317.