Whole-Genome Identification of Regulatory Function of CDPK Gene Families in Cold Stress Response for Prunus mume and Prunus mume var. Tortuosa

Calcium-dependent protein kinases (CDPK) are known to mediate plant growth and development and respond to various environmental changes. Here, we performed whole-genome identification of CDPK families in cultivated and wild mei (Prunus mume). We identified 14 and 17 CDPK genes in P. mume and P. mume var. Tortuosa genomes, respectively. All 270 CPDK proteins were classified into four clade, displaying frequent homologies between these two genomes and those of other Rosaceae species. Exon/intron structure, motif and synteny blocks were conserved between P. mume and P. mume var. Tortuosa. The interaction network revealed all PmCDPK and PmvCDPK proteins is interacted with respiratory burst oxidase homologs (RBOHs) and mitogen-activated protein kinase (MAPK). RNA-seq data analysis of cold experiments show that cis-acting elements in the PmCDPK genes, especially PmCDPK14, are associated with cold hardiness. Our results provide and broad insights into CDPK gene families in mei and their role in modulating cold stress response in plants.


Introduction
Calcium ion (Ca 2+ ) is a universal secondary messenger of signal transduction in many plant physiological processes, especially in plant response to environmental stimuli and plant growth and development [1]. When plants are subjected to various environmental stresses, the changes of transient fluctuations in cytoplasmic Ca 2+ concentration are perceived by different calcineurin sensors, including calcium-dependent protein kinases (CDPKs) [2], calmodulins (CaMs), CaM-like proteins and calcineurin B-like proteins (CBLs) [3][4][5]. Accordingly, Ca 2+ signals can be activated and convert to downstream targets. As a vital sensor-transducer molecule, a typical CDPK harbors four domains, consisting of a variable N-terminal domain with myristoylation or palmitoylation sites [6], an inhibitory-junction domain [7], a Ser/Thr protein kinase domain and a calmodulin-like domain constructed by four EF-hands [8]. Therefore, these special functional domains allow CDPKs to have a dual function of a Ca 2+ sensor and responder that can directly convert upstream Ca 2+ signals to target genes through phosphorylation [9].

Structural Characterizations and Protein Domain Conservation Analysis of PmCDPK and PmvCDPK Gene Families
The evolutionary analysis further classified the 31 PmCDPK and PmvCDPK genes into four groups according to conserved domains among the proteins (Figure 2A). MEME software was further used to predict conserved domains of PmCDPK and PmvCDPK proteins. A total of 18 motif were identified and all PmCDPK and PmvCDPK proteins structural domains were conserved and contain motifs 1-4 and 7. PmCDPK13 and PmvCDPK16 (in clade I), PmCDPK5 and PmvCDPK6 (in clade II) specifically contained motifs 14 and 17 in the N-terminal of the proteins. Motif 18 was specifically in the proteins of clade I expect for PmCDPK9, PmvCDPK11, PmCDPK13 and PmvCDPK16. Motif 16 was specifically in the PmCDPK10 and PmvCDPK12 proteins of clade I. Motif 13 was specifically in the proteins of clade III expect for PmCDPK12 ( Figure 2B and Figure S1). The exon-introns organization can indicate the evolutionary relationships within PmCDPK and PmvCDPK families. Phylogenetic tree analysis showed that the PmCDPK and PmvCDPK genes of each subfamily had similar gene structures. In Clade I, a total of 6 PmCDPK and PmvCDPK genes contained 7 exons, PmCDPK8 and PmvCDPK9 contained one exons each, while PmCDPK9 and PmvCDPK16 contained 8 exons each. There were 8 PmCDPK and PmvCDPK genes with 8 exons in clades II. In Clade III, 3 CDPK genes contained 9 exons, 4 CDPK genes had 8 exons, PmCDPK1 and PmvCDPK1 had 7 exon each, while PmCDPK12 only contained 6 exon. PmvCDPK10 was identified with 7 exon in clade IV. However, PmCDPK14 and PmvCDPK15 (in clade IV) had the largest number of exons as 12 ( Figure 2C). The majority of PmCDPK and PmvCDPK genes have comparable numbers of exons and introns, according to these findings, potentially indicating common genetic evolution.

