Expression Profiling of Mitogen-Activated Protein Kinase Genes Reveals Their Evolutionary and Functional Diversity in Different Rubber Tree (Hevea brasiliensis) Cultivars

Rubber tree (Hevea brasiliensis) is the only commercially cultivated plant for producing natural rubber, one of the most essential industrial raw materials. Knowledge of the evolutionary and functional characteristics of kinases in H. brasiliensis is limited because of the long growth period and lack of well annotated genome information. Here, we reported mitogen-activated protein kinases in H. brasiliensis (HbMPKs) by manually checking and correcting the rubber tree genome. Of the 20 identified HbMPKs, four members were validated by proteomic data. Protein motif and phylogenetic analyses classified these members into four known groups comprising Thr-Glu-Tyr (TEY) and Thr-Asp-Tyr (TDY) domains, respectively. Evolutionary and syntenic analyses suggested four duplication events: HbMPK3/HbMPK6, HbMPK8/HbMPK9/HbMPK15, HbMPK10/HbMPK12 and HbMPK11/HbMPK16/HbMPK19. Expression profiling of the identified HbMPKs in roots, stems, leaves and latex obtained from three cultivars with different latex yield ability revealed tissue- and variety-expression specificity of HbMPK paralogues. Gene expression patterns under osmotic, oxidative, salt and cold stresses, combined with cis-element distribution analyses, indicated different regulation patterns of HbMPK paralogues. Further, Ka/Ks and Tajima analyses suggested an accelerated evolutionary rate in paralogues HbMPK10/12. These results revealed HbMPKs have diverse functions in natural rubber biosynthesis, and highlighted the potential possibility of using MPKs to improve stress tolerance in future rubber tree breeding.


Introduction
The rubber tree (Hevea brasiliensis) is the most important commercially cultivated natural rubber-bearing plant. It is a member of the spurge family (Euphorbiaceae) and is closely related to cassava (Manihot esculenta) and castor bean (Ricinus communis). Natural rubber (cis-1,4-polyisoprene) is strategically important, and no other synthetic alternatives possess the exact same physicochemical properties. Therefore, rubber tree has substantial economic value as the only commercially cultivated

Identification and Genomic Location of HbMPK Genes
Previously reported MPKs of Arabidopsis [15] and rice [16] were retrieved and subjected to a Basic Local Alignment Search Tool for protein query (BLASTP) against the Hevea genome database, in order to comprehensively annotate MPKs [4]. Candidate HbMPK gene sequences were then submitted to InterProScan [30] to assess the MPK conserved motif (IPR003527) and protein kinase domain (IPR000719). A careful manual review of the resultant HbMPK candidate genes was performed to correct any potential mistakes in the genome database.
The genomic locations of these identified HbMPK genes were further determined based on the results of a BLAST nucleotide (BLASTN) query against the Hevea genome sequence. MapInspect software was used to draw the locations of the HbMPK genes.

Phylogenetic Tree, Motif Distribution and Alignment of the Mitogen-Activated Protein Kinase Family
The amino acids that compose MPKs from Hevea, two model plants Arabidopsis and rice, and a closely related cassava species [31] were identified, after which a phylogenetic tree was constructed by MEGA5.0 [32] using neighbour-joining (NJ) method. Bootstrap tests were performed and consisted of 1000 replicates.
Multiple sequence alignment was performed using ClustalW [33]. The molecular weight and isoelectric points of predicted HbMPK proteins were estimated using the ExPASy proteomics server [34]. The subcellular localization of HbMPKs was predicted using Softberry [35]. Conserved motifs of the HbMPK proteins were analysed using the Multiple Expectation Maximization for Motif Elicitation (MEME) program [36].

Gene Structure, Duplication Events and Syntenic Analyses
Gene Structure Display Server (GSDS) software [37] was used to analyse the exon-intron distribution based on the comparison of the open reading frame (ORF) sequence and genomic coding region of each HbMPK gene.
Gene duplication events were determined by multiple sequence alignment and collinear analyses. Homologous HbMPK genes were considered as paralogues only when their nucleotide identities were >90% [38,39]. The Ka (nonsynonymous substitution rate) and Ks (synonymous substitution rate) were calculated using DnaSP 5.0 software [40]. The Ka/Ks ratios for the HbMPK genes were calculated to assess the selection pressure on duplicated genes; a Ka/Ks ratio >1, <1, or =1 indicates positive, negative, or neutral evolution, respectively [41]. Tajima relative rate tests [42] were performed using the amino acid sequences of the duplicated HbMPK pairs. The syntenic relationships of paralogues and/or orthologues among Arabidopsis, cassava, and rubber tree were analysed using the Circos program [43].

