The pure line PY4-7 of P. yezoensis, established by our laboratory at Ocean University of China, was used in this study. Fresh leafy gametophytes were cultured in sterilized seawater with Provasoli’s enriched seawater medium at 12°C and 60 μmol photons ⋅ m^–2.s^–1 with a 12:12 h light/dark cycle. The seawater was bubbled continuously with filter-sterilized air and renewed every 3 days. For water stress, the thalli were desiccated until the total water content lost 30, 50, and 80% of full hydration. And then after losing 80% water content, the samples were recovered in normal seawater for 30 min. Three biological replicates were set for each treatment. These samples were harvested and placed in liquid nitrogen for gene expression analysis.

The P. yezoensis MAPKKK gene family was identified following the methods described by Rao et al., with some modifications (Rao et al., 2010; Li et al., 2018). Our laboratory has completed the sequencing of P. yezoensis genome and the sequencing data were deposited under DDBJ/ENA/GenBank WMLA00000000. Firstly, the local protein database predicted by our laboratory was blasted with 155 known MAPKKK gene sequences including A. thaliana (80) and O. sativa (75). The threshold was e-value of le^–5. The alignments of all the putative MAPKKK sequences were further used to construct hidden Markov model (HMM) profile using the hmmbuild tool in HMMER3.0¹, and then the HMM profile was searched against the local protein database. HMMER and BLAST hits were compared and the common sequences were selected as the putative PyMAPKKK proteins, which were submitted to the NCBI Batch CD-search database² and PFAM database³ to confirm the presence of the kinase domain. Finally, 33 proteins remained and were considered to be PyMAPKKKs. The prediction methods of MAPKKKs for other red algae (Cyanidioschyzon merolae, P. haitanensis, Porphyra umbilicalis, and Porphyridium purpureum) were the same as above. Their genomes have been sequenced and the accession numbers are GCA_000091205.1, GCA_008729055.1, GCA_002049455.2, and GCA_000397085.1, respectively. The theoretical pI and Mw of the putative PyMAPKKKs were calculated using the Compute pI/Mw tool in the ExPASy database⁴. Subcellular localization of each PyMAPKKK cascade kinase was predicted by the WoLF PSORT tool⁵.

The kinase domain sequences of the predicted MAPKKKs in P. yezoensis, the four red algae species mentioned above, and A. thaliana and O. sativa were multi-aligned using the MUSCLE (Multiple Sequence Comparison by Log-Expectation) program. A phylogenetic tree was constructed using MrBayes version 3.2.5 (Ronquist et al., 2012) with the following parameter settings: prset aamodelpr = mixed, mcmc ngen = 10 million samplefreq = 10. The operation terminated when the standard deviation of the split frequencies was ≤ 0.01, and the operation results were visualized using sump burnin = 250 and sumt burnin = 250.

The DNA sequences of PyMAPKKK gene containing exon sequences were submitted to Gene Structure Display Server⁶ to display the exon/intron structures. To assess the regulatory elements of PyMAPKK genes, 2-kb upstream sequences of gene ORF were selected. All the sequences were submitted to the PlantCARE database⁷ to identify the cis-acting regulatory elements (Lescot et al., 2002).

In order to gain insight into the expression characteristics of the PyMAPKKK gene family when the alga is exposed to different stressors, including different dehydration and temperature treatments, transcriptome sequencing data were obtained to investigate the differential expression of PyMAPKKKs. The RNA-sequencing (RNA-Seq) data were obtained from a previous study by our laboratory under the accession numbers PRJNA236059 and PRJNA401507 (Sun et al., 2015). The treatments of water and temperature stress were set according to RNA-sequencing. Three biological replicates were tested for each condition. TopHat and Cufflinks were used to analyze gene expression based on the RNA-Seq data (Trapnell et al., 2012). The FPKM value (Fragments Per Kilobase of transcript per Million fragments mapped) of each PyMAPKKK gene was calculated, and the log2-transformed values of the 33 PyMAPKKK genes were used to generate heatmaps. The heatmaps and dendrograms were created in R × 86_64-w64-mingw32/ × 64 platform (64-bit) (Maindonald, 2008).