Structural Characterizations and Protein Domain Conservation Analysis of PmCDPK and PmvCDPK Gene Families
The evolutionary analysis further classified the 31 PmCDPK and PmvCDPK genes into four groups according to conserved domains among the proteins (Figure 2A). MEME software was further used to predict conserved domains of PmCDPK and PmvCDPK proteins. A total of 18 motif were identified and all PmCDPK and PmvCDPK proteins structural domains were conserved and contain motifs 1-4 and 7. PmCDPK13 and PmvCDPK16 (in clade I), PmCDPK5 and PmvCDPK6 (in clade II) specifically contained motifs 14 and 17 in the N-terminal of the proteins. Motif 18 was specifically in the proteins of clade I expect for PmCDPK9, PmvCDPK11, PmCDPK13 and PmvCDPK16. Motif 16 was specifically in the PmCDPK10 and PmvCDPK12 proteins of clade I. Motif 13 was specifically in the proteins of clade III expect for PmCDPK12 ( Figures 2B and S1). The exonintrons organization can indicate the evolutionary relationships within PmCDPK and PmvCDPK families. Phylogenetic tree analysis showed that the PmCDPK and PmvCDPK genes of each subfamily had similar gene structures. In Clade I, a total of 6 PmCDPK and PmvCDPK genes contained 7 exons, PmCDPK8 and PmvCDPK9 contained one exons each while PmCDPK9 and PmvCDPK16 contained 8 exons each. There were 8 PmCDPK and PmvCDPK genes with 8 exons in clades II. In Clade III, 3 CDPK genes contained 9 exons 4 CDPK genes had 8 exons, PmCDPK1 and PmvCDPK1 had 7 exon each, while PmCDPK12 only contained 6 exon. PmvCDPK10 was identified with 7 exon in clade IV. However PmCDPK14 and PmvCDPK15 (in clade IV) had the largest number of exons as 12 ( Figure  2C). The majority of PmCDPK and PmvCDPK genes have comparable numbers of exons and introns, according to these findings, potentially indicating common genetic evolution In order to gain insight into the three-dimensional (3D) protein structure of CDPKs in P. mumu and P. mume var. tortuosa, homology modeling was performed on all PmCDPKs and PmvCDPKs. All protein structures were modeled with 100% confidence using template proteins. 31 proteins were modeled with coverage ranging from 66% (PmvCDPK16) to 98% (PmCDPK4) ( Table S3). The α-helix contributed to 37-52% of the protein structure, whereas the β-strands varied from 7 to 24% (Table S3). All CDPK proteins in P. mumu and P. mume var. tortuosa contained catalytic kinase domains, an inhibitory junction domain (JD) and the N-lobe and C-lobe, respectively with each lobe containing two EF hand motifs ( Figure 3). By transmembrane structures analysis, all the PmCDPK and PmvCDPK proteins did not have the transmembrane domain except for PmCDPK3 ( Figure S2).   In order to gain insight into the three-dimensional (3D) protein structure of CDPKs in P. mumu and P. mume var. tortuosa, homology modeling was performed on all PmCDPKs and PmvCDPKs. All protein structures were modeled with 100% confidence using template proteins. 31 proteins were modeled with coverage ranging from 66% (PmvCDPK16) to 98% (PmCDPK4) ( Table S3). The α-helix contributed to 37-52% of the protein structure, whereas the β-strands varied from 7 to 24% (Table S3). All CDPK proteins in P. mumu and P. mume var. tortuosa contained catalytic kinase domains, an inhibitory junction domain (JD) and the N-lobe and C-lobe, respectively with each lobe containing two EF hand motifs ( Figure 3). By transmembrane structures analysis, all the PmCDPK and PmvCDPK proteins did not have the transmembrane domain except for PmCDPK3 ( Figure S2).

Chromosome Localization of PmCDPK and PmvCDPK Gene Families
Chromosomal location with gene density analyses revealed that 14 PmCDPK and 17 PmvCDPK genes were unevenly distributed across chromosomes ( Figure 4). In P. mume, CDPK genes were distributed on only seven chromosomes. Chr1, Chr4, Chr5 each carried one CDPK gene, Chr2 and Chr6 possessed two genes each, whereas Chr8 and Chr3 possessed three and four genes, respectively. In P. mume var. Tortuosa, 17 CDPK genes were distributed on only 8 chromosomes. Chr1, Chr4, Chr7 each carried one CDPK gene, Chr5 and Chr6 possessed two genes each and three genes were mapped on Chr2 and Chr8, respectively. Four CDPKs were located on Chr3.

Chromosome Localization of PmCDPK and PmvCDPK Gene Families
Chromosomal location with gene density analyses revealed that 14 PmCDPK and 17 PmvCDPK genes were unevenly distributed across chromosomes ( Figure 4). In P. mume, CDPK genes were distributed on only seven chromosomes. Chr1, Chr4, Chr5 each carried one CDPK gene, Chr2 and Chr6 possessed two genes each, whereas Chr8 and Chr3 possessed three and four genes, respectively. In P. mume var. Tortuosa, 17 CDPK genes were distributed on only 8 chromosomes. Chr1, Chr4, Chr7 each carried one CDPK gene, Chr5 and Chr6 possessed two genes each and three genes were mapped on Chr2 and Chr8, respectively. Four CDPKs were located on Chr3.

