MAPK analysis

Of the genomes under study, the angiosperm clade has the highest number of MAPK members per family, with A. thaliana having the highest number of genes encoding MAPKs (20) followed by rice (O. sativa) and grape (Vitis vinifera) with 15 and 14 genes, respectively. In addition, species A. filiculoides and S. cucullata in the pteridophytes clade have 15 and 14 genes, respectively. The large number of MAPK genes in Azolla and Salvinia might be a result of whole genome duplication events in the Salviniales [60]. On the other hand, O. tauri in the algae group has the lowest number of genes of the analyzed species (3), while C. reindhardtti has 6. Bryophytes have an increase in the number of MAPK genes with respect to that of algae. The members of this clade, P. patens, S. fallax, and M. polymorpha contain 11, 8, and 6 genes, respectively. The MAPK gene family expansion in bryophytes is congruent with the expansion of gene families in land plants upon terrestrialization and might be associates with plant facing new types of stress, such as dehydration, gravity, exposure to ultraviolet light, etc. after land colonization [61, 62] (Fig 1; Table 1).

Among the 162 MAPK proteins identified, those belonging to group D are the longest. VvMPK7 (GSVIVT01018883001) of V. vinifera is the longest protein with 770 amino acids. In general, MAPKs proteins that belong to organisms from higher taxa such as angiosperms and gymnosperms have longer sequences than organisms of lower taxa except CrMPK18 (Cre17.g745447) of the alga C. reinhardtii with a length of 765 amino acids. On the other hand, the PaMPK2 protein from Picea abies has the shortest length (172 amino acids). Regarding the structure of the MAPK gene, there seems to be no specific pattern of intron organization in terms of the taxonomic level to which each species belongs including genes with zero (6 sequences), one (18), two (8), three (5), four (5), five (46), six (8), seven (10), eight (8), nine (33), ten (20), eleven (2), twelve (2), fourteen (2), fifteen (7) and sixteen (1) introns. The MpMPK3 gene of M. polymorpha contains the maximum number of introns (16), and none of the three genes identified in the unicellular algae O. tauri have introns. Furthermore, 18 of the sequences listed above contain at least one intron in the 3’ or 5’ UTR (UnTranslated Region) according to their annotated gene models (S1 Table). Regarding the organization of the groups, however, those proteins that belong to groups D and E have the highest number of introns while the members of group C have the least amount. Notably, gymnosperm and angiosperms organisms analyzed contain at least one intron. The predicted molecular weights of each of the MAP kinases reported here range from 20,559 Da (PaMPK2) to 87,480 Da (VvMPK7); while their isoelectric points (Ip) range from 4.66 (PaMPK2) to 9.72 (AfMPK16a) (S1 Table). Interestingly, group D and E comprise the largest proteins, with an average peptide length of 558 and 461, respectively (63,358 and 52,287 Da); these groups contain the most basic theorical isoelectric points averaging 8.45 and 8.10, respectively (S3 Fig).

Analysis of motifs and domains conserved in the MAPK family

MAPKs are characterized by the presence of multiple domains. Eight domains, named I-V, VIa, VIb, and VII are distributed towards the N-terminal preceding the activation domain. Towards the C-terminal end of the activation domain, there are five other domains previously described as typical of MAPKs. Domain VIII, IX, and X precede the CD domain; finally, the XI domain is localized around 30 amino acids after the CD domain (Table 2) [33, 35]. We used the MEME’s bioinformatics software to address domain conservation among the identified MAPK sequences [44].

Table 2

Consensus sequences of the conserved domains in MAPKs.