Plant Materials and Abiotic Stress Treatments
Three rubber tree cultivated varieties, Reyan8-79 (high yielding), Reyan7-33-97 (medium yielding), and PR107 (low yielding), were used in this study. Different tissues (roots, stems, leaves and latex) were collected from the ten-year-old mature trees, and immediately frozen in liquid nitrogen. One-year-old Reyan7-33-97 seedlings were used for abiotic stress treatments. For osmotic and salt stress treatments, the seedlings were irrigated with 20% polyethylene glycol (PEG) 6000 and 1 M NaCl, respectively. For oxidative stress, the leaves were sprayed with 2% (v/v) H 2 O 2 . For cold stress, the seedlings were placed in a growth chamber at 4 • C. For RNA extraction, the leaves were collected at 3 h, 12 h, and 24 h after stress treatment. Untreated rubber tree leaves were used as 0 h control.

RNA Extraction and Expression Profiling of HbMPKs
The total RNA from different tissues (roots, stems, leaves and latex) of three cultivated varieties and the leaves of Reyan7-33-97 following abiotic stress treatments were extracted using an RNAprep Pure Plant Kit (Tiangen, Beijing, China). First-strand complementary DNA (cDNA) was generated using a reverse transcription system (Takara, Kusatsu, Japan) in accordance with the manufacturer's instruction.
The qRT-PCR was conducted using the gene specific primers listed in Supplementary Table S1. The reactions were performed using a Stratagene Mx3005P Real-Time Thermal Cycler (Agilent, Santa Clara, CA, USA) and SYBR Green Master Mix Reagent (Takara, Kusatsu, Japan), in accordance with previously described PCR conditions [44]. The specificity of the reaction was verified by melting curve analysis. The 2 −∆∆CT method was used to calculate the relative expression level of the target genes. The HbActin (GenBank Acc. HQ260674.1) served as a reference gene. Three independent biological replications were performed for each gene. Heat maps were generated using Multi Experiment Viewer (MeV) software as described [45]. Semi-quantitative PCR assays were performed using the same cDNA templates and primers used for qRT-PCR. The template concentrates were normalised using the HbActin gene as a reference. The expression levels of each HbMPK gene were visualised using agarose gel electrophoresis of the corresponding PCR products.

Cis-Element Distributions in HbMPK Promoter Regions
To investigate the cis-regulatory elements in the promotor regions of HbMPKs, a 1500 bp region of the genomic sequence up stream of the start codon of each gene was retrieved from the Hevea genome database. These sequences were then submitted to online PlantCARE software [46] in order to predict putative cis-elements.

Identification and Characterization of MPK Members in Hevea
A total of 45 MPK candidates were identified by submitting the A. thaliana MPK (AtMPK) and O. sativa MPK (OsMPK) sequences to BLASTP queries against the Hevea genome database. Furthermore, 20 candidates that had both a protein kinase domain and a conserved MPK active site were determined to be HbMPKs after sequence analysis ( Table 1). The 20 putative HbMPKs were carefully reviewed, after which the coding sequences (CDS) of HbMPK15 (scaffold0949_294805) were manually corrected by investigating the downstream genomic sequences. The ORF lengths of the 20 HbMPKs ranged from 894 bp (HbMPK17) to 1848 bp (HbMPK12) and their encoded proteins consisting of 298 to 616 amino acids. The calculated molecular weights of these proteins ranged from 34.85 to 69.95 kDa, and their isoelectric points ranged from 4.81 to 9.28. These proteins were predicted to localise in the nucleus (16 members) or cytoplasm (4 members).
Previously obtained isobaric tag for relative and absolute quantitation (iTRAQ) based proteomic data of the Reyan7-33-97 rubber tree latex [47] were submitted to the Hevea genome database to determine the expression patterns of HbMPK proteins. As a result, 4 out of the 20 HbMPKs were validated by high-throughput tandem mass spectrometry (MS/MS) data using a threshold of 95% confident peptides ≥2; two of them are induced by ethylene stimulation (Table 2).  The raw spectra of isobaric tag reporters ( Figure 1) and identified peptides that mapped to the HbMPKs are provided (Supplementary Figure S1). These data validated the above 20 putative HbMPKs at the proteomic level. Previously obtained isobaric tag for relative and absolute quantitation (iTRAQ) based proteomic data of the Reyan7-33-97 rubber tree latex [47] were submitted to the Hevea genome database to determine the expression patterns of HbMPK proteins. As a result, 4 out of the 20 HbMPKs were validated by high-throughput tandem mass spectrometry (MS/MS) data using a threshold of 95% confident peptides ≥2; two of them are induced by ethylene stimulation (Table 2). The raw spectra of isobaric tag reporters ( Figure 1) and identified peptides that mapped to the HbMPKs are provided (Supplementary Figure S1). These data validated the above 20 putative HbMPKs at the proteomic level. Figure 1. Proteomic validation of Hevea brasiliensis mitogen-activated protein kinases (HbMPK) expression patterns in rubber tree latex before and after ethephon treatment. Representative mass spectra indicate the signal intensities of isobaric tags for the validated MPKs. The spectra correspond to the peptides of "APELCGSFFSK", "DYVHQLR", "DLKPSNLLLNANCDLK" and "SFDFEQNALTEEQMK" for HbMPK12, HbMPK14, HbMPK16 and HbMPK19, respectively.
Genomic location analysis showed that HbMPK family members are distributed onto 20 individual scaffolds (Figure 2). The rubber tree was considered to have undergone two rounds of genome duplication events [4]. The expansion of the MPK gene family in Hevea was further investigated by analysing duplication events in MPK genes. Four duplication events involving 10 paralogues were identified (HbMPK3/HbMPK6, HbMPK8/HbMPK9/HbMPK15, HbMPK10/HbMPK12, and HbMPK11/HbMPK16/HbMPK19), which indicates that segmental duplication events play significant roles in MPK gene expansion in the Hevea genome. Genomic location analysis showed that HbMPK family members are distributed onto 20 individual scaffolds ( Figure 2). The rubber tree was considered to have undergone two rounds of genome duplication events [4]. The expansion of the MPK gene family in Hevea was further investigated by analysing duplication events in MPK genes. Four duplication events involving 10 paralogues were identified (HbMPK3/HbMPK6, HbMPK8/HbMPK9/HbMPK15, HbMPK10/HbMPK12, and HbMPK11/HbMPK16/HbMPK19), which indicates that segmental duplication events play significant roles in MPK gene expansion in the Hevea genome.