Gene Duplication and Synteny Analysis PmCDPK and PmvCDPK Gene Families
To examine the expansion patterns of the CDPK gene family, we analyzed the duplicated events within the P. mume and P. mume var. Tortuosa genomes ( Figure 5). We identified 4 gene pairs of CDPKs derived from segmental gene duplications in P. mume (PmCDPK3/PmCDPK9, PmCDPK4/PmCDPK5, PmCDPK6/PmCDPK7 and PmCDPK7/PmCDPK9) and 3 pairs in P. mume var. Tortuosa (PmvCDPK5/PmvCDPK6, PmvCDPK8/PmvCDPK10, PmvCDPK11/PmvCDPK13), but did not find tandem duplication events. Surprisingly, 17 pairs of collinearity relationships between PmCDPKs and PmvCDPKs were detected. There was a high collinearity relationship among most of the P. mume and P. mume var. Tortuosa chromosomes, suggesting that most of the CDPK genes in both species had a similar origin and evolutionary process. To further explore the specific retention of PmCDPK and PmvCDPK genes, their collinearity analyses included AtCDPKs, PpCDPKs, MdCDPKs and RcCDPKs were performed using the MCScanX. ( Figure 6, Table S4). A synteny of 16 P. mume chromosomes with A. thaliana while 19 of P. mume var. Tortuosa chromosomes with A. thaliana. Similarly, 16 pairs of homologous genes between P. mume and P. persica, while 17 between P. mume var. Tortuosa and P. persica. A total of 21 and 25 homologous genes pairs were found between M. domestica and P. mume and P. mume var. Tortuosa, respectively. 12 gene pairs were detected in P. mume and R. chinensis. The result was the same in P. mume var. Tortuosa. The collinear complexity of P. mume and P. mume var. Tortuosa with M. domestica was much higher than that other species, indicating that P. mume and P. mume var. Tortuosa was more closely related to M. domestica.

Gene Duplication and Synteny Analysis PmCDPK and PmvCDPK Gene Families
To examine the expansion patterns of the CDPK gene family, we analyzed the duplicated events within the P. mume and P. mume var. Tortuosa genomes ( Figure 5). We identified 4 gene pairs of CDPKs derived from segmental gene duplications in P. mume (Pm-CDPK3/PmCDPK9, PmCDPK4/PmCDPK5, PmCDPK6/PmCDPK7 and PmCDPK7/PmCDPK9) and 3 pairs in P. mume var. Tortuosa (PmvCDPK5/PmvCDPK6, PmvCDPK8/PmvCDPK10, PmvCDPK11/PmvCDPK13), but did not find tandem duplication events. Surprisingly, 17 pairs of collinearity relationships between PmCDPKs and PmvCDPKs were detected. There was a high collinearity relationship among most of the P. mume and P. mume var. Tortuosa chromosomes, suggesting that most of the CDPK genes in both species had a similar origin and evolutionary process. To further explore the specific retention of PmCDPK and PmvCDPK genes, their collinearity analyses included AtCDPKs, PpCDPKs, MdCDPKs and RcCDPKs were performed using the MCScanX. ( Figure 6, Table S4). A synteny of 16 P. mume chromosomes with A. thaliana while 19 of P. mume var. Tortuosa chromosomes with A. thaliana. Similarly, 16 pairs of homologous genes between P. mume and P. persica, while 17 between P. mume var. Tortuosa and P. persica. A total of 21 and 25 homologous genes pairs were found between M. domestica and P. mume and P. mume var. Tortuosa, respectively. 12 gene pairs were detected in P. mume and R. chinensis. The result was the same in P. mume var. Tortuosa. The collinear complexity of P. mume and P. mume var. Tortuosa with M. domestica was much higher than that other species, indicating that P. mume and P. mume var. Tortuosa was more closely related to M. domestica.

Promoter Analysis of PmCDPK and PmvCDPK Gene Families
Identification of promoter binding sites helps to understand the function of genes. In this study, the key components of the PmCDPK and PmvCDPK genes include core promoter elements ( Figure 7A) and plant-inducible promoter elements including light response elements, stress response elements and hormone response elements. Transcription-associated cis element (CAAT-box and TATA-box) were conserved in all PmCDPK and PmvCDPK genes. Among them, the types and numbers of light-responsive elements were the largest, for example, Box 4 was present in the promoter regions of each CDPK gene, except for PmCDPK11 and PmvCDPK13. Some PmCDPK and PmvCDPK gene family had unique light-responsive elements, such as GA-motif, AT1-motif, LAMP-element, Gap-box, 3-AF3 binding site, Box II, chs-CMA1a, which were found only in one CDPK gene in P. mume and P.mume var. Tortuosa respectively ( Figure 7B). Notably, the PmCDPK and PmvCDPK promoters contained 11 elements involved in phytohormone response and stress response ( Figure 7C,D, Table S5). Among these five types of hormone-responsive elements, 86% of the PmCDPKs and 82% of the PmvCDPKs contained MeJA response elements (CGTCA-motif, TGACG-motif), making it the most frequent hormone-responsive element, the number of abscisic response element is second. 93% of the PmCDPKs and 94% of the PmvCDPKs contained anaerobically inducible elements (AREs). In addition, 36% of the PmCDPKs and 29% of the PmvCDPKs contained low-temperature responsive binding site (LTR): PmCDPK1, PmvCDPK1 and PmvCDPK16 had 2 LTRs each, PmCDPK3/4/12/13 and PmvCDPK3/4/17 had only one LTR each. All CDPKs contained phytohormone response and stress response elements, suggesting that PmCDPK and PmvCDPK genes containing these elements may be induced to be express by stresses and plant hormones.