Domain	Consensus sequence
I	[V/P]-[I/V]-G-[K/R]-G-[S/A]-Y-G-[V/I]-V-C-S-A
II	E-X-V-A-I-K-K-I-X-[N/D]-[A/V/I]-F-[E/D]-[N/H]-X2-D-A
III	R-[T/I]-L-R-E-[I/L]-K-L-L-R-[H/L]-[L/M]-[R/D]
IV	P-X-[R/K]-X2-F-X-D-[V/I]-Y
V	V-[F/Y]-E-L-M-[E/D]-[T/S]-D-L-H-Q-[V/I]-I-[K/R]
VIa	[F/Y]-F-L-Y-Q-[L/I/M]-L-R-[G/A]-L-K-Y
VIb	H-[S/T]-A-N-[V/I]-[L/F/Y]-H-R-D-K-L-P-[K/S]-N-[L/I]-L-[A/L]-N
VII	C-[D/K]-L-K-I-[C/A]-D-F-G-L-A-R-[V/T]
VIII	[V/A]-T-R-W-Y-R-A-P-E-L-[L/C]-[L/G]-[S/N]
IX	A-I-D-[I/V/M]-W-S-[V/I]-G-C-I-F-[A/M]-E-[L/I/M]-[L/M]
X	P-[L/I]-F-P-G-X3-[V/L]-X-Q-L-X-L-[I/M]-T-[D/E]
XI	F-D-P-X2-R-[I/P]-[T/S]-[A/V]-X-[E/D]-A-L-X-[H/D]-P-Y-[F/L]
CD	L-H-D-X2-D-E-P

Of the 90 sequences analyzed that belong to angiosperms, 83 present all the canonical domains. The remaining seven (2 sequences from A. trichopoda and 5 from V. vinifera) lack at least one domain: VvMPK1 and VvMPK7 (domains I-IV), VvMK5a and VvMPK5b (III-V and VIa), VvMPK20 (I -III), AmtrMPK5 (IV), and AmtrMPK16a (XI). The fact that the vast majority of sequences belonging to angiosperms contain all the canonical domains could reflect a bias because this clade is the most studied (S4 Fig; S2 Table). Most of the sequences of the analyzed angiosperm species lack multiple domains. In P. taeda: PtMPK5a (III-V, and VIa), PtMPK5b (III-V, VIa, and IX), PtMPK5c (III-V, VIa, and IX-X), PtMPK11 (I -III), PtMPK16 (II-V, and VIa), and PtMPK20 (IX-X). In P. abies: we note PaMPK2 (I-VII), PaMPK4a (V), PaMPK7a (VIb and VII-XI), PaMPK7b and PaMPK10 (III, V, and VIa), PaMPK11 (I-III, and V), and PaMPK16 (IX and XI). Eight of these sequences belong to group E MAPKs, which include MHKs and whose inclusion as MAPKs, as already mentioned, has been debated due to their sequence similarity to CDKs [11]. Whether group E proteins function as MAPKs has been poorly studied, although an ortholog of these proteins in the fungus Ustilago maydis is activated by phosphorylation by a MAPKK and contains the activation motif TXY [63]. This observation supports the hypothesis that group E MAPKs should be regarded as a subfamily of MAPKs. In fern species, the AfMPK6b and AfMPK16b sequences of A. filiculoides lack domain I and domain IX, respectively. The AfMPK4d sequence lacks domains IV and IX. On the other hand, in S. cucullata, the ScMPK6b sequence lacks the V, VIa and VIb domains; the ScMPK16c sequence lacks domain XI. Domains III, IV and VII are absent in the ScMPK2 sequence. Finally, domain I is absent in the ScMPK3, ScMPK4a, and ScMPK4b sequences. In the case of I. echinospora, the sequence IeMPK4 lacks I-V domains while IeMPK7 and IeMPK11 lacks domain XI (S4 Fig; S2 Table).

The lack of several domains in the angiosperm and fern sequences is due to the fact that the genomes of these species are still fragmented, in the databases. Some of the ORFs that are found as hits might come from fragmented sequences and therefore the sequences corresponding to certain domains are missing. In the bryophytes M. polymorpha, S. fallax, and P. patens, we see MpMPK3 and MpMPK5 (III-V, and VIa), MpMPK4b (III and IV), PpMPK4b (II-IV, VIa, and IX), PpMPK4c (III-V, VIa, and IX), PpMPK15 (III, IV, VIa, and IX), PpMPK16a, and PpMPK16b (IX). For their part, the three sequences identified here from the alga O. tauri lack domain IV. Thus, the results of this analysis allow us to suggest that the presence of the activation motif (T-X-Y) together with five domains (II, III, V, and VII) located towards the N-terminal end of the activation domain and present in all the analyzed sequences are sufficient to identify MAP kinases from any plant species. In addition to this, the presence of domains VIII-XI in most Viridiplantae MAPK sequences suggests that these domains are of critical importance for maintaining either the function or structure of these proteins (S4 Fig; S2 Table).

