Introduction
============

Based on their exposure to environmental stress, plants have constructed complex signaling networks to adapt to different stress conditions.^[@bib1]^ Protein kinases (PKs) play important roles in both growth and reproduction processes.^[@bib1],[@bib2]^ As a highly conserved sub-family of PKs, the mitogen-activated protein kinases (MAPKs) consist of the following three classes: (1) MAPKs, (2) MAPK kinases (MAPKKs), and (3) MAPKK kinases (MAPKKKs), which act as a signal transmission cascade in the cellular metabolic and transcriptional response centers.^[@bib1],[@bib2]^ MAPKs are activated by MAPKKs via phosphorylation of conserved threonine and tyrosine residues in the Thr--X--Tyr (T--X--Y) motif. MAPKKs, in turn, are activated by MAPKKKs when serine and threonine residues in the S/TXXXXXS/T motif are phosphorylated.^[@bib3; @bib4; @bib5; @bib6; @bib7; @bib8; @bib9]^

In eukaryotic organisms, MAPK pathways are relatively well-characterized cascades.^[@bib10]^ In plants, genome-wide analyses of MAPK cascades (MAPKKK---MAPKK-MAPK) have been reported in several species.^[@bib11; @bib12; @bib13; @bib14; @bib15; @bib16; @bib17; @bib18]^ In *Arabidopsis*, 23 MAPKs, 10 MAPKKs, and 80 putative MAPKKKs have been identified. MAPK and MAPKK genes in plants are each divided into four distinct classes (classes A, B, C, and D), whereas MAPKKK genes are classified into the following three subfamilies: (1) MEKK, (2) Raf and (3) ZIK.

Plant MAPKs are involved in both biotic stresses, such as AtMAPK4,^[@bib18],[@bib19]^ MsMMK2, and MsMMK3,^[@bib1]^ and in abiotic stresses, such as AtMAPK4,^[@bib20]^ OsMAPK4 ^[@bib21]^ and AtMAPK7.^[@bib22]^ MKK9-MPK3/MPK6 regulates ethylene signaling and may also play a role in leaf senescence.^[@bib9],[@bib23],[@bib24]^ Plant MKKs(MAPKKs) contain a phosphorylation site (S/T-X5-S/T) and a putative MAPK docking site in their N-terminal. Currently, only A and C type MAPKKs, such as AtMKK1, have been demonstrated to be activated by multiple abiotic stresses.^[@bib25]^ Preivous studies have examined that MEKK1, MKK2, MPK4, and/or MPK6 respond to salt, drought, and cold stress in *Arabidopsis*.^[@bib26]^ CsNMAPK regulated NO~3~^−^ stress,^[@bib27]^ ROS, and osmotic adjustment under salt stress.^[@bib28]^ Several reports revealed that MAPKKK play roles in drought resistance. A study in rice showed that DSM1, a RAF-Like MAPKKK gene has function in drought and oxidative stresses signaling.^[@bib29]^ A tobacco MAPKKK (NPK1) has been shown to enhance drought tolerance in transgenic maize.^[@bib30]^ The plant immune response is negatively controlled by the MEKK1(AtMEKK1)-MKK1(AtMKK1) and MKK2---MPK4-MKS1/WRKY33 cascades.^[@bib31],[@bib32]^ AtMKK1-AtMPK6 are involved in ABA---induced CAT1 expression,^[@bib33]^ MEKK1-MKK4/5 have an important function in plant innate immunity,^[@bib34]^ and MKK6 can regulate cytokinesis.^[@bib35],[@bib36]^

Studies on the evolution would be another interesting topic because their evolutionary history can shed light on functions. In land plants, the number of MAPKs and MAPKKs has expanded.^[@bib37],[@bib38]^ Before the divergence of the rice and Brachypodium, MAPKKK genes experienced the segment duplication.^[@bib38]^ The conservation of MAPKKKs have also implied by the strong purifying selection during the evolution.^[@bib39]^ During the long evolutionary history, most of higher plant lineages have undergone polyploidization.^[@bib40]^ In parallel with the whole-genome duplication (WGD) in angiosperm genomes, the three large families (MAPK, MAPKK and MAPKKK) experienced expansion and functional divergence. If the genes were involved in networks or function in a dose-sensitive manner, which should be retained that gene dosage hypothesis predicts.^[@bib41]^ *Brassica rapa* has undergone two duplication events (WGD α and β) and one whole-genome triplication event (WGT γ) within the Brassicaceae lineage.^[@bib42]^ Specifically, the WGT caused extensive fractionation in the *B. rapa* genome, which has the chance to study the molecular evolution of the MAPK cascade. Comparative analyses of many other plant species can yield the origin and evolution of these genes. Therefore, in addition to the detailed study of the MAPK cascade in *Brassica rapa,* we comparatively analyzed MAPK, MAPKK, and MAPKKK genes from eight representative plant species, including an algae (*Volvox carteri*), the ancestor of land plants,^[@bib43]^ a bryophyte (*Physcomitrella patens*), the earliest sequenced land plant,^[@bib44]^ a lycophyte (*Selaginella moellendorffii)*, an early vascular plant,^[@bib45]^ and five angiosperm plants: *Amborella trichopoda,* one of the earliest angiosperm, *Vitis vinifera,* the basal eurosid, and *Populus trichocarpa, Carica papaya* and *A. thaliana*.