Promoter Analysis of PmCDPK and PmvCDPK Gene Families
Identification of promoter binding sites helps to understand the function of genes. In this study, the key components of the PmCDPK and PmvCDPK genes include core promoter elements ( Figure 7A) and plant-inducible promoter elements including light response elements, stress response elements and hormone response elements. Transcriptionassociated cis element (CAAT-box and TATA-box) were conserved in all PmCDPK and PmvCDPK genes. Among them, the types and numbers of light-responsive elements were the largest, for example, Box 4 was present in the promoter regions of each CDPK gene, except for PmCDPK11 and PmvCDPK13. Some PmCDPK and PmvCDPK gene family had unique light-responsive elements, such as GA-motif, AT1-motif, LAMP-element, Gap-box, 3-AF3 binding site, Box II, chs-CMA1a, which were found only in one CDPK gene in P. mume and P. mume var. Tortuosa respectively ( Figure 7B). Notably, the PmCDPK and PmvCDPK promoters contained 11 elements involved in phytohormone response and stress response ( Figure 7C,D, Table S5). Among these five types of hormone-responsive elements, 86% of the PmCDPKs and 82% of the PmvCDPKs contained MeJA response elements (CGTCA-motif, TGACG-motif), making it the most frequent hormone-responsive element, the number of abscisic response element is second. 93% of the PmCDPKs and 94% of the PmvCDPKs contained anaerobically inducible elements (AREs). In addition, 36% of the PmCDPKs and 29% of the PmvCDPKs contained low-temperature responsive binding site (LTR): PmCDPK1, PmvCDPK1 and PmvCDPK16 had 2 LTRs each, PmCDPK3/4/12/13 and PmvCDPK3/4/17 had only one LTR each. All CDPKs contained phytohormone response and stress response elements, suggesting that PmCDPK and PmvCDPK genes containing these elements may be induced to be express by stresses and plant hormones.

Expression Pattern of PmCDPKs in Different Underground and Aerial Tissues
As illustrated in Figure 9A, all the PmCDPKs genes were differentially expressed in flower buds, fruits, leaves, roots and stems. Among them, PmCDPK5/9/12 were expressed in roots, PmCDPK7 showed higher expression levels in leaves, PmCDPK3/6 presented relatively higher expression levels in bud, but their expression levels were low in other tissues. All PmCDPKs were expressed during the flower bud dormancy period and also at specific stages of development ( Figure 9A, Table S6). Six genes (PmCDPK1/7/8/10/12/14) showed the high level of expression from EDI to EDIII, then abruptly decreased in the NF stage. The expressions of six genes (PmCDPK2/3/4/5/6/9) were exhibited specifically high level in the NF stage. PmCDPK11 and PmCDPK13 showed the highest level of expression in the EDI and EDII, respectively ( Figure 9B, Table S7).
Expression profiles in stems of cold-insensitive cultivar 'Songchun' were analyzed under three low temperatures at three test sites: Beijing (BJ), Chifeng (CF), and Gongzhuling (GZL) (Figure 10, Table S8). The expression of PmCDPK6 was not detected and PmCDPK3 expression was low in RNA-seq. At the Beijing site, a total of seven genes (PmCDPK2/12/14/7/9/4/10) showed up-regulation in autumn (7.3 °C) and winter (−5.4 °C) but then decreased in spring (3.2 °C). PmCDPK1 and PmCDPK8 showed higher expression in spring, the expression of PmCDPK5, PmCDPK11 and PmCDPK13 was up-regulated in autumn but down-regulated in winter and spring. The results showed that the gene expression trends at CZL and CF sites were relatively consistent with that in BJ. In general, PmCDPK12/14/7/9/4/10 genes were up-regulated in autumn and winter, then down-regulated in the early-spring in three sites, while PmCDPK1 and PmCDPK8 expression showed the opposite trend, the expression levels of the PmCDPK2/5/11/13 genes showed different trends with the different test sites, indicating that the expression patterns of these genes were different in response to low temperature ( Figure 10A). As shown in Figure 10B, for the same season, the expression of most PmCDPKs at different test sites was slightly different, among which PmCDPK3 was highly expressed at Chifeng site in spring, while PmCDPK9 was highly expressed at GZL site in winter. More than a half of PmCDPKs genes showed a relatively low expression in expression in spring, while PmCDPK1 and PmCDPK8 was highly expressed at Beijing site in spring.