MAPK proteins contain a T-X-Y motif whose phosphorylation leads to MAPK activation [14]. The Viridiplantae multiple sequence alignment shows that all analyzed species contain the canonical activation domain in the form of T-E-Y (98 sequences), T-D-Y (77 sequences), and T-S-Y (1 sequence). The alignment also shows the presence of non-canonical activation domains including M-E-Y (2 sequences), T-E-M (2 sequences), T-H-E (4 sequences), T-H-L (1 sequence), T-H-Q (1 sequence), T-K-T (2 sequences), T-Q-M (1 sequence), and T-S-Y (1 sequence). The distribution of these activation domains on the MAPK phylogeny shows that the T-E-Y domain is present in MAPKs from groups A, B, C, and E; while the T-D-Y is present on the D and E groups (S2 Fig; S2 Table). The activation domain contained in group E sequences exhibits the highest variability, which is especially evident on Bryophyta including T-D-Y, T-H-E, T-K-T, and T-Q-M variants (S2 Fig).

In addition to the 11 characteristic domains, MAP kinases also contain the sequence (L/H)-D-X₂-D-E-P known as the common docking domain (CD domain) [19]. Our results show that the CD domain is present in most proteins from the A, B, and C MAPK group; although it is also present in a few sequences from the groups D and E, it is generally absent in proteins from these clades (S4 Fig; S2 Table). The CD includes a motif involved in recognition and binding to MAPK substrates, as well as in protein-protein interaction with MEKs. The adjacent aspartate (D) and glutamate (E) residues in the CD, are essential for interaction with the basic residues lysine (K) and arginine (R), located in the binding site of MEK proteins [64]. It remains to be experimentally demonstrated whether the absence of the CD domain has any impact on the function of these group D and E MAPKs either for interaction and phosphorylation by MEKs or in MAPKs substrate specificity.

Identification of novel and distinctive domains for each group of MAPKs

The domain analysis of the retrieved sequences suggests the existence of six novel domains, named 12–17 (Fig 3; S3 Table). The presence of these domains might predict the group to which a MAPK protein belongs facilitating MAPK identification during genome or transcriptome annotation. The consensus sequences of such domains consist of 12–17 amino acid residues (Table 3).

Fig 3

Phylogenetic analysis and identification of novel domains in MAPKs.

The analysis of novel domains was performed with the MEME suite [44]. The distribution of the novel domains is group specific and could, therefore, facilitate MAPK identification and annotation. Sequence nomenclature follows a two-letter code to indicate genus and species (see text).

Table 3

Consensus sequences of the novel domains discovered in MAPKs.

Domain	Consensus sequence
12	G-N-X-F-E-V-[T/S]-X-K-Y
13	Y-X-[M/L]-W-[Q/R]-[T/S]-X-F-E-I-D-T-K-Y
14	H-[P/K]-D-I-V-E-[I/V/K]-[K/L]-[I/H/N]-[I/K]-[M/L]-L-P
15	M-L-X-F-[D/N]-P-X2-R-I-[T/S]
16	E-[L/V]-[I/L]-G-[T/S]-P-X-[E/D]-X-D-L-X-F-[L/I/V]
17	A-[R/K]-[R/K]-Y-[L/I/V]-X2-[L/M]-[R/P]-X3-[P/R/K]-X-[P/S]

Domain 12 is located approximately 15 residues before domain I from group A and B proteins. It is also present in OtMPK13 from O. tauri. Domain 13 is also located towards the N-terminus end of domain I but is only present in the proteins from group C therefore providing a signature for group C MAPKs. Domain 14 is located between domains II and III and typifies proteins from group D. Domains 15–17 are distributed towards the C-terminal end and, if present, they are always located before the CD domain between the canonical domains X and XI. Domain 15 is present in proteins from the B and E group with the sole exception of PtMPK16b from group E. Except VvMPK3 and VvMPK6 that belong to group A, domain 16 is only located in group B MAPKs. Domain 17 is present in most proteins except those in groups C and E. Interestingly, domain 17 is also present in the Fus3 kinase of S. cerevisiae this prompts the hypothesis that this domain was present in the ancestral MAPK sequences and was lost in groups C and E (Fig 3; S3 Table).