Multiple Sequence Alignment and Motif Distribution of HbMPKs
To investigate the sequence conservation of HbMPK proteins, multiple sequence alignment and conserved motif distribution analyses were performed. The results of the multiple sequence alignment showed that all identified HbMPKs contain the previously reported MPK-specific conserved domains: an N-terminal ATP-binding domain (P-loop), a phosphorylation catalytic loop (C-loop), and an activation loop (TEY or TDY). Notably, the HbMPKs that contained the TEY-type ( Figure 3A) activation loop also have an additional evolutionarily conserved C-terminal common docking domain (CD-domain), whereas the TDY-type ( Figure 3B) subfamily proteins do not, indicating a different evolutionary origin and functional diversity among these proteins.
Furthermore, 12 conserved motifs were recognised by the MEME program [36]. The distributions of the 12 conserved motifs for all 20 HbMPKs were highlighted ( Figure 4A). The MPKs were phylogenetically ordered to better understand their evolutionary relationships. Motifs 2, 5, 6 and 9 are conserved in each HbMPK and are organised in the exact same order. Motif 8 is specific to the MPKs of groups A and B, which have the conserved CD-domain. The conserved amino acid sequence information of these 12 motifs was also shown ( Figure 4B). Motifs 1, 3, 6, and 8 correspond to the conserved C-loop, P-loop, activation loop, and CD-domain, and were identified in the multiple sequence alignment (black boxes in Figure 4B).

Multiple Sequence Alignment and Motif Distribution of HbMPKs
To investigate the sequence conservation of HbMPK proteins, multiple sequence alignment and conserved motif distribution analyses were performed. The results of the multiple sequence alignment showed that all identified HbMPKs contain the previously reported MPK-specific conserved domains: an N-terminal ATP-binding domain (P-loop), a phosphorylation catalytic loop (C-loop), and an activation loop (TEY or TDY). Notably, the HbMPKs that contained the TEY-type ( Figure 3A) activation loop also have an additional evolutionarily conserved C-terminal common docking domain (CD-domain), whereas the TDY-type ( Figure 3B) subfamily proteins do not, indicating a different evolutionary origin and functional diversity among these proteins.
Furthermore, 12 conserved motifs were recognised by the MEME program [36]. The distributions of the 12 conserved motifs for all 20 HbMPKs were highlighted ( Figure 4A). The MPKs were phylogenetically ordered to better understand their evolutionary relationships. Motifs 2, 5, 6 and 9 are conserved in each HbMPK and are organised in the exact same order. Motif 8 is specific to the MPKs of groups A and B, which have the conserved CD-domain. The conserved amino acid sequence information of these 12 motifs was also shown ( Figure 4B). Motifs 1, 3, 6, and 8 correspond to the conserved C-loop, P-loop, activation loop, and CD-domain, and were identified in the multiple sequence alignment (black boxes in Figure 4B).

