Genome-wide identification, expression profiles and regulatory network of MAPK cascade gene family in barley
BMC Genomics volume 20, Article number: 750 (2019)
Mitogen-activated protein kinase (MAPK) cascade is a conserved and universal signal transduction module in organisms. Although it has been well characterized in many plants, no systematic analysis has been conducted in barley.
Here, we identified 20 MAPKs, 6 MAPKKs and 156 MAPKKKs in barley through a genome-wide search against the updated reference genome. Then, phylogenetic relationship, gene structure and conserved protein motifs organization of them were systematically analyzed and results supported the predictions. Gene duplication analysis revealed that segmental and tandem duplication events contributed to the expansion of barley MAPK cascade genes and the duplicated gene pairs were found to undergone strong purifying selection. Expression profiles of them were further investigated in different organs and under diverse abiotic stresses using the available 173 RNA-seq datasets, and then the tissue-specific and stress-responsive candidates were found. Finally, co-expression regulatory network of MAPK cascade genes was constructed by WGCNA tool, resulting in a complicated network composed of a total of 72 branches containing 46 HvMAPK cascade genes and 46 miRNAs.
This study provides the targets for further functional study and also contribute to better understand the MAPK cascade regulatory network in barley and beyond.
To coordinate the biotic and abiotic stresses during growth and development, plants have evolved to form the complex mechanisms to perceive and transmit environmental stimuli by inducing or repressing a series of genes to express . The Mitogen-activated protein kinase (MAPK) cascades are characterized as evolutionarily conserved and fundamentally universal signaling transduction pathways, playing the vital roles as diverse receptors/sensors from the extracellular environment to intracellular transcriptional and metabolic centers in eukaryotes . The canonical MAPK cascade is composed of three specific kinases, namely MAPK, MAPK kinase (MAPKK) as well as MAPKK kinase (MAPKKK), which was activated sequentially by phosphorylation at certain activation sites [3, 4]. In general, MAPKs are phosphorylated at their conserved threonine and tyrosine residues in the activation loop (T-loop) by MAPK kinase, and in turn, MAPKK are activated by MAPKKKs when development or environmental signals incurred as their serine and serine/threonine residues located in the S/T activation site are phosphorylated [1, 2].
In plants, extensive studies have revealed that the MAPK cascades widely involved in regulating many biological processes, including cell division, plant development, growth and hormonal response as well as in response to diverse biotic and abiotic stresses, such as drought, salt, heat and pathogen infection [5,6,7]. In light of their importance, a large number of MAPK genes have been functionally identified in several plants, including Arabidopsis , rice [9,10,11], Brachypdoium [12, 13] and maize [14, 15]. At the same time, a series of plant MAPK signaling cascades have also been well constructed and studied. The AtMEKK1-MKK4/5-MPK3/6 cascades is the first identified MAPK signaling module in plants, which was involved in plant innate immunity of flg22 signal transmission [16, 17]. The complete MAPK signaling cascade of ANP3-MKK6-MPK4 and YDA-MKK4/5-MPK3/6 is determined to control the stomatal development and patterning in Arabidopsis . MEKK1-MKK1/2-MPK4 module was found to play the important role in the defenses against abiotic stresses and contributed to the freezing tolerance in Arabidopsis [19,20,21]. The ABA(abscisic acid)-activated MEKK17/18-MKK3-MPK1/2/7/14 module displayed stress signaling to ABA and regulated the expression of a series of ABA-dependent genes . In tobacco, the NPK1-MEK1-Ntf6 cascade was identified to confer the resistance to tobacco mosaic virus via mediating the resistant protein N . Additionally, the NPK1-NQK1/NtMEK1-NRK1 module is found to be a positive regulator of tobacco cytokinesis during meiosis as well as mitosis . Barley (Hordeum vulgare L.) is one of the earliest domesticated and also one of the most important staple crops, which holds the significance for agriculture drawn and human civilization [25, 26]. Furthermore, barley is also well-studied in terms of cytology, genetics and genomics and thus qualifies as the model for Triticeae research . The survey of MAPK family in barley has also been conducted and a total of 16 HvMAPKs were identified based on the full-length cDNA, EST(expressed sequence tag) and genomic survey database . However, the incomplete data used by Krenek et al.  might cause the incomplete prediction and identification of MAPKK and MAPKKK family is not performed in barley up to now. The recently published reference-quality barley genome  makes it possible to conduct a comprehensive identification of its MAPK cascade gene families at whole genome scale and then construct the MAPK signal transduction pathway.
In this study, we systematically identified the MAPK, MAPKK and MAPKKK gene family based on a genome-wide search against barley reference genome. Then, the gene structures, chromosomal locations, gene duplication events and evolutionary dynamics were investigated. Furthermore, the expression patterns at diverse development stages and under different abiotic stresses were also analyzed. Finally, we constructed the regulatory networks of MAPK-MAPKK-MAPKKK signal pathway based on the co-expression patterns from a total of 173 RNA-seq datasets. This study reported the genomic organization, expression and phylogenetic relationships of the MAPK, MAPKK and MAPKKK gene families in barley, which could provide the candidates for further functional analysis and also contribute to illuminate the MAPK signal cascade-mediated pathway of barley and beyond.
Results and discussion
Genome-wide identification of MAPK cascade genes in barley
Availability of the reference-quality barley genome  made it possible for the first time to systematically identify all the MAPK cascade genes in this model crop species. Using the methods as described below, a total of 20 HvMAPK, 6 HvMAPKKs and 156 HvMAPKKKs were obtained, respectively (Table 1). The conserved domain analysis showed that all of them have the serine/threonine-protein kinase-like domain (PFAM accession No. PF00069) (Additional file 7: Table S1). We further validated the identified genes using the public ESTs to provide the expression support. Results showed that majority (19 out of 20 HvMAPKs, 5 out of 6 HvMAPKKs and 103 out of 156 HvMAPKKKs) of the predicted genes had the existing EST hit supports (Table 1). Given the limit of available ESTs, the non-supported HvMAPK cascade gene might not be detected under specific conditions or low levels of expression that can’t be investigated experimentally. Compared to previous study that only 16 HvMAPKs were identified by Krenek et al , this study found 20 HvMAPKs, which covered the 16 previous predicted ones, suggesting the whole genome-search could provide more comprehensive prediction of barley MAPK family.
Furthermore, the physical and chemical properties of these genes were investigated and compared. The length of MAPK cascade related proteins varied from 100 to 1332 amino acids, with an average of 596 in length. The putative molecular mass ranged from 11.2 kDa to 147.1 kDa, and the isoelectric points varied from 4.22 to 9.73, respectively (Table 1), which is similar to that of wheat and Brachypodium [29, 30]. The significance difference of physical and chemistry properties between the members of barley MAPK genes suggested that the subfunctionalization and neofunctionalization may have occurred among the MAPK cascade genes in barley . Analysis of subcellular location showed that 52 (30%) HvMAPK cascade genes were predicted to be located in nuclear, followed by PlasmaMembrane (45) and Cytoplasmic (43), while the remaining ones were predicted to be located in chloroplast, mitochondrial and extra-cellular.
These 182 barley MAPK cascade genes can be classified into three major clades in coordination to MAPK, MAPKK and MAPKKK with the specific conserved signature motifs, respectively (Fig. 1). Among them, 20 genes harboring the specific conserved signature motifs of T(E/D)YVxTRWYRAPE(L/V), and 6 genes possessing the VGTxxYMSPER conserved signature, which were categorized into MAPK and MAPKK subfamilies, respectively [3, 31]. Consistent with the other species [3, 10], these HvMAPKs could be assigned into the 10 TDY- and 10 TEY-subtype members (Fig. 2a and Additional file 1: Figure S1). We also investigated the docking site CD (Common docking) domain in HvMAPKs. Results showed that the TDY-subtype HvMAPKs lacked this domain (Fig. 2c and Additional file 2: Figure S2), which was the same as that of Arabidopsis . All of MAPKK members contained the VGTxxYMSPER motif and the putative MAPK docking sites [K/R][K/R][K/R]x(1–5)[L/I]x[L/I] (Additional file 3: Figure S3). The remaining 156 genes belonged to MAPKKK subfamily. The barley MAPKKK genes could be further divided into three groups, which owned the conserved motifs of G(T/S)Px(W/Y/F)MAPEV, GTxx(W/Y)MAPE and GTPEFMAPE(L/V)Y for MEKK, Raf-like and ZIK subfamilies, respectively (Additional file 4: Figure S4). Remarkably, the Raf-like subfamily had 124 members, ranking the largest group of MAPKKK in barley, whereas the ZIK subfamilies possessed only 4 members as the smallest group, which was consistent with the abundance and composition of MAPKKK genes in other species, especially in wheat [29, 30] (Table 2).
Phylogenetic relationship, gene structure and motifs analysis
To further support the subfamily grouping, phylogenetic analysis were performed using the full-length protein sequences of these barley MAPK cascade genes (Fig. 3). Consistent with specific conserved signature motifs , the MEKK, Raf-like and ZIK subfamilies belonging to MAPKKK family were also clustered into independent sub-clade, respectively. For MAPK, it could be further divided into TDY and TEY two sub-clades, and TEY sub-clade was further assigned into A to C subgroups. We further performed phylogenetic analysis of these HvMAPK and the reported rice and Arabidopsis MAPKs. Results found they could clustered into different groups and the orthology pairs of them were obtained depending on phylogenetic relationship (Additional file 5: Figure S5). These results could provide some clues for candidate selection for further functional study as some orthologous genes in rice and Arabidopsis has been extensively functionally characterized [16, 18].
Gene structure played vital roles in the evolution of gene families and provided extra evidence to estimate the functional diversifications . Thus, the exon-intron organization of these barley MAPK cascade genes was further analyzed (Fig. 2b). Result found that there were significant intron abundance variations between these genes. It is reported that C- and D-group of MAPKKs tend to have no introns in Arabidopsis . The C-group of HvMAPKKs also showed intron-less while D-group have abundant introns. For instance, HvMAPKK3 and HvMAPKK4, which assigned into D subgroup, possessed 7 and 9 introns, respectively. Furthermore, the intron count of HvMAPKKK gene family ranged from 1 to 24, showing obviously variations even in the same subgroup. For the MEKK subfamily, more than half (54.2%) of the genes possessed no or one intron, while the other MEKK members had 6 to 24 introns. The intron number of the ZIK subfamily varied from 2 to 5, whereas the RAF genes with the intron number ranged from 1 to 20 and presented the highest level of variation among them.
Additionally, the conserved protein domains in the barley MAPK cascade genes were identified and compared. A total of 32 conserved motifs were detected (Fig. 2c). The protein kinase domain was found in each member of the MAPK cascade proteins. A certain degree of conservation could be observed in the HvMAPK and HvMAPKK genes that almost all of them harbored the ATP (Adenosine triphosphate) binding site and serine/threonine-protein kinase active site. Similar to the intron/exon structure, the composition of conserved motifs was also highly variable in HvMAPKKK family. Apart from the protein kinase and its related domains, a series of other functional motifs was widely distributed, such as Bulb-type lectin domain, S-locus glycoprotein domain and PAN/Apple domain, suggested they are widely involved in growth and development as well as signaling transduction . The PAS domain, S-locus glycoprotein domain and Concanavalin A-like lectin/glucanase domain were possessed by 4, 1 and 3 Raf subfamily members. The EF-hand domain pair, EF-Hand 1, calcium-binding site and EF-hand domain were uniquely found in MEKK subfamily, whereas no domains were specific to the ZIK subfamily. On the whole, the MAPK cascade proteins clustered into the same group phylogenetically tended to share similar motifs composition.
Finally, the 1.5 kb genomic sequences upstream of the transcriptional start sites of HvMAPK genes were extracted and used to identify the cis-regulatory elements. Totally, 27 cis-elements were obtained, of which SARE(salicylic acid responsiveness) domain and the TGA(auxin-responsive) domain were found to be present only in 3 and 7 genes respectively, whereas the Skn-1 motif was shared by 159 genes, which ranked the least and most abundant motifs (Additional file 7: Table S2). Skn-1 motif is reported to be a cis-acting regulatory element required for endosperm expression and oxidative stress response in eukaryotes , suggesting the MAPK cascades played the important role in regulating the barley development and stress response. In addition, a large amount of plant growth and development (including circadian, meristem and endosperm), hormone-related (e.g., abscisic acid, auxin, MeJA, ethylene, gibberellin) cis-elements were found in these promoter regions, suggesting that MAPK cascade genes widely involved in regulating the signal transduction network of diverse developmental processes. Meanwhile, the cis-element related to biotic (e.g. fungal and wound) and abiotic stress response (e.g. salt, extreme temperature, dehydration) were also identified in the promoter region of the HvMAPK cascade genes, which suggested that these MAPK cascade genes might have potential functions in stress adaptation and signaling pathways .
Gene duplication and synteny analysis
In order to investigate the mechanism of expansion of the MAPK cascade genes in barley, we further investigated the segmental and tandem duplication events by genome synteny analysis. Results showed that 13 paralogs composed of 26 HvMAPK cascade genes were identified, of which 5 were segmental duplications and 8 were tandem duplication events (Fig. 4 and Additional file 7: Table S3). In detail, 3 and 2 segmental events were found in HvMAPKs and HvMAPKKKs, as well as 8 tandem repeats events in HvMAPKKKs, suggesting that segmental duplication played important roles in the expansion of MAPKs while tandem repeat duplication was the driven force for HvMAPKKK gene family expansion. It is noteworthy that the segmental events mainly occurred at chromosome 1 and chromosome 3, whereas the tandem duplication blocks distributed throughout the whole genome, of which 1, 1, 4, 1, 1 paralogous pairs were mapped to chromosome 1, 2, 3, 4 and 5, respectively (Fig. 4). In order to detect the selection effect during gene divergence after duplication, the Ka/Ks substitution ratio of the duplicated pairs were further calculated. Result showed that Ka/Ks ratios of MAPK cascade genes ranged from 0.001 to 0.4727, with an average of 0.1964, suggesting that they have undergone purifying selection pressure during the process of evolution in barley .
Furthermore, the comparative analysis between barley with other six species (Brachypodium, sorghum, maize, rice, soybean and grape) was performed to determine the origin and evolutionary relationships of MAPK cascade genes (Fig. 5). Through whole genome-wide syntenic analysis, a total of 84, 80, 77, 67, 5 and 7 barely MAPK cascade genes were identified to have orthologous counterpart in Brachypodium, rice, sorghum, maize, grape and soybean (Additional file 7: Table S4 to S9). The average Ka/Ks value was maximum between barley and Brachypodium (0.1641), followed by rice and sorghum (0.1544) as well as maize (0.43), suggesting the genes pairs between barley and those species appeared to have undergone extensive intense purifying selection. Besides, we found that most of MAPK cascade genes showed syntenic bias towards particular chromosomes of sorghum, maize, rice, which indicated that the chromosomal rearrangement events like duplication and inversion may predominantly shape the distribution and organization of MAPK genes in these genomes .
Comprehensive analysis of the expression profiles of barley MAPK cascade genes
To preliminarily predict the biological function of these barley MAPK cascade genes, gene ontology (GO) analysis was firstly performed (Additional file 6: Figure S6) and they could be annotated into 40 GO terms, including 9 terms of molecular function, 19 of biological processes and 11 of cellular components, respectively. In the cellular components category, cell and cell part were main annotation terms, whereas binding, catalytic nucleoside and transferase were the most presented function in the molecular function category. In the biological process category, cellular metabolic, cellular, metabolic and macromolecule metabolic process occupied most of the proportion. By employing the fisher statistical test method, a total of 17 terms were significant enriched (P < 0.05 and Q < 0.05) when taking the whole barley genome as customized backgrounds, including 5 biological process categories, 6 molecular function categories and 6 cellular component categories (Additional file 7: Table S10). These results revealed that the MAPK cascade genes played diverse roles in diverse development and stress response pathways in barley.
Furthermore, the expression profiles of MAPK cascade genes at 16 developmental stages were investigated using RNA-Seq data. A total of 75 genes were found to be expressed in at least one organ or stage (Fig. 6). A high variance in the expression levels among these MPAK cascade genes was observed, of which a series of them showed relatively high expression in all the tested tissues, such as HvMAPK1, HvMAPK4, HvRaf-like63, HvRaf-like87 and HvZIK2, The ortholog of HvZIK2 in Arabidopsis is AtZIK4(WNK1), which is found to regulating internal circadian rhythm and flowering time . It highly expressed in different organs, suggesting it also played the indispensable role in organ formation and development. Additionally, the tissue- and stage-specific MAPK cascade genes were also identified. HvRaf-like103 and HvRaf-like49 were found to be predominantly expressed in senescing leaf, whereas HvRaf-like66, HvRaf-like47, HvRaf-like93 and HvMAPK7 showed preferential expression in the root, lemma, seedling root and epidermis, respectively, suggesting that these genes may mainly involve into organ- or tissue-specific development in barley.
To get insight into the roles of MAPK cascade genes in response to abiotic stresses, the expression profiles of them under drought, heat, salt were investigated to discover the abiotic stress-responsive candidates. Results showed that a total of 123 genes were detected to be expressed under drought stress (Fig. 7a). Among them, 10 and 24 genes were significantly up-regulated, whereas 5 and 19 MAPK cascade genes were significantly down-regulated in flowers and leaves when subjecting to drought. Meanwhile, 114 MAPK cascade genes were found to express under heat stress (Fig. 7b). Remarkably, HvRaf-like124 and HvMAPKK5 presented about 62 and 21 times higher expression level under heat stress compared to control. Previous study found the MPK20 have the defense function in cotton, while its ortholog HvMAPKK5 involved in regulating heat stress adaptation in barley, suggesting it might have divergent function in different species . The expression patterns of MAPK cascades genes under salt stress were also examined (Fig. 7c). Totally, 5, 7 and 9 genes showed up-regulated in the root Z1, Z2 and Z3 respectively, of which the expression level of HvRaf-like28 and Hv-Raf-like113 were up-regulated with more than 10 fold at the Z1 zone and HvMAPKK1 showed 34-fold change at the Z2 zone. Besides, a total of 7, 11 and 4 genes were identified to be down-regulated at root Z1, Z2 and Z3 zone respectively. HvZIK4 and HvRaf-like56 was 862 and 558 time lower expression at Z1 and Z2 zone of root under salt stress than that of control.
Finally, the expression profiles of these genes under zinc metal poisoning and iron were investigated (Fig. 7d). When in response to iron stress, 9 genes showed up-regulated and 7 showed down-regulated after 6 h treatment. Furthermore, 8 up-regulated and 11 down-regulated genes were found after 24 h treatment. Among them, HvMAPK17, HvRaf-like4, HvRaf-like70, HvRaf-like109 HvZIK3and HvRafZIK4 all presented up-regulated under iron stress after both 6 h and 24 h treatment, whereas HvMAPK2 and HvRaf-like41 showed down-regulated. Under zinc stress, a total of 13 and 12 up-regulated genes as well as 14 and 16 down-regulated genes were found after 6 h and 24 h treatment, respectively. Among them HvMAPKK5, HvMEKK7, HvMEKK26, HvRaf-like28 and HvRaf-like58 were all down-regulated at all treatment, whereas HvZIK3, HvRaf-like65, HvRaf-like4, HvRaf-like108, HvMEKK14, HvMEKK10 and HvRaf-like108 displayed up-regulated after both 6 h and 24 h treatment. Obviously, HvZIK3, HvRaf-like4, HvRaf-108 showed up-regulated expression under both iron and zinc treatment, which might play the important roles in regulating signal transduction process for metal poisoning response and detoxification.
Network construction of HvMAPK cascade genes
To get the network of miRNA targeting on MAPK cascade genes, the putative miRNAs targeted HvMAPK cascade genes were analyzed. Results found that 26 MAPK cascade genes including 3 MAPKs and 23 MAPKKK genes were predicted to be targeted by 11 miRNAs, while no miRNA target was found for HvMAPKK genes, which might be due to the limited barley miRNA reported at present (Additional file 7: Table S11). Totally, 36 miRNA-MAPK interactions were constructed based on the target relationship. The barley cascade genes were mainly inhibited by miRNAs through transcript cleavage (94.44%), while HvRaf-like12 and HvRaf-like12 and HvRaf-like76 were inhibited to translation by miRNAs. Additionally, miRNAs mainly targeted on the CDS region but behind the protein kinase domain of these MAPK cascade genes to function gene silence.
The co-expression regulatory network was further constructed to detect the interaction among these barley MAPK cascade genes based on weighted correlation of their expressions using a big datasets of 173 RNA-seq data. Only the relations between MAPKKK and MAPKK as well as MAPKK and MAPK were presented. A total of 40 interactions composed of 25 genes were constructed, including 7 MAPK, 3 MAPKK and 15 MAPKKK genes respectively (Fig. 8). Among them, some MAPK cascade modules has been verified in model plants, such as MKK3-MPK6 in Arabidopsis and MAPK18-MAPKK2-MEKK4 in Brachypodium . Furthermore, a total of 18 genes including 2 MAPK, 10 MEKK, 2 HvRaf-like and one ZIK gene were predicted to be interacted with HvMAPKK3, suggesting that it may be the hub gene of the co-expression regulatory network, playing the key role in barley MAPK cascade signaling pathway. In Arabidopsis, MAPKK3 is found to be expressed in all organs, and plays a vital role in photomorphogenesis to regulate gene expression under various light conditions, as well as involved in cell expansion, pathogen signaling and jasmonate signaling pathway, indicating it is critical for development and signaling transduction [39, 40]. Thus, the barley ortholog HvMAPKK3 might also play the hub role in co-expression network in barley response to development and stresses. In addition, there was 10 MAPK-MAPKK, 30 MAPKK-MAPKKK interactions were also obtained to use to subsequently experimental validation. Combined with miRNA-target interaction mentioned above, the regulatory network containing a total of 46 HvMAPK cascade genes and 46 miRNAs were constructed and 72 branches were linked for each other, which provided the indispensable resource to facilitate the MAPK pathway and signal transduction mechanism studies in barley and beyond.
This is the first study to identify the MAPK cascade genes in barley at genomic level. Totally, 20 HvMAPKs, 6 HvMAPKKs and 156 HvMAPKKKs were obtained, which was further supported the existence by EST or full-length cDNA sequences. The phylogenetic relationships, intron-exon structure as well as conserved motif analysis all strongly supported the prediction. Furthermore, both segmental and tandem duplication events contributed to the expansion of the MAPK cascade genes in barley. The expression profiles of these MAPK cascade genes during development and under abiotic stresses were investigated and the tissue-specific or stress-responsive genes were identified, which could be considered as the candidates for further functional studies. Finally, the co-expression regulatory network of the MAPK cascade genes was constructed using WGCNA tool based on a total of 174 RNA-seq data. A total of 30 MAPKKK-MAPKK, 10 MAPKK-MAPK potential interactions were identified, which contributed to better understanding the MAPK signal transduction pathway in barely.
Identification of MAPK cascade genes in barley
The protein sequences of the latest updated barley genome Morex v2.0  were retrieved from the IPK website (http://webblast.ipk-gatersleben.de/barley_ibsc/). Then, the MAPK cascade proteins of Arabidopsis from the TAIR database, were used as queries to search against the barley proteins using BLASTP program with an e-value of 1e-5 and identity of 50% as the threshold. The HMMER 3.0 program was employed to conduct for Hidden Markov Model (HMM) algorithm search using the serine/threonine-protein kinase-like domain (PF00069) as the query with the threshold of E < 1e− 5. The HMMER hits were further integrated with the BLASTP results and parsed by manual editing to remove redundant. Those genes displayed the consensus sequences as Jonak et al described were considered as the potential MAPK cascade genes . The candidates were subsequently submitted to SMART and PFAM web tool to verify the kinase domain. Additionally, the putative MAPK cascade genes were further verified through searching against the barely ESTs by BLASTN tool. The theoretical isoelectric point (pI), molecular weight (MW) and gravy of the identified barley MAPK cascade genes were evaluated using ProtParam tool (http://web.expasy.org/protparam/) integrated in ExPASy database. The cello online server (http://cello.life.nctu.edu.tw/) was used to detect the subcellular localization and protein solubility was predicted by PROSOII tool (http://mips.helmholtz-muenchen.de/prosoII).
Phylogenetic relationship and conserved motif analysis
Multiple sequence alignment were performed using ClustalX v2.0 with default parameter . A neighbor–joining (NJ) phylogenetic tree was constructed based on the full-length protein sequences using the MEGA software with a bootstrap of 1000 replications . The gene structures were obtained from the GTF annotation file of barley genome and then were displayed by Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/index.php). Furthermore, the protein domain and conserved motifs of barley MAPK cascade genes were predicted using InterProScan tool. Finally, the upstream 1.5 kb genomic DNA sequences of each gene were extracted from barley genome, and then submitted to PlantCARE database to detect the putative cis-regulatory elements .
Gene duplication and molecular selection analysis
Gene duplication events were defined based on the following three criteria: 1) the alignment should cover more than 70% of the longer gene; (b) the identity of the aligned region should be more than 70%; 3) for tightly linked genes only one duplication event was counted . The gene synteny between barley and other species, including Brachypodium distachyon, Sorghum bicolor, Zea mays, Oryza sativa, Vitis vinifera and Glycine max was conducted using the MCScanX toolkit . The linked genes pairs were displayed using the Circus tool. The rate of Ka (non-synonymous substitution)/Ks (synonymous substitution) was employed to assess the codon evolutionary rate between the synteny genes using the codeml program embedded in the PAML package . The formula T = Ks/2λ was employed to calculate the duplication and divergence time, where λ referred to the mutation rate, was considered as 6.5 × 10− 9 synonymous substitutions per site per year.
Expression profiles and co-expression networks construction
The MAPK cascade genes were firstly searched against the NR protein database using the local BLASTx with an E-value cut off of 10–5. Based on the Nr annotation, Blast2GO  program was used to retrieved the GO (gene ontology) annotation. AgriGO v2 (http://systemsbiology.cau.edu.cn/agriGOv2/index.php) was applied to conduct the singular enrichment analysis. Furthermore, a total of 172 public available RNA-seqs (Additional file 7: Table S12) including multiple tissues and developmental stages as well as biotic and abiotic stresses were downloaded from the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) database to investigate the expression profiles of these genes. The FPKM (fragments per kilobase of transcript per million fragments mapped reads) value were calculated by Hisat2 and Stringtie software . Then, differentially expressed genes were identified with the following threshold values: fold change≥2, FDR(false discovery rate) ≤ 0.01, and the absolute ratio of log2 ≥ 1. All FPKM data was finally reported by log2 counts and the heat map was visualized using pheatmap package in R. WGCNA was used to construct the co-expression network based on all of the downloaded transcriptome data . Besides, all the identified MAPK cascade transcripts were submitted to the psRNATarget tool  to search the barley miRNAs targets in the miRBase. The regulatory network of Hvu-miRNA and HvMAPK cascade genes were visualized using cytoscape tool (http://www.cytoscape.org/).
Availability of data and materials
The datasets supporting the conclusions of this article are included within the article and its additional files.
Expressed sequence tag
False discovery rate
Fragments per kilobase of transcript per million fragments mapped reads
Hidden Markov Model
Mitogen-activated protein kinase
Mitogen-activated protein kinase kinase
Mitogen-activated protein kinase kinase
Salicylic acid responsiveness
Putative activation motif in MAPK gene TDY1
Tena G, Asai T, Chiu WL, Sheen J. Plant mitogen-activated protein kinase signaling cascades. Curr Opin Plant Biol. 2001;4(5):392–400.
MAPK Group (Ichimura K, et al.) (2002) Mitogen-activated protein kinase cascades in plants: a new nomenclature. Trends Plant Sci, 7(7), 301–308.
Jonak C, Okrész L, Bögre L, Hirt H. Complexity, cross talk and integration of plant MAP kinase signalling. Curr Opin Plant Biol. 2002;5(5):415.
Zhang T, Liu Y, Yang T, Zhang L, Xu S, Xue L, An L. Diverse signals converge at MAPK cascades in plant. Plant Physiol Biochem. 2006;44(5):274–83.
Hamel LP, Nicole MC, Sritubtim S, Morency MJ, Ellis M, et al. Ancient signals: comparative genomics of plant MAPK and MAPKK gene families. Trends Plant Sci. 2006;11:192–8.
Pitzschke A, Schikora A, Hirt H. MAPK cascade signalling networks in plant defence. Curr Opin Plant Biol. 2009;12(4):421–6.
Hao L, Wen Y, Zhao Y, Lu W, Xiao K. Wheat mitogen-activated protein kinase gene TaMPK4 improves plant tolerance to multiple stresses through modifying root growth, ROS metabolism, and nutrient acquisitions. Plant Cell Rep. 2015;34:2081–97.
Menzel W, Stenzel I, Helbig LM, Krishnamoorthy P, Neumann S, Eschen-Lippold L, Heilmann M, Lee J, Heilmann I. A PAMP-triggered MAPK cascade inhibits phosphatidylinositol 4,5-bisphosphate production by PIP5K6 in Arabidopsis thaliana. New Phytol. 2019;224(2):833-47.
Rohila JS, Yang Y. Rice mitogen-activated protein kinase gene family and its role in biotic and abiotic stress response. J Integr Plant Biol. 2007;49(6):751–9.
Rao KP, Richa T, Kumar K, Raghuram B, Sinha AK. In silico analysis reveals 75 members of mitogen-activated protein kinase kinase kinase gene family in rice. DNA Res. 2010;17(3):139–53.
Wankhede DP, Misra M, Singh P, Sinha AK. Rice mitogen activated protein kinase kinase and mitogen activated protein kinase interaction network revealed by in-silico docking and yeast two-hybrid approaches. PLoS One. 2013;8(5):e65011.
Chen L, Hu W, Tan S, Wang M, Ma Z, Zhou S, Deng X, Zhang Y, Huang C, Yang G. Genome-wide identification and analysis of MAPK and MAPKK gene families in Brachypodium distachyon. PLoS One. 2012;7(10):e46744.
Feng K, Liu F, Zou J, Xing G, Deng P, Song W, Tong W, Nie X. Genome-wide identification, evolution and co-expression network analysis of mitogen-activated protein kinase kinase kinases in Brachypodium distachyon. Front Plant Sci. 2016;7(556):1400.
Liu Y, Zhang D, Wang L, Li D. Genome-wide analysis of mitogen-activated protein kinase gene family in maize. Plant Mol Biol Report. 2013;31(6):1446–60.
Kong X, Pan J, Zhang D, Jiang S, Cai G, Wang L, Li D. Identification of mitogen-activated protein kinase kinase gene family and MKK–MAPK interaction network in maize. Biochem Biophys Res Commun. 2013;441(4):964–9.
Asai T, Tena G, Plotnikova J, Willmann MR, Chiu WL, Gomezgomez L, Boller T, Ausubel FM, Sheen J. MAP kinase signaling cascade in Arabidopsis innate immunity. Nature. 2002;415(6875):977.
Galletti R, Ferrari S, De LG. Arabidopsis MPK3 and MPK6 play different roles in basal and oligogalacturonide- or flagellin-induced resistance against Botrytis cinerea. Plant Physiol. 2011;157(2):804–14.
Eckardt NA. A complete MAPK signaling cascade that functions in stomatal development and patterning in Arabidopsis. Plant Cell. 2007;19(1):7.
Petersen M, Brodersen P, Naested H, Andreasson E, Lindhart U, Johansen B, Nielsen HB, Lacy M, Austin MJ, Parker JE. Arabidopsis MAP kinase 4 negatively regulates systemic acquired resistance. Cell. 2000;103(7):1111–20.
Qiu JL, Zhou L, Yun BW, Nielsen HB, Fiil BK, Petersen K, Mackinlay J, Loake GJ, Mundy J, Morris PC. Arabidopsis mtogen-activated protein kinase kinases MKK1 and MKK2 have overlapping functions in defense signaling mediated by MEKK1, MPK4, and MKS1. Plant Physiol. 2008;148(1):212–22.
Furuya T, Matsuoka D, Nanmori T. Membrane rigidification functions upstream of the MEKK1-MKK2-MPK4 cascade during cold acclimation in Arabidopsis thaliana. FEBS Lett. 2014;588(11):2025–30.
Danquah A, De ZA, Boudsocq M, Neubauer J, Frei DFN, Leonhardt N, Pateyron S, Gwinner F, Tamby JP, Ortizmasia D. Identification and characterization of an ABA-activated MAP kinase cascade in Arabidopsis thaliana. Plant J Cell Mol Biol. 2015;82(2):232–44.
Liu Y, Schiff M, Dineshkumar SP. Involvement of MEK1 MAPKK, NTF6 MAPK, WRKY/MYB transcription factors, COI1 and CTR1 in N-mediated resistance to tobacco mosaic virus. Plant J Cell Mol Biol. 2004;38(5):800–9.
Soyano T, Nishihama R, Morikiyo K, Ishikawa M, Machida Y. NQK1/NtMEK1 is a MAPKK that acts in the NPK1 MAPKKK-mediated MAPK cascade and is required for plant cytokinesis. Genes Dev. 2003;17(8):1055.
Mayer KFX, Waugh R, Langridge P, Close TJ, Wise RP, Graner A, Matsumoto T, Sato K, Schulman A, Muehlbauer GJ. A physical, genetic and functional sequence assembly of the barley genome. Nature. 2012;491(7426):711–6.
Mascher M, Gundlach H, Himmelbach A, Beier S, Twardziok SO, Wicker T, Radchuk V, Dockter C, Hedley PE, Russell J. A chromosome conformation capture ordered sequence of the barley genome. Nature. 2017;544(7651):427–33.
Lv SZ, Nie XJ, Wang L, Biradar SB, Jia XO, Weining S. Identification and characterization of microRNAs from barley (Hordeum vulgare L.) by Solexa sequencing. Int J Mol Sci. 2012;13:2973–84.
Křenek P, Niks RE, Vels A, Vyplelová P, Šamaj J. Genome-wide analysis of the barley MAPK gene family and its expression patterns in relation to Pucciniahordei infection. Acta Physiol Plant. 2015;37(11):254.
Wang M, Yue H, Feng K, Deng P, Song W, Nie X. Genome-wide identification, phylogeny and expressional profiles of mitogen activated protein kinase kinase kinase (MAPKKK) gene family in bread wheat (Triticum aestivum L.). BMC Genomics. 2016;17:668.
Jiang M, Wen F, Cao J, Li P, She J, Chu Z. Genome-wide exploration of the molecular evolution and regulatory network of mitogen-activated protein kinase cascades upon multiple stresses in Brachypodium distachyon. BMC Genomics. 2015;16:228.
Lehti-Shiu MD, Shiu SH. Diversity, classification and function of the plant protein kinase superfamily. Philos Trans R Soc Lond Ser B Biol Sci. 2012;367(1602):2619–39.
Ye J, Yang H, Shi H, Wei Y, Tie W, Ding Z, Yan Y, Luo Y, Xia Z, Wang W, Peng M, Li K, Zhang H, Hu W. The MAPKKK gene family in cassava: genome-wide identification and expression analysis against drought stress. Sci Rep. 2017;7(1):14939.
Colcombet J, Hirt H. Arabidopsis MAPKs: a complex signalling network involved in multiple biological processes. Biochem J. 2008;413(2):217.
Wilson MA, Iser WB, Son TG, Logie A, Cabral-Costa JV, Mattson MP, Camandola S. Skn-1 is required for interneuron sensory integration and foraging behavior in Caenorhabditis elegans. PLoS ONE. 2017;12(5):e0176798.
Akhunov ED, Sunish S, Hanquan L, Shichen W, Akhunova AR, Gaganpreet K, Wanlong L, Forrest KL, Deven S, Hana S. Comparative analysis of syntenic genes in grass genomes reveals accelerated rates of gene structure and coding sequence evolution in polyploid wheat. Plant Physiol. 2013;161(1):9090–100.
Wang Y, Liu K, Liao H, Zhuang C, Ma H, Yan X. The plant WNK gene family and regulation of flowering time in Arabidopsis. Plant Biol (Stuttg). 2008;10(5):548–62.
Wang C, He X, Li Y, Wang L, Guo X, Guo X. The cotton MAPK kinase GhMPK20 negatively regulates resistance to Fusarium oxysporum by mediating the MKK4-MPK20-WRKY40 cascade. Mol Plant Pathol. 2018;19(7):1624–38.
Sethi V, Raghuram B, Sinha AK, Chattopadhyay S. A mitogen-activated protein kinase cascade module, MKK3-MPK6 and MYC2, is involved in blue light-mediated seedling development in Arabidopsis. Plant Cell. 2014;26(8):3343–57.
Schikora A, Carreri A, Charpentier E, Hirt H. The dark side of the salad: Salmonella typhimurium overcomes the innate immune response of Arabidopsis thaliana and shows an endopathogenic lifestyle. PLoS One. 2008;3(5):e2279.
Dóczi R, Brader G, Pettkó-Szandtner A, Rajh I, Djamei A, Pitzschke A, Teige M, Hirt H. The Arabidopsis mitogen-activated protein kinase kinase MKK3 is upstream of group C mitogen-activated protein kinases and participates in pathogen signaling. Plant Cell. 2007;19(10):3266–79.
Larkin MA, Blackshields G, Brown N, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R. ClustalW and Clustal X version 2.0. Bioinformatics. 2007;23(21):2947–8.
Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9.
Lescot M, Déhais P, Thijs G, Marchal K, Moreau Y, Van de Peer Y, Rouzé P, Rombauts S. PlantCARE, a database of plant cis-acting regulatory elements and a portal to tools for in silico analysis of promoter sequences. Nucleic Acids Res. 2002;30(1):325–7.
Gu Z, Cavalcanti A, Chen FC, Bouman P, Li WH. Extent of gene duplication in the genomes of Drosophila, nematode, and yeast. Mol Biol Evol. 2002;19(3):256–62.
Wang Y, Tang H, DeBarry JD, Tan X, Li J, Wang X, Lee T-h, Jin H, Marler B, Guo H. MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Res. 2012;40(7):e49.
Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.
Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.
Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11:1650–67.
Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4(1):17.
Dai X, Zhuang Z, Zhao PX. psRNATarget: a plant small RNA target analysis server (2017 release). Nucleic Acids Res. 2018;46(W1):W49–54.
We thank the High Performance Computing center of Northwest A&F University for providing computational resources in this work.
This work was mainly supported by the Open Project Program of the State Key Laboratory of Crop Stress Biology in Arid Areas, Northwest A&F University (Grant No. CSBAA2019001) and National Natural Science Foundation of China (Grant No. 31771778). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Multiple sequence alignment of the partial sequences of 20 HvMAPK proteins to identify the TDY and TEY motif. The red color marked sequence is the TDY or TEY motif.
Multiple sequence alignment of the full length sequence of 20 HvMAPK proteins to identify the conserved kinase motifs. The color marked indicated the conserved motifs found.
Multiple sequence alignment of the HvMAPKK to identify the conserved kinase motifs. The red color marked are the signature motif of MAPKK proteins.
Multiple sequence alignment of the HvMAPKKK to identify the conserved kinase motifs. The red color marked are the signature motif of MEKK, Raf and ZIK three sub family.
Evolutionary relationships and grouping among barley, rice and Arabidopsis MAPKs.
GO annotation of these identified barley MAPK cascade genes.
Motif identification based on PFAM database. Table S2. Characteristics of cis-acting regulatory elements in the promoter region of these identified barley MAPK cascade genes. Table S3. Duplicated MAPK cascade gene pairs identified in barley. Table S4. The Ka/Ks ratios for orthologous MAPK cascade proteins between barley and brachypodium. Table S5. The Ka/Ks ratios for orthologous HvMAPK cascade proteins between barley and rice sorghum. Table S6. The Ka/Ks ratios for orthologous MAPK cascade proteins between barley and maize. Table S7. The Ka/Ks ratios for orthologous MAPK cascade proteins between barley and sorghum. Table S8. The Ka/Ks ratios for orthologous MAPK cascade proteins between barley and soybean. Table S9. The Ka/Ks ratios for orthologous MAPK cascade proteins between barley and grape. Table S10. GO annotation of the identified barley MAPK cascade genes. Table S11. List of the putative miRNAs targeted on HvMAPK cascade genes identified by psRNATarget online tool. Table S12. Accession number and sample information of RNA-seq data using in this study.
About this article
Cite this article
Cui, L., Yang, G., Yan, J. et al. Genome-wide identification, expression profiles and regulatory network of MAPK cascade gene family in barley. BMC Genomics 20, 750 (2019). https://doi.org/10.1186/s12864-019-6144-9
- Gene family
- MAPK cascade
- Regulatory network