RNA-Seq Analysis Reveals MAPKKK Family Members Related to Drought Tolerance in Maize

The mitogen-activated protein kinase (MAPK) cascade is an evolutionarily conserved signal transduction pathway that is involved in plant development and stress responses. As the first component of this phosphorelay cascade, mitogen-activated protein kinase kinase kinases (MAPKKKs) act as adaptors linking upstream signaling steps to the core MAPK cascade to promote the appropriate cellular responses; however, the functions of MAPKKKs in maize are unclear. Here, we identified 71 MAPKKK genes, of which 14 were novel, based on a computational analysis of the maize (Zea mays L.) genome. Using an RNA-seq analysis in the leaf, stem and root of maize under well-watered and drought-stress conditions, we identified 5,866 differentially expressed genes (DEGs), including 8 MAPKKK genes responsive to drought stress. Many of the DEGs were enriched in processes such as drought stress, abiotic stimulus, oxidation-reduction, and metabolic processes. The other way round, DEGs involved in processes such as oxidation, photosynthesis, and starch, proline, ethylene, and salicylic acid metabolism were clearly co-expressed with the MAPKKK genes. Furthermore, a quantitative real-time PCR (qRT-PCR) analysis was performed to assess the relative expression levels of MAPKKKs. Correlation analysis revealed that there was a significant correlation between expression levels of two MAPKKKs and relative biomass responsive to drought in 8 inbred lines. Our results indicate that MAPKKKs may have important regulatory functions in drought tolerance in maize.


Introduction
The world population is increasing at an alarming rate, while abiotic stresses play crucial roles in crop failures and reductions in field crop productivity. Drought is the most important environmental stress affecting agricultural production. Maize (Zea mays L.), which is frequently exposed to drought stress conditions, is one of the most important cereal crops in the world, together with rice and wheat. In recent years, some studies have found that several genes encoding protein kinases can activate responses to abiotic stresses, such as drought [1,2]. Drought stress is one of the major limiting factors for maize production. Some studies have revealed that drought stress can affect maize, especially during the reproductive stage. Their working hypothesis involves signaling events associated with increased ABA levels, decreased glucose levels, the disruption of ABA/sugar signaling, and the activation of programmed cell death/senescence through the repression of a phospholipase C-mediated signaling pathway [3]. A total of 524 non-synonymous single nucleotide polymorphisms (nsSNPs) that were associated with 271 candidate genes for drought tolerance involved in plant hormone regulation, carbohydrate and sugar metabolism, signaling molecules regulation, redox reaction and acclimation of photosynthesis to environment were detected by common variants (CV) and cluster analyses with the availability of maize B73 reference genome and whole-genome resequencing of 15 maize inbred lines [4]. These selected genes will not only facilitate our understanding of the genetic basis of the drought stress response but will also accelerate genetic improvement through marker-assisted selection in maize.
The mechanisms of drought tolerance in plants are complex. Effects of drought on plant hormone signal transduction, protein modification and photosynthesis have been reported [1][2][3]5]. It has been demonstrated that the accumulation of proline can increase tolerance to water stress in plants. Transgenic plants with increased proline content produce higher biomass under drought [6]. An analysis of edr1 mutation determined that the EDR1 gene functions in negative regulation after salicylic acid treatment [7]. The stress-induced increase in reactive oxygen species (ROS) in plant cells results from an imbalance between generation and degradation. Drought stress can successively induce stomatal closure, moderate increases in ROS and decreases in photosynthesis. As described above, the response of plants to drought stress involves numerous genes and pathways related to diverse mechanisms.
Some molecular and biochemical studies have revealed that drought stress can elicit defense responses through MAP kinase pathways [8]. For example, a functional analysis of Arabidopsis demonstrated that the expression of AtMPK3, which shows high sequence similarity to WIPK, was increased when plants underwent drought [9]. The expression of OsMAPK5, OsMSRMK2, OsMAPKK44 and OsMKK1 was inducible by drought treatment in rice [10][11][12][13]. An analysis of transgenic rice plants strongly revealed that OsMAPK5 could regulate drought tolerance. Moreover, in maize, the transcript level of ZmMPK3 was increased by drought stress [14]. An investigation in alfalfa indicated that MKK4 kinase is transiently induced by drought treatment [15]. In addition, it has been reported that the tomato MAPK gene SIMPK4, which shows homology with Arabidopsis AtMPK4, can improve tolerance to drought stress [16]. In this study, we investigated the relationship between MAPKKK genes and drought stress in maize. Several previous studies have indicated that MAPKKK genes can be activated by drought treatment. In Arabidopsis, the mRNA level of the MAPKK kinase AtMEKK1, which can activate its downstream factors ATMKK2 and MEK1, can be increased by drought [9,17]. A MAPKKK gene that responds to drought has also been identified in rice. Furthermore, NPK1 isolated from tobacco [18] has been reported to positively regulate drought tolerance. NPK1 transgenic maize and rice both have higher yields under drought conditions compared with their nontransgenic counterparts [19][20][21].
In addition, it is well documented that MAPKKKs are related to many stress-response pathways, such as plant hormone signal transduction and oxidative signaling. The Raf-like MAPKKK gene in rice, DSM1, functions as an early signaling component to regulate the response to drought stress through ROS (reactive oxygen species) scavenging [51]. The results were identified by examining dsm1 mutants. DSM1-RNA interference lines were hypersensitive to drought stress and were more sensitive to oxidative stress. DSM1 may also function in the drought stress signaling response through an ABA-independent pathway that involves oxidative stress signaling [51]. Moreover, ANP1 has been shown to be involved in auxin signaling transduction and the oxidative stress signaling pathway [52,53]. In turn, oxidative stress can activate MAPK cascades to mediate the induction of oxidative stress-responsive genes. The overexpression of ZmMKK1 in Arabidopsis can enhance drought tolerance and enhance the expression of ROS-scavenging enzyme-related genes [54]. Analyses in maize have revealed the involvement of ZmMPK3 and ZmMPK5 in ROS signaling pathways [14,55]. Data from recent studies provide evidence that ZmMKK1, ZmMKK3 and ZmMKK4 can enhance drought tolerance through ROS scavenging [38,54,56]. Additionally, antioxidant enzyme activities can be enhanced in ZmMKK4-overexpressing plants compared with wild-type plants [56]. In brief, MAPK cascades can respond to drought stress through different pathways. However, the involvement of MAPKKKs in the drought response in maize is not yet clear.
In maize, 74 MAPKKKs had been reported by searching against the maize genome database and NCBI [57]. In this work, we will search against the maize B73 genome in the MaizeGDB to identify MAPKKK genes in maize using the BLAT program with a rigorous limiting condition. To further explore whether the identified maize MAPKKKs are related to drought stress, a variety of maize grown widely in SouthWest China, ZD619, and eight maize inbred lines were used in the RNA-seq, qRT-PCR and biomass analyses respectively. Additionally, important metabolic processes and regulatory process pathways related to the identified MAPKKK genes were analyzed to investigate the MAPKKKs involved in the drought stress response in maize.