Expression Pattern of PmCDPKs in Different Underground and Aerial Tissues
As illustrated in Figure 9A, all the PmCDPKs genes were differentially expressed in flower buds, fruits, leaves, roots and stems. Among them, PmCDPK5/9/12 were expressed in roots, PmCDPK7 showed higher expression levels in leaves, PmCDPK3/6 presented relatively higher expression levels in bud, but their expression levels were low in other tissues. All PmCDPKs were expressed during the flower bud dormancy period and also at specific stages of development ( Figure 9A, Table S6). Six genes (PmCDPK1/7/8/10/12/14) showed the high level of expression from EDI to EDIII, then abruptly decreased in the NF stage. The expressions of six genes (PmCDPK2/3/4/5/6/9) were exhibited specifically high level in the NF stage. PmCDPK11 and PmCDPK13 showed the highest level of expression in the EDI and EDII, respectively ( Figure 9B, Table S7).  Expression profiles in stems of cold-insensitive cultivar 'Songchun' were analyzed under three low temperatures at three test sites: Beijing (BJ), Chifeng (CF), and Gongzhuling (GZL) (Figure 10, Table S8). The expression of PmCDPK6 was not detected and Pm-CDPK3 expression was low in RNA-seq. At the Beijing site, a total of seven genes (Pm-CDPK2/12/14/7/9/4/10) showed up-regulation in autumn (7.3 • C) and winter (−5.4 • C) but then decreased in spring (3.2 • C). PmCDPK1 and PmCDPK8 showed higher expression in spring, the expression of PmCDPK5, PmCDPK11 and PmCDPK13 was up-regulated in autumn but down-regulated in winter and spring. The results showed that the gene expression trends at CZL and CF sites were relatively consistent with that in BJ. In general, PmCDPK12/14/7/9/4/10 genes were up-regulated in autumn and winter, then downregulated in the early-spring in three sites, while PmCDPK1 and PmCDPK8 expression showed the opposite trend, the expression levels of the PmCDPK2/5/11/13 genes showed different trends with the different test sites, indicating that the expression patterns of these genes were different in response to low temperature ( Figure 10A). As shown in Figure 10B, for the same season, the expression of most PmCDPKs at different test sites was slightly different, among which PmCDPK3 was highly expressed at Chifeng site in spring, while PmCDPK9 was highly expressed at GZL site in winter. More than a half of PmCDPKs genes showed a relatively low expression in expression in spring, while PmCDPK1 and PmCDPK8 was highly expressed at Beijing site in spring.

RT-PCR of PmCDPKs under Cold Treatment
To further investigate the key P. mume CDPK genes involved in cold stress, all identified PmCDPK genes that showed a higher expression level in 4 • C low temperature treatment, were selected for further RT-qPCR analysis in cold-tolerant cultivar 'Meirenmei' and coldsensitive cultivar 'Jinsheng'. It is well known that lysis curves can test the specificity of amplification reactions. Therefore, we determined the suitable annealing temperature for the primers based on the shape of the single peak of the solubility curve and the Tm value. Results from gel electrophoresis plots of PCR amplification products showed that each primer showed a single band, indicating good specificity for use in subsequent qRT-PCR experiments ( Figure S3). There were the PmCDPK3, PmCDPK6 and PmCDPK13 that were not detected, a finding that is consistent with the transcriptome data. A total of 12 PmCDPK members showed significant expression profiles differed during the cold stress treatment period in the two cultivars. PmCDPK1/8/11 changed only slightly in the two cultivars, except that the expression level changed significantly at some point in the cold stress. PmCDPK2 reached its peak at 12 h ('Jinsheng') and 24 h ('Meirenmei'), respectively. The expression level of PmCDPK5 in 'Jinsheng' was higher than that in 'Meirenmei', whereas PmCDPK4 the genes were highly expressed in the 'Meirenmei', by contrast, repressed in the 'Jinsheng'. The expression levels of these genes (PmCDPK9/10/12/14) were significantly increased in the early stage of stress and decreased in the late stage of stress in 'Meirenmei', but the changes were not significant in 'Jinsheng'.