The conservation and distribution of these novel MAPK domains in a group-specific manner could suggest that they are involved in either maintaining the structure of MAPKs or in fine-tuning the selectivity towards substrates or protein interactors. This hypothesis is supported by homology structure prediction of selected MAPKs, which shows that domains 12, 13, 15, 15, and 17 are exposed to the solvent and could therefore provide an interaction surface. Meanwhile domain 14 is in a pocket (S5 Fig). Moreover, a computational prediction of phosphorylation sites suggests that domains 15 and 16 contain amino acid residues susceptible to phosphorylation (S4 Table). Whether these residues are indeed phosphorylated or if the identified novel domains are of relevance for the MAPKs structure and function remains to be experimentally demonstrated and extends beyond the scope of this study.

Phylogenetic analysis

A multiple alignment of the retrieved MEK sequences was generated this showed that MEK proteins are highly similar and that the 8 characteristic domains of MEK proteins are highly conserved as well as the ATP-binding and MAPK-binding domains (S8 Fig). Ichimura et al. classified MEK proteins into four groups, but this classification system was generated from angiosperm sequences only [14]. To avoid bias in the phylogeny reconstruction sequences retrieved from gymnosperms, seedless vascular plants, and non-vascular plants were included as well as the MEK protein Ste7 from S. cerevisiae as an outgroup to root the tree. The obtained ML topology shows that MEK proteins partition into five groups named A-E (Fig 4).

Fig 4

Phylogenetic tree of MEK proteins in Viridiplantae.

Maximum-likelihood tree reconstructed from amino acid sequences. The tree was rooted with the Ste7 protein from Saccharomyces cerevisiae. Bootstrap values from 1,000 replicates are indicated at nodes with support values >90%. Each MEK group is indicated with a distinctive color. Different species are indicated with different shapes and colors.

MEK proteins are involved in a plethora of biological processes; for example, some members of group A, such as AtMEK1 and AtMEK2 are involved in responses to environmental stimuli [35, 65] while AtMEK6 participates in cell division [66]. The MEK proteins included in group B have an extended C-terminal with the conserved sequence E-[K/R]-L-V-H-V-V-E-[K/N]-L-[H/Q]-C-X-A-X_1-4-G-[I/V]-X-I-R-V. This group had been reported earlier [11] and includes AtMEK3 whose C-terminal extension is similar to the NTF2 domain, which is involved in the nuclear import of protein cargos [11, 12, 49, 67]. Although the NTF2 domain is found in most eukaryotic lineages, its fusion to a MEK protein has only been observed in Viridiplantae [11, 49]. AtMEK3 is the only member of group B MEK proteins that has been characterized; it plays important roles in jasmonate [68] and blue light responses [69]. To the best of our knowledge, no biological function has yet been ascribed to any other member of the group not even for the NTF2 domain per se; therefore, the biological processes in which group B MEKs are involved remain obscured. The presence of the NTF2 domain in all group B MEKs including sequences from angiosperms, gymnosperms, pteridophytes and bryophytes but no algae strongly suggests that the NTF2–MEK domain combination occurred after the divergence of embryophytes this is estimated to have ocurred around 500 MYA [12]. Moreover, its conservation in taxa might reflect an evolutionary advantage in photosynthetic eukaryotes [49] although biochemical or functional characterization of group B MEKs are needed to corroborate this hypothesis.

The available information about groups C and D MEKs is scarce, but some members from this clade are involved in hormone signaling pathways. For example, AtMEK4/5 belongs to group C and mediates the ABA signaling pathways in response to biotic stress [34]. AtMEK7 and AtMEK9, belong to group D. They participate in the polar transport of auxins [70], and respond to ethylene [71], respectively. The novel group E proposed here contains sequences from bryophytes (M. polymorpha, P. patens and S. fallax), lycophytes (S. moellendorffii) pteridophytes (A. filiculoides and S. cucullata), gymnosperms (P. abies), and the early divergent angiosperm A. trichopoda. Interestingly, no group E MEKs were identified in monocot or dicotyledonous species, which suggests that these genes were lost prior to the Mesangiospermae origin, which diverged ~180 MYA [72]. Nonetheless, group E MEKs appears to be a sister clade to group C and D, which is well supported from bootstrap values (Fig 4).