Total RNA was extracted using the RNeasy Plant Mini Kit (OMEGA) according to the manufacturer’s procedure. qRT-PCR analysis was followed by Kong et al. (2017). Then 1 μg total RNA was used to synthesize the first-strand cDNA using a First Strand cDNA Synthesis Kit (Roche, Germany). qRT-PCR was performed with Light Cycler^®480 (Roche, Germany). Sequence-specific primers of the genes studied in this research were designed using Primer 5.0 software (Premier Biosoft International, Palo Alto, CA, United States) (Table 1). The PyTub1 and PyeIF4A were selected as internal controls. The 2^–ΔΔCt method was used to assess the expression level of genes (Livak and Schmittgen, 2001). The expression levels were compared using one-way analysis of variance (ANOVA) (P < 0.05) in SPSS 19.0 software. The qPCR results shown were the average (± SD) of three biological repeats.

Protein.links.v10.5.txt.gz and protein.sequences.v10.5.fa.gz were download from STRING⁸ (Szklarczyk et al., 2010). Protein.sequences.v10.5.fa.gz was used as a standard database, and the PyMAPKKK and P. yezoensis protein databases were compared against the standard database, and then the links file (Protein.links.v10.5.txt.gz) was used to establish the relevance of the two databases (PyMAPKKK sequence database and P. yezoensis protein database). Furthermore, sequences with an interaction score ≥ 400 were selected, and the functional annotations were searched in the non-redundant (Nr) database. Finally, the network was visualized using Cytoscape version 3.5.1. Additionally, Gene Ontology (GO) analysis was also performed using OmicShare Tools⁹.

PyMAPKKK gene family members were identified based on the complete genome sequence of P. yezoensis. A total of 155 known MAPKKK protein sequences were used as a query to search the P. yezoensis protein database by local BLAST search. 56 hits were obtained as candidate sequences. Out of these hits, only 33 MAPKKK protein sequences contained a kinase domain and their genes contained the full-length ORF. The filtered MAPKKK genes in P. yezoensis were designated as PyMAPKKK1 to PyMAPKKK33.

The length of the cDNA of the PyMAPKKK genes ranged from 642 bp (PyMAPKKK1) to 6,762 bp (PyMAPKKK5). The encoded polypeptides ranged from 213 to 2,253 amino acids (aa). The predicted molecular weight (Mw) ranged from 22.74 kDa (PyMAPKKK1) to 219.15 kDa (PyMAPKKK5), and the isoelectric point (pI) values ranged from 5.45 (PyMAPKKK11) to 11.14 (PyMAPKKK32). The subcellular location prediction indicated that most of the MAPKKK proteins were localized in the chloroplast (57.6%), followed by the nucleus (27.3%) and the cytoplasm (12.1%). Only PyMAPKKK2 was predicted to localize in the plasma membrane. The diverse subcellular localization of the PyMAPKKK is similar with that of higher plants. The accession number, the length of the cDNA and protein, the Mw, and the pI and subcellular localization of the PyMAPKKKs are summarized in Table 1.

Furthermore, we also analyzed the MAPKKK family in other red algae, including single-cellular and multi-cellular species for which complete genome sequences were available. The results showed that there are 19 MAPKKKs in C. merolae, 43 MAPKKKs in P. purpureum, 31 MAPKKKs in P. umbilicalis, and 28 MAPKKKs in P. haitanensis (Table 2). As mentioned above, there are over 60 putative MAPKKK genes in higher plants. This indicated that the number of MAPKKK genes from P. yezoensis and red algae is significantly lower than in higher plants.

In plants, the MAPKKK family can be divided into MEKK, Raf, and ZIK subfamilies based on their conserved domain. The MEKK-specific signature is G (T/S) PX (F/Y/W) MAPEV; ZIK members have a GTPEFMAPE (L/V) (Y/F) conserved signature motif; and Raf members possess the GTXX (W/Y) MAPE motif (Feng et al., 2016). To subcategorize PyMAPKKKs, we further analyzed the conserved signature motif in PyMAPKKKs. In our research, 33 PyMAPKKK members contained three types of conserved motifs. Twenty-five PyMAPKKKs contained G (T/S) xxxxAPE, which is similar to the Raf type. PyMAPKKK1, PyMAPKKK2, PyMAPKKK10, PyMAPKKK13, PyMAPKKK17, and PyMAPKKK19 contained G (T/S) P (L/F/Y) WMAPE (L/V), which is similar to the ZIK type. PyMAPKKK1 and PyMAPKKK21 contained the signature of GTPHxxAPE, which is similar to MEKK (Figure 1). Based on the results, the Raf subfamily was the largest subfamily, whereas the MEKK subfamily was the smallest subfamily in P. yezoensis. In comparison to the conserved signature motifs of higher plant MAPKKKs, some minor variations existed in the PyMAPKKKs, yet the similar classification results indicated that the divergence of the MAPKKK gene family occurred prior to the separation of higher plants and rhodophytes.

