Genome-wide identification of MAPKKK genes and their response to phytoplasma stress in Chinese

to various biotic and abiotic stresses. In the previous study we have identified 10 ZjMAPKs and 5 ZjMAPKKs in Chinese jujube genome and found some crucial members of ZjMAPKs and ZjMAPKKs might function importantly in the process of phytoplasma infection. But how these ZjMAPKKs were modulated by ZjMAPKKKs during this process is still elusive and little information is known about the MAPKKKs in Chinese jujube.


Conclusions
The identification and classification analysis of ZjMAPKKKs was firstly reported and some key individual ZjMAPKKKs genes might play essential roles in response to phytoplasma infection. This could provide initial understanding for the mechanism that how the ZjMAPKKKs were involved in jujube -phytoplasma infection.

Background
Mitogen activated protein kinase (MAPK) cascades, consisting of three specific kinase families: MAP kinase kinase kinases (MAPKKKs), MAP kinase kinases (MAPKKs) and MAP kinases (MAPKs), are essential middle signal modules from signal sensing to the activation of related transcript factors in response to biotic and abiotic stresses such as drought, salt, cold and pathogen attack [1][2][3]. Among the signalling transduction, the conserved serine/threonine MAPKKKs could be activated by the plasma membrane receptors, then phosphorylate the MAPKKs which further activate the MAPKs by sequentially phosphorylation. Finally, the MAPKs could regulate other kinases or related transcription factors in response to various stresses [4,5]. Furthermore, each MAPK cascade family consists of a series of members and the number of them is significantly different. For instance, , the MAPKKKs family composes greater members and shows more complex sequence diversity comparing to the other two families. The members belong to this family could be classified into Raf, ZIK and MEKK subfamilies, according to characteristic sequence motifs [6]. The structural disparity of MAPKKKs in each subfamily could be observed, the Raf subfamily has a C-terminal kinase domain and a long Nterminal regulatory domain, while ZIK proteins only have the N-terminal kinase domains and MEKK subfamily has less conserved kinase domain. In addition, thethe long N-terminal regulatory domain is the backbone for Raf and ZIK subfamilies [1,6].
MAPK cascades have been implicated in the signal transduction in distinct plants innate immunity [7,8]. In Arabidopsis, the MEKK1-MKK4/5-MPK3/6-WRKY22/WRKY29/FRK1 cascade was involved in the innate immunity signalling transduction and the MEKK1-MKK1/MKK2-MPK4 kinase cascade could negatively activate MEKK2 which further leads to SUMM2-mediated immune responses [9,10]. In tobacco, the NPK1-MEK1-Ntf6 could regulate WRKY/MYB transcription factors to participate in the tobacco mosaic virus infection pathway [11]. In addition, the MAPKKKα-MKK2/MKK4-MPK2/MPK3 cascades took part in the Pto-mediated effect or triggered immunity (ETI) pathway by regulating the transcription factor TGA in tomato [12]. Hence, the MAPKs pathway is indeed involved in the pathogen attack and might also play essential roles in response to phytoplasma infection in Chinese jujube.
Hitherto, the MAPKKKs family has been characterized in plant kingdom. 80 MAPKKKs was firstly identified in Arabidopsis in 2002 [6,13]. From then on, an array of MAPKKKs has been identified from different plants, including rice (75 members), Zea mays (71 members), Vitis vinifera (45 members), Malus domestic (72 members), Musa nana (77 members), and so on [14][15][16][17]. However, little is known about the biological information and function of MAPKKK gene family in Chinese jujube, even though the detailed information of ZjMAPKKs and ZjMAPKs has been reported by our group [18].
Jujube witches' broom disease (JWB) caused by 'Candidatus Phytoplasma ziziphi' has become a devastating disease in Asia [19]. This disease shows three typical symptoms, including witches' broom, phyllody and yellowing [20]. The physiological and biochemical behaviors of jujube infected by phytoplasma has been widely studied [21], but the molecular mechanism behind is still elusive.
Recently, the MAPKs in response to phytoplasma infection have been reported with analysis of RNAseq transcriptomics in Chinese jujube [22], this study provides the insight that MAPKs must play important roles during this process. In addition, ZjMPK2, ZjMKK2 and ZjMKK4 have been demonstrated to be the main and key genes involved in Chinese jujube-phytoplasma interaction and ZjMKK2 could interact with ZjMPK2 with the yeast two hybrid analysis [18,23]. All these results demonstrate the important function of MAPK cascades in response to phytoplasma infection and it is necessary to do the identification and initial functional analysis of ZjMAPKKKs to build the complete MAPK cascade signalling transduction pathway. Thus, in this study, the ZjMAPKKKs were identified by the genome wide analysis and the phylogenetic analysis, gene structure and conserved motifs of ZjMAPKKKs were also predicted. Furthermore, the expression profiles of these ZjMAPKKKs in response to phytoplasma infection were investigated by qPCR. These investigation and identification of the ZjMAPKKKs and the detection of their initial functions to phytoplasma infection could facilitate the understanding of the mechanism how the ZjMAPKs cascades involved in JWB defense response.