Discussion
A type of serine/threonine-protein kinase known as CDPK, which plays a pivotal role in controlling plant growth and maturation. However, to date, detailed or complete information of CDPKs in mei had not been obtained. Increasing effort has been made on the genomic research of mei, high quality assemblies of the P. mume and P. mume var. Tortuosa genomes provided powerful genomic information to explore and identify CDPKs for cultivated and wild mei [40,41]. Among the long genetic improvement history of mei, cold is one of the main constraints that seriously affected the transplantation of mei from the south to the north. When plants are exposed to cold, it leads to various changes in physiology and gene expression patterns [42]. CDPKs have been known for many years to participate in Ca 2+ -related signal transduction induced by cold stress, particularly, isoforms were implicated in ABA-mediated signaling [43]. CDPK-dependent changes in ion flux, metabolic changes, or alterations in gene expression [44]. It has been reported that plants overexpressing some CDPK genes of A.thaliana [7], rice [27], and other plant species were shown to exhibit enhanced cold tolerance [45]. Some CDPK substrates have been directely linked to abiotic stress response [46], In view of the research progress on CDPK members involved in plant cold stress responses, the identification of CDPK members in cultivated and wild mei has important implications of the selection of cold tolerant mei plants.
On this basis, we identified 14 and 17 CDPKs in P. mume and P. mume var. Tortuosa at the genome level, respectively, which contained the PF00069 and PF13499 domains. Comparative analysis of both species showed that some CDPK genes are present in P. mume var. Tortuosa but absent in P. mume (Table 1). For example, P. mume var. Tortuosa had one CDPK gene each on Ch2/5/7, compared with P. mume. In addition, this study also identified a total of 205 CDPK genes that have been identified in 11 other species, besides the number of CDPKs in these plants varied considerably. Such a wide range of occurrence may be due to the high degree of variation at duplicated genome and ploidy level in plants. Phylogenetic analysis of Rosaceae plants suggested similarity in CDPK genes. The CDPK proteins of mei were distributed into four groups through protein sequence analyses (Figure 1), which is consistent with reports of A. thaliana, rice, maize [47]. Furthermore, PmCDPK13/PmvCDPK16, PmCDPK8/PmvCDPK9, PmCDPK11/PmvCDPK13, PmCDPK10/PmvCDPK12, PmCDPK6/PmvCDPK8, PmCDPK2/PmvCDPK2, PmCDPK1/ PmvCDPK1 and PmCDPK4/PmvCDPK15 showed close associations with one another (Figure 1). Most CDPKs have acylation sites, such as N-myristoylation sites and Spalmitoylation sites, which are thought to be key biological processes affecting a variety of cellular functions through modulation of membrane targeting [48]. The variable N-terminal domain of some CDPKs contain N-myristoylation and palmitoylation, five PmCDPKs and seven PmvCDPKs have N-myristoylation, ten PmCDPK and thirteen PmvCDPKs, which have palmitoylation motifs in P. mume and P. mume var. Tortuosa (Table 1). The N-myristoylation is the binding site of protein and membrane, while palmitoylation plays an important role in membrane binding stability.
Phylogenetic tree analysis showed that CDPK genes of P. mume and P. mume var. Tortuosa were close to each other in every branch compared to other species, suggesting a high degree of similarity between cultivated and wild mei, as previously observed in sweet potato [49]. The phylogenetic group members of PmCDPKs and PmvCDPKs had relatively conserved gene structure with only a few exceptions, which are similar to those reports for other gene families in mei, such as GASA [50], SWEET [51], AUXIN/INDOLE ACETIC ACIDs (Aux/IAAs) [52]. In P. mume and P. mume var. Tortuosa, there are overall differences in the number, length and position of introns, as well as CDS in CDPK genes. These variations may result in different lengths of CDPKs genes in cultivated and wild mei. There was little variation in the motif composition of CDPKs in both species. Similar protein 3D structural features were obtained in PmCDPK and PmvCDPK (Figure 3), with the typical EF-hand and protein kinase domain motifs (Figure 2), which were also observed in other species, such as A. thaliana [11], Solanum habrochaites [53], rice, maize, sorghum and ginger (Zingiber officinale) [54,55]. This suggests that the CDPK protein structure had conserved function and structural mechanism. The uneven distribution of CDPKs on chromosomes is a common feature in Rosaceae (Figure 4) and is similar in strawberry [19] and peach [20], which may be related to species evolution and genetic variation. These findings indicate a close relationship and high level of sequence similarity between P. mume and P. mume var. Tortuosa, suggesting that despite significant differences, they are quite similar at the genomic level, with sequences conserved from wild to cultivated species.
During evolution, duplication of genes is considered to be the main factor in the evolution and diversification of angiosperms. Differences between species become more pronounced especially when species are subjected to selection pressure and restrictive growth conditions [56]. No tandem duplication involving PmCDPK and PmvCDPK genes was discovered, but 4 segmental duplications and 3 segmental duplications were found in P. mume and P. mume var. Tortuosa, respectively ( Figure 5). Therefore, segmental duplication is the main driving force of duplication for CDPK gene families in two species. It has been reported that segmental duplication as a major factor leading to the expansion of the CDPKs in other species, i.e., grapevine [57], tomato [53]. These duplicated genes may decrease the probability of extinction, as well as contribute to improved stress tolerance in plants [58]. In addition, 18 pairs of collinear genes were detected between PmCDPK and PmvCDPK, and there was a high collinearity relationship among P. mume and P. mume var. Tortuosa chromosomes, suggesting that most of the CDPK genes in both species had a similar origin and evolutionary process ( Figure 5). Compared to that with A. thaliana and other Rosaceae plants, the collinearity between apple and mei showed higher complexity in CDPK gene family, which is a logical reason for the relatively high number of MdCDPKs (Figure 6), indicating that as an ancient tetraploid plant recent whole-genome duplication events and the increase of chromosomes have led to a significant expansion of the apple CDPKs family.
PlantCARE analysis revealed the presence of various cis-acting elements of stress, growth, light, phytohormones and development-related elements in the promoters of PmCDPK and PmvCDPK genes. These cis-acting elements can help analyze gene function by predicting the binding sites of transcription factors in the promoter regions of genes. In this study, CDPK genes family possesses two key promoter elements including TATA box and CAAT box ( Figure 7A). TATA box connects transcription initiation with RNA polymerase and CAAT box regulates gene transcription efficiency. The results indicated that PmCDPK and PmvCDPK can be transcribed normally. Many CDPKs have been identified in different species in relation to plant cold resistance. A total of 10 PmCDPK and PmvCDPK genes contained one or more low-temperature responsive cis-elements, suggesting that these genes may be involved in the regulation of low-temperature stress. Moreover, the promoter regions of CDPK members contain different stress-responsive cis-acting elements, suggesting that different PmCDPK and PmvCDPK memebers may play important roles in different stresses ( Figure 7B-D).
Our interaction network predicted that PmCDPKs and PmvCDPKs were associated with RBOHB, RBOHC and RBOHD. RBOHs are integral membrane proteins that produce superoxide anions (·O 2− ) and subsequently promote the production of ROS, which play a virtual role in the cold response of plants [35]. FvRBOHD can be rapidly induced to be expressed in response to cold stress in strawberry [59]. CDPKs have been shown to phosphorylate the N-terminal fragment of RBOHs to participate in signal transmission in defense responses. Other predicted interacting proteins from PmCDPKs and PmvCDPKs network analysis include PbrMAPKs and PtrXTH proteins increased the cold resistance in pears (Pyrus × bretschneideri) and poplar (Populus simonii × Populus nigra), respectively [36,37].
To further validate the function of CDPK genes, we investigated the expression pattern of CDPK genes under cold stress conditions using RNA-seq data (previously published) and RT-qPCR analysis. In this study, the expression pattern at different tissues in P. mume showed tissue-specific expression of CDPKs genes, whereas some PmCDPK genes (Pm-CDPK5/10/12/14) showed enhanced expression in all tissues ( Figure 8A), indicating that they might have a wide range of regulatory functions. Similar tissue-specific expression patterns of CDPK genes were also observed in grape and strawberry [57,60]. PmCDPK3 and PmCDPK6 were only expressed in the floral buds, not detected in any other tissue, and they may be involved in regulating bud development, dormancy and so on ( Figure 9B). The expression levels of the PmCDPK gene family members were up-regulated or down-regulated in different seasons and locations ( Figure 10), suggesting that CDPK gene members play different roles in response to cold stress. Many CDPK genes have been identified to be closely associated with cold stress. In rice, OsCDPK13 and OsCDPK24 are positive regulators of cold stress tolerance, whereas OsCPK17 gene expression reduces cold tolerance but dose not affect the expression of key cold stress-inducible genes [27,61,62]. In grapevine, overexpressing the VaCPK20 gene increases the expression of stress-responsive genes such as COR47, NHX1, KIN1, or ABF3 to improve cold resistance [57,63]. In ripe strawberry fruit, FaCDPK1 transcript levels were increased in response to low temperature [60]. In cucumber and melon [4], almost all CsCDPKs and CmCDPKs were up-regulated under cold stress [4,64]. In this study, the PmCDPK14, which was a homologue of CDPK28 in A. thaliana, was highly expressed in both flower buds and stems at different sites under cold stress. In addition, qRT-PCR analysis showed that PmCDPK14 gene reached its peak value at 4h of cold tolerance cultivar 'Meirenmei', which was 8.06 times that of the control, and reached its peak value at 12 h of cold-sensitive cultivar 'Jinsheng', which was 4.7 times that of the control (Figure 11). AtCDPK28 is activated rapidly upon cold shock and then phosphorylates, leading to promote the nuclear translocation of NIN-LIKE PROTEIN 7 (NLP7), which specifies transcriptional reprogramming of cold-responsive gene sets in response to Ca 2+ [65]. These results suggest that PmCDPK14 involved in regulating cold-related genes promotes cold tolerance in 'Meirenmei'. Previous studies have shown that PeCPK10 was rapidly induced when the transgenic lines were subjected to −4/−8 • C for 8 h and showed enhanced freezing tolerance [45]. The observed gene expression patterns showed that some PmCDPKs (PmCDPK7/9/10/14) were more highly expressed in winter than in autumn and spring. Through gene expression patterns analysis, we speculated that these genes played a role in cold stress and freezing stress.
qRT-PCR analysis indicated that three PmCDPK genes (PmCDPK3/6/13) were considered not to be expressed, which was consistent with the transcription data. PmCDPK genes was transcriptionally activated at 4°C low temperature and increased or decreased in expression with the extension of treatment time. PmCDPK2/5/7/9/10/12/14 were up-regulated, PmCDPK1 was negatively regulated under cold stress, suggesting that these genes might be positively or negatively regulated by cold stress. The discrepancy in expression patterns between PmCDPK4/8/11 is potentially due to genetic differences. qRT-PCR analysis indicated that three PmCDPK genes (PmCDPK3/6/13) were considered not to be expressed, which was consistent with the transcription data. PmCDPK genes was transcriptionally activated at 4℃ low temperature and increased or decreased in expression with the extension of treatment time. PmCDPK2/5/7/9/10/12/14 were up-regulated, PmCDPK1 was negatively regulated under cold stress, suggesting that these genes might be positively or negatively regulated by cold stress. The discrepancy in expression patterns between PmCDPK4/8/11 is potentially due to genetic differences.