Analysis of motifs and domains conserved in the MEK family

MEK proteins contain several conserved sequences to maintain the proper structure and function of these enzymes [73]. The MAPK-binding domain or docking-site (Table 5; S5 Table) is typically located towards the N-terminal end and its basic residues (lysine and arginine) interact with the acidic residues from the binding domain of MAPKs [64]. In addition, MEKs also contain an ATP-binding domain and the activation domain (Table 5; S5 Table) where serine/threonine residues are prone to phosphorylation and are involved in signal transduction [29]. Some proteins such as OsMEK7 or BdMEK7a have insertions within its activation domain, with one of serine/threonine residue absent or located around 3 to 5 residues upstream of the canonical position (S5 Table). Other proteins, such as ScMEK2, completely lack the MAPK-binding domain (S8 Fig; S5 Table). It is unclear whether these genes are expressed or if their gene products would exhibit MEK activity in MAP kinase cascades. If expressed, these proteins could represent neo or sub-functionalized versions arising from gene duplication events or even act as a noncatalytic scaffolding protein [12, 74].

Table 5

Consensus sequences of the conserved domains in MEKs.

Domain	Consensus sequence
I	G-X-[S/A/N]-[G/S]-G-X-V-X-[K/L]-[V/A]-X-H-[K/R]
II	P-X-[V/L]-V-X-[C/F]-[H/Y]-X2-[F/Y]
III	[L/M]-E-[Y/F]-M-D-X-G-S-L-[A/E]
IV	[A/S]-X6-L-X-G-L-X-Y-L-H
V	V-G-T-X2-Y-M-S-P-E-R-[I/F]
VI	[G/S/A]-D-[I/V]-W-S-[L/F]-G-[L/V]-X2-L-E
VII	[F/A]-S-X-E-[F/L]-[R/C]-X-F-[I/V]-X2-C
VIII	[S/T]-[A/V]-X2-L-L-X-H-P-F-[I/V/L]
ATP-binding	H-K-X3-8-A-L-K-X4-[N/D]-X-[D/E/Q]
MAPK-binding	[V/I]-H-R-D-[I/L]-K-P-[S/A]-N-L-L
Activation	V-S-X5-[S/T]-[M/L]-[D/G/A]
NTF2-binding	E-[K/R]-[L/I]-[V/I]-H-V-V-E-[K/N]-L-[Q/H]-C

A scan of the retrieved MEK sequences with the MEME suite confirmed the existence of eight domains named I-VIII (Table 5; S5 Table). Domains I-IV are located prior the activation domain and are present in the vast majority of the MEK proteins analyzed except for PpMEK3a, PaMEK6, and IeMEK10a, which lack domain I and VvMEK9 that lacks domain II. VvMEK5 lacks domains I and II. BdMEK9c and CrMEK6 do not contain domain III. PaMEK10a lacks domains I-III, and AhMEK3 and ScMEK2 lack domain IV. In contrast, domains V-VIII are located towards the C-terminal end after the activation site. Domains V and VI are absent only in AfMEK6b; domain VII is not present in VvMEK2, VvMEK9, IeMK10b, and IeMEKc; and domain VIII is not found in VvMEK2, VvMEK9, AmtrMEK6, PtMEK6b, AtMPEK8, and PaMEK6b. Additionally, six other sequences lack at least two domains: PaMEK3 (V, VI, and VIII); PaMEK4a (VII and VIII); PAMEK4B (V and VII); PaMEK6a and IeMEK6 (VII and VIII); and IeMEK10 (VII and VIII) (S8 Fig; S5 Table). Nevertheless, the presence of these signature sequences could be used as a criterion to refine MEK homology searches or annotation.

Identification of novel and distinctive domains for each group of MEKs

MEK proteins include a MAPK binding domain with the consensus sequence [V/I]-H-R-D-[I/L]-K-P-[S/A]-N-L-L. This domain is involved in protein-protein interaction with its MAPK partner, and therefore is essential for MAPK signaling cascades [75]. A closer examination of this sequence reveals group-specific variations of the MAPK binding domain, which suggests specificity, at least to some level, in MEK selection towards its MAPK interactors prompting the hypothesis that each MEK group might be involved in specific biological processes (Table 6).

Table 6

Predicted specific MAPK binding domain of plant MEKs.