In this study, we comprehensively analyzed the MAPK cascade genes in the *B. rapa* genome, as well as representative genomes in the tree of plant life, to understand the cascade retention and expansion patterns following the WGT event. In addition, divergent tissue-specific expression patterns and response to stresses of these genes will help us to understand the function of the MAPK cascade genes.

Materials and methods
=====================

Retrieval of sequences
----------------------

To identify MAPK cascade genes in *B. rapa*, *Arabidopsis*, and *O.sativa*, the sequences were retrieved from BRAD (<http://brassicadb.org/brad/>),^[@bib42]^ TAIR ([http://www.*arabidopsis*.org/](http://www.arabidopsis.org/)) and RGAP (<http://rice.plantbiology.msu.edu/>),^[@bib46]^ respectively. The sequences of others species were downloaded from Phytozome v9.1 (<http://www.phytozome.net/>).^[@bib47]^ MAPK cascade genes of *Arabidopsis* were performed against the proteins from other species by BLASTP (*e*-value \<1×10^−10^, identity \>70%). The serine/threonine-protein kinase-like domain (PF00069), was downloaded from Pfam (<http://pfam.sanger.ac.uk/>), and used to search the *B. rapa* protein database by HMMER v3.0 (<http://hmmer.janelia.org/>). Potential sequences were detected by SMART (<http://smart.embl-heidelberg.de/>).

Identification of gene synteny and duplicated MAPK cascade genes
----------------------------------------------------------------

According to previous reports, gene synteny and duplicated types were detected by the MCScanX (Multiple Collinearity Scan toolkit).^[@bib48]^ The positions of syntenic MAPK cascade genes in the *B. rapa* on the ancestral genomic blocks were searched from the BRAD website (<http://brassicadb.org/brad/searchSynteny.php>).^[@bib42]^ Circos software was used to draw syntenic diagram.^[@bib49]^

Phylogenetic analysis and characterization of the MAPK cascade genes
--------------------------------------------------------------------

MEGA v6.0 was used to construct the maximum-likelihood (ML) trees with a bootstrap value of 1000 replications. MEGA v6.0 was used to calculate genetic distances. The MEMEv.4.10.1 tool (<http://meme.sdsc.edu/meme/>) was used to identify conserved motifs in the MAPK cascade proteins. The MAPK cascade gene structures were drawn by the online tool Gene Structure Display (GSDS, <http://gsds.cbi.pku.edu.cn/>). The STRING tool was used to construct interaction network of MAPK cascade proteins in *B. rapa* (<http://string-db.org/>).^[@bib50]^

The MAPK cascade genes Ka/Ks analysed
-------------------------------------

The MAPK cascade genes from *B. rapa* were aligned using ClustalX 2.0 with the default parameters.^[@bib51]^ The KaKs_calculator was used to calculate the Ks and Ka values.^[@bib52]^ The Ks values for all the MAPK cascade homologous genes in *B. rapa* and *A. thaliana* were shown by boxplot using the R program (boxplot package).^[@bib53]^ The following formula: *T*=Ks/2*r* was used to calculate divergence time (*r*=of 1.5×10^−8^, synonymous substitutions per site per year).^[@bib54]^ Nucleotide distance was used to analyze the relationships among Classes I, II, III and IV.

Expression pattern analysis for MAPK cascade genes
--------------------------------------------------

Illumina RNA-seq data^[@bib55]^ were used to characterize the expression of the MAPK cascade genes in *B. rapa*. The expression level of MAPK cascade genes were analysed by the fragments per kilobase of exon model per million mapped reads (FPKM). The AtGenExpress Visualization Tool (AVT; <http://jsp.weigelworld.org/expviz/expviz.jsp>) was used to analyse the expression level^[@bib56]^ in *A. thaliana*. The R program (VennDiagram package) was used to draw Venn diagrams.

Promoter sequence analysis
--------------------------

Search for *cis*-regulatory elements in promoter of MAPK cascade genes sequences was performed using the PLACE Web Signal ScanPLACE (<http://bioinformatics.psb.ugent.be/webtools/plantcare/html/>).

Results
=======

Different retention of MAPK cascade genes following the WGT event in *Brassica rapa*
------------------------------------------------------------------------------------

BLAST software was used to identify the number of MAPK, MAPKK, and MAPKKK genes in *A. thaliana* and *B. rapa*. Then, the MAPK cascade genes containing the serine/threonine-protein kinase-like domain was searched by an HMM (PF00069) in *B. rapa*. Finally, all the homogeneous candidate MAPK cascade sequences were subjected to Pfam, SMART and NCBI database analyses. After multiple screening and validation steps, we finally identified 34 BrMAPK, 14 BrMAPKK, and 112 BrMAPKKK genes ([Supplementary Table S1](#xob1){ref-type="supplementary-material"}).

To investigate different syntenic genes in the MAPK cascade during the *B. rapa* WGT event, we completed a syntenic analysis with MCScanX. In total, 31 MAPKs, 14 MAPKKs, and 100 MAPKKKs were located in the syntenic regions ([Figure 1](#fig1){ref-type="fig"}). About 91% (31/34) of the MAPK genes, 93% (13/14) of the MAPKK genes, and 89% (100/112) of the MAPKKK were retained in the syntenic regions ([Figure 1](#fig1){ref-type="fig"}). We also counted the number of genes and analyzed the distribution of the three subgenomes by comparing the retention of MAPKs, MAPKKs and MAPKKKs ([Figures 1d and e](#fig1){ref-type="fig"}). The proportion of retained MAPK cascade genes was higher in the least fractionated (LF) subgenome than in the medium fractionated (MF1) and most fractionated (MF2) subgenome.This finding is consistent with a previous report in which the degree of retained genes in these three sub-genomes (LF, MF1, and MF2) decreased gradually.^[@bib42]^ More MAPKKKs were retained in the MF2 sub-genome than in the MF1 sub-genome. MAPK (59%), MAPKK (40%), and MAPKKK (46%) genes retained more than two homologous genes, and the MAPK genes had higher retention than the two other gene families (MAPKK and MAPKKK). However, some MAPK (18%), MAPKK (10%) and MAPKKK (16%) genes were lost ([Figure 1e](#fig1){ref-type="fig"}). Therefore, MAPK genes were preferentially retained compared to the other two gene families. The MAPKK genes were similarly retained to MAPKKK genes following the WGT in *B. rapa* ([Figure 1](#fig1){ref-type="fig"}, [Supplementary Table S2](#xob1){ref-type="supplementary-material"}).

Characteristic structure analysis of BrMAPK cascade proteins
------------------------------------------------------------

To study the characteristic of MAPK cascade genes in *B. rapa*, we used the Maximum Likelihood method to construct phylogenetic trees ([Figure 2](#fig2){ref-type="fig"}). All the BrMAPKs were classified into four classes, I, II, III and IV ([Figure 2](#fig2){ref-type="fig"}). We searched using the MEME tool to identify conserved motifs of BrMAPK proteins ([Figure 2a](#fig2){ref-type="fig"}). Eight of the 12 motifs were conserved in all the BrMAPK proteins, but the remaining 4 motifs (9--12) were not. In addition, each subfamily had its specific motif. For example, aside from the conserved motifs, MAPK proteins in classes I and II had specific motif 10 in their N-terminal region, whereas those in class IV contained motifs 8 and 9 in their C-terminal region ([Figure 2](#fig2){ref-type="fig"}). The annotated motifs from MEME revealed that eight of the 10 motifs (1--8) corresponded with the kinase domains in MAPKs ([Supplementary Figure S1](#xob1){ref-type="supplementary-material"}). Class II BrMAPK genes had no introns and all the BrMAPKK geenes contained a protein kinase domain ([Supplementary Figure S1](#xob1){ref-type="supplementary-material"}).

The BrMAPKK genes were clustered into four distinct gene classes (I, II, III and IV) ([Figure 2b](#fig2){ref-type="fig"}), and varied in kinase domain presence or absence, adenosine triphosphate (ATP) binding site and serine/threonine protein kinase active site according to phylogenetic analysis ([Supplementary Figure S1](#xob1){ref-type="supplementary-material"}). The BrMAPKK class III and IV subfamily genes have no introns, except for BrMAPKK9, which had one intron ([Supplementary Figure S1](#xob1){ref-type="supplementary-material"}). The BrMAPKKK gene family in *B. rapa* was clustered into three subfamilies (MEKK, RAF and ZIK) ([Figure 2c](#fig2){ref-type="fig"}). All the BrMAPKKK proteins also had a protein kinase domain ([Supplementary Figure S1c](#xob1){ref-type="supplementary-material"}).

Duplication and Ks analysis for MAPK cascade genes in *A. thaliana* and *B. rapa*
---------------------------------------------------------------------------------

We reconstructed 24 conserved ancestral genomic blocks (GBs) according to a previous report.^[@bib57]^ The color-coding of these blocks was based on their positions in a proposed ancestral karyotype (AK1--8).^[@bib57],[@bib58]^ To investigate the BrMAPK cascade genes of chromosomal location, all of them were mapped onto ten chromosomes in *B. rapa*, except *BrZIK4*, which is located on Scaffold000203 ([Supplementary Figure S2](#xob1){ref-type="supplementary-material"}). Among the ten sets of chromosomes in *B. rapa*, chromosome 09 contained more BrMAPK cascade genes (30 genes) than any of the other chromosomes. Among the LF, MF1 and MF2 subgenomes, LF contained more MAPK cascade genes than the MF1 and MF2 subgenomes. The most BrMAPKKK genes (18%) clustered in the AK8 region, whereas for BrMAPK and BrMAPKK genes, the most (20%) clustered in AK1 region ([Supplementary Figure S2](#xob1){ref-type="supplementary-material"}).

Gene duplication was the driving force for gene families evolution.^[@bib59]^ We used the MCScanX software to detect the MAPK cascade genes duplication type. All of the BrMAPK cascade genes were traced to one of five duplication types (singleton, WGD, tandem, proximal or dispersed duplication types) ([Figure 3](#fig3){ref-type="fig"}). Twenty-seven BrMAPK, 8 BrMAPKK and 91 BrMAPKKK genes in *B. rapa* were duplicated from WGD or segmental events compared to only 1 BrMAPK, 1 BrMAPKK and 1 BrMAPKK genes from tandem duplication and 6 BrMAPK, 5 BrMAPKK and 27 BrMAPKK genes from dispersed duplication ([Figure 3](#fig3){ref-type="fig"}). The expansion of the BrMAPK cascade gene family was derived by the WGD or segmental duplication and dispersed gene duplication in *B. rapa*. The BrMAPKKK gene Ks values ranged from 0.2 to 0.5 and focused on about 0.25 (\~8.5 Myr), whereas the BrMAPKK genes ranged from 0.4 to 0.7 and had a mean of about 0.40 (\~12.5 Myr; [Figure 3a](#fig3){ref-type="fig"} and [Supplementary Table S3](#xob1){ref-type="supplementary-material"}). The BrMAPK gene Ks values ranged from 0.25 to 0.5 and centered \~0.3 (\~10.5 Myr).The divergence time for the BrMAPKKK duplicated gene pairs was 8.49 MYA (Million years ago), which indicates that their divergence occurred during the *Brassica* triplication event (5--9 MYA). The BrMAPK and BrMAPKK duplicated gene happened 10 to 16.6 MYA and occurred during the divergence of *B. rapa* from *Arabidopsis* (9.6--16.1 MYA).

Expansion and evolution of MAPK genes in plants
-----------------------------------------------

We performed a phylogenetic tree of MAPK genes with additional eight key species in the tree of plant life to survey the evolutionary relationships of the MAPK gene family in plants. We selected four angiosperms *Vitis vinifera, P. trichocarpa* and *C. papaya* because these species did not undergo α and β duplications, and *A. trichopoda* is a basal angiosperm that did not undergo the γ duplication event.^[@bib60; @bib61; @bib62]^ A total of 204 MAPK genes were identified from the nine plant species ([Supplementary Table S4](#xob1){ref-type="supplementary-material"}; [Figure 4a](#fig4){ref-type="fig"}). All the MAPK genes segregated into four distinct classes (classes I, II, III and IV) with 100% bootstrap support ([Supplementary Figure S3](#xob1){ref-type="supplementary-material"}). The nucleotide distance between class I subfamily with class II, III and IV subfamilies was higher than those between MAPK classes II and III subfamilies, class II and IV subfamily, and classes III and IV subfamilies ([Supplementary Figure S4](#xob1){ref-type="supplementary-material"}). These findings indicated that class I had a distant phylogenetic relationship with the other three subfamilies, which is consistent with the results from the phylogenetic trees ([Figure 4](#fig4){ref-type="fig"}).

The gene numbers for the four classes in the MAPK gene family were marked in the nine plant species ([Figure 4b](#fig4){ref-type="fig"}). Only one class I and three class II subfamilies of MAPK genes were found in *Volvox carteri*, suggesting that MAPK gene expansion happened after the divergence of green algae. All these 9 species had conserved numbers (1--5) of class II MAPK genes, except *P. trichocarpa*, which had 10 MAPK genes. Based on the phylogenetic relationships, the class II subfamily of MAPK genes originated from green algae, whereas classes III and IV subfamilies originated from land plants ([Figure 4c](#fig4){ref-type="fig"}). Class II subfamily of MAPK genes diverged from Bryophyta class III and IV subfamilies, and that the rapid expansion of classes II, III and IV subfamilies may have primarily led to MAPK gene expansion in angiosperm plants. The evolution of MAPK genes was re-constructed and suggests that classes III and IV subfamilies originated from the class II subfamily, which was the initial group, and that the class III subfamily separated from the class IV subfamily in plants classified in the Bryophyta ([Figure 4d](#fig4){ref-type="fig"}).

Expansion and evolution of MAPKK genes in plants
------------------------------------------------

Based on the phylogenetic trees for all the MAPKKs, we tried to re-construct the expansion history of the MAPKK families. A total of 84 MAPKK genes were identified in the nine species, although only two MAPKK genes were detected in *V. carteri*. According to the phylogenetic tree topologies, the MAPKK genes could be grouped into four classes (classes I--IV) ([Supplementary Figure S5](#xob1){ref-type="supplementary-material"}). The MAPKK genes from *P. patens* and *S. moellendorffii* were clustered into only three MAPKK subfamilies included classes I, II and III, whereas all MAPKK subfamilies were present in the angiosperms, which expanded rapidly, and was accompanied by several WGD ([Supplementary Figure S6](#xob1){ref-type="supplementary-material"}). The number of MAPKK genes in most angiosperms was greater than that in the other non-angiosperm species (*V. carteri* and *P. patens*), and was most likely due to WGDs, which leads to gene expansion ([Figure 5](#fig5){ref-type="fig"}). The number of MAPKK genes in class II was a highly conserved in the all species except *A. trichopoda* (2), *P. pteridophyta* (3) and *V. carteri* (no class II MAPKK genes were identified; [Figures 5a and b](#fig5){ref-type="fig"}). One class I MAPKK gene subfamily was found in *V. carteri* and suggested that MAPKK gene expansion occurred in land plants. However, no class II MAPKK genes were detected in *V. carteri* that indicated the class II gene subfamily may be unique to land plants.

The four MAPKK classes have similar domain architectures and motif compositions ([Supplementary Figure S7](#xob1){ref-type="supplementary-material"}). In addition to the domain compositions identified by pfam, 10 common motifs embedded in MAPKK gene domains were identified by MEME for all of the MAPKK genes. Class I subfamily of MAPKK genes had motif 8. Although the two MAPKK genes in *V. carteri* do not have motif 8, all the other species have similar motif compositions. From algae to angiosperm, the MAPKK gene family has highly conserved domains and motifs.

Phylogenetic trees were constructed in the same nine species to investigate the evolution of the MAPKK genes ([Supplementary Figure S6](#xob1){ref-type="supplementary-material"}). Genetic distances were analyzed to examine the relationships among the four classes ([Figure 5c](#fig5){ref-type="fig"}). The genetic distances for class III subfamily versus class IV subfamily were smaller than class III versus class II subfamily, and class III subfamily versus class I subfamily, which indicated that class III subfamily had a closer relationship with class IV subfamily. Moreover, the genetic distances for classes I, II, and III subfamilies were similar, whereas that of subfamily IV was the smallest, which suggested that the degree of sequence divergence for class IV subfamily was lower than the other subfamilies ([Figure 5c](#fig5){ref-type="fig"}). The evolutionary history of MAPKKs was surmised that class I subfamily was the initial class from which MAPKKs originated. Additionally, the class IV subfamily separated from class III subfamily in angiosperm plants.

Expansion and evolution of MAPKKK genes in plants
-------------------------------------------------

We constructed a phylogenetic tree of the MAPKKK genes that clustered into three subfamilies (MEKK, ZIK and Raf) ([Figure 6a](#fig6){ref-type="fig"}, [Supplementary Figure S8](#xob1){ref-type="supplementary-material"}). The MEKK, ZIK and Raf subfamily were detected in *V. carteri* ([Figure 6b](#fig6){ref-type="fig"}) that suggested they all originated from the green algae. Genetic distances were analyzed to further determine the relationships among the three subfamilies ([Figure 6c](#fig6){ref-type="fig"}). The genetic distances between MKKK subfamily and the Raf subfamily were smaller than MEKK subfamily versus other subfamilies ([Figure 6c](#fig6){ref-type="fig"}). These results indicated that MEKK subfamily has a close relationship with Raf subfamily and that MEKK and Raf subfamilies may share a common evolutionary origin.

Analysis of MAPKKKs family size in nine plants indicated that it underwent rapid expansion during evolution, and further expansion in the Brassicaceae ([Figure 6b](#fig6){ref-type="fig"}). Phylogenetic trees were constructed for each species and the MAPKKK family divided into three classes (MEKK, ZIK and Raf) ([Supplementary Figure S8](#xob1){ref-type="supplementary-material"}). The WGD had important role in the expansion and evolution of gene families in plant genomes, the rapidly increased gene numbers of the MAPKKK family in *B. rapa* (18) and *A. thaliana* (11) ([Figure 6b](#fig6){ref-type="fig"}) gave a good example for that. Notably, the number of MAPKKKs in *P. trichocarpa* and *B. rapa* is more than other species, which was caused by salicoid duplication and Brassica-specific WGD events. During expansion of the MAPKKK gene family, the Raf subfamily played a major role in this process.

We compared gene structure characteristics between algae and higher plants to better understand the evolutionary characteristics of MAPKKK genes in the plant kingdom ([Figure 6d](#fig6){ref-type="fig"}). The gene structure of Raf was quite simple in lower species with most genes having just five introns, whereas the number of introns increased to about 10 in the angiosperms ([Figure 6d](#fig6){ref-type="fig"}). The length of the protein has been stable at about 386 amino acids. MEKKs have no introns in *V. carteri* and *P. pteridophyta* ([Figure 6d](#fig6){ref-type="fig"}).

Comparative expression patterns of the MAPK cascade genes in different tissues from *B. rapa* and *A. thaliana*
---------------------------------------------------------------------------------------------------------------

The MAPK cascade gene expression patterns between *A. thaliana* and *B. rapa* were investigated in roots, stems, leaves, flowers and siliques (callus only was studied in *B. rapa* and mature pollen studied in *A. thaliana*; [Figure 7](#fig7){ref-type="fig"} and [Supplementary Tables S5 and S6](#xob1){ref-type="supplementary-material"}). Gene expression of BrMAPKs was detected in the roots, stems, leaves and flowers (28 BrMAPKs, 82%) and siliques (31BrMAPKs, 91%; [Figure 7a](#fig7){ref-type="fig"}). All AtMAPKs showed expression (mean---normalized value \>1) in all tissues ([Figure 7a](#fig7){ref-type="fig"}). Gene expression of BrMAPKK and AtMAPKK was detected in at least one of the tissues, except *BrMAPKK13* and *BrMAPKK14* in which no expression was found in any of the tissues. Among the 196 BrMAPKKKs (including 84 AtMAPKKKs and 112 BrMAPKKKs), 5 (*BrRaf36*, *BrRaf70*, *BrRaf63*, *BrZIK11*, and *BrRaf71*) were not expressed, and 2 (*BrMEKK7* and *BrRaf21*) were expressed only in the flower tissue ([Figure 7c](#fig7){ref-type="fig"} and [Supplementary Figure S9](#xob1){ref-type="supplementary-material"}). Additionally, a total of 174 MAPKKKs had high expression levels (FPKM value \>10) in the five tissues. We also analyzed MAPK cascade genes for tissue-specific expression profiles in five tissues; some BrMAPK cascade genes were highly expressed (FPKM value \>10) in flowers, including *BrMAPK32*, *BrMAPKK12,* and *BrMEKK7*, whereas they had no or low expression in the other tissues ([Figures 7a and c](#fig7){ref-type="fig"}).

MAPK cascade genes in the five tissues on the phylogenetic tree were selected to survey divergent functions of the homologous genes ([Figure 7](#fig7){ref-type="fig"}). AtMAPKKs Class IV subfamily showed high expression (mean-normalized value \>1) in the five tissues, except *BrMAPKK13* and *BrMAPKK14* were not expressed. Class III AtMAPKs also showed expression in all five tissues, except *BrMAPK30* and *BrMAPK34* were not expressed, and *BrMAPK31*, *BrMAPK32* and *BrMAPK33* were expressed only in the roots. These results indicated that class IV MAPKK and class III MAPK genes may have lost some of their functions after duplication events according to the expression patterns.

*Cis*-elements and interaction network analysis among MAPK cascade proteins in *B. rapa*
----------------------------------------------------------------------------------------

To analyze the *cis*-regulatory elements, we used the Plant *Cis*-acting Regulatory DNA Elements (PLACE) online tool to identify MAPK cascade genes in *B. rapa*. A total of 11 common *cis*-regulatory elements were detected in the promoter regions of the MAPK cascade genes and were highly conserved among all of the MAPK cascade genes studied in *B. rapa* ([Figure 8](#fig8){ref-type="fig"}). Three common *cis*-regulatory elements, ARBE, the CGTCA-motif, and the GARE-motif, were responsive to plant hormones, including ABA, JA and GA, and suggested that they could affect MAPK cascade gene expression levels in *B. rapa*. Some common *cis*-regulatory elements responsive to both abiotic and biotic stresses were fungal elicitor-responsive elements (W-box and Box-W1), a light-responsive element (G-Box), low-temperature responsiveness (LTR), defense and stress responsiveness (TC-rich repeats), and a drought-responsive element (MBS). The diverse regulatory elements indicated the importance of MAPK cascade genes in stress tolerance. There were two more elements including a regulatory element essential for anaerobic induction (ARE) and an elicitor-responsive element (EIRE). In addition, we statistics the number of genes that contained *cis*-regulatory in each subfamily ([Figure 8d](#fig8){ref-type="fig"}). In the MAPK, most genes contianed more MBS *cis*-regulatory elements than other *cis*-regulatory elements. Only the class II subfamily of MAPKK had six *cis*-regulatory elements. In MAPKKK, most genes also contained *cis*-regulatory elements of MBS, which is similar to MAPK genes. These results indicated that MAPK cascade genes were responsive to various stresses, which may be due to the binding of corresponding *cis*-elements that regulate the expression of MAPK cascade genes.

The interaction network of BrMAPK cascade proteins was constructed by STRING software ([Supplementary Figure S10](#xob1){ref-type="supplementary-material"}). The lines stand for the positive correlations that were retrieved from STRING database corresponding to the protein interactions. The interaction netwok contained 42 nodes and 167 pairs of interacting genes. It showed a very complicated correlation among the BrMAPK cascade proteins ([Supplementary Figure S10](#xob1){ref-type="supplementary-material"}). This phenomenon indicated that MAPK cascade genes are involved in many fundamental mechanisms and regulate many downstream factors and/or are regulated by many upstream genes.

Discussion
==========

MAPK cascade genes have been shown to be involved in growth, development and responses to a variety of stress stimuli. The MAPK signaling pathway is one of the most important and conserved ubiquitous modules for signal transduction in plants ^[@bib14],[@bib20],[@bib21],[@bib63],[@bib64]^ ([Supplementary Figure S11](#xob1){ref-type="supplementary-material"}).

If genes products are involved in complex regulatory networks, via different metabolic pathways or transcriptional regulation, they are likely to be retained.^[@bib41],[@bib65]^ In this study, we identified 34 BrMAPK, 15BrMAPKK and 112 BrMAPKKK genes in the *B. rapa* genome that displayed high retention after the WGD. More MAPK (59%) and MAPKKK (46%) genes retained two or three copies in *B. rapa* than MAPKK (40%) genes. By analyzing the duplicated types, we found more segmental duplication in BrMAPKs (79.4%) and BrMAPKKKs (81.3%) than in BrMAPKKs (53.3%). BrMAPKs and BrMAPKKKs that participate in development, stress and defense, and play particularly important roles in signaling pathways. This result is agreement with the gene dosage hypothesis.^[@bib66]^ Comparison of nucleotide distances showed that BrMAPKs and BrMAPKKKs diverged approximately 8.5 MYA, earlier than BrMAPKKs that separated at 12.5 MYA.

During evolutionary, angiosperm genomes have experienced polyploidization events, which have led to the MAPK cascade gene expansion, followed by divergence by (1) subfunctionalization, (2) neofunctionalization, or (3) nonfunctionalization.^[@bib62],[@bib67; @bib68; @bib69]^ These fates have provided choices for duplicated genes to gain functional diversification.

We determined the evolutionary histories of the MAPK cascade genes through comparison analyses that examined expansion, phylogenetic relationships, characteristic structures and nucleotide distances. The MAPK class II subfamily had a close relationship with the classes III and IV subfamilies. We constructed a MAPK evolutionary model whereby classes III and IV subfamilies might have originated from class II, which was the ancestral group, followed by class III subfamily separating from class IV subfamily in Bryophyta ([Figure 4](#fig4){ref-type="fig"}). At present, MAPK cascade gene family members have been reported from a few plant species, including rice (*Oryza sativa*),^[@bib70]^ Arabidopsis (*A. thaliana*),^[@bib71]^ maize (*Zea mays*),^[@bib72]^ Brachypodium distachyon,^[@bib73]^ canola (*B. napus*),^[@bib74]^ apple (*Malus domestica*)^[@bib75]^ and grape (*V. vinifera*).^[@bib76]^ MAPK genes had conserved amino acids sequences that reflects their common evolutionary ancestry.^[@bib77]^ Before the divergence of the monocots and eudicots, the plant MPK gene family structure already was established from an ancient pattern of diversification.^[@bib77]^ MAPK cascades play central roles in the eukaryotic signal transduction networks, and they have highly conserved component kinases during evolution.^[@bib78]^

In flowering plants, MAPK genes number was twice that of MAPKKs.^[@bib37]^ The MAPKK gene family has retained highly conserved domains and motifs from algae to angiosperms ([Supplementary Figure S7](#xob1){ref-type="supplementary-material"}). In their evolutionary, class II and III subfamilies of MAPKK genes might have originated from class I subfamily, which was the ancestral class, followed by class IV subfamily separating from class III subfamily in angiosperms. The MAPKK genes in *Arabidopsis* and *B. rapa* possessed similar gene complements and experienced several recent gene duplication events. After comparing the MAPKK genes in *V. vinifera* and *B. rapa*, we found that *B. rapa* had a relatively large number of MAPKK genes compared to *V. vinifera*. The reason might be that *V. vinifera* had not undergone α and β duplication events. The MAPKK gene family formed four distinct subfamilies in the phylogenetic tree, which is consistent with the result for *B. rapa* and *A. thaliana*.^[@bib37]^

In the three subfamilies in the MAPKKK family, MEKK has a closer relationship to Raf than MEKK to ZIK, and plant MEKKs and Rafs may share a common evolutionary origin. Moreover, the increase of Rafs played a major role in the expansion of the MAPKKK gene family. An evolutionary history of MAPK cascade gene family was inferred in the plant kingdom; all three gene families exist in the genomes from Bryophyta to Angiospermae ([Figure 9](#fig9){ref-type="fig"}, [Supplementary Figure S12](#xob1){ref-type="supplementary-material"}).

Many of the MAPK pathway members have been characterized according to their response to various environmental stimuli in plants.^[@bib63]^ The MEKK1-MPK4 cascade plays an important role in ROS (reactive oxygen species) metabolism,^[@bib79]^ and MKK1--MPK6 cascade is involved in H~2~O~2~ metabolism.^[@bib33]^ A TIPK (*Trichoderma*-induced MAPK) was reported to response to pathogen,^[@bib80]^ whereas CsNMAPK regulates NO~3~^−^ stress,^[@bib27]^ and salt stress.^[@bib28]^ Plant MAPKs play roles in different hormones pathways,^[@bib81]^ such as MKK3/MPK6 participate in jasmonic acid (JA) signaling.^[@bib82]^ MAPKs were involved in abscisic acid (ABA) signaling pathway.^[@bib33]^ To better understand the possible transcriptional regulation functions of MAPK cascade genes in *B. rapa*, we identified common conserved *cis*-regulatory elements in the promoter regions of MAPK cascade genes. We also found that 11 common *cis*-elements in all the MAPK cascade genes of *B. rapa* were responsive to plant hormones and environmental stresses.

Conclusions
===========

In summary, genome-wide analysis of MAPK cascade gene families was first performed in *B. rapa*, and 34 BrMAPK, 14 BrMAPKK and 112 BrMAPKKK genes were identified. We analyzed the characteristic structures, evolutionary patterns, gene duplication, expression divergence and *cis*-regulatory elements of MAPK cascade genes in plants. Comparative evolutionary analysis from eight other plants was also conducted in the tree of plant life. Our study provides new insights into the evolution of the MAPK cascade gene family in *B. rapa* and the plant kingdom in general. We further confirmed the significant regulation of the MAPK cascade genes in *B. rapa* during growth, development and stress tolerance. Our study provides useful information regarding the MAPK cascade genes functional divergence and conservation, and also aids in understanding the effects of polyploidy during evolution with regards to MAPK cascade genes.

This work was supported by the key project of National Natural Science Foundation of China (No. 31330067), the National Natural Science Foundation of China (No. 31471886), and National vegetable industry technology system (CARS-23-A-06), and Jiangsu modern agriculture (vegetable) industrial technology system (SXGC\[2017\]273).

[Supplementary Information](#xob1){ref-type="supplementary-material"} for this article can be found on the *Horticulture Research* website (http://www.nature.com/hortres)

The authors declare no conflict of interest.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

![The syntenic MAPK cascade genes between *Brassica rapa* and *Arabidopsis thaliana* and their different retentions. (**a**--**c**) MAPK, MAPKK, MAPKKK syntenic genes lines are shown among the 10 *B. rapa* chromosomes (Br01-Br10) and the five *A. Thaliana* chromosomes (At01--At05). The red lines represent the syntenic genes pairs between *B. rapa* and *Arabidopsis*, and the blue lines represent the syntenic genes in *B. rapa*. (**d**) Retention of homoeologs of MAPK cascade genes in the three subgenomes of LF (least fractionized subgenome), MF1(moderately fractionized subgenome), and MF2 (most fractionized subgenome) in *B. rapa*. (**e**) Numbers of MAPK cascade genes after genome triplication and fractionation in *B. rapa*.](hortres201779-f1){#fig1}

![Schematic diagram of amino-acid motifs of MAPK cascade genes in *Arabidopsis* and *B. rapa*. (**a**) Phylogenetic relationships and protein motifs of MAPK genes. (**b**) Phylogenetic relationships and protein motifs of MAPKK genes. (**c**) Phylogenetic relationships and protein motifs of MAPKKK genes.](hortres201779-f2){#fig2}

![Pairwise comparison of the Ks values for MAPK cascade homologous genes in *B. rapa* and *A. thaliana* and their different duplicated types. (**a**) The distribution of Ks values for MAPK cascade genes between *B. rapa* and *A. thaliana*. (**b**--**d**) Singleton, dispersed, proximal, tandem and segmental duplicated types of MAPK cascade genes in *B. rapa*.](hortres201779-f3){#fig3}

![The duplication origin of MAPK genes in plants. (**a**) Phylogenetic relationships among MAPK genes. (**b**) Comparison of numbers of MAPK genes in representative species. (**c**) The phylogenetic relationships of MAPK genes. The green stars denote the inferred whole-genome duplication in angiosperms. (**d**) Schematic representation of the duplication history of MAPKs in plants. The red, yellow, blue and violet lines indicate class I, II, III and IV MAPK genes, respectively.](hortres201779-f4){#fig4}

![The analysis of MAPKK gene evolution. (**a**) Comparison of numbers of MAPKK genes in the phylogenetic tree representing plant species. (**b**) Venn diagram shows the number of common gene families and genes in 10 plants. (**c**) Genetic distance among the different classes of MAPKK genes.](hortres201779-f5){#fig5}

![The analysis of MAPKKK genes evolution. (**a**) Phylogenetic relationships among MAPKKK genes. (**b**) The percentage and numbers of MAPKKK genes in representative species. (**c**) Genetic distance among the different clades of MAPKKK genes. (**d**) Distribution of structural characteristics of MAPKKK genes in plant kingdom.](hortres201779-f6){#fig6}

![Expression pattern analysis of MAPK cascade genes in *A. thaliana* and *B. rapa.* Heat map representation and hierarchical clustering of MAPK cascade genes ((**a**) MAPK, (**b**) MAPKK, (**c**) MAPKKK) in root, stem, leaf, flower and silique. Venn diagram depicting the distribution of shared expression of the MAPK cascade genes among five *Brassica rapa* tissues.](hortres201779-f7){#fig7}

![The *cis*-elements of the promoter regions of *B. rapa* MAPK cascade genes. The micro-segments in different colors were the putative elements sequences. The description of the ten *cis*- elements are in brackets.](hortres201779-f8){#fig8}

![Evolutionary history of the MAPK cascade genes in the plant kingdom.](hortres201779-f9){#fig9}

[^1]: PW conceived and designed the experimental design. YL and WW contributed to the experimental work. XH contributed to the interpretation of the results and coordinated the study. All the authors read and approved the final manuscript.
