Transcriptome-wide identification of MAPKKK genes in bermudagrass (Cynodon dactylon L.) and their potential roles in low temperature stress responses

Identification of MAPKKK genes in bermudagrass

To identify candidate CdMAPKKKs, transcriptome data of the low temperature-tolerant bermudagrass cultivar WBD-128 was produced. Plants were grown in a greenhouse for two months under natural sunlight and at temperatures of 30/20 °C (day/night) and were then transferred to controlled-environment growth chambers for a 48-h treatment at 4 °C, after which they were exposed to −5 °C for 4 h. We used previously published clean RNA sequencing reads of leaf samples collected at 0 h (termed CdR_0), 24 h (CdRCA_24) and at 48 h (CdRCA_48) of a 4 °C treatment and after 4 h of a −5 °C treatment following the cold acclimation (CA) treatment (CdRCA_4) (Chen et al., 2015a), and a transcriptome read assembly was produced using Trinity software (Haas et al., 2013) with min_kmer_cov set to 2 and all other parameters at default settings. 151,222 ORFs with the minimum length of 100 amino acids were identified among 501,735 Trinity transcript sequences using TransDecoder v2.0.1 (available at http://transdecoder.github.io). CD-HIT (version: cd-hit-v4.6.6) clustered the remaining genes with a 90% identify threshold, the longest representative sequence in each cluster sequence were selected and generated a final set of 93,182 potential protein coding unigenes (Fu et al., 2012). After this, a local protein database was constructed, and a BLASTP search was performed using 228 known MAPKKK protein sequences of A. thaliana (80) (Jonak et al., 2002), Oryza sativa (75) (Rao et al., 2010), and Brachypodium distachyon (73) (Feng et al., 2016), which were downloaded from the Phytozome v11.0 database (Goodstein et al., 2012) with an e-value of 1e−10 and minimum amino acid identity of 50%. A hidden Markov model (HMM) profile of the MAPKKK family domain PF00069 (https://pfam.xfam.org/family/PF00069) (Wang et al., 2015; Finn et al., 2016) was used as a query to search the local protein database using HMMER 3.0 software (Prakash et al., 2017). After this step, HMMER hits were compared with BLASTP results, and a self-blast program was performed to remove redundancy among these obtained sequences. Finally, the candidate sequences were submitted to the databases SMART (http://smart.embl-heidelberg.de/) (Schultz et al., 1998), Pfam (http://pfam.xfam.org/) (Finn et al., 2016), and the Conserved Domain Database (http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) (Marchler-Bauer et al., 2015) to confirm presence and integrity of the serine/threonine-protein kinase domain. The putative protein sequences of 55 CdMAPKKKs were listed in Table S1. Molecular weights and theoretical isoelectric point of each CdMAPKKK protein was estimated using the ExPASy Compute pI/Mw tool (http://www.expasy.ch/tools/pi_tool.html) (Wilkins et al., 1999), and the subcellular localization of CdMAPKKKs were predicted using WoLF PSORT (http://www.genscript.com/wolf-psort.html) (Horton et al., 2007). In order to investigate distribution patterns of MAPKKK genes in Poaceae, candidate MAPKKKs of seven species including Aegilops tauschii, Triticum urartu, Oropetium thomaeum, Sorghum bicolor, Echinochloa crus-galli, Panicum virgatum, and Setaria italica were identified using the same identification strategy, and data on members of the MAPKKK family of five species including O. sativa, B, distachyon, Zea mays, Hordeum vulgare, and Triticum aestivum were obtained from previous studies (Rao et al., 2010; Kong et al., 2013; Feng et al., 2016; Wang et al., 2016; Cui et al., 2019), respectively.

Multiple sequence alignment and phylogenetic analysis of MAPKKK proteins

Amino acid sequences of 228 known MAPKKK genes of A. thaliana (80) (Jonak et al., 2002), O. sativa (75) (Rao et al., 2010), and B. distachyon (73) (Feng et al., 2016), as well as those of CdMAPKKKs identified in the present study were aligned using the ClustalW tool of BioEditV7.0.5.3 (Larkin et al., 2007). MEGA 5.2 software (Tamura et al., 2011) was used to construct a phylogenetic tree using the neighbor-joining method (Saitou & Nei, 1987) with 1,000 bootstraps and the partial-deletion option.

Co-functional gene network analysis of CdMAPKKK proteins

Potential co-functional CdMAPKKKs genes were identified based on orthologous genes of bermudagrass and Arabidopsis using the AraNet V2 tool (https://www.inetbio.org/aranet/), which is a probabilistic functional gene network database of Arabidopsis (Lee et al., 2010; Lee et al., 2015). And then, a Gene Ontology (GO) functional analysis for putative co-functional genes of CdMAPKKKs was carried out by agriGO tool (http://bioinfo.cau.edu.cn/agriGO/analysis.php) using the singular enrichment analysis method (Du et al., 2010). The latest Arabidopsis genome Athaliana_167_TAIR10 was used as a reference and significantly enriched GO terms were identified with all parameters at previously settings (Wang et al., 2017b).

Analysis of putative cis-elements in promoter regions of CdMAPKKKs

Upstream sequences of 12 CdMAPKKKs of up to 2 kb from the translation start codon were retrieved manually (Table S2). Putative stressors or hormone responsive cis-regulatory elements in these sequences were then predicted using Plant Cis-acting Regulatory DNA Elements (PLACE,  http://www.dna.affrc.go.jp/PLACE/) (Higo et al., 1999). Predicted cis-elements were mapped on the sense and anti-sense strand of each CdMAPKKKs using the TBtools software (Chen et al., 2020).

Interaction network analysis of CdMAPKKK proteins

Interaction networks of all CdMAPKKKs were constructed. In detail, homologs of each CdMAPKKK in rice were identified using blastp software. After this, the corresponding rice homologs of each CdMAPKKK gene were submitted to the String databases (https://string-db.org/cgi/input.pl) (Szklarczyk et al., 2017) to predict and construct interaction networks using Markov clustering with an inflation factor of 8.5.

Analysis of expression levels of CdMAPKKKs using RNA sequencing data

For in-silico analysis of CdMAPKKKs in response to low temperature stress, transcriptome obtained using the Trinity software provided a reference sequence. RSEM software was applied to estimate the expression level of transcripts (Li & Dewey, 2011). Initially, Bowtie 2.0 was used to map the clean data of each of the cDNA libraries to the reference transcriptome sequence (Langmead & Salzberg, 2012), and then, the expression level of each transcript was normalized using RSEM method and reported in Fragments per kilobase of exon model per million mapped reads (FPKM)(Li & Dewey, 2011). Finally, the edgeR package (Robinson, McCarthy & Smyth, 2010) was used to identify differentially expressed genes to compare the experimental treatments with the control (4 °C for 0 h treatment, termed CdR_0) based on the counts matrix obtained from the results of normalization step. Fold change values of the 55 CdMAPKKKs were used to analyze their expression patterns in bermudagrass plants subjected to low temperature stress. In addition, a heatmap was produced and cluster analysis of the candidate CdMAPKKKs was performed using MEV (version 4.9) and the hierarchical clustering method (Howe et al., 2011).

Plant material and treatments

To investigate expression of MAPKKKs in bermudagrass subjected to low temperatures, the low temperature-tolerant bermudagrass cultivar WBD-128 was used. Plants were grown in plastic pots containing a mixture of sand and peat soil (v/v 1/1) in a greenhouse under natural sunlight and at temperatures of 30/20 °C (day/night). All plants were ramets of one clone, thus the genetic background of these plants was uniform. In order to induce low temperature stress, two-months-old WBD-128 plants were transferred to controlled-environment growth chambers and were exposed to 4 °C for 72 h under 60% relative humidity and with a 14/10 h (light/dark) photoperiod. Leaf samples of three biological replicates were collected from control and treatment group after low temperature treatment for 1, 3, 6, 12, 24, 48 and 72 h, and stored immediately at −80 °C for qRT-PCR analysis.

RNA isolation and qRT-PCR analysis

The total RNA was extracted using Plant RNA Kit reagent (Omega Bio-Tek, USA), and first-strand cDNA was synthesized using HifairTM II 1st Strand cDNA Synthesis SuperMix (Yeasen, China) according to the manufacturer’s instructions. The quantitative real time-PCR was performed on the StepOnePlus Real-Time PCR Systems (Applied Biosystems, United States) using SYBR Premix Ex Taq (TaKaRa), and the relative expression levels of CdMAPKKKs at 1, 3, 6, 12, 24, 48 and 72 h after exposure to low temperature treatment were determined by qRT-PCR in comparison with control conditions for each time point with CdACT2 was used as the reference gene for normalization. Quantitative primers for CdMAPKKKs were designed using GenScript Real-time PCR (TaqMan) primer design tools (Table S3), and the 2^−ΔΔCT method was used to calculate the expression level of each CdMAPKKK gene (Livak & Schmittgen, 2001).

Identification and characteristics of the MAPKKK gene family

The identification procedure of MAPKKK gene family members in bermudagrass is shown in Fig. S1. In brief, BLASTP and HMMER searches resulted in primary identification of 392 potential MAPKKK protein sequences in the bermudagrass transcriptome data. Redundant sequences were removed from the candidate gene list using a self-blast method, and the kinase domains were analyzed using SMART and the NCBI CDD tool. A total of 55 putative MAPKKKs with ORFs longer than 100 bp were identified in bermudagrass. As there is no standard nomenclature for MAPKKKs, identified MAPKKKs were designated as CdMAPKKK1 to CdMAPKKK55 based on the blast scores obtained from the HMM search output; a similar nomenclature system is used in rice (Rao et al., 2010) and wheat (Wang et al., 2016). These CdMAPKKKs encoded proteins ranged from 126 (CdRAF38) to 1,247 (CdRAF2) amino acids in length, with molecular weights ranging from 14.03 (CdRAF37) to 134.04 kDa (CdRAF1) and predicted isoelectric points in a range of 4.95 (CdZIK3) to 9.74 (CdMEKK2). The prediction of subcellular localizations using WoLF PSORT revealed that a total of 21 CdMAPKKKs were localized in the cytoplasm, 14 in the chloroplast, and 13 in the nucleus, whereas the remaining CdMAPKKKs were predicted to be located in mitochondria, in extracellular space, and in the cytoskeleton (Table 1).

Multiple alignment and conserved motif analysis

As reported in Arabidopsis and other plants, MAPKKKs can be classified into MEKK, Raf, and ZIK subfamilies according to their specific and conserved catalytic kinase domains; MEKK subfamily members are characterized by catalytic domains with G (T/S) PX (F/Y/W) MAPEV, Raf by catalytic domains with GTXX (W/Y) MAPE, and ZIK by catalytic domains with GTPEFMAPE (L/V/M) (Y/F/L) (Tena et al., 2001; Group, 2002). For predicting and subcategorizing the identified CdMAPKKKs, we examined the conserved catalytic domains of these CdMAPKKKs (Fig. 1). Our results showed that 55 CdMAPKKKs possessed at least one of the three conserved catalytic domains. In detail, 13, 38, and 4 CdMAPKKKs shared the conserved catalytic domain of the MEKK-like, Raf-like, and ZIK subfamily, respectively (Figs. 1A–1C). In addition, multiple sequence alignment of MEKK, Raf, and ZIK members confirmed that these MAPKKKs occurred in bermudagrass.

To investigate conserved domains of MAPKKK proteins in bermudagrass, multiple sequence alignments and conserved motifs were analyzed using the MEME website (http://meme-suite.org/tools/meme) and were schematically illustrated based on their evolutionary relationships (Fig. 2, Fig. S2). A total of 10 predicted conserved motifs was identified, referred to as motifs 1-10, and the number of motifs in each protein sequence varied from 2 to 10 (Figs. 2A–2B). Among these motifs, motif 2 was annotated as STKc_MAP3K-like domain, and all CdMAPKKK proteins contained this motif apart from CdMEKK11 and CdMEKK13. Further observation of the protein sequences of CdMEKK11 and CdMEKK13 showed that CdMEKK11 and CdMEKK13 also contained the conserved core amino acid residue MAPEV which is specific of the MAPKKK gene family (Fig. 1A). Additionally, CdMAPKKKs in the same subfamily frequently display high similarity motif patterns. For example, the member of the ZIK subfamily contained the same motifs including motifs 1, 2, 3, 4, 6, and 8, however, CdZIK1, CdZIK2, CdZIK3 and CdZIK4 also contained motif 9, 5, 5 and 7, respectively. Motif analyses thus provided additional evidence supporting the close evolutionary relationship of CdMAPKKKs in the same subfamily (Fig. 2B).

Comparative phylogenetic analysis of MAPKKKs among different organisms

Using protein sequences, we examined phylogenetic relationships of MAPKKK genes using bermudagrass and the three angiosperms A. thaliana, O. sativa, and B. distachyon. The resulting dendrogram illustrated that the 55 CdMAPKKKs were grouped into three distinct clades which belonged to the MEKK, Raf, and ZIK subfamily, respectively, and contained 13, 38, and 4 CdMAPKKKs, respectively, in accordance with the orthologous MAPKKKs in B. distachyon, A. thaliana, and O. sativa (Fig. 3). In addition, the phylogenetic analysis showed that each MAPKKK subfamily existed in both monocotyledons and dicotyledons, and some closely related orthologous MAPKKKs of bermudagrass occurred in A. thaliana, O. sativa, and B. distachyon, suggesting that an ancestral set of MAPKKK genes existed before the divergence of monocotyledons and dicotyledons. Overall, MAPKKKs of bermudagrass are more closely related with MAPKKKs of B. distachyon and O. sativa than with those of A. thaliana, which is in line with the APG taxonomic system.

In order to investigate the distribution of the MAPKKK family, members of the MAPKKK family of seven species were identified in the present study, and those of five species were compiled from previous reports (Rao et al., 2010; Kong et al., 2013; Feng et al., 2016; Wang et al., 2016; Cui et al., 2019) (Figs. 4A–4B). The members of each MAPKKK subfamily showed no significant variation in eight diploid Poaceae species apart from barley which contained 156 MAPKKK genes, and the Raf subfamily was the largest subfamily, in comparison with other diploid species (Fig. 4B). In addition, polyploid species including Panicum virgatum, Triticum aestivum, and Echinochloa crus-galli showed more MAPKKK genes than any diploid species, apart from barley. This result suggests that MAPKKK family members were conserved among almost of all these diploid species, and species-specific polyploidy and higher duplication ratios of barley in which approximately 84% of the genome is comprised of mobile elements or other repeat structures resulted in an expansion of the MAPKKK family in Poaceae (International Barley Genome Sequencing et al., 2012).

Potential co-functional genes of CdMAPKKKs

To detect potential functions of CdMAPKKKs, a co-functional gene network was constructed using the orthology-based method and the AraNet v2 database (Lee et al., 2015). A total of 714 gene pairs of co-functional links were identified, with an average of 12.98 co-functional genes per MAPKKK gene in bermudagrass (Table S4). Substantial variation in the number of co-functional genes among CdMAPKKKs was observed. For example, CdMEKK9 and CdMEKK10 had 123 putative co-functional genes, whereas CdRAF24 and CdRAF31 had only two putative co-functional genes. In addition, the GO enrichment analysis suggested that the top 30 enriched biological processes of co-functional genes were predominantly associated with protein phosphorylation (GO:0006468, FDR = 7.89e−79), signal transduction (GO:0007165, FDR = 9.60e−40), response to stimuli (GO:0050896, FDR = 1.07e−34), with the abscisic acid-activated signaling pathway (GO:0009738, FDR = 2.66e−27), and with other biological processes (Fig. 5A; Fig. S3). Co-functional genes of CdMAPKKKs were enriched in the following molecular processes: kinase activity (GO:0016301, FDR = 8.40e−106), transferase activity (GO:0016740, FDR = 2.10e−58), calmodulin-dependent protein kinase activity (GO:0004683, FDR = 1.40e−18), signal transducer activity (GO:0004871, FDR = 1.32e−18), and other processes (Fig. 5B and Fig. S4). In addition, 30 and 19 co-functional genes associated with responses to cold and heat stress, respectively, were identified, suggesting the potential roles of these genes and their co-functional CdMAPKKKs in response to temperature stress (Table S4).

Characterization of putative cis-regulatory elements in promoter regions of CdMAPKKKs

To gain insights into transcriptional regulation of CdMAPKKK genes, putative stress- or hormone-responsive cis-regulatory elements in the upstream region of 12 CdMAPKKK genes were predicted using the PLACE software (Fig. 6). A total of 30 types of cis-elements were identified, including 17 stress response and 13 hormone response elements. Among them, MYBCORE (S000176) cis-element, which is involved in the response to dehydration, was identified in the promoter regions of all 12 CdMAPKKKs, and other cis-elements involved in responses to dehydration including DRECRTCOREAT (S000418), CBFHV (S000497), MYCCONSENSUSAT (S000407), MYB1AT (S000408), MYB2CONSENSUSAT (S000409), and MYCATRD22 (S000174) were found in the promoter regions of different CdMAPKKK genes. Some cis-elements that are responsive to low temperatures were also identified, such as LTRECOREATCOR15 (S000153) in eight CdMAPKKKs, LTRE1HVBLT49 (S000250) in four MAPKKKs, LTREATLTI78 (S000157) in three CdMAPKKKs, and CRTDREHVCBF2 (S000411) in CdZIK2. Some cis-elements responsive to ABA (ABRERATCAL, ABRELATERD1, ABREOSRAB21, ABREBZMRAB28, ABREATRD22, and ABREATCONSENSUS), gibberellin (GARE1OSREP1 and GAREAT), and auxin (ASF1MOTIFCAMV, ARFAT, and CATATGGMSAUR) were found in the examined sequences. Furthermore, cis-element GCCCORE (S000430) which is associated with responses to pathogens was also included in the promoter regions of nine CdMAPKKKs. These results suggest a regulatory role of CdMAPKKKs in bermudagrass development and in responses to various stressors.

Interaction network analysis of CdMAPKKK proteins

To investigate potential biological functions of CdMAPKKK proteins, interaction networks among CdMAPKKK proteins were constructed based on experimentally validated interactions using STRING software and Markow clustering (MCL) with an inflation factor of 8.5 (Fig. 7). A total of 34 protein pairs were predicted to interact between 29 CdMAPKKK proteins, among them CdRAF3 showed interaction with seven CdMAPKKK proteins (CdMEKK3, CdMEKK5, CdRAF6, CdRAF12, CdRAF19, CdRAF23, and CdZIK3). In addition, the interaction network analysis suggested that CdMAPKKKs were involved in various biological processes including regulation of plant immune responses (CdMEKK1/YODA), responses to various environmental stressors (CdMEKK8/ANP1), stress responses and ethylene signaling (CdRAF6/EDR1), seed development (CdRAF8/ SIS8), abiotic stress, development, and defense (CdRAF12/CTR1), and flowering time (CdZIK1/WNK1). These results suggest that CdMAPKKKs may play important roles in various biological processes and provide some useful clues for further functional studies.

Expression analysis of CdMAPKKKs in response to low temperature stress

To reveal temporal and spatial expression patterns of MAPKKK genes in silico with respect to their responses to low temperature stress, the expression level of 55 MAPKKK genes were retrieved from RNA sequencing of leaves of the low temperature-tolerant bermudagrass cultivar WBD-128 (Chen et al., 2015a). The criteria of P < 0.05 and —log₂ fold change— ≥ 1 were used to identify CdMAPKKKs which were differentially expressed between low temperature stressed and control plants, and a heat map of the 55 CdMAPKKKs was produced (Fig. 8A and Table S5). The results indicated that expression levels of 25 (45.4%) CdMAPKKKs including 12 (21.8%) and 13 (23.6%) were significantly up- or downregulated due to low temperature stress at least at one point in time. In addition, expression levels of 9 CdMAPKKKs changed significantly at all three time points, with 4 upregulated and 5 downregulated CdMAPKKKs (Fig. 8B; Table S5). For example, CdMEKK7 was significantly upregulated in response to exposure to 4 °C for 24 h (CdRCA_24) and for 48 h (CdRCA_48), and it was also substantially upregulated after exposure to −5 °C for 4 h after the 48-h treatment at 4 °C (CdRCA_4); expression levels of CdMEKK3 decreased rapidly after the same treatment.

To validate the results determined by RNA sequencing, qRT-PCR was performed with ten randomly chosen differentially expressed CdMAPKKKs. The results of qRT-PCR were generally in line with the transcriptome data (R ² = 0.80) (Fig. S5). Furthermore, to identify low temperature stress-associated MAPKKK candidate genes, six out of 12 upregulated CdMAPKKKs exhibited a significant log ₂fold change ≥ 3.0 at least at one point in time were selected for qRT-PCR analysis regarding responses to low temperatures (Table S5). Overall, the relative expression levels of all 6 genes were upregulated in response to the 4 °C treatment, which were generally also consistent with the RNA-Seq data apart from CdMEKK11, CdMEKK12 and CdRAF38 exhibited decreased tendency at 4 °C for 24 h or for 48 h compared to control plants by qRT-PCR (Figs. 8C–8H). The expression levels of CdMEKK7 and CdMEKK8 were significantly upregulated at multiple time points after 4 °C treatment, and reached the highest expression level at 24 h. Transcription levels of CdMEKK11 and CdMEKK12 in bermudagrass leaves were rapidly upregulated at earlier time points and reached the highest expression levels at 3 h but were decreased at subsequent time points. The expression level of CdMEKK13 increased with the duration of the low temperature treatment, whereas the expression level of CdMEKK38 was highest at 1 h and subsequently decreased.