To further analyze the phylogenetic relationships of the MAPKKK genes in P. yezoensis, the kinase domain sequences of MAPKKKs from P. yezoensis and four other Rhodophyta species including three macroalgae and two primitive algae were aligned and a Bayesian phylogenetic tree was constructed. Based on the phylogenetic analysis, the PyMAPKKK family could be classified into three major groups with strong bootstrap values, each including members of the Raf, ZIK, and MEKKs subfamily, respectively (Figure 2). For PyMAPKKKs, 25 members clustered together, forming the largest group (Raf group), whereas PyMAPKKK11 and PyMAPKKK21 clustered together and formed the smallest group (MEKK group). Others, including PyMAPKKK1, PyMAPKKK2, PyMAPKKK10, PyMAPKKK17, and PyMAPKKK19 formed the ZIK group (Figure 2). In the Raf group, 25 members could be further divided into two subgroups with 10 and 15 members, respectively. This may indicate that obvious gene expansion or divergence occurred in the Raf group of PyMAPKKK during evolution, forming two sister groups. The phylogenetic tree from five red algae species also showed three groups with strong bootstrap values (Figure 3 and Table 2). Similar with PyMAPKKK classification, the Raf subfamily was the largest, followed by the ZIK subfamily, and MEKK subfamily was the least. Yet in higher plants, Raf subfamily was also the largest, while ZIK subfamily was the least subfamily. In general, most members containing the typical signal domain were classified into the major group, especially for Raf and MEKK family members. In ZIK family, CmMAPKKK6 did not cluster together with other members and formed a single clade. It was inferred that it was an independent member of CmMAPKKK. Furthermore, the MAPKKK orthologs of the multicellular red algae clearly clustered together, and the PyMAPKKKs almost clustered together with the PuMAPKKKs and PhMAPKKKs, which was in consistent with the evolutionary relationship of species. Compared to the numbers of MAPKKKs in unicellular algae, it maybe that the MAPKKK gene family of multicellular red algae experienced significant evolutionary expansion. Gene expansion event also existed in MAPKKK gene family of higher plants, especially in Raf and MEKK subfamily (Feng et al., 2016).

To analyze the evolutionary relationships of MAPKKK proteins between rhodophytes and higher plants, a phylogenetic tree was built based on an alignment of the amino acid sequences of all 309 MAPKKK kinase domains (33 in P. yezoensis, 19 in C. merolae, 43 in P. purpureum, 31 in P. umbilicalis, 28 MAPKKKs in P. haitanensis, 80 in Arabidopsis, and 75 in rice). As indicated in Figure 4, most MAPKKKs still formed three large groups with strong bootstrap values, suggesting that the divergence of the three MAPKKKs subfamilies occurred prior to the divergence of rhodophytes and green plants. It is worth noting that CmMAPKKK10, PyMAPKKK4 and PuMAPKKK13 belonging to Raf subfamily, were clustered with ZIK subfamily members. It was inferred that although they contained typical kinase domain, yet more variations occurred in other domain of their protein. Furthermore, in each group, a few MAPKKKs clustered together with those of higher plants, while most red algae MAPKKKs clustered together, indicating that some rhodophyte MAPKKKs evolved independently. ML phylogenetic tree was also constructed to show the same topological tree (Supplementary Figure 1), which indicated the constructing tree method with Bayesian was reliable.