Evolution and Exon-Intron Organization of HbMPKs
To investigate the evolutionary relationships between HbMPKs and the MPKs from other plants, a phylogenetic tree was constructed using the protein sequences of 17 MPKs from the monocot model plant rice (O. sativa), 20 MPKs from the dicot model plant A. thaliana, 16 MPKs from the rubber tree relative cassava (M. esculenta), and 20 from Hevea in this study. These MPKs were divided into four groups according to phylogenetic analysis: TEY-type MPKs compose groups A, B and C, and TDY-type MPKs compose group D ( Figure 5). No species-specific clades were identified; however, groups A and B contain many more dicot MPK members than monocot MPK members. This finding indicates that duplication events may have occurred after the divergence of monocots and dicots. At the same time, group D contains 9 OsMPKs, 6 AtMPKs, 6 M. esculenta MPKs (MeMPKs), and 8 HbMPKs, which suggests multiple evolutionary and functional diversities for these MPK members. An A. thaliana-specific duplication event was identified in group C, and two

Evolution and Exon-Intron Organization of HbMPKs
To investigate the evolutionary relationships between HbMPKs and the MPKs from other plants, a phylogenetic tree was constructed using the protein sequences of 17 MPKs from the monocot model plant rice (O. sativa), 20 MPKs from the dicot model plant A. thaliana, 16 MPKs from the rubber tree relative cassava (M. esculenta), and 20 from Hevea in this study. These MPKs were divided into four groups according to phylogenetic analysis: TEY-type MPKs compose groups A, B and C, and TDY-type MPKs compose group D ( Figure 5). No species-specific clades were identified; however, groups A and B contain many more dicot MPK members than monocot MPK members. This finding indicates that duplication events may have occurred after the divergence of monocots and dicots. At the same time, group D contains 9 OsMPKs, 6 AtMPKs, 6 M. esculenta MPKs (MeMPKs), and 8 HbMPKs, which suggests multiple evolutionary and functional diversities for these MPK members. An A. thaliana-specific duplication event was identified in group C, and two potential Hevea-specific duplication events were identified in groups B (HbMPK11, 16 and 19) and D (HbMPK8, 9 and 15) by phylogenetic analysis (Figure 5). potential Hevea-specific duplication events were identified in groups B (HbMPK11, 16 and 19) and D (HbMPK8, 9 and 15) by phylogenetic analysis ( Figure 5). Gene structure divergence plays important roles in the evolution of gene families and can be used to assess phylogenetic relationships [48]. Analyses of the phylogenetic relationships and exon-intron distributions of both AtMPKs and HbMPKs revealed that most members in the same group shared very similar exon-intron structures ( Figure 6). This similarity provides evidence of MPK classification and evolutionary relationships between Arabidopsis and Hevea. Group A and B MPKs contain 4-7 exons of similar length, and 3-6 introns of variable size, whereas group C members have 2 parallel exons and only 1 intron. Group D shows a complex distribution of exons and introns. We also observed that the gene structures are conserved between HbMPK and AtMPK homologues. Additionally, similar gene structures were observed not only among paralogues, but also among homologues, which indicates the possible origin of these MPKs and evolutionary duplication events (AtMPK6 compared with HbMPK11/HbMPK16/HbMPK19, AtMPK16 compared with HbMPK8/HbMPK9/HbMPK15, and AtMPK20 compared with HbMPK10/HbMPK12). Moreover, we determined the Hevea-specific duplication relationships in the abovementioned paralogues by comparing the gene structure of MPKs to those of the related spurge family plant M. esculenta, and the result indicates that these duplicated MPKs are potentially involved in rubber tree-specific Gene structure divergence plays important roles in the evolution of gene families and can be used to assess phylogenetic relationships [48]. Analyses of the phylogenetic relationships and exon-intron distributions of both AtMPKs and HbMPKs revealed that most members in the same group shared very similar exon-intron structures ( Figure 6). This similarity provides evidence of MPK classification and evolutionary relationships between Arabidopsis and Hevea. Group A and B MPKs contain 4-7 exons of similar length, and 3-6 introns of variable size, whereas group C members have 2 parallel exons and only 1 intron. Group D shows a complex distribution of exons and introns. We also observed that the gene structures are conserved between HbMPK and AtMPK homologues. Additionally, similar gene structures were observed not only among paralogues, but also among homologues, which indicates the possible origin of these MPKs and evolutionary duplication events (AtMPK6 compared with HbMPK11/HbMPK16/HbMPK19, AtMPK16 compared with HbMPK8/HbMPK9/HbMPK15, and AtMPK20 compared with HbMPK10/HbMPK12). Moreover, we determined the Hevea-specific duplication relationships in the abovementioned paralogues by comparing the gene structure of MPKs to those of the related spurge family plant M. esculenta, and the result indicates that these duplicated MPKs are potentially involved in rubber tree-specific tissues and development processes (Supplementary Figure S2). However, the last duplication events, involving HbMPK3/HbMPK6, were observed in all of the investigated plant species investigated in this study (AtMPK1/AtMPK2 and MeMPK1/MeMPK2). We are interested in the potential functional diversity of Hevea-specific duplicated MPKs. events, involving HbMPK3/HbMPK6, were observed in all of the investigated plant species investigated in this study (AtMPK1/AtMPK2 and MeMPK1/MeMPK2). We are interested in the potential functional diversity of Hevea-specific duplicated MPKs.
To further investigate the expansion of the MPK gene family in Hevea, syntenic analysis of Hevea, M. esculenta, and A. thaliana MPKs was performed. The results were visualised using Circos software. We identified eight segmental duplication pairs in Hevea MPKs, but no tandem duplication events were detected, as our results were partly limited by the lack of chromosomal assembly information for rubber tree genome (Figure 7). Again, the duplication events described above were determined to be syntenic genes: HbMPK3/HbMPK6-MeMPK2, HbMPK10/HbMPK12-MeMPK20, HbMPK11/HbMPK16/HbMPK19-MeMPK6, and HbMPK8/HbMPK9/HbMPK15-MeMPK16-2. Modes of evolutionary selection can be estimated by the Ka/Ks ratio [41]. A Ka/Ks ratio > 1 indicates a positive selection, a Ka/Ks ratio < 1 indicates a purifying selection, and a Ka/Ks ratio = 1 indicates a neutral selection. The Ka/Ks ratios of the duplicated HbMPKs indicated that they all were subjected to purifying selection (Table 3). Furthermore, Tajima relative rate tests were conducted to investigate whether the HbMPK duplicates evolved at an accelerated rate following the duplication events. A statistically significant increase in evolutionary rate occurred between the HbMPK10/HbMPK12 duplicated pairs (Table 4), which indicates a potential functional divergence of these duplicated paralogues. To further investigate the expansion of the MPK gene family in Hevea, syntenic analysis of Hevea, M. esculenta, and A. thaliana MPKs was performed. The results were visualised using Circos software. We identified eight segmental duplication pairs in Hevea MPKs, but no tandem duplication events were detected, as our results were partly limited by the lack of chromosomal assembly information for rubber tree genome (Figure 7). Again, the duplication events described above were determined to be syntenic genes: HbMPK3/HbMPK6-MeMPK2, HbMPK10/HbMPK12-MeMPK20, HbMPK11/HbMPK16/HbMPK19-MeMPK6, and HbMPK8/HbMPK9/HbMPK15-MeMPK16-2.
Modes of evolutionary selection can be estimated by the Ka/Ks ratio [41]. A Ka/Ks ratio > 1 indicates a positive selection, a Ka/Ks ratio < 1 indicates a purifying selection, and a Ka/Ks ratio = 1 indicates a neutral selection. The Ka/Ks ratios of the duplicated HbMPKs indicated that they all were subjected to purifying selection (Table 3). Furthermore, Tajima relative rate tests were conducted to investigate whether the HbMPK duplicates evolved at an accelerated rate following the duplication events. A statistically significant increase in evolutionary rate occurred between the HbMPK10/HbMPK12 duplicated pairs (Table 4), which indicates a potential functional divergence of these duplicated paralogues.   a The Tajima relative rate test was used to examine the equality of the evolutionary rate between rubber tree paralogues; b Mt is the sum of the identical sites in all three sequences tested; c M1 is the number of unique differences in the first paralogue; d M2 is the number of unique differences in the second paralogue; e If p < 0.05, the test rejects the equal substitution rates between the two duplicates and infers that one of the two duplicates has an accelerated evolutionary rate.   paralogues; b Mt is the sum of the identical sites in all three sequences tested; c M1 is the number of unique differences in the first paralogue; d M2 is the number of unique differences in the second paralogue; e If p < 0.05, the test rejects the equal substitution rates between the two duplicates and infers that one of the two duplicates has an accelerated evolutionary rate.