Plant materials and stress treatments
In this experiment, the RNA-seq used ZD619, a variety of growing widely in SouthWest China, which showed strong drought resistance in the field. Eight inbred lines consisting of J24, J853, X178, Q319, B73, E28, C8605-2 and 200B as well as ZD619 were used to further analyze on their expression pattern by qRT-PCR and phynotype discrimination responsive to drought stress. The maize seeds were pre-germinated in the plant incubator. The uniformly germinated seeds were chosen and sown in plastic pots (23cm×16 cm) which filled with a 1:1:1 mix of soil: vermiculite:nutrient soil in greenhouse under 16h of light (25°C) and 8h of dark (20°C). Four uniformly strong plants finally retained at the 2-3 leaf stage.
The experiment included control (80% of relative soil water content, RSWC) and drought stress (35% of RSWC). Each treatment was comprised of two replications. The RSWC for the control plants exceeded to 80%. The drought treatment plants did not receive water from three-week-old seedling until the RSWC decreased to 35%, the drought treatment takes about one week, then kept stable for 2-3 days.

Plant biomass determination
Each plant except for qRT-PCR analysis was rapidly collected and weighed for the fresh biomass of two replicates in the greenhouse. Then the entire plant was put into a paper packet and was dried at approximately 75°C in a drying cabinet. The dry biomass was determined until 48 hours.

RNA isolation and qRT-PCR
The samples were collected and were immediately frozen in liquid N 2 for further use. Each total RNA sample of three tissues, leaf (the top three leaves), stem and root, was extracted from four uniformly plants using an RNeasy Plant Mini Kit followed by RNase-Free DNase Set (Qiagen, Hilden, Germany). The extracted RNA was assessed by a NanoDrop 2000 Spectrophotometer, 1% agarose gel electrophoresis. The poly(A) mRNA was isolated from purified total RNA using oligo(dT) magnetic beads. The mRNA was fragmented into short pieces by adding the fragmentation buffer. The first-strand cDNA was synthesized by reverse transcriptase and random primers using mRNA fragments as templates. Second-strand cDNA synthesis was produced by RNase H and DNA polymerase I. Double-stranded cDNA fragments went through end repair, 3'dA tailing and adapter ligation. The required fragments were then purified and enriched by PCR to create the final cDNA library. The sequences were generated by Illumina HiSeq 2000 sequencer.
In qRT-PCR three biological replicates were conducted and each biological replicate was technically repeated three times. The qRT-PCR was carried out using Maxima SYBR Green/ ROX qPCR Master Mix (2×) (Thermo Scientific) performed with a 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA) according to the supplier's protocols. Each PCR reaction mixture contained 12.5μl of 2×real-time PCR mix, 1μl of gene-specific primers, 1μl reverse transcribed cDNA product and water. The thermal cycle applied was as follows: 95°C for 10min followed by 40 cycles of 95°C for 15s, 60°C for 30s, 72°C for 30s, 82°C for 5s. At the end of the PCR cycles, melting curve analysis was performed using a single cycle consisting of 95°C for 15s and 60°C for 1min followed by a slow temperature increase to 95°C at the rate of 0.3°C/s. Relative gene expression was calculated according to the delta-delta Ct method of the system.