Protein kinases play roles by protein phosphorylation (Jiang et al., 2015). In this research, we found that a kinase domain (IPR00719) existed in all the PyMAPKKK proteins, and most of proteins contained the Ser/Thr kinase active site (IPR008271) in the central part of the catalytic domain, except for five PyMAPKKKs (PyMAPKKK1, PyMAPKKK2, PyMAPKKK10, PyMAPKKK12, and PyMAPKKK13). The existence of these typical features reflected the MAPKKK conservation of structure and function, which is the same in red algae and green higher plants. The ATP-binding site (IPR017441) is also located on the catalytic domain of most PyMAPKKKs, which is a conserved sequence in the kinase family (Wang et al., 2015), suggesting that ATP as a substrate is involved in the PyMAPKKK phosphorylation cascade, activating the activity of downstream protein kinases to complete the signal transduction process. In addition, the PyMAPKKKs also possessed some other conserved domains. For example, an ankyrin repeat-containing domain (IPR020683), which is related to protein-protein interaction and has been found in different kinds of proteins, such as cell cycle regulators, ion transporters, transcriptional initiators, the cytoskeleton, and signal transducers (Gorina and Pavletich, 1996; Sedgwick and Smerdon, 1999), also existed in the domain of PyMAPKKK4 and PyMAPKKK5. PyMAPKKK23 possesses the TyrA (ACT) domain, which is involved in metabolic enzymes by responding to amino acid concentrations changes (Aravind and Koonin, 1999). PyMAPKKK11 and PyMAPKKK21 possess the armadillo-like helical (IPR011989), which has been implicated in protein-protein interactions (Han et al., 2012). PyMAPKKK9 and PyMAPKKK18 contain the domain of the toll/interleukin-1 receptor homology (IPR000157), which has been demonstrated to play a role in immune responses for pathogen recognition (Wang et al., 2016b), and thus PyMAPKKK9 and PyMAPKKK18 may be involved in the regulation of innate immunity in P. yezoensis. Identification and analysis of these conserved domains demonstrated that the PyMAPKKKs possess a variable functional domain, which may be helpful in the adaptation to harsh abiotic and biotic intertidal zone environments. These conserved domains were also found in the MAPKKK proteins of higher plants, indicating that the functions of MAPKKKs were conserved during in evolution.

Gene structure analysis facilitates a better understanding of gene evolution, function, and regulation (Wu et al., 2014). The gene structure of the PyMAPKKKs was relatively variable. Generally, 15 PyMAPKKK gene members contained no intron, 14 members contained one intron, three members contained two introns, and only one gene contained three introns (Figure 2). Although the number of introns among the three subfamilies did not show too much variation, the MEKK group contained the highest variation in structure among the three groups. Among the Raf subfamily, two out of 25 genes had two introns, 10 genes had one intron, and 13 gens had no intron, while two of these genes (PyMAPKKK5 and PyMAPKKK22) had longer introns. Overall, the number of MAPKKK introns in P. yezoensis was less than that in higher plants, which usually contain seven to eight introns per gene, and up to 20 introns in total (Feng et al., 2016; Wang et al., 2016a). Intron gain or loss are considered to be the result of selection pressures during evolution, and thus it can be inferred that the speed of intron gain in P. yezoensis was much slower than that of green plants during gene evolution. The genes of P. yezoensis have few introns with 0.63 of the number of introns per gene (Nakamura et al., 2013). This phenomenon is also very common in the P. yezoensis genome. A low number of introns has also found in other red algae genes and is possibly a reason resulting in a reduction in genome size of red algae (Collén et al., 2013; Nakamura et al., 2013).