Expression Profiling of MPKs in Different Tissues from Three Cultivated Varieties of Rubber Tree
To gain insights on the tissue-specific expression patterns of the 20 HbMPKs, the expression patterns of HbMPKs were visualised as a heat map and separated into three clusters ( Figure 8A). Members of cluster TI, including HbMPK11/19 (these paralogues could not be distinguished by gene-specific primers), HbMPK13, HbMPK6, HbMPK16, and HbMPK20, showed the highest expression level in latex, and the members in cluster TI were distributed in all groups. Members of cluster TII showed universal expression patterns in all four tissues; the highest expression levels usually occurred in photosynthetic tissues. Most members in cluster TII were distributed in Group D, only HbMPK3 was distributed in Group A. The remaining HbMPKs were barely expressed. The cluster TIII did not have the members in Group A. The semi-quantitative PCR results were also shown by the electrophoresis analysis of the PCR products ( Figure 8B). Notably, the MPKs in cluster TI (HbMPK11/19, HbMPK13,  HbMPK6, HbMPK16, and HbMPK20), which are mostly expressed in the latex, also demonstrated different expression patterns among different rubber tree varieties, indicating that these genes might be involved in latex production processes in Hevea.

Expression Profiling of MPKs in Different Tissues from Three Cultivated Varieties of Rubber Tree
To gain insights on the tissue-specific expression patterns of the 20 HbMPKs, the expression patterns of HbMPKs were visualised as a heat map and separated into three clusters ( Figure 8A). Members of cluster TI, including HbMPK11/19 (these paralogues could not be distinguished by gene-specific primers) , HbMPK13, HbMPK6, HbMPK16, and HbMPK20, showed the highest expression level in latex, and the members in cluster TI were distributed in all groups. Members of cluster TII showed universal expression patterns in all four tissues; the highest expression levels usually occurred in photosynthetic tissues. Most members in cluster TII were distributed in Group D, only HbMPK3 was distributed in Group A. The remaining HbMPKs were barely expressed. The cluster TIII did not have the members in Group A. The semi-quantitative PCR results were also shown by the electrophoresis analysis of the PCR products ( Figure 8B). Notably, the MPKs in cluster TI (HbMPK11/19, HbMPK13, HbMPK6, HbMPK16, and HbMPK20), which are mostly expressed in the latex, also demonstrated different expression patterns among different rubber tree varieties, indicating that these genes might be involved in latex production processes in Hevea.