Group	Sequence
A	H-X2-[R/K]-[H/R]-[I/V]-I-H-R-D-[I/L]-K-P-S-N-L-L
B	H-X-V-R-H-L-V-H-R-D-I-K-P-A-N-[L/M]-L
C	H-X2-R-X-I-V-H-R-D-X-K-P-[S/A]-N-L-L
D	H-X3-K-I-V-H-R-D-I-K-P-X-N-L-L
E	H-K-X2-[H/N]-K-I-V-H-R-D-I-K-P-S-N-L-L

Six new motifs spanning 10–18 amino acid residues were identified in addition to the canonical MEK domains. Remarkably, these new motifs are group-specific and therefore can be used to quickly assign new MEK sequences to their respective ortholog group. These domains were named with consecutive cardinal numbers [76–81] (Fig 5; Table 7). With the exception of PaMEK6b, IeMEK6, PtMEK6a, and OtMEK2, all proteins from group A contain domain 9 that is located towards the N-ter. Proteins from group B contain domain 10 towards the N-terminal and domains 12 and 13 are towards the C-terminal. MEKs from the C and E group contain domains 11 and 14, respectively, located towards the C-terminal of the protein. MEK proteins from the D group do not contain any domains in addition to the canonical ones (Fig 5; S6 Table). A computational analysis of post-translational modifications showed that all the putative novel domains contain potential phosphorylation sites, which might suggest that these domains are functionally relevant and possibly regulate MEK activity; this hypothesis, however, requires experimental validation (S4 Table).

Fig 5

Phylogenetic analysis and identification of novel MEK domains.

The analysis of novel domains was performed with the MEME suite [44]. The combinatorial distribution of the novel domains is group-specific and could facilitate future MEK identification and annotation. Sequence nomenclature follows a two-letter code to indicate genus and species (see text).

Table 7

Consensus sequences of the novel domains discovered in MEKs.

Domain	Consensus sequence
9	[A/Q]-S-G-T-F-X-D-G-D-[L/I]-X-[L/V]-N-X2-G
10	N-L-L-S-R-S-X3-Y-N-[I/F]-N-E-X-G-[L/F]
11	L-P-[L/M]-P-X-R-X2-D-X2-4-S-L-A-V-P-L
12	M-L-[T/A]-[V/I]-H-Y-Y-[L/M]-L-F-[D/N]-G-X-D
13	G-[V/I]-X-I-R-S-G-S-F-[I/V]-V-G-X2-F
14	[D/E]-P-P-X-P-P-X3-[S/T]-P-X-F-X2-F-I

Current knowledge on MEKK family

Plant MEKK are the least studied component of MAPK signaling cascades; therefore, there is limited information for these proteins compared to the information available for MAPKs and MEKs [10, 16, 35, 82]. The Arabidopsis genome encodes at least 80 genes distributed in three kinase subfamilies whose orthologs in animals have been shown to have MEKK activity: MAPKKK (MAPK/ERK kinase kinase, 21 members), ZIK (ZR1-interacting kinase, 11 members) and RAF (rapid acceleration fibrosarcoma, 48 members) [16]. Kinase activity has been confirmed for several members of the MAPKKK subfamily in plants, but no evidence has been reported for members of the ZIK subfamily [16, 34, 36, 54]. CTR1 (Constitutive Triple Response 1) and EDR1 (Enhanced Disease Resistance 1), which belongs to the RAF subfamily and acts primarily as a negative regulators of MAP kinase modules [71, 83–85], but no functional characterization has been made for other RAF members.

Plant MPKKKs have a variable MAP kinase domain [14]. Most proteins from the MEKK group have a regulatory domain and a serine/threonine residue at the activation site [86]. The proteins belonging to the MEKK group have the consensus sequence G-[T/S]-P-X-[F/Y/W]-M-A-P-E-V in the kinase domain [19, 21]. Those from the ZIK group contain the sequence G-T-P-E-F-M-A-P-E-[L/V/M]-[Y/F/L] and are located towards the N-terminal end. Meanwhile, members of the RAF group have the sequence G-T-X-X-[W/Y]-M-A-P-E. In addition, most proteins from the RAF group have the kinase domain located towards the C-terminal end, while the N-terminal end contains a regulatory domain.