The promoter is the DNA region of the transcription factor recognizing and binding and plays a key role in the regulation of gene spatial and temporal expression (Kong et al., 2012). Prediction analysis of the 2 kb upstream region of the start codon ATG of PyMAPKKKs was used to screen for cis-regulatory elements (Supplementary Table 1). The results showed that the cis-acting element could be divided into two groups: motifs related to plant growth and development and motifs related to the stress response, which is similar to results obtained in Brachypodium, tomato, and cucumber (Wu et al., 2014; Jiang et al., 2015; Wang et al., 2015). Plant growth and development-related motifs include five categories: light response elements (G-box, sp1, ACE, 3-AF1 binding site), meristem expression elements (CAT-box and CCGTCC-motif), an endosperm expression element (sk-1_motif), a metabolism related element (O₂-site), and a circadian clock control element (circadian motif). Abiotic stress-related motifs include five categories: a cis-acting regulatory element essential for anaerobic induction (ARE), a cold-stress related element (LTR), drought related elements (TC-repeat/DRE and MBS), wounding related elements (Box-W1, W-box, and box S), and hormonal regulation related elements (ABRE, CE3, AuxRR-core, P-box, CARE-motif, TGA-element, TCA-motif, CGTCA-motif, and TGACG-motif), indicating that PyMAPKKKs may be involved in regulating various stress responses and hormone signaling transduction processes. It was noted that nine PyMAPKKKs contained LTR, and 26 PyMAPKKK genes possessed the drought stress related elements TC-repeat/DRE and MBS, which are involved in an ABA-independent pathway in response to dehydration and salt stress (Dolan-O’Keefe and Ferl, 2000). This indicates that PyMAPKKKs play an important role in the adaptation of P. yezoensis to the inter-tidal environment where it is exposed to low temperatures and desiccation. Various regulatory elements related to hormones, including a TCA-element related to salicylic acid (SA) expression, P-motif and CARE-motif related to gibberellic acid (GA)-responsive elements, a TGA-element related to an auxin-responsive element, and CGTCA-motif and TGACG-motifs related to methyl jasmonate (MeJA)-responsiveness, were also identified in the PyMAPKKKs (Zheng et al., 2017). In particular, 31 PyMAPKKK contained an ABA responsive element (ABRE element), which is involved in ABA responsiveness. ABRE is a key cis-regulatory element in ABA signaling and mediates plant responses to many biotic and abiotic stresses, such as water stress, high salinity, temperature extremes, and exposure to UV light (Watanabe et al., 2017). For example, OsMAPK5 was shown to be activated by ABA and was found to further increase tolerance to abiotic (pathogen infection) and biotic (wounding, drought, salt, and cold) stresses in rice (Xiong and Yang, 2003), and thus it was inferred that PyMAPKKK genes may be involved in the mediation of ABA in the response of P. yezoensis to abiotic stress.

To gain further insight into the spatial and temporal patterns of different MAPKKK family members in P. yezoensis in response to stress factors, their expression profiles under dehydration and temperature stress were compared. The 33 PyMAPKKK genes were investigated using available RNA-Seq data for two different stress conditions. The heatmaps in Figures 5, 6 represent the abundance of each transcript at hydration and temperature stress. The expression levels of the PyMAPKKKs varied significantly in the different conditions based on the log2-transformed (FPKM) values. As indicated in studies in higher plants, drought can seriously limit the productivity and distribution of plants. As the core regulatory components of the MAPK cascade, MAPKKKs participate widely in the response of plants to abiotic stress (Shou et al., 2004; Shitamichi et al., 2013; Virk et al., 2015; Jia et al., 2016). To assess the expression profiles of PyMAPKKK genes under dehydration stress, the samples were treated with different degrees of dehydration and rehydration stress (Figure 5 and Supplementary Table 2). Generally, the expression profiles of the 33 PyMAPKKKs under dehydrated conditions showed three different patterns. Cluster A showed that, with an increase in water loss, the expression of MAPKKKs was down-regulated. Following rehydration, yet their expression showed a different trendency, especially PyMAPKKK1, PyMAPKKK13, PyMAPKKK117, and PyMAPKKK24 with up-regulation expression pattern. The set of MAPKKK transcripts in Cluster B first increased in abundance at mild dehydration (50%) and then decreased, and then further increased following rehydration. PyMAPKKK transcripts in Cluster C showed that gene expression increased with an increase in water loss, with the highest expression level observed at 80% dehydration. After rehydration, expression level further decreased. In summary, the different gene family members exhibited great disparities under different dehydration and rehydration stresses. This phenomenon also exists in many higher plants (Wang et al., 2014; Hu et al., 2016). In addition, many cis-acting regulatory elements were identified in response to drought stress, clearly indicating that the PyMAPKKK gene family is involved in dehydration and rehydration, not only at a transcript level but also at a regulatory level.