Definition of the maize MAPKKKs and phylogenetic analysis
Maize protein sequences were downloaded from MaizeGDB database (http://www.maizegdb. org/). The maize MAPKKKs was identified through BLAT searches against Arabidopsis, rice MAPKKK query sequences [25,35,37]. 25% identity was taken as the threshold, and the query sequences which aligned 5 BLAST alignments passed the threshold in MAPKKK gene family in Arabidopsis and rice were collected as MAPKKKs in maize.
Multiple sequence alignments were conducted on the amino acid sequences of MAPKKK proteins in maize, Arabidopsis and rice genomes using ClustalW with default settings. Phylogenetic tree was constructed by MEGA5.0 software based on alignments using Maximum Likelihood method with 1000 bootstraps to investigate the evolutionary relationship among MAPKKK proteins in three species. Analyze the conserved domain of the full protein sequences of MAPKKK genes by BioEdit.

Subcellular localization and chromosomal locations
The predicted subcellular localization of maize MAPKKKs was retrieved using the CELLO v2.5 server (http://cello.life.nctu.edu.tw/). Information about the physical locations of maize MAPKKK genes on chromosomes was obtained from MaizeGDB database. The chromosomal locations of the MAPKKK genes in maize were mapped on chromosomes using PERL script and Adobe Illustrator Artwork software.

RNA-seq analysis and identification of DEGs
The RNA-seq reads for each tissue were mapped to maize reference genome B73 using the Tophat [58]. Then the differentially expressed genes between drought processed sample and well-water processed sample in three tissues were extracted after processing Cufflinks and Cuffdiff (http://cufflinks.cbcb.umd.edu/manual.html).
We identified differentially expressed genes (DEGs) by fold change greater than two and pvalue 0.05.
Our raw data and the processed RNA-seq data have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus (GSE71377).

GO Analysis and Functional classification
We used an online tool Venn (http://bioinfogp.cnb.csic.es/tools/venny/index.html) to compare the DEGs in three tissues. The cluster was done by cluster3.0 using the Complete linkage method. The heatmap was drew by treeview and we chose the log2 (fold change) absolute value one as scale. Functional classification was performed by using the tool agriGO (http://bioinfo. cau.edu.cn/agriGO/). The functional clusters enrichment analysis was calculated by comparing with the whole maize genome V5a by Singular Enrichment Analysis (SEA) method, and the highest classification stringency was chosen for clustering by FDR 0.05.
We used an interactive ontology tool PageMan to generate overview graphs for profiling experiments with maize database and a user-driven tool MAPMAN to display genomics data sets onto diagrams of metabolic pathways and other biological processes. The Wilcoxon test likes a t-test was used to test the hypothesis that objects within one functional class behave differently from the rest of the objects. The significance of the change is reflected in the intensity and the direction by the color. The number used to do the heatmap was calculated by p-value. All p-values above 0.05 are set to a z-score of 0 to avoid misinterpretation. A highly saturated color indicates a high absolute value, whereas smaller values are indicated by lower color saturation. The two different colors can be selected to distinguish between categories where the average of the signals for all the genes in a category increases (red) or decreases (blue).

Co-expression analysis
We carried out pathway enrichment analysis by KAAS (http://www.genome.jp/tools/kaas/) for our DEGs in three tissues with Arabidopsis database. 14 common regulated pathways and 5 tissue specific regulated pathways were used for the heat map show. An R package WGCNA was used for our sequenced maize transcripts. We kept the transcripts by at least FPKM>0.1 in six samples. Firstly, we obtained the DEGs in each selected pathways and the co-expressed MAPKKKs and the co-expression score between each two genes. Secondly, the cutoff of 0.2 was used to select co-expression level in all the three tissues and a standard Z-score for each term is given by the following formula: where Sm is the mean coexpression value for genes selected with the cutoff, μ is the mean co-expression value of the entire gene list from WGCNA, m is the number of MAPKKK genes co-expressed with the DEGs, and δ is the SD of the co-expression value in the entire gene list [3]. A cutoff value of 0.2 was used to select highly co-expression events.

Genes specific primers
First-strand cDNA was used as a template in a qRT-PCR synthesized from 5μg total RNA using the SuperScript III First-Strand Synthesis System for RT-PCR (Invitrogen Corp.). The endogenous reference zSSIIb from zea mays starch synthase II was used as an internal control to normalize the data. The eight genes specific primers which designed based on the cDNA sequences were used to detect the gene expression level. The qRT-PCR primers were check the specification by NCBI primer-BLAST and validated by PCR(S4 File).

Physical location and phylogenetic analysis of MAPKKKs in maize
To identify MAPKKK genes in maize, the Arabidopsis and rice MAPKKK protein sequences were employed to perform a BLAST search against the protein sequences of maize B73. In maize, 71 MAPKKKs were identified, including 57 reported MAPKKKs [57] and 14 novel ones ( Table A in S1 File), which were basically consistent with previous research. These reported MAPKKKs were named according to the foregoing study [57], while since there was no standard nomenclature followed for MAPKKKs neither in Arabidopsis nor in rice, the 14 putative MAPKKKs were designated according to their group (Table A in S1 File).We obtained the chromosomal locations of these 71 MAPKKKs, and the predicted subcellular localization was determined from the CELLO v2.5 server (Table A in S1 File). The physical location data showed that 70 MAPKKKs were distributed on all 10 maize chromosomes (Fig 1D), except the GRMZM2G011070 (ZmMAPKKK29) gene, which was situated on an undetermined chromosome (chromosome unknown). Among the 14 novel MAPKKKs, GRMZM2G459824 (Zm MAPKKK75) was located on chromosome 1; GRMZM2G404078 (ZmMAPKKK76) and GRMZM2G335826 (ZmMAPKKK77) were located on chromosome 3; and GRMZM2G158860 (ZmMAPKKK78), GRMZM2G032619 (ZmMAPKKK81) GRMZM2G002531 (ZmMAPKKK85), GRMZM2G028709 (ZmMAPKKK86) and GRMZM5G852329 (ZmMAPKKK87) were located on chromosome 4. GRMZM2G127632 (ZmMAPKKK88) was located on chromosome 5. Chromosome 6 contained AC204050.4 (ZmMAPKKK79). Chromosome 7 contained GRMZM 2G034779 (ZmMAPKKK82) and GRMZM2G023444 (ZmMAPKKK83). GRMZM2G072395 (ZmMAPKKK84) was located on chromosome 8. GRMZM2G378852 (ZmMAPKKK80) was located on chromosome 9. The phylogenetic analysis was performed using the protein sequences in Arabidopsis, rice and maize. We aligned the full protein sequences of all the MAPKKKs in Arabidopsis, rice and maize using ClustalW, and we built a phylogenetic tree using MEGA5.0 (Fig 1). We identified 80 MAPKKKs in Arabidopsis and 75 MAPKKKs in rice [25,35,37], which could be subdivided into three major subtypes, Raf, MEKK and ZIK [30]. Based on the subtypes subdivided in a recent analysis, the MAPKKKs in maize could also be divided into three major groups, includ-

RNA-seq analysis of well-watered and drought-stressed maize transcriptomes in three tissues
The maize variety ZD619 was selected to analyze differentially expressed genes in response to drought stress (Fig 2A). RNA samples from leaf, root and stem tissues of ZD619 under wellwatered and drought stress conditions were isolated for next-generation sequencing using the Illumina HiSeq 2000 platform. The raw data included~40 million reads from each sample. The quality of the sequencing data is very well. The Q20 which represents for the percentage of the base whose quality is greater than 20 is approximately 98%. And the Q30 can also be around 93%. (Table B in S1 File) After trimming the adaptor and low-quality reads, approximately 85% of the processed data could be aligned to the B73 reference using TopHat [58] (Table C in S1 File). Then, using Cuffdiff, 5,866 differentially expressed genes (DEGs) were identified in leaf, stem and root (S2 File). Among them, 2,319 DEGs were identified in leaf, including 1,206 that were up-regulated and 1,113 that were down-regulated. 2,371 DEGs were identified in stem, including 1,057 that were up-regulated and 1,314 that were down-regulated. 2,181 DEGs were also identified in root, including 1,544 that were up-regulated and 637 that were down-regulated (Table 1, S3 File, Fig C in S5 File).
Through comparing the DEGs in leaf, stem and root, we found that only 74 genes were commonly regulated in all three tissues in response to drought, which suggested that these genes may play key roles during drought stress and that most of these genes were tissue-specific ( Fig 2B). The heat map of the DEGs also demonstrated significant differences among the three tissues ( Fig 2C). Out of those genes, there are 74 genes commonly regulated in the three tissues, approximately half of them were oppositely regulated in different tissues, which also indicated a large difference in the manner of regulation in different tissues (Fig 2D). Among the 74 highly significant co-regulated genes, GRMZM2G476477 (ZmMAPKKK20) is a member of the MAPKKKs, which showed significant up-regulation in all three tissues. In addition, it has been reported that GRMZM053669, GRMZM011598 and GRMZM059836 were up-regulated under dehydration conditions [59]. GRMZM2G051943 is a gene related to plant defense [60]. GRMZM2G162359 is a gene related to Aspergillus flavus pathogenesis [61]. GRMZM2G137108 and GRMZM2G455124 are required for leaf development [62]. GRMZM2G419239 and GRMZM2G341410 are involved in cell wall-related biogenesis [63].
Then, a GO analysis was conducted with respective DEG terms in the three tissues using the agriGO tool (see methods). Four GO classifications were commonly enriched after comparison among the three tissues, including response to stimulus, oxidation-reduction, polysaccharide metabolic process and metabolic process. The numbers of genes enriched in the three tissues for each commonly enriched GO classification were similar, and only a fraction of these overlapped among the three tissues. For example, there were 190, 167, and 167 genes enriched in leaf, stem and root, respectively, in the "response to stimulus" classification, but only 8 genes were the same, which accounted for approximately 5% of the number in each gene set. Correspondingly, there were 5.8%, 6.3%, and 4.4% commonly enriched genes in leaf, stem and root for oxidation-reduction, polysaccharide metabolic process and metabolic process, respectively. Seven GO classifications were enriched in two tissues. For example, response to chemical stimulus, response to stress and carbohydrate metabolic process were enriched in leaf and root; and cellulose biosynthetic process, cellulose metabolic process, multicellular organismal process and response to abiotic stimulus (which contains drought stress) were enriched in leaf and stem (Fig 3). Furthermore, 74 GO classifications were tissue specific, especially metabolismrelated classifications (Fig B in S5 File). Our results suggested that different tissues of plants may have similar stress signal reception and conduction pathways, yet the tissue-specific genes had a greater impact on plant growth and development.

MAPKKK genes revealed apparently different expression profiles in different tissues
Among all the MAPKKK genes, four genes (ZmMAPKKK18, ZmMAPKKK19, ZmMAPKKK20 and ZmMAPKKK56) in leaf, three genes (ZmMAPKKK19, ZmMAPKKK20 and ZmMAPKKK 21) in root, and seven genes (ZmMAPKKK18, ZmMAPKKK19, ZmMAPKKK20, ZmMAPKKK 21, ZmMAPKKK22, ZmMAPKKK26 and ZmMAPKKK73) in stem revealed significantly differential expression in response to drought stress. A comprehensive analysis of these three tissues revealed that eight MAPKKK genes may be involved in the response to drought stress. We conducted a cluster analysis of the three MAPKKK sub-families' expression profiles, including the ZIK family genes (Fig 4A), the MEKK family genes (Fig 4B), and the Raf family genes (Fig 4C). The data suggested that most of the differentially expressed MAPKKKs belonged to the MEKK family. We detected an important gene, GRMZM2G165099 (ZmMAPKKK19), which was predicted to function similar to NPK1 (Nicotiana protein kinase), in the MEKK sub-family. The expression of Nicotiana protein kinase (NPK1) has been shown to enhance drought tolerance in transgenic maize [19,20]. In our results, ZmMAPKKK19 showed an up-regulation pattern in all three tissues, especially in root, which is in accordance with the previous research. According to the subcellular localization results (Table A in S1 File) and the expression patterns shown in Fig 4B, the highly up-regulated MAPKKKs in leaf, ZmMAPKKK17, ZmMAP KKK18, ZmMAPKKK19, ZmMAPKKK20 and ZmMAPKKK79, were all located in the chloroplast, which is similar to the genetic mode of their expression in leaf. The results indicated that MAPKKKs may play an important role in drought-related stress. In addition, the pattern of the functional MAPKKKs displays clearly high expression under drought. Some of the differentially expressed genes in leaf, root, and stem of ZD619 that were identified using RNA-seq were confirmed by qRT-PCR to validate whether these genes were connected with the drought response ( Fig 4D). In leaf, the expression levels of ZmMAPKKK18 and ZmMAPKKK56 were markedly up-regulated, with greater than 4-fold increases, whereas the changes in ZmMAPKKK19 and ZmMAPKKK20 were not significant and were less than twofold. The expression of genes in stem showed a slight increase in ZmMAPKKK20. The expression levels of ZmMAPKKK18, ZmMAPKKK19, ZmMAPKKK21 and ZmMAPKKK22 were up-regulated in stem.

DEGs were enriched in different pathways in response to drought stress
As the above results showed, different tissues showed different responses to drought stress, and the DEGs in the three tissues represented distinct functional categories. Several MAPKKKs were including among the DEGs affected by drought. We investigated how drought regulated the differentially expressed genes, particularly the MAPKKKs, using PageMan to categorize the drought-mediated gene expression data for leaf, stem and root tissues into known metabolic pathways and different regulatory processes (Fig 5).
There were many DEGs enriched in hormone signaling, photosynthesis, CHO metabolism, secondary metabolic processes and regulation of transcription (Fig 5). Genes enriched in hormone signaling included more up-regulated genes in leaf and more down-regulated genes in root. Genes enriched in photosynthesis included many down-regulated genes in leaf. Genes enriched in CHO metabolism included many up-regulated genes in root. Genes enriched in secondary metabolic processes included more down-regulated genes in the three tissues. Genes enriched in regulation of transcription included more up-regulated genes in leaf and stem. To identify important enriched pathways in our DEGs, we analyzed an overview of enrichment in the metabolic processes TCA, starch synthesis and decomposition; proline metabolic pathways (Fig D in S5 File); and hormone signaling pathways (Fig E in S5 File). Genes that were significantly involved in metabolism were more down-regulated in leaf but more up-regulated in root. Genes that were significantly involved in regulation were more down-regulated in stem but up-regulated in leaf and root. These results suggested that the differential regulation of genes in different tissues may lead to particular patterns of regulation.
Furthermore, the overview of regulation revealed three regulated MAPKKKs (Zm MAPKKK19, ZmMAPKKK20 and ZmMAPKKK56) in leaf, four regulated MAPKKKs (Zm MAPKKK20, ZmMAPKKK21, ZmMAPKKK22 and ZmMAPKKK26) in stem, and four regulated MAPKKKs (ZmMAPKKK17, ZmMAPKKK19, ZmMAPKKK20 and ZmMAPKKK21) in root. ZmMAPKKK19 and ZmMAPKKK20 were common in at least two tissues. They exhibited distinctly differential expression after drought stress (Fig 4B). According to the phylogenetic tree (Fig 1B), ZmMAPKKK17, ZmMAPKKK19, ZmMAPKKK20, ZmMAPKKK21 and ZmMAPKKK22 belong to the MEKK subfamily. These results indicated that MAPKKKs may have important regulatory functions in drought tolerance in maize.

Co-expression analysis of DEGs and MAPKKKs
To examine patterns of correlation between the differentially expressed genes and MAPKKKs, we investigated their co-expression using WGCNA [64][65][66]. All 26,901 filtered transcripts were clustered using WGCNA to obtain an overview of their expression relationships. Networks and models for our data sets were constructed to determine the co-expression value between each pair of genes from our filtered data sets.
We focused on the MAPKKKs that were especially co-expressed with DEGs that were enriched in important pathways. After filtering with at least 10 genes, the most significant pathways could be identified, and most of these pathways coexisted in all three tissues. Then, we chose fourteen common and five tissue-specific pathways for further analysis. The heatmap in the left of Fig 6 shows the different co-expression patterns in the three tissues based on their Z-scores. The plant hormone signaling pathway categories oxidation stress, photosynthesis, starch and sucrose metabolism, arginine and proline metabolism, TCA cycle, auxin, cytokinin, gibberellin, ABA, ethylene, BR, jasmonic and salicylic acid presented high Z-scores in the coexpression analysis.
To better understand the relationships between the DEGs and MAPKKKs in the selected pathways, a heat map was constructed to study the expression patterns of single DEGs that were co-expressed with MAPKKKs corresponding to the nineteen enriched pathways (right in Fig 6). The data revealed that the DEGs were highly related to the MAPKKKs (Table D in S1 File). The DEGs that were co-expressed with MAPKKKs in these pathways also showed a broad differential expression pattern. It was found that most of the DEGs were up-regulated in root. In the high Z-score pathways, the photosynthesis pathway contained most of the downregulated pathway-enriched DEGs, whereas the proline, ethylene and salicylic acids pathways contained many up-regulated pathway-enriched DEGs, especially in root.
Furthermore, thirteen co-expressed MAPKKKs were co-expressed with DEGs in the above nineteen enriched pathways; these MAPKKKs were ZmMAPKKK11, ZmMAPKKK12, ZmMAPKKK13, ZmMAPKKK22, ZmMAPKKK34, ZmMAPKKK38, ZmMAPKKK39, ZmMAPKKK51, ZmMAPKKK56, ZmMAPKKK65, ZmMAPKKK72, ZmMAPKKK74 and ZmMAPKKK83 (Table E in S1 File). Although some MAPKKKs were not found to be coexpressed with DEGs in any of the nineteen pathways, they may play roles in minor pathways under drought stress; these genes included ZmMAPKKK19 and ZmMAPKKK21, among other. ZmMAPKKK26 and ZmMAPKKK73 play roles in several pathways. Our results confirmed that ZmMAPKKK19, ZmMAPKKK21, ZmMAPKKK22 and ZmMAPKKK56 were up-regulated in different tissues under drought conditions. ZmMAPKKK56 in leaf and ZmMAPKKK22 in stem revealed significantly differential expression in response to drought stress. ZmMAPKKK21 in stem and root, ZmMAPKKK19 in all three tissues showed significantly differential expression under drought conditions. These MAPKKKs may be key regulators of the tolerance of plants to drought. The co-expressed DEGs in these nineteen pathways were very peculiar. Only three DEGs were co-expressed with MAPKKKs both in the oxidative and photosynthesis pathways, only five DEGs were co-expressed with MAPKKKs in both the oxidative and TCA pathways. The oxidative pathway, especially ROS, may be the most core element in drought stress. All the results revealed that MAPKKKs provide a link between gene response expression and drought tolerance. The drought stress-responsive MAPKKK genes may be crucial to maize under drought conditions.

Differential expression levels of the eight MAPKKK genes in response to drought stress in maize
To understand the differences of eight MAPKKK genes in response to drought stress among various maize varieties, a hybrid ZD619 and eight inbred lines of maize were used to further qRT-PCR analysis and biomass determination ( Table 2). Plant biomass and kernel yield are usually used as integrative indicators to judge the drought-tolerant potential of plants [67]. In this study, seedling biomass was regarded as a criterion to identify drought tolerance capability. The maize seedling biomass was measured as fresh weight and dry weight and subsequently calculated as drought tolerance index (DTI) to compare the drought response performance of various inbred lines. Results indicated that drought stress exhibited enough influence to elicit differential phenotypes. Each inbred line showed different drought resistance capability, among them X178 was the most tolerant one to drought stress and E28 was the most sensitive one to drought stress.
Eight MAPKKK genes for expression analysis by qRT-PCR were shown in Fig 7. The expression levels of the eight genes differed among the nine varieties. The expression levels of ZmMAPKKK18 and ZmMAPKKK56 were up-regulated in response to drought stress in all varieties. The expression levels of ZmMAPKKK19, ZmMAPKKK20 and ZmMAPKKK73 were up-regulated under drought conditions in seven inbred lines except E28. In addition, ZmMAPKKK21 was up-regulated in seven inbred lines except B73. The expression levels in various inbred lines had significantly different. The inbred line X178, in expression level increased by 58-fold in ZmMAPKKK18, 20-fold in ZmMAPKKK19, 11-fold in ZmMAPKKK56, 8-fold in ZmMAPKKK20, 7-fold in ZmMAPKKK22, and 4-fold in ZmMAPKKK21. In the inbred lines E28 and B73, the expression levels were not significantly changed, and some genes were even down-regulated. In E28, the expression levels of ZmMAPKKK19, ZmMAPKKK20 and ZmMAPKKK73 were markedly down-regulated. The results indicated that different expression levels of differential expression MAPKKK genes had certain relation with the variety's drought resistance ability.
We then performed a correlation analysis of the different expression levels of eight MAPKKK genes in eight inbred lines. Pearson's correlation coefficients were calculated between each gene's expression pattern and the plant biomass DTI values in the eight inbred maize lines ( Table 2). The results suggested that the expression levels of the genes GRMZM 2G305066 (ZmMAPKKK18) and GRMZM2G063069 (ZmMAPKKK56) were significantly correlated with changes in plant biomass in response to drought (P<0.05). The expression levels Fig 5. Pathway enrichment of differentially expressed genes involved in different regulatory processes under drought stress. A q-value cutoff of 0.05 was used to select enriched gene sets in all three tissues. The heat map represents the Z-scores obtained from a parametric analysis of gene set enrichment q-values for term enrichment. Red represents enriched genes in the treatment group that were over-represented compared with the control set. Blue represents the enriched genes in the treatment group that were under-represented compared with the control set. The absolute values represent the enrichment level. The bar represents the Zscore region from -3 to 3. of GRMZM2G165099 (ZmMAPKKK19) and GRMZM2G476477 (ZmMAPKKK20) were slightly correlated with plant biomass (P<0.1). This analysis suggested that the eight MAPKKK genes, especially ZmMAPKKK18 and ZmMAPKKK56, were obviously related to biological characteristics of maize under drought stress. Among the eight genes, according to this analysis, ZmMAPKKK18 and ZmMAPKKK56 may have crucial functions in the response of maize to drought stress.

MAPKKK genes are related to drought stress in maize
Plants have developed complicated signaling pathways to adapt to environmental stress. In some species, it has been verified that MAPKKK can be rapidly activated during suffering to drought stress. In Arabidopsis, the mRNA levels of AtMEKK1, which is structurally related to MAPKKK, and AtMPK3, which is structurally related to MAPK, has increased in response to dehydration stress [9]. In rice, a Raf-like MAPKKK gene, DSM1, can respond to drought stress through ROS scavenging [51]. The NPK1 (MAPKKK) isolated from tobacco is conserved among different organisms [18]. Transgenic rice and maize with NPK1 showed significantly improved tolerance to drought [19][20][21].
In this study, seventy-one MAPKKK genes have been obtained by a computational analysis of the entire maize genome and eight MAPKKK genes that exhibited significantly differential expression in response to drought stress are identified by an RNA-seq analysis. Among these Fig 6. Co-expression analysis of tissue-specific DEGs and MAPKKK genes. The weight value obtained from the WGCNA package was used as a parameter for the parametric analysis of gene co-expression levels. A cutoff of 0.2 was used to select highly co-expressed genes in all three tissues. The left side of the heat map represents the Z-scores obtained from a parametric analysis of gene co-expression. The lower left bar represents the degrees of the Zscore scale. The right side of the heat map represents the expression patterns of the DEGs co-expressed with MAPKKKs that were enriched in the corresponding pathways on the left. The lower right bar represents the log2 of the drought/control ratio. genes, ZmMAPKKK26 and ZmMAPKKK73, which show down-regulated expression in stem under drought stress, are located in the nucleus. ZmMAPKKK18, ZmMAPKKK19, ZmMAPKKK20 and ZmMAPKKK22, which are up-regulated under drought stress, are located in the chloroplast. Four MAPKKKs are predicted with a serine/threonine-protein kinase function. Notably, ZmMAPKKK19 is predicted to have NPK1-related protein kinase function. In addition, ZmMAPKKK19 is highly homologous with ZmMAPKKK18. Indeed, the expression levels of these two genes are up-regulated after drought treatment. The results indicate that ZmMAPKKK18 and ZmMAPKKK19 may be involved in the drought response. These proteins have high sequence similarity, which may means the correlation with similar expression patterns. Furthermore, qRT-PCR is performed with eight MAPKKK genes to compare their expression levels in different varieties between drought conditions and well-watered conditions. The results show that the expression levels of these eight MAPKKKs in the nine maize varieties are different. In all the varieties, the relative expression levels of ZmMAPKKK18 are markedly upregulated when maize seedlings are under drought stress condition. Additionally, the relative expression levels of ZmMAPKKK19, which is predicted to have NPK1-related protein kinase function, are up-regulated in seven inbred lines. The genes ZmMAPKKK18 and ZmMAPKKK 19 may play a crucial role in response to drought stress. However, to elucidate the roles of ZmMAPKKK18 and ZmMAPKKK19 in drought tolerance, many details need to be clarified by further research. A comparison of these genes' expression among the different maize lines show that the abilities of genes to tolerate drought stress probably differ among the different lines.
In our study, the analysis of drought tolerance capability which identified by DTI values indicates that drought stress exhibit enormous influence for the different phenotypes. Among the eight inbred lines, X178 was the most tolerant line to drought stress and E28 was the most sensitive line to drought stress. The correlation analysis reveals that the eight MAPKKK genes, especially ZmMAPKKK18 and ZmMAPKKK56, are significantly correlated with the biological characteristics under drought stress in maize. The data suggests that, among the eight genes, ZmMAPKKK18 and ZmMAPKKK56 may have crucial functions in the response of drought stress in maize.
In maize, the expression levels of eight MAPKKK genes are regulated under drought stress, and the results are confirmed in different maize varieties. These observations suggests that the MAPKKK genes may be involved in the drought stress response in maize.

MAPKKK genes are involved in multiple pathways under drought stress
Plants tend to adapt to drought stress by serial changes. The mechanisms of drought tolerance in plants are complex. The data in our study suggests that drought stress can lead to expression changes in some genes, including MAPKKKs. Under drought stress, several MAPKKKs, especially those that co-expressed with DEGs, are enriched in some important pathways, such as the photosynthesis pathway, hormone signal transduction pathway, oxidation pathway, and protein modifications(including proline) pathway (Fig 8). In the photosynthesis pathway, To determine whether the relative expression levels of the drought stress-responsive MAPKKK genes differed among varieties induced by drought, ZD619 and the 6 inbred maize lines J24, J853, X178, E28, C8605-2, 200B, Q319 and B73 were used. Lines X178, J24 are drought-resistant lines; 200B and E28 have poor drought tolerance. The eight MAPKKK genes are GRMZM2G305066 (ZmMAPKKK18), GRMZM2G165099 (ZmMAPKKK19), GRMZM2G476477 (ZmMAPKKK20), GRMZM2G173965 (ZmMAPKKK21), GRMZM2G041774 (ZmMAPKKK22), GRMZM2G021416 (ZmMAPKKK26), GRMZM2G063069 (ZmMAPKKK56), GRMZM2G474546 (ZmMAPKKK73). * indicates significant differences in comparison with the control at P < 0.05 respectively. Error bars indicate standard deviation for three replicates. there are more down-regulated DEGs in leaf. There are many up-regulated DEGs in the three tissues, especially in root. There are obviously more up-regulated DEGs in oxidation and protein modification (including proline) pathways, primarily in root. It is well known that drought stress encumbers photosynthetic carbon fixation by limiting the entry of CO 2 into leaves via down-regulating the photosynthetic metabolism or accelerating stomatal closure. Furthermore, proline synthesis moderately increases under drought stress, which can help plants to resist drought stress. The accumulation of proline in transgenic plants is associated with higher biomass under drought stress [5,6]. In our study, the enriched DEGs are down-regulated in the photosynthesis pathway in leaf and are up-regulated in proline modification, which are in accordance with previous studies. Some MAPKKKs are found to be co-expressed with DEGs in these pathways, which suggests that MAPKKKs are involved in the photosynthesis pathway and proline modification pathway under drought stress condition. Among them, the downregulated ZmMAPKKK26 and ZmMAPKKK73 are co-expressed with DEGs in the photosynthesis pathway, while the up-regulated ZmMAPKKK19, ZmMAPKKK21, ZmMAPKKK22 and ZmMAPKKK56 are co-expressed with DEGs in proline modification pathway.
Moreover, the accumulation of ethylene and salicylic acid may be caused by water stress. Several earlier investigations have found that some ethylene response factor genes and endogenous salicylic acid levels can induce water stress tolerance in plants [68][69][70][71]. Seed which is pretreated with SA promotes plant vigor index under drought stress condition [60]. A study in Medicago and Arabidopsis plants showed that a MAPK pathway was involved in ethylene signaling [61]. Drought inhibits maize growing by causing oxidative damage to biological membranes and disturbing the water status in tissues. Reactive oxygen species (ROS) are central signaling molecules in response to drought. To reduce the toxicity of ROS, plant cells have developed an antioxidative system. The increasing ROS to toxic levels can activate PCD to remove damaged cells. Ethylene and salicylic acid response factor proteins may be connected to the ROS pathway through defending cells against ROS attacks. Salicylic acid can enhance the antioxidant status and the tolerance to drought stress [62]. Ethylene plays an important role in lowering ROS generation and protecting plants under drought stress [63]. Ethylene can modulate the programmed cell death induced by O 3 exposure [72]. Many studies have reported that MAPK cascades play crucial roles in ROS signaling pathway under water stress. MAPK, which is induced by oxidative stress, can promote ROS scavenging. Arabidopsis MKK4 increases plant drought tolerance and decreases the production of ROS [73]. Tolerance to drought stress, which is enhanced by NPK1, may activate the oxidative signaling pathway [19]. In tobacco, GbMPK3 increases ROS scavenging under drought stress [74]. The MAPKKKs that are activated under drought stress may promote ROS scavenging to control the concentrations of ROS. In our study, DEGs are enriched in the ethylene and SA hormone signal transduction pathways. The expression levels of most of the genes related to salicylic acid and ethylene are up-regulated. The up-regulation of proteins may be involved in redox homeostasis, ROS scavenging and its damaging effects on proteins and DNA. In addition, some MAPKKKs are coexpressed with DEGs in ethylene and SA plant hormone signal transduction pathways, such as ZmMAPKKK19, ZmMAPKKK21, ZmMAPKKK22 and ZmMAPKKK56. The data supports the theory that MAPKKKs may promote ROS scavenging by the ethylene and SA pathways in response to drought. The ethylene and SA hormone signal transduction pathway-related genes, especially in root, show an up-regulation pattern, which reveals that roots may play a key role in response to drought stress. Under drought condition, MAPKKKs are activated by ROS, which leads to the differential expression of some genes. The expression levels of eight MAPKKK genes are increased under drought treatment. The genes ZmMAPKKK18 and ZmMAPKKK19, which are highly homologous to NPK1, are up-regulated in leaf, stem and root, which may play roles in regulating the response to drought stress. The differentially expressed MAPKKKs, ZmMAPKKK19, ZmMAPKKK21, ZmMAPKKK22 and ZmMAPKKK56 are co-expressed with some DEGs in most pathways.
In conclusion, seventy-one MAPKKK genes in maize are identified in our study. Among them, 8 differentially expressed MAPKKK genes responsive to drought stress are obtained. The results suggest that MAPKKKs may react to drought stress through various pathways. To our knowledge, this is the first study of the expression of maize MAPKKK genes in relation to drought stress. The results of our study should provide further information about the characteristics of MAPKKKs in maize and their roles in the response to drought stress. Further biological functional analyses of MAPKKKs and their downstream MAPKK and MAPK targets will help elucidate the mechanisms underlying the response of MAPK cascades to drought in maize.
Supporting Information S1 File. Supplementary tables. S1 File contains five Supplementary tables: Table A. 71 predictive MAPKKKs in maize. Table B. The quality of the RNA-seq data for Q20 and Q30. Table C. RNA-seq data and mapping rates.  Green represents root. The enrichment figure was constructed of significantly enriched GO terms only. The X-axis is the -log 10 (p-value), which represents the level of enrichment, and the p-value cutoff was 0.05. Fig C. Venn analysis for the DEGs with up-regulated and downregulated in different organs. Each circle represents an organ with different regulated pattern. Fig D. Differentially expressed transcripts involved in different metabolic processes under drought stress. (a, b and c) represent drought-mediated expression changes in different metabolic processes in the leaf, stem and root meristem, respectively. The images were obtained using MapMan and show the different functional categories that fulfilled the criteria for differential expression (a p-value less than 0.05 and greater than 2-fold change). Fig E. Differentially expressed transcripts involved in different regulatory processes under drought stress. (a, b and c) represent drought-mediated expression changes in different regulatory processes in the leaf, stem and root meristem, respectively. The images were obtained using MapMan and show the different functional categories that fulfilled the criteria for differential expression (a p-value less than 0.05 and greater than 2-fold change). (ZIP)