Identification and comparative analysis of the CIPK gene family and characterization of the cold stress response in the woody plant Prunus mume

Identification of CIPK Genes in P. mume and other species

Based on the CIPK protein sequences reviewed (SWISS-PROT, https://www.uniprot.org/) (The UniProt Consortium, 2015), including those from A. thaliana and O. sativa, a CIPK model was built using HMMScan software (Finn et al., 2015). The CIPK model along with the protein kinase (PF00069) and NAF domains (PF03822), which were downloaded from the Pfam database (https://pfam.xfam.org/) (Finn et al., 2008), were used as queries to search the genome sequences of P. mume (Zhang et al., 2012), Prunus persica (The International Peach Genome Initiative et al., 2013), Fragaria vesca (Shulaev et al., 2011), Rosa chinensis (Raymond et al., 2018), Prunus avium (Shirasawa et al., 2017), Malus × domestica (Velasco et al., 2010), Pyrus bretschneideri (Wu et al., 2013), Pyrus communis (Chagné et al., 2014), Rubus occidentalis (Chagné et al., 2014), and Prunus yedoensis (VanBuren et al., 2016). The P. mume genome was obtained from the P. mume genome project, while other Rosaceae genomes were downloaded from the Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html#) and GDR databases (https://www.rosaceae.org/). Then, a HMMER search (E-value < 1e-6) was employed to identify the CIPK gene family members. Based on positional information in the P. mume genome project, the physical chromosomal locations of PmCIPK genes were plotted using MapChart software (Voorrips, 2002).

Phylogenetic tree construction and calculation of Ka/Ks ratios

To study the phylogenetic relationships between CIPK genes in P. mume and other species, CIPK proteins from three species (P. mume, A. thaliana, and O. sativa) and CIPK proteins from nine Rosaceae species were used in a multiple sequence alignment with ClustalX 2.0.11 software (Larkin et al., 2007). Subsequently, a phylogenetic tree based on the sequences was constructed via the maximum likelihood (ML) method using 1,000 replicate bootstrap values and the Jones-Taylor-Thornton model using MEGA X (Kumar et al., 2018). Phylogenetic trees from CIPK sequences were annotated and embellished using the online Evolview tool (http://www.evolgenius.info/evolview/#login) (He et al., 2016).

Duplicate PmCIPK sequences were found using a plant genome duplication database (http://chibba.agtec.uga.edu/duplication/) (Lee et al., 2013). The nonsynonymous (Ka) and synonymous (Ks) substitution of PmCIPK sequences were estimated using DnaSP software (Lee et al., 2013) to predict the divergence time (t) and evolutionary strain (Ka/Ks ratio). Based on two common rates (λ) of 1.5 × 10⁻⁸ or 6.1 × 10⁻⁹ substitutions per site per year (Blanc & Wolfe, 2004; Lynch & Conery, 2000), the divergence time was calculated using the formula t = Ks/2λ × 10⁻⁶ Mya.

Gene structures, protein tertiary structures and motif prediction

The exon/intron structures of PmCIPK genes were obtained with Gene Structure Display Server 2.0 (Hu et al., 2015) using genomic sequences and structural information. The PmCIPK protein sequences were submitted for multiple expectation maximization for motif elicitation (http://meme-suite.org/index.html) (Brown et al., 2013) analysis to identify conserved motifs and structural divergence. The PmCIPK proteins were submitted to Pfam and SMART (http://smart.embl-heidelberg.de/) for analysis and confirmation of the CIPK-specific functional domains. The tertiary structures and homologs of PmCIPKs were predicted using the online server Phyre2 (http://www.sbg.bio.ic.ac.uk/phyre2/html) (Kelley & Sternberg, 2009).

Promoter cis-element and miR167/miR172 target site analysis

A total of 1.5 kb-upstream sequences of each PmCIPK gene were predicted using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002). All of the PmCIPK gene family members were predicted for Pmu-miR167 and Pmu-miR172 by psRNATarget (http://plantgrn.noble.org/psRNATarget/) (Dai, Zhuang & Zhao, 2018). The Mfold web server (http://unafold.rna.albany.edu/) (Bindewald, Kluth & Shapiro, 2010) was used to predict miRNA secondary structure.

Expression analysis of PmCIPK genes

To investigate the PmCIPK expression patterns involved in tissue development and cold stress tolerance, raw data from the RNA sequencing of various tissues (i.e., bud, leaf, root, stem, fruit, and flower buds sampled in November, December, January, and March) in P. mume were obtained (Zhang et al., 2018b). Transcriptome data from three periods (autumn, October; winter, January; and spring, March) and three different places (Beijing (BJ, N39°54′, E116°28′), Gongzhuling (GZL, N43°42′, E124°47′), and Chifeng (CF, N42°17′, E118°58′)) were analysed under low temperature growth conditions. To compare the gene expression patterns of CIPK genes across P. mume, A. thaliana, and O. sative during cold stress, the publicly-available RNA-seq data were downloaded from the National Center for Biotechnology Information GEO DataSets (GSE112225 and GSE67373). The HeatMapper online tool (Babicki et al., 2016) was used to generate the heat map. Molecular interactions of PmCIPKs were analysed using KEGG PATHWAY (https://www.kegg.jp/kegg/pathway.html) (Kanehisa et al., 2017) and STRING (https://string-db.org/) (Szklarczyk et al., 2015).

Plant material and qRT-PCR analysis

A total of 6-month-old seedlings at 24 °C under long-day conditions (16-h light/8-h dark) were used for examining the effect of PmCIPK genes on the cold response. We incubated seedlings in soil at 4 °C at approximately 65% humidity. Leaves from treated seedlings were sampled at 0, 1, 4, 6, 12, and 24 h for total RNA isolation. First strand cDNA synthesis was performed using a TIANScript First Strand cDNA Synthesis Kit (Tiangen, Beijing, China) according to the manufacturer’s instructions. qRT-PCR was carried out using a PikoReal real-time PCR system (Thermo Fisher Scientific, CA, USA) with SYBR® Premix ExTaq TM (TaKaRa, Dalian, China). The reactions were performed in a 10 μL volume containing 5 μL of SYBR® Premix Ex Taq II, 0.25 μL each of forward and reverse primers (Table S1), 0.5 μL of cDNA and 3 μL of ddH₂O. The reactions were completed with the following conditions: 95 °C for 30 s, 40 cycles of 95 °C for 5 s and 60 °C for 40 s, 60 °C for 30 s, and an end step at 20 °C. The analyses were confirmed in triplicate. The relative expression levels of PmCIPK genes were calculated using the 2^−ΔΔCt method with the protein phosphatase 2A gene from P. mume as the reference gene. The final data were subjected to an analysis of variance test.

Identification of CIPK genes in P. mume

Based on the HMMER search using the CIPK model, 16 non-redundant PmCIPKs were identified in the P. mume genome, and 194 CIPKs were identified in the other 10 species from the Rosaceae genome. The putative CIPKs were named based on the hmmsearch E-value of NAF domain (Table S2). Among PmCIPK proteins, there were sequences with a high level of similarity (i.e., more than 51% identical on average). The predicted sizes of the 16 PmCIPKs ranged from 420 (PmCIPK15) to 508 (PmCIPK16) amino acids (aa), and average molecular weight (MW) was 51.24 kDa. The predicted isoelectric points (pI) varied from 7.01 (PmCIPK15) to 9.59 (PmCIPK3), and all of these proteins were alkaline (pI > 7) (Table 1). PmCIPK genes were unevenly distributed in the P. mume genome. The PmCIPK genes were located in all of the chromosomes and were densely distributed on Chr1, Chr2, and Chr3, which contained five, five and three genes, respectively. However, Chr4, Chr6, and Chr7 contained no PmCIPK genes (Fig. S1). A similar phenomenon of unbalanced chromosomal distribution of CIPK genes was also shown in the apple (Niu et al., 2018). The unbalanced distribution of genes may be related to species evolution and genetic variation.

Phylogenetic tree analysis and calculation of Ka/Ks ratios

According to multiple sequence alignments, we used the ML method to construct a phylogenetic tree of all CIPK sequences from P. mume, A. thaliana, and O. sativa. Based on the reviewed (SWISS-PROT) AtCIPKs and OsCIPKs, 16 PmCIPKs were divided into three groups (i.e., Group I, Group II, and Group III) (Fig. S2). To better understand the evolutionary relationships between CIPKs, a ML phylogenetic tree was constructed using the CIPKs from 12 species, including 10 Rosaceae species (Fig. 1). Similar to the tree structure described above, the CIPKs were divided into three groups, and every group included similar clades (I-a, I-b; II-a, II-b; and III-a, III-b). Group I was the largest group, containing 100 CIPKs, whereas group III was the smallest group, consisting of 77 members, indicating that CIPKs were distributed unevenly in the different groups. The CIPKs from the Rosaceae genus were distributed uniformly in every small clade, whereas CIPKs from O. sativa tended to cluster together. Notably, all intron-rich PmCIPKs were clustered in Group I, similar to the clustering of CIPKs from Zea mays (Chen et al., 2011). The PmCIPKs, PpCIPKs, and PaCIPKs were clustered together and had similar distributions in the phylogenetic tree.

In the present study, three PmCIPK gene pairs were identified (PmCIPK11/3, PmCIPK16/3, and PmCIPK10/9), indicating that whole genome duplication/segmental duplications (SDs) may have been involved in the expansion of the CIPK gene family in P. mume (Table S3). During their evolution, a Ka/Ks ratio > 1 implies positive selection (adaptive evolution), a ratio = 1 implies neutral evolution (drift), and a ratio <1 implies negative selection (conservation). The Ka/Ks ratio of only one PmCIPK gene pair (PmCIPK10/9), which was 0.28, has meaning and shows there has been synonymous change. The divergence period of the PmCIPK10/9 gene pair ranged from 41.24 to 76.38 Mya.

Motif prediction, gene structure, and protein tertiary structures

The serine/threonine protein kinases are likely found in all eukaryotes and have roles in a multitude of cellular processes. The kinase, NAF, and PPI domains are the basic characteristics of the CIPK family in plants. The conserved domain showed that all 16 PmCIPK contained the kinase, NAF, and PPI domains (Fig. 2B). The results of MEME analysis showed that all PmCIPK proteins shared six motifs that were contained within the typical domains of plant CIPKs (protein kinase domain and NAF domain). A total of 14 distinct motifs were identified, and motif 1, motif 2, motif 3, motif 4, motif 5, and motif 9 were found in all of the PmCIPKs (Fig. S3). Furthermore, most of the PmCIPKs contained similar motifs, and these motifs are in similar places with similar sizes in the PmCIPK sequences. Similar permutations of motifs exhibited close phylogenetic relationships between the PmCIPK genes (Fig. 2).

The exon-intron diversity is not only an important part of gene family evolution but also may influence their expression. To understand the annotated features of the PmCIPK genes, we used GSDS 2.0 to construct the coordinates of the exon-intron. The results showed that out of 16 PmCIPKs, 10 genes (62.5%) comprised the largest/smallest numbers of exons, which were mainly concentrated in the range of one to three, while six genes (37.5%) contained multiple exons. We found that closely clustered PmCIPK genes in the same clades have similar exon number, UTR exon number, and even the number of exons with significant differences in the PmCIPK genes. Interestingly, the evolutionary relationships of the PmCIPKs were divided into two groups based on different numbers of exons (Fig. 2C).

The three-dimensional structures of 16 PmCIPK proteins were predicted (Fig. S4). Using a template model, all CIPK sequences were modelled with 100.0% confidence by the single highest scoring template. Five identical top-scoring proteins for all PmCIPK sequences were found. Hypothetical protein c6c9dB belonged to the serine/threonine-protein kinase MARK1, while the other four belonged to the protein kinase family. Serine/threonine-protein kinase MARK1 (c6c9dB) shared 28–45% sequence identity with the PmCIPKs, which was anticipated as it is a homologue of CIPK, with 100% probability. On the other hand, the other four protein kinases (c5ebzF, c4wnkA, c3pfqA, and c4cfh) shared 24–44% sequence identity with the PmCIPKs (Table 2; Table S4).

Cis-acting regulatory elements and miR167/miR172 target site analysis

The structure of a promoter affects the affinity of the promoter and RNA polymerase, thereby affecting the level of gene expression. Cis-acting regulatory elements related stress responses, including MYB, MYC, DRE, ABRE, and TC-rich repeats were identified in the promoter of PmCIPKs (Fig. 3). These elements were unevenly distributed in the promoter regions and the number of MYBs was relatively abundant. MYBs are considered crucial TFs in response to water, salt, drought, and cold (Li, Ng & Fan, 2015; Urao et al., 1993). Most of the PmCIPKs, except PmCIPK4, PmCIPK10, PmCIPK14, and PmCIPK15, contained ABRE (ACGTG) cis-promoter elements. ABRE is involved in abscisic acid (ABA) and VP1 responsiveness to abiotic stresses (Hattori, Terada & Hamasuna, 1995; Hobo et al., 1999). Meanwhile, DRE and TC-rich repeats were found in PmCIPK promoter regions, which respond to dehydration and cold stress or high-salt stress (Germain et al., 2012; Narusaka et al., 2003).

MicroRNAs are important regulators of gene expression associated with abiotic stress in plants. Recent studies showed that miR167 and miR172 play a key role in the low temperature stress response in plants (Koc, Filiz & Tombuloglu, 2015; Zhang et al., 2008). In P. mume, we found that four miR172 (Pmu-miR172a-d) and four miR167 (Pmu-miR167a-d) members had been identified (Wang et al., 2014). Notably, two PmCIPKs (PmCIPK5 and PmCIPK6) were targeted by Pmu-miR172s, and only PmCIPK13 was targeted by Pmu-miR167s in the psRNATarget online software (Fig. 4). We drew the energy dot plot for Pmu-miR172a, Pmu-miR172c, and Pmu-miR167b, and the δG values were 3.7, 5.9, and 3.9 kcal/mol at 37 °C in the plot file, respectively (Fig. S5). The energy required for the microRNA target is considered an important determinant of the response of mRNAs to miRNAs (Kertesz et al., 2007). miR167 and miR172 are highly conserved miRNA family members among plants (Jones-Rhoades & Bartel, 2004). Wang et al. (2014) have shown that Pmu-miR167 and Pmu-miR172 family members negatively regulate their targets during the flower bud development process.

Expression pattern analysis of PmCIPKs

In different tissues (i.e., bud, leaf, root, stem, fruit), 16 PmCIPK gene heat maps showed differential expression, suggesting specific functions at particular stages of development (Fig. 5A; Table S5). For instance, seven genes (PmCIPK1, PmCIPK2, PmCIPK5, PmCIPK6, PmCIPK9, PmCIPK12, and PmCIPK13) exhibited high expression in the leaf. All of the PmCIPK genes were expressed during the bud dormancy process, and most of them had specific functions at particular stages of development (Fig. 5B; Table S6). We found that 6 PmCIPKs in subset I showed high expression levels at the NF stage, whereas subset II genes were significantly down-regulated except for PmCIPK10. Some of the genes were significantly up-regulated during the dormancy-breaking process. For example, four PmCIPK genes (i.e., PmCIPK1, PmCIPK2, PmCIPK6, and PmCIPK10) were simultaneously up-regulated at the EDII stage, while others, including PmCIPK5 and PmCIPK14, were induced at the EDIII stage.

To further analyse the expression patterns of PmCIPK genes under cold response, PmCIPK genes were clustered based on different growth regions and periods. We found that the expression profiles of PmCIPK genes were significantly clustered in three different periods (autumn, winter, and spring). The expression levels of most PmCIPK genes was significantly up-regulated in autumn and winter, showing a response to cold acclimation and low-temperature stress. However, PmCIPK14 and PmCIPK15 showed high expression levels in the spring. According to the clustering analysis of PmCIPK genes in different regions, the results are similar to those obtained at different periods (Fig. 5C; Table S7). Most up-regulated PmCIPK genes were clustered in winter. The clustering results showed that the expression levels of most PmCIPK genes was up-regulated at low temperature. Among the up-regulated genes, we identified a SnRK (SNF1-related protein kinase, Pm004487) that could respond to the ABA complex and MAPK pathways of stress signalling in plants.

To examine the signalling pathways among PmCIPK proteins, KOALA (KEGG Orthology And Links Annotation) analysis based on the K number assignment of KEGG GENES using SSEARCH computation was carried out. The PmCIPK13 sequence was identified in the scoring scheme for K number, and the K number was K07198. The K07198 number is associated with the AMP-activated protein kinase (AMPK) signalling pathway (ko04152). The major stress signalling pathways in plants are associated with the mammalian AMPK kinase, suggesting that these pathways may have evolved from an energy-aware pathway (Zhu, 2016). The catalytic domain of PmCIPK13 was similar to that of the mammalian AMPK. To further elucidate the putative function of PmCIPKs, the predicted interactions of PmCIPKs were established based on homologous proteins in A. thaliana (Fig. 6). The interaction network showed there is no direct interaction between the CIPKs, but the interaction network links between the CIPKs through CBLs interaction and has been reported in many species. PmCIPK13 may play an important role in the interaction network, which are associated with abiotic stress tolerance.

Quantification of gene expression

Cold acclimation was found to contribute to the maximum freezing tolerance of plants. To investigate the function of PmCIPKs, the expression levels of PmCIPKs under low temperature treatment were examined by qRT-PCR experiments. The PmCIPK genes were up-regulated or down-regulated in the time course expression levels during 4 °C treatment (Fig. 7). We have compared the gene expression patterns of these CIPK genes across P. mume, A. thaliana, and O. sativa during cold stress and some CIPK genes expression levels are similar in the plants (Fig. S6). The majority of PmCIPKs were up-regulated within 4 or 6 h except for PmCIPK1, PmCIPK2, PmCIPK4, and PmCIPK7, and the highest expression levels were found at 6 h before they were down-regulated with the prolongation of cold treatment time. We found that the target genes of Pmu-miR172/167 (PmCIPK5, PmCIPK6, PmCIPK13) were up-regulated following cold treatment compared to those of the control while the expression levels of other genes (PmCIPK3, PmCIPK16) were relatively stable. Only PmCIPK7 was rapidly up-regulated at 24 h under 4 °C treatment. The expression levels of PmCIPK5, PmCIPK6, PmCIPK8, and PmCIPK13 were up-regulated and had similar changes as the PmCBF gene that has been proven to be related to cold response (Zhang et al., 2013b).