Expression Profiling of HbMPKs in Reyan7-33-97 under Different Abiotic Stresses
The genes of the MPK family are crucial abiotic stress responsive genes [49]. The expression patterns of HbMPKs under different abiotic stresses were analysed to further understand the functional diversity of HbMPK members. As shown in Figure 8, the 20 HbMPKs were grouped into three clusters according to their response to different abiotic stresses. The MPKs in cluster AI, especially HbMPK6, showed a weak or slightly reduced response to the abiotic stress treatments. Cluster AII included MPK family members that were barely expressed in leaves, either before or after abiotic stress. Cluster AIII contained the most HbMPK members, and showed different degrees of increased expression levels after abiotic stress ( Figure 8C,D). Most members in cluster AIII were distributed in Group D, which means the members in Group D may play roles in diverse abiotic stresses.

Expression Profiling of HbMPKs in Reyan7-33-97 under Different Abiotic Stresses
The genes of the MPK family are crucial abiotic stress responsive genes [49]. The expression patterns of HbMPKs under different abiotic stresses were analysed to further understand the functional diversity of HbMPK members. As shown in Figure 8, the 20 HbMPKs were grouped into three clusters according to their response to different abiotic stresses. The MPKs in cluster AI, especially HbMPK6, showed a weak or slightly reduced response to the abiotic stress treatments. Cluster AII included MPK family members that were barely expressed in leaves, either before or after abiotic stress. Cluster AIII contained the most HbMPK members, and showed different degrees of increased expression levels after abiotic stress ( Figure 8C,D). Most members in cluster AIII were distributed in Group D, which means the members in Group D may play roles in diverse abiotic stresses.
These expression pattern analyses of HbMPKs in different tissues and under different abiotic stresses suggested that MPK genes might have multiple functions in Hevea to adapt to various environmental changes (Figure 8). In addition, the duplicated paralogues showed very different expression profiles, indicating functional and evolutionary divergence emerged following the duplication events.

Analysis of cis-Elements in the Promotor Regions of MPKs in Hevea
The distribution of cis-elements in the promoter regions of HbMPKs was further investigated to demonstrate the regulatory patterns of the MPK family members, especially duplicated paralogues ( Figure 9). To facilitate this understanding, cis-elements were symbolised by capital letters in different colour; detailed classification and sequence information are listed (Supplementary Table S2). No significant positive correlations were observed between the distribution patterns of cis-elements and the expression patterns of HbMPKs, partly because of the lack of promoter sequence information for several genes. However, we did observe different cis-element organisational patterns in paralogues that exhibited different expression patterns in various tissues and/or under different abiotic stresses (HbMPK3-HbMPK6, HbMPK10-HbMPK12, and HbMPK11/19-HbMPK16). These results suggest the existence of different regulatory patterns and potential evolutionary fates in these paralogues.  (Figure 8). In addition, the duplicated paralogues showed very different expression profiles, indicating functional and evolutionary divergence emerged following the duplication events.

Analysis of cis-Elements in the Promotor Regions of MPKs in Hevea
The distribution of cis-elements in the promoter regions of HbMPKs was further investigated to demonstrate the regulatory patterns of the MPK family members, especially duplicated paralogues ( Figure 9). To facilitate this understanding, cis-elements were symbolised by capital letters in different colour; detailed classification and sequence information are listed (Supplementary Table S2). No significant positive correlations were observed between the distribution patterns of cis-elements and the expression patterns of HbMPKs, partly because of the lack of promoter sequence information for several genes. However, we did observe different cis-element organisational patterns in paralogues that exhibited different expression patterns in various tissues and/or under different abiotic stresses (HbMPK3-HbMPK6, HbMPK10-HbMPK12, and HbMPK11/19-HbMPK16). These results suggest the existence of different regulatory patterns and potential evolutionary fates in these paralogues.   Table S2).

Identification and Characteristics of HbMPK Genes and TheirAssociated Proteins
MPKs in various higher plants, including Arabidopsis, rice, maize, tobacco, poplar, and cassava, etc., have been systematically investigated [15][16][17][18][19][20][21][22][23][24]. However, for the first time, we verified 20 MPK genes in the most important natural rubber-bearing plant. Rubber tree MPK family members share three conserved domains and could be separated into four known groups; similar results were observed in other plants. Three of these groups have a TEY-type activation loop, and 1 has a TDY-type activation loop [14]; the TEY-type MPKs also contain additional C-terminal CD-domains (Figures 3-5). The genome of rubber tree is believed to have undergone two rounds of duplication events, resulting in similar numbers and characteristics of HbMPKs compared with those of Arabidopsis and M. esculenta.
To further verify the identified MPKs in rubber tree, previous iTRAQ data [47] obtained from rubber tree latex were introduced (Table 2; Figure 1). This is the first time that proteomic data were combined with genome-wide gene family studies for rubber-bearing plant research. The results support the gene family identification data and provide insights into the proteome patterns of MPK members in Hevea. Similar to genome research projects, proteomic data should be increasingly introduced in future gene-wide gene family studies. Although the proteomic data shown here did not provide the expression information about MPK proteins under abiotic stresses, we still believe that it is necessary to add proteomic evidence in the gene family investigations.

Evolutionary Divergence of Rubber Tree MPK Genes
The expansion of a gene family within a genome mainly occurs by three kinds of duplication events: tandem duplication of individual genes, segmental duplication of multiple genes, and background duplications that cannot easily be classified [50]. Given the limited chromosome assembly information of Hevea, the present study revealed that there is no or little contribution from tandem duplication to the expansion of the rubber tree MPK family (Figures 2 and 7).
One important purpose of the evolutionary study of gene families is to determine whether there are species-specific family members. This information can aid in understanding the functional divergence of a specific gene family [51]. Although no rubber tree-specific groups were determined by phylogenetic, gene structure, or gene syntenic analyses of the HbMPK family, we identified rubber tree-specific duplication events (Figures 5-7). Three Hevea-specific duplication events could be determined, but they were not observed for Arabidopsis or the spurge family plant cassava. These duplication events are HbMPK11/HbMPK16/HbMPK19-AtMPK6, HbMPK8/HbMPK9/HbMPK15-AtMPK16, and HbMPK10/HbMPK12-AtMPK20. The remaining paralogues demonstrate no species specificities, and HbMPK3/HbMPK6 duplication was observed in all of the investigated plants. These results indicate that all these duplication events occurred at different time points of evolution ( Figure 6 and Supplementary Figure S2).
Functional diversity caused by gene duplication might result in altered expression profiles and protein properties, and gene duplication is a major evolutionary driver for increasing the fitness of plants to new environment [52]. Furthermore, these paralogous pairs show tissue-and variety-specific expression patterns (Figure 8) as well as abiotic stress-specific patterns (Figure 8). One of the duplicated genes of three paralogous pairs shows global expression patterns in most tissues (HbMPK6, HbMPK8/15, HbMPK12, HbMPK11/19, and HbMPK16), whereas the other paralogue copies are barely expressed in root, stem, leaf, and latex tissues (HbMPK3, HbMPK9, and HbMPK10). Notably, the paraloguesHbMPK11/HbMPK16/HbMPK19 show universal expression patterns in different tissues (roots, stems, leaves, and latex) and under different abiotic stresses (osmotic, oxidative, salt and cold); however, HbMPK11/19 are mainly expressed in rubber latex (Figure 8). These results suggest that HbMPK16 might play similar role as its AtMPK6 homologue in Arabidopsis (which is universally expressed in various tissues, and is strongly associated with numerous abiotic stresses) [53] and that HbMPK11/19 somehow acquired new functions during the evolution of Hevea, and thus, are involved in nature rubber biosynthesis.
The selection pressure on HbMPK paralogues was estimated by calculating the ratio of Ka/Ks of rubber tree MPK paralogues. All Ka/Ks ratios were less than 0.1, which suggests that these paralogues were subjected to purifying selection, and mutations that alter amino acid sequences were not accepted during evolution (Table 3). These results are consistent with MPK family members playing very important roles during plant development, which requires highly conserved amino acid sequences to maintain protein functionality. Furthermore, the Tajima relative rates of Hevea gene pairs were calculated to evaluate the evolutionary rate between paralogues (Table 4). Notably, the duplicated gene pair HbMPK10/HbMPK12 has a p-value of 0.00555, which means that one of the two duplicates has a significant accelerated evolutionary rate. After combining and analysing the different expression patterns between HbMPK10 and HbMPK12 both in various tissues and under abiotic stresses, we assumed that the two genes are undergoing an accelerated separation in Euphorbiaceae plants, which suggests that they potentially play specific roles in Hevea.

Expression Characteristics of HbMPKs in Various Tissues and Cultivated Varieties under Different Stresses
MPK genes in groups A and B are more like the mainly expressed MPKs [54]. We observed different expression patterns of MPKs in vegetative tissues (roots, stems, and leaves) and in specialised rubber tree tissue (latex); these MPKs mainly belong to groups A and B.
Five MPK genes (HbMPK6, 11/19, 13, 16 and 20) were expressed in latex, according to the qRT-PCR and semi-quantitative PCR results ( Figure 8). Notably, HbMPK11/19, HbMPK13, and HbMPK6 showed significantly different expression levels across the high (Reyan8-79), medium (Reyan7-33-97), and low (PR107) yielding varieties of rubber tree. In particular, the expression levels of HbMPK11/19 and HbMPK13 in latex showed similar change patterns across the cultivated varieties, whereas HbMPK6 showed opposite patterns, which indicates that these genes may participate in the induction or repression of the biogenesis of natural rubber. The phosphorylation of functional proteins, including rubber elongation factor (REF), small rubber particle protein (SRPP), and osmotin, plays a very important role in natural rubber biosynthesis [55][56][57], and thus provides potential targets for HbMPK6, HbMPK13, and HbMPK11/19. In addition, two other MPKs, named HbMPK20 and HbMPK5, have shown very different expression patterns across the three cultivated varieties, the reason for which remains unclear.
The majority of the HbMPKs can respond to various abiotic stresses in this study ( Figure 8). The expression patterns under abiotic stress do not match the cis-element distributions in their promoter regions; however, the duplicated paralogues that have various cis-element distributions really exhibit completely different expression patterns both in different tissues and in different cultivated varieties under various abiotic stresses (Figures 8 and 9). This phenomenon supports the hypothesis that gene loss occurs very rapidly following whole-genome duplication, and that functional divergence can explain why duplicated genes escape extinction [58].

Conclusions and Prospects
In summary, based on our proteomic data and the released Hevea genome, a total of 20 Hevea MPK family members were characterised. Phylogenetic and syntenic analyses determined that 10 of these HbMPKs are paralogues. Evolutionary analysis showed that all the paralogues were subjected to purifying selection, and that only one duplicated gene pair (HbMPK10/HbMPK12) has an accelerated evolutionary rate. Furthermore, expression profiling of the 20 HbMPKs in the roots, stems, leaves, and latex of three rubber tree cultivated varieties showed very different expression patterns among the paralogues. The paralogues of HbMPK11/19 and HbMPK16 showed very different expression levels in the latex of various rubber tree varieties, which indicates functional divergence between these gene pairs. Significantly different cis-element distributions were observed in the promoter regions of the duplicated paralogues. Our data provided the first systematic and evolutionary aspects of the MPK gene family of the rubber-bearing plant Hevea, and suggest that members of duplicated gene pairs, such as HbMPK11/19 and HbMPK13, might play crucial roles in the biogenesis of natural rubber.
However, the limitation of transgenic technology in rubber tree and the lack of Hevea MPK-specific antibodies lead to a deficiency in functional validation in the present study. Moreover, additional investigations of several HbMPK paralogues (HbMPK8/HbMPK9/HbMPK15 and HbMPK10/HbMPK12) that showed interesting expression patterns should be conducted in the near future.