As an intertidal zone species, the thalli of P. yezoensis are exposed to air during tidal changes and endure dramatic changes in temperature, which fluctuates from 10 to 20°C. To evaluate the role of PyMAPKKK genes under temperature stress, their expression was analyzed under different temperature stress treatments. The expression profile of PyMAPKKK genes under temperature conditions is provided in Figure 6 and Supplementary Table 3. Overall, it is evident that the different gene members differ in their sensitivity to different temperature stresses. The expression profiles could be grouped into four main clusters: cluster A showed that high temperature could inhibit some PyMAPKKK gene expression, while low temperature induced PyMAPKKK gene expression. Furthermore, with the decrease in temperature, the expression patterns of these genes also showed variable sensitivity; for example, the expression levels of PyMAPKKK17 and PyMAPKKK19 increased significantly at −8°C. In cluster B, extreme low temperature and high temperature dramatically inhibited the expression of some genes, particularly PyMAPKKK33, PyMAPKKK15, and PyMAPKKK10. In cluster C, high temperature (24°C) significantly induced gene expression, while low temperature inhibited gene expression to some extent. In contrast, in cluster D, only at an extremely low temperature (−8°C) did the expression of this set of genes increase significantly, while no significant changes were observed under the other temperature stresses. Many studies in higher plants have indicated the role of MAPKs in temperature stress, including low temperature and high temperature (Teige et al., 2004; Sinha et al., 2011; Wu et al., 2014; Jia et al., 2016). Furthermore, as observed with drought-related cis-regulatory elements, some PyMAPKKK genes also contained elements related to temperature, such as the LTR motif, which further confirmed that PyMAPKKs respond to low and high temperature at transcript and regulatory levels, exhibiting diverse expression patterns.

In addition to the differential expression patterns of gene members under a single stress factor, we also found that the same gene member could respond to more than one stress. For example, PyMAPKKK12, PyMAPKKK30, and PyMAPKKK31 were up-regulated under drought stress and were also specifically expressed under extreme cold stress. In higher plants, paralogous genes show similar or different expression profiles responding to the same stressors; for example, among the MAPKKK genes belonging to the Raf family in bread wheat, three genes showed specific expression under salt stress, while another two genes were specifically expressed under drought stress (Wang et al., 2014, 2016a). This phenomenon was also observed in the present study. PyMAPKKK4/PyMAPKKK11/PyMAPKKK13 and PyMAPKKK1/PyMAPKKK2/PyMAPKKK19 showed similar expression patterns under cold stress and drought stress treatments, while PyMAPKKK10/PyMAPKKK17 and PyMAPKKK10/PyMAPKKK19 exhibited different expression patterns under the same stress treatments. This suggested that different paralogous genes in the same gene family experienced different selective pressures during evolution, and their functions show convergence or differentiation in stress factor adaptation. Regardless, these stress-induced MAPKKK genes still are helpful for further revealing the roles of PyMAPKKKs in the response of P. yezoensis to diverse stress processes.

To understand the interaction between PyMAPKKKs and other P. yezoensis genes, their regulatory networks were predicted using the STRING database and local BLAST (Figure 7). A total of 28 PyMAPKKKs interacted with 965 proteins, while five PyMAPKKKs (PyMAPKKK4, PyMAPKKK17, PyMAPKKK19, PyMAPKKK26, and PyMAPKKK28) did not interact with any proteins, suggesting that PyMAPKKKs are widely involved in regulatory networks and metabolic processes in P. yezoensis.

Among these interacting proteins, very few interacted with MEKKK subfamily members (PyMAPKKK11/PyMAPKKK21). PyMAPKKK11 only interacted with an extracellular signal-regulated protein kinase, and PyMAPKKK21 was predicted to interact with a hypothetical protein. In the ZIK subfamily, 66 proteins were predicted to interact with three ZIK members (PyMAPKKK1, PyMAPKKK2, PyMAPKKK10 and PyMAPKKK13), among which 52 proteins were Ser/Thr protein kinases. Other predicted proteins included MAPK, ADP-ribosylation factors, and heat shock proteins. Another two ZIK members (PyMAPKKK17 and PyMAPKKK19) were predicted to have no interaction. In the largest Raf subfamily, among the interacted proteins, the majority of proteins were protein kinases (Ser/Thr protein kinase, MAPK, and protein phosphatase 2C), followed by ubiquitin related enzymes or proteins (ubiquitin-conjugating enzyme, ubiquitin-like protein, ubiquitin-protein ligase E3, and ubiquitin-specific protease).