Plants Genome Resources and Identification of CDPK Genes
The reliable genome assemblies of the P. mume were obtained from the Prunus mume genome database (https://github.com/lileiting/prunusmumegenome, accessed on 5 September 2022 The sequences of 34 AtCDPKs and 31 OsCDPKs were searched from the Arabidopsis Information Resource (TAIR, https://www.arabidopsis.org, accessed on 25 September 2022) and the rice genome annotation database (http://rice.uga.edu/, accessed on 28 September 2022). Theses initial sequences was used as query sequences to identify candidate CDPK genes by using the local library software BLAST-P [66] with an E-value of 1 × 10 −7 . Then, all non redundant sequences were verified the sequence by HMMER software (version 3.0, http://hmmer.org/, accessed on 30 September 2022) and the E-value was set up less than 1 × 10 −7 . Finally, NCBI-CDD database (https://www.ncbi.nlm.nih.gov/cdd, accessed on 2 October 2022) and SMART database (http://www.smart.embl-heidelberg. de/, accessed on 10 October 2022) was used to scan and delete the proteins without domains. ExPASy (https://www.expasy.org/, accessed on 15 October 2022) were used to predict the MW, PI and the N-terminal myristoylation sites of PmCDPKs and PmvCDPKs. GPS-Palm (http://gpspalm.biocuckoo.cn/, accessed on 17 October 2022) and WoLF PSORT (https: //wolfpsort.hgc.jp/, accessed on 23 October 2022) was used to predict palmitoylation sites and subcellular localization, respectively [48]. The 3D structures of all PmCDPK and PmvCDPK proteins were predicted by homology modeling using the PHYRE2 web portal (http://www.sbg.bio.ic.ac.uk/phyre2, accessed on 28 October 2022). PHYRE2 uses advanced remote homology detection methods to construct 3D models of protein sequences, with a single highest score template model for all proteins and 100% model confidence.

Protein Interaction Network of CDPKs
The protein interaction network of PmCDPK and PmvCDPK was predicted by STRING (v 11.5, https://cn.string-db.org/, accessed on 15 November 2022) and was visualized using Cytoscape [70].

Transcriptome Data Analysis
RNA-seq datasets of PmCDPKs in five different tissue (flower buds, fruits, leaves, roots, and stems) [40] and flower bud dormancy of P. mume ('Zaolve') from November to February [71] were downloaded. The expression profiles of PmCDPK genes in the stem of P. mume cultivar 'Songchun' were analyzed in different geographical locations including Beijing (BJ, 39 •

qRT-PCR Analysis
The annual branches of the cold-tolerated cultivar 'Meirenmei' and the cold-sensitive cultivar 'Jinsheng' were used as experimental materials and were maintained in water overnight at 22 • C. Plants were treated in a light incubator at 4 • C (16-h light/8-h dark) and sampled for RNA extraction after 1, 3, 6 12, 24, 36, 48, and 72 h.
Total RNA isolation and cDNA synthesis were performed using RNAprep Pure Plant Plus Kit (Tiangen, Beijing, China) and SYBR ® Green Premix Pro Taq HS qPCR Kit (Accurate-Biology, Hunan, China). The reaction system consisted of 20 µL with a 10 µL SYBR ® Green Premix Pro Taq HS qPCR Kit, 0.4 µL forward and reverse primers mix, 7.2 µL ddH 2 O and 2 µL 10 × cDNA samples. The PCR amplification conditions was as follows: (1) 95 • C for 30 s; (2) 95 • C for 5 s, 55-60 • C for 30 s, 72 • C for 30 s for 40 cycles; (3) 72 • C for 30 s. The relative expression levels of the genes were calculated using the 2 −∆∆Ct method [73]. The primers used in this study were listed in Table S9, the specificity tests of the primers were shown in Figure S3.

Conclusions
The present research is the first in-depth and methodical report of genome-wide characterization of CDPK gene families in cultivated and wild mei. We identified 31 high confidence CDPK genes in the genomes of P. mume and P. mume var. Tortuosa, which were divided into four subgroups based on a phylogenetic analysis. Gene motifs, structure, chromosomal location and cis-acting elements were analyzed to explain the various and traits of the identified CDPK genes. Duplication events occurred in both P. mume and P. mume var. Tortuosa genomes were identified by collinearity analysis. Importantly, RNAseq data in five tissues and geographic locations revealed some tissue-specific expression and significant up-regulation of cold responsive CDPK genes. Different expression trends were also observed following RT-qPCR under cold stress, indicating the important role of PmCDPKs in the cold stress signaling pathway. Ultimately, the knowledge gained from this research will contribute to breed cold-tolerance cultivars of mei.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants12132548/s1, Figure S1: Sequences of the 18 motifs; Figure S2: The transmembrane structures of PmCDPK and PmvCDPK proteins. Table S1: Information for the proteins used in the present study. Table S2: The specific number of genes in the Clades used in the present study. Table S3: Details of various structural features in 3-dimensional structure of PmCDPK and PmvCDPK Proteins. Table S4. Duplication events between P. mume and P. mume var. Tortuosa and A. thaliana, P. persica, M. domestica and R. chinensis. respectively. Table S5: Analysis of cis acting elements of CDPK gene family in P. mume and P. mume var. Tortuosa. Table S6: Expression profiles of PmCDPK genes in five different tissues(root, stem, leaf, bud and fruit) (RPKM). Table S7: Expression profiles of PmCDPK genes during the process of flower bud dormancy release (RPKM). Table S8: Expression profiles of PmCDPK genes in different regions and seasons (FPKM). Table S9 The primer list used in this study.