Identification of ZjMAPKKKs in Chinese jujube
The ZjMAPKKKs gene family was identified according to our previous study on the identification of ZjMAPKKs and ZjMAPKs with some modifications [18]. Firstly, the whole protein sequences of MAPKKKs in Arabidopsis were retrieved from TAIR databases (Additional file 1). These sequences were as queries to search against in the whole jujube genome database (accession JREP00000000) [24]. In addition, the alignments of all Arabidopsis MAPKKK sequences were used to construct a HMM profile (http://hmmer.org/download.html) to search the other potential MAPKKKs members in jujube genome database. All the potential ZjMAPKKKs genes were confirmed by HMMER tools which contain protein kinase domain (PF00069) [25], while the redundancy was removed, and the remaining sequences were identified as ZjMAPKKKs family. Furthermore, the ExPASy Proteomics Server (http://expasy.org/) was used to calculate the theoretical pI (isoelectricpoint) and Mw (molecular weight) of the putative ZjMAPKKKs [26].

Gene structure, protein domain and motif analysis of ZjMAPKKKs genes
The open reading frames (ORF) of 56 ZjMAPKKKS were analyzed through the National Center for Biotechnology Information (NCBI, http://www.ncbi.nlm.nih.gov) ORF finder and the main protein domains of ZjMAPKKKs were recorded by blastp of NCBI as well (Additional file 2). GSDS (Gene Structure Display Server, http://gsds.cbi.pku.edu.cn/) was exploited to determine the exon/intron structures of individual ZjMAPKKKs by aligning the cDNA sequences to their corresponding genomic DNA sequences [27]. MEME database was exploited to identify the conserved motifs of ZjMAPKKKs [28]. The ideal motif widths were set to be between six and fifty, and each protein sequence of the ZjMAPKKKs subfamily was used to find the higher number of conserved domains [29].

Sequence alignment, phylogenetic and gene duplication analysis
All the protein sequences of ZjMAPKKKs were aligned by ClusterX software with default parameters.
Then the alignment of the protein sequences of 56 ZjMAPKKKs and 80 AtMAPKKKs were performed to do the phylogenetic analysis using MEGA 6.06 alignment explorer. The parameters of alignment as follows: gap opening penalty, 10.00; gap extension penalty, 0.50 (both in pairwise alignment and multiple alignment); protein weight matrix, gonnet; residue-specific penalties, on; hydrophilic penalties, on; gap separation distance, 0; end-gap separation, on; use negative matrix, off; and delay divergent cutoff (%), 30. Phylogenetic trees were constructed by the neighbor-joining (NJ) method.
The parameters of the constructed trees were: statistical method, neighbor-joining; scope, all selected taxa; test of phylogeny, bootstrap method; number of bootstrap replications, 1000; substitutions type, amino acid; model/method, poisson model; rates among sites, uniform rates; pattern among lineages, same (homogeneous); and gaps/missing data treatment: complete deletion. Multiple Collinearity Scan toolkit (MCScanX) was adopted to analyze the gene duplication events, with the default parameters and the homologous relationships were drawn using the Circos software [30,31].

Plant materials and treatments
Ziziphus jujuba Mill. 'Dongzao' which was cultivated in the Experimental Station of Chinese Jujube, Hebei Agricultural University, were used as the experimental materials. The leaf samples were collected from at least three healthy jujube trees and three trees infected with JWB for each time point. All these experimental trees were cultivated in the naturally environmental conditions [21]. The visual JWB symptoms, such as phyllody leaves (floral organs becoming leaf-like) and apparent normal leaves (infected but none symptom), were firstly observed earlier in June, and then witches' broom leaves (shoot with small leaves) were observed at the middle of Junein each year. Finally, the mature leaves of witches' broom leaves (shoot with small leaves), phyllody leaves (floral organs becoming leaf-like) and apparent normal leaves (infected but none symptom), were collected from the same branch of diseased trees in June, July, August, and September (on the 15th day of each month). The leaves from healthy trees were used as control for each time point. The different type materials were shown in Additional file 3. The detection of JWB phytoplasmas in the infected materials in June was detected by DAPI staining at histological level and quantitative real-time PCR analysis (qRT-PCR) at the molecular level from June to September [21]. Each material was sampled with three replicates from independent trees. In addition, the sterile cultivated JWB diseased plantlets of Ziziphus jujuba Mill. 'Goutouzao' were used as another parallel testing group and the healthy plantlets were used as control (As shown in Additional file 4). The detection of JWB phytoplasmas in the diseased and healthy sterile cultivated plantlets was detected by DAPI staining at histological level (Additional file 5) [23]. Every ten plantlets were pooled as one sample, three independent biological replications were sampled separately. All samples were rapidly frozen in liquid nitrogen and kept at -80 °C for RNA isolation and qPCR analysis.

RNA extraction and qRT-PCR analysis
Total RNA from the leaves was extracted by the TIANGEN RNA Extraction Kit. The genomic DNA contamination from the samples was removed by digesting with DNase. The cDNA synthesis was performed by TaKaRa RNA PCR Kit (AMV) Ver.3.0 (TaKaRa) according to the protocol of the manufacturer protocol using 1 µg of RNA template.
The qRT-PCR was carried out on the Bio-Rad iQ™ 5 using TransStart Top Green qPCR SuperMix AQ131 (TransGen Biotech, China). The 20 µL reaction system contained 10 µL of 2×SYBR Premix ExTaq™, 0.4µL each of 10 µM primers, 1µL diluted cDNA and 8.2 µL ddH 2 O. The thermal profile was preincubation for 3 min at 94 °C, followed by 40 cycles of 5 s at 94 °C, 15 s at 55~63 °C and 15 s at 72°C.Relative expression levels of ZjMAPKKKs were calculated with the 2 -ΔΔCt method [30] using ZjActin as endogenous control for normalization [31]. The primer sequences of ZjMAPKKKs for qPCR were shown in additional file 6. The gene expression with more than or less than two-fold changes was considered as significant difference [23].

Heatmap construction
The expression profiles of all ZjMAPKKKs in different samples were illustrated by a color gradient heatmap. The heatmap was constructed by heatmap software Heml 1.0 using Log2 based expression fold-changes.  Table 1). The ZjMAPKKKs distributed all over the 12 pseudo-chromosomes, excepted

Genome-wide identification of ZjMAPKKKs
ZjMAPKKK44-56 could not match to corresponding chromosome.
Specific information of each CDS and amino acid sequences of ZjMAPKKKs were listed in the Additional file 1 and 7. In addition, according to the specific conserved signature motif, all the ZjMAPKKKs could be divided into two subfamilies (Raf and MEKK); no ZIK subfamily members could be identified. Furthermore, as shown in Table 1, the length of the CDS sequence ranged from 762 bp (ZjMAPKKK36) to 4455 bp (ZjMAPKKK7), with an average length of 1804 bp. The amino acid sequence length of ZjMAPKKKs varied from 253 (ZjMAPKKK36) to 1484 (ZjMAPKKK7) amino acids (aa); average length was 600 aa. The predicted molecular weight (Mw) of these proteins ranged from 28.29 (ZjMAPKKK36) to 160.66 (ZjMAPKKK7) and the theoretical isoelectric point (pI) ranged from 4.78 to 9.34, respectively.

Phylogenetic analysis of ZjMAPKKKs genes
In order to assess the phylogenetic relationships between Chinese jujube and Arobidopsis, the phylogenetic tree was constructed with the total 136 protein sequences (56 ZjMAPKKKs and 80 AtMAPKKKs). As illustrated in Fig. 1, the members of AtMAPKKKs could be clustered into three categories, Raf, ZIK and MEKK, indicating that the method for the phylogenetic tree was reliable.
However, the 56 members of ZjMAPKKKs could only be clustered into two subfamilies, Raf and MEKK.
In addition, the largest Raf subfamily consisted of 41 members and the other 15 members of ZjMAPKKKs belonged to MEKK subfamily, and none of them could be defined into ZIK subfamily.
Morever, some ZjMAPKKKs located on the same chromosome, shown little divergence, but clustered into the same group, such as ZjMAPKKK36, 37 and 40, ZjMAPKKK15 and 16, ZjMAPKKK38 and 39.
These results indicated that some duplication of ZjMAPKKKs took place during the evolutionary process of jujube.

Conserved domains and gene structure analysis of ZjMAPKKKs
Within the analysis of MEME software, five main conserved motifs were identified in all 56 ZjMAPKKKs (Fig.2). The motif 1, 3 and 4 were founded in all the ZjMAPKKKs, while the other two motifs were observed in all the members of Raf subfamily. Additionally, the members of MEKK subfamily could be divided into two groups, one group contains motifs 1-4, including ZjMAPKKK21, 56, 6, 31, 10 and 25.
The remaining members only consisted of motifs 1, 3 and 4. These results illustrated that the ZjMAPKKKs share the same conserved motifs which further indicate that the protein structures for each subfamily were highly conserved.
For the analysis of the contents of exon/introns, the difference among ZjMAPKKKs was significant. As shown in Fig. 3 and Additional file 8, the number of exons in ZjMAPKKKs ranged from 1 (ZjMAPKKK9, 12,29,35,36,44 and 54) to 19 (ZjMAPKKK42). Interestingly, the members of

Synteny analysis of ZjMAPKKK genes
Furthermore, a tandem duplication event was firstly analyzed according to the principle that two or more genes located on a chromosomal region within 200 kb [34]. As shown in fig.5, One pair of ZjMAPKKKs (15/16) was the only tandem duplication event on LG5. In addition, 13 segmental duplication events with 22 ZjMAPKKKs were also identified. These results indicated that some ZjMAPKKKs were possibly generated by gene duplication and the segmental duplication events played a major driving force for ZjMAPKKKs evolution.

Phytoplasma detection in different phytoplasma infected tissues
In order to get insight to understand the function of ZjMAPKKKs involved in phytoplasma infection, the expression level of individual ZjMAPKKKs was detected by qPCR in two kind of plant materials one material was infected by phytoplasma with three different symptoms, including witches' broom leaves, phyllody leaves and apparent normal leaves from diseased plants (in vivo) and the other plant material was sterile cultivated tissues of JWB plantlets (in vitro). Firstly, the phytoplasma concentration in in the first plant material with three sympotems has been detected by Xue et al.
(2018) [21], and the phytoplasma determination in the sterile cultivated tissues of JWB plantlets shown that the fluorescent spots formed a large bright circle in the petiole phloem (Additional file 5). These results guaranteed the following test on the function of ZjMAPKKKs in response to phytoplasma infection.

Expression analysis of ZjMAPKKKs in witches' broom leaves
As shown in Additional file 9 and Fig. 6 (A), the heat map showed the expression levels of ZjMAPKKKs with significantly different patterns in witches' broom leaves from June to September. 42 candidates expressing level could be detectable, but the expression levels of the other 14 ZjMAPKKKs were very low and could not be detectable. Finally these undetectable expressing ZjMAPKKKs genes were considered as the redundant candidates in our research and were not selected for further calculation and analysis. Among them, the most significant transcript induction took place at the early stage (June or July) when the witches' broom began to grow. For instance, the ZjMAPKKK13, 14, 15,23,34,42,44,47 and 56 were highly induced in June or July (above 2-folds), afterwards, began to decrease from August to September as shown in Fig. 6 (B). However, the ZjMAPKKK3, 43 and 50 were down-expressed from June to September. In addition, two members of ZjMAPKKKs (26 and 45) remained high level expression constantly from June to September, indicating these two candidates might be the potentially key MAPKKKs in response to phytoplasma infection. However, it should be kept in mind that the clustering of the expression profiles of ZjMAPKKKs was not consistent with the gene similarities, illustrating the gene function did not only rely on the gene structure.

Expression analysis of ZjMAPKKKs in phyllody leaves
As depicted above, the transcript abundance of ZjMAPKKKs were as well as investigated in the tissue of phyllody leaves. The heat map of the expressing ZjMAPKKKs was indicated in Fig. 7 (A). Several of the ZjMAPKKKs were highly expressed in June or July and then the expression levels of them decreased from August to September, while most of the ZjMAPKKKs showed no significant changes or down regulated. For the details of each single ZjMAPKKKs expressing level could be seen in the Fig. 7 (B), the ZjMAPKKK10, 14, 15 34, 44 and 56 were significantly up regulated at the early stage (June or July). However ten of the ZjMAPKKKs showed significantly down regulated, including ZjMAPKKK3, 16,18,41,43,50,51,52,53 and 55. As same as the expressing pattern of ZjMAPKKK26 and 45 in the tissue of witches' broom leaves, these two members were also highly up regulated from June to September in the tissue of phyllody leaves.

Expression analysis of ZjMAPKKKs in apparent normal leaves
Furthermore, the apparent normal leaves that were infected by phytoplasma but shown no apparent symptoms was used to test which ZjMAPKKKs really play roles in response to phytoplasma infection in the non-strong disease leaves . Interestingly, the heat map figure showed different expressing pattern of ZjMAPKKKs in the tissue of apparent normal leaves (Fig. 8B). Seldom genes were highly up regulated and most of them showed down regulated pattern. For example, the ZjMAPKKK1, 3,7,16,17,18,19,41,43,50,51 and 53 were down regulated from June to September, while the ZjMAPKKK28, 34 and 47 were significantly up regulated in June or July, and the ZjMAPKKK27 and 54 were up regulated from August or September. However, ZjMAPKKK26 and 45 showed the same high expressed pattern in the tissue of apparent normal leaves from June to September (Fig. 8 A).
To sum up, in the above-mentioned three phytoplasma infected tissues from four different time series, ZjMAPKKK26 was significantly up regulated (8 times) and the ZjMAPKKK45 were highly induced 10 times. Moreover, with the infecting development of phytoplasmas, the symptoms of which observed on jujube leaves were gradually stronger from the apparent normal leaves, phyllody leaves to withches' broom leaves [21]. For the analysis of this aspect, the expression level of ZjMAPKKK26 was highly induced in phyllody leaves (~6 fold) in June, but not within other two symptom leaves.
Then with time and infecting development, ZjMAPKKK26 was up regulated in witches' broom leaves (~3 fold) in July and down regulated in phyllody leaves (~2 fold) in the meantime. However, ZjMAPKKK45 was highly induced in apparent normal leaves and phyllody leaves (~3 and ~6 fold, respectively) in June and began to be down regulated from July to September, but it was also induced in witches' broom leaves constantly (~2 fold) (Fig. 6, 7 and 8). These results demonstrated within the infection of phytoplasmas, the ZjMAPKKK26 was quickly response in phyllody leaves and then highly induced in witches' broom leaves, however, ZjMAPKKK45 response faster than ZjMAPKKK26 because of its high expression in apparent leaves in June. In contrast to these two up regulated ZjMAPKKK genes, theZjMAPKKK3, 43 and 50 were down regulated 9, 9, 10 times, respectively.

Expression analysis of ZjMAPKKKs in the sterile cultivated tissues of JWB plantlets
After investigating the expression profiles of ZjMAPKKKs in field tissues, we further detected their expression levels in the sterile cultivated (in vitro) tissues of JWB plantlets and the healthy plantlets were used as control. As shown in Fig. 9, the expressing profiles of ZjMAPKKKs indicated significant disparities compared with the above results. Only four members of ZjMAPKKKs were significantly induced in the disease plants, i.e., 4, 10, 25 and 44. While the ZjMAPKKK7, 30,35,37,40,41,43 and 46 were significantly down regulated. The others showed none significantly changes.

Discussion
MAPK cascades, as conserved signal transducing modules in eukaryotes, are widely studied in the plant kingdom, including Arabidopsis, rice, maize, apple and so on [6,10,13,14,16]. In Chinese jujube, 10 MAPKs and 5 MAPKKs have been identified according to our previous study [18]. The structures of these genes were mostly as similar as to other plants, such Arabidopsis, poplar or apple [18]. In the current study, it is the first time to report the MAPKKKs in Chinese jujube genome, 56 ZjMAPKKKs were identified and the number of ZjMAPKKKs was approximately larger than that of VviMAPKKKs [15], but significantly smaller (nearly 45%) compared with the number of MdMAPKKKs [16]. Moreover, all the ZjMAPKKKs could be classified into two main subfamilies (Raf and MEKK), and interestingly none of them belong to ZIK subfamily. The reason for the disappearance of ZIK subfamily may due to the loss function of them during the evolutionary process in Chinese jujube, because the Raf subfamily was considered as the original members of MAPKKKs family [32]. In addition, the rate of intron loss is faster than the rate of intron gain after segmental duplication [33] and in our study the number of the conserved motifs and exons were all higher in Raf subfamily compared with that in MEKK subfamily, and 13 segmental duplication events including 22 ZjMAPKKKs were also identified, these data may further illustrate that the Raf subfamily contains original genes and the segmental duplication occurred during the long evolutionary history. This is consistent with what have been found in maize [14]. Furthermore, only 43 members of ZjMAPKKKs were sparsely located on the 12 chromosomes and the ZjMAPKKK44-56 could not find the corresponding position in chromosomes ( Table 1). All these findings may demonstrate that the evolutionary duplication process of ZjMAPKKKs took place and the unknown location of ZjMAPKKKs may confer a number of paralogous genes and play critical roles during the biological process, such as ZjMAPKKK44, 45, 46, and 50, which might involve in the process of phytoplasma infection. This is different with our previous study on ZjMAPKs and ZjMAPKKs that they did not experience genome duplication in the evolutionary process [18].
MAPK cascades have been demonstrated as key important signalling modules in response to biotic and abiotic stresses, particularly the apparently wide implication in pathogen attack [34]. Among the MAPKKKs, the MEKK subfamily has been widely studied, while the biological function of Raf subfamily is still elusive at present. In Arabidopsis, the MAPK cascades signaling module MEKK1-MKK4/MKK5-MPK3/MPK6 was proposed to be activated in the interaction with flg22 treatment [35]. EDR1, as member of Raf subfamily, was indicated to be responsible for the salicylic acid induced powdery mildew attack [36]. In tomato, MAPKKKε and MAPKKKα play important roles in the cell death signaling associated with plant immunity [37,38]. In wheat, a MAPKKK named TaFLR could be activated by leaf rust pathogen Puccinia triticina [39]. While in grapevine, VqMAPKKK38 could be highly induced by powdery mildew infection [40]. All the evidence suggested that the member of MAPKKKs is essential in the pathogen attack signal transduction. In this study, we demonstrate that the ZjMAPKKK26 and  30,35,37,40,41,43 and 46 were significantly down regulated. The potential candidate genes responses to phytoplasmas were different in vivo and in vitro. This was consistent with our previous work that ZjMPK1 was considered as the potential genes related with phytoplasmas infection in vitro [18] which was also different with our late work that ZjMPK2 was the main genes which could be regulate with ZjMEK2 in vivo [23]. The reason behind this might be that more environmental factors (light, temperature, etc.) affect the phytoplasma growth, evolution and infection in vivo. Even though, all these candidates which have been depicted above might be the potential genes which involved in phytoplasma infected signal transduction by recruiting the MAPKKs and MAPK families.
Moreover, the ZjMAPKKK43 which is homologous to AtRaf1 in Arabidopsis was down regulated in all the phytoplasma infected tissues, and the ZjMAPKKK10 was highly induced in the sterile cultivated tissues of JWB plantlets, this gene is homologous to AtMEKK1 in Arabidopsis. Thus, these two genes might be the key gene in response to phytoplasma infection, because the AtMEKK1 has already been demonstrated function importantly in the innate immunity in Arabidopsis by activation of MKK4/MKK5-MPK3/MPK6 [35]. Moreover, the ZjMKK2 could activate the ZjMPK2 to play essential roles in JWB defense response [23], Ye et al. [22] have illustrated that after the phytoplasma infection, the MAPKs could also be activated and furthermore, the transcription factor WRKY33 was regulated. Taking all together, the potential ZjMAPKKKs mediates ZjMKK2-ZjMPK2 to WRKY transcription factors in response to phytoplasma infection need to be illustrated in the future study. In addition, the scaffold RACK 1 (Receptor for Activated C Kinase 1) was identified in Arabidopsis which tethers the MAPKKKs to the plasma membrane and associates with Gb subunit to involve in immune responses [41]. Therefore, it is worth noting to demonstrate the function of MAPKKKs candidates and the relationship to RACK1in response to phytoplasma infection in Chinese jujube in the further study.

Consent for publication
Not applicable.

Ethical standards
This research does not contain any studies with human participants or animals.

Availability of data and materials
All data and materials are presented in the main paper and additional file. In addition, the whole protein sequences of MAPKKKs in Arabidopsis were retrieved from TAIR databases. The CDS and genome sequences of MAPKKKs in jujube were retrieved from the whole jujube genome database (accession JREP00000000) in NCBI.

Figure 1
The phylogenetic analysis of ZjMAPKKKs (Ziziphus jujuba Mill.) and AtMAPKKKs (Arabidopsis thaliana) with total 136 protein sequences. The MEGA 6.0 was used to construct the phylogenetic tree with the neighbor-joining (NJ) method, and 1000 bootstrap replications were performed to indicate the reliability. All the ZjMAPKKKs were clustered into two groups, called Raf and MEKK subfamilies. Identifing the conserved motifs of ZjMAPKKKs in corresponding to the phylogenetic tree. The MEME database was used to do the motifs' identification in according to the protein sequences.