Of the total predicted interacting proteins, Ser/Thr kinases were the most abundant, followed by MAPKs. Ser/Thr protein kinases, believed to be the predominant protein kinases in eukaryotes, function in signal transduction related to various stress responses by catalyzing protein phosphorylation (Pereira et al., 2011; Pan et al., 2017). Therefore, Ser/Thr kinases may act as an upstream switch of PyMAPKKK to activate the MAPK signaling pathway of P. yezoensis by phosphorylating PyMAPKKK. In this research, 12 PyMAPKKKs were predicted to interact with MAPK, among which five PyMAPKKKs (PyMAPKKK2, PyMAPKKK3, PyMAPKKK5, and PyMAPKKK9) each interacted with seven PyMAPKs. A previous experiment demonstrated that OMTK1 in alfalfa interacts with MMK3 to form a complex involved in hydrogen peroxide-induced cell death (Nakagami et al., 2004). Furthermore, we found that four MAPKKKs (PyMAPKKK9, PyMAPKKK21, PyMAPKKK23, and PyMAPKKK24) interacted with both the 14-3-3 protein and protein phosphatase 2C (PP2C), and thus the activity of the four PyMAPKKKs may be regulated by 14-3-3 and PP2C. A point mutation experiment in HEK293 Epstein–Bar nuclear antigen (EBNA) cells showed that the 14-3-3 is associated with MEKK3 phosphorylation and can prevent dephosphorylation by PP2A (Fritz et al., 2006). In addition, PyMAPKKK2 and PyMAPKKK14 were predicted to interact with Hsp90 and E3 ubiquitin ligase, and the transcriptome expression profile also showed that their expression was up-regulated at 24°C. Further experiments are thus required to evaluate whether they interact to establish a pathway in response to high temperature stress. In A. thaliana, ZEITLUPE E3 ubiquitin ligase can interact with HSP90 to lead to proteasomal degradation, and further enhance the plants resistance to high temperature (Gil et al., 2017).

Interestingly, we also found some new proteins that interacted with PyMAPKKKs, such as histone acetyltransferase (HAT) and histone deacetylase (HDAC), which interacted with PyMAPKKK14 and PyMAPKKK24. Studies have demonstrated that HAT and HDAC contain nuclear receptor complexes that can modify the chromatin structure through binding the promoter of target genes and thus regulate their transcription through the MEK/ERK1/2 pathway (Zassadowski et al., 2012). The manner in which PyMAPKKK14 and PyMAPKKK24 phosphorylate HAT and HDAC to target the substrate and regulate gene transcription requires further experimental validation. Furthermore, PyMAPKKK2 interacted with a calcium dependent protein kinase, and 10 PyMAPKKKs interacted with casein kinase. As is well-known, calcium dependent protein kinases and casein kinases are Ser/Thr-selective enzymes that function as regulators of signal transduction pathways in most eukaryotes (Knippschild et al., 2005; Valmonte et al., 2014). This indicated that a cross-link signaling pathway existed in P. yezoensis in response to environment stress.

GO functional enrichment analysis was performed to evaluate the potential functions of the predicted genes. The results showed that a total of 606 interacting proteins were enriched to GO terms, involved in diverse modules related to biological processes (43%), molecular functions (30%), and cellular components (27%) (Supplementary Figure 2 and Supplementary Table 4). In the biological process category, the ZIK genes interacting with proteins were significantly enriched in metabolic processes, cellular processes, and single-organism processes, and the MEKK and Raf genes were mostly enriched in cellular processes and metabolic processes. In the cellular component category, the MEKK and Raf genes interacting with proteins were significantly enriched in cell part, while none within the ZIK subfamily were enriched in this term. In molecular function category, the MEKK, Raf, and ZIK genes interacting with proteins were all enriched in binding protein and catalytic activity. It has been shown that most MAPKKK sequences consist of long C-terminal regions with a catalytic domain that interact with other signaling molecules to regulate their function and direct their activity in MAPK pathways (Wu et al., 1999).

Gene expression profiles provide the important clue for understanding the function of PyMAPKKKs genes responding to abiotic stress. To further verify the expression patterns of these PyMAPKKKs, 16 genes with significant differences were selected to detect their expression patterns through qRT-PCR analysis (Figure 8). Among them, 15 genes showed highly consistent with those using RNA-seq except from PyMAPKKK4. These results further suggested PyMAKKK were involved in dehydration/rehydration-stress and can be as the candidates for further studying their function in Pyropia abiotic response.