In this study, genome-wide analysis of CDPK-SnRK gene family has been performed on the basis of the completed B. rapa genome sequence (Wang et al., 2011). Based on previously reported methods (Harmon et al., 2001; Hrabak et al., 2003), the homogeneous candidate CDPK-SnRK genes between Brassica rapa and other species were identified by BLASTP (Supplementary Table 1). Subsequently, all candidate protein sequences were subjected to Pfam and SMART analyses. Finally, we identified 49 BrCPKs, 14 BrCRKs, 3 BrPPCKs, 5 BrPEPRKs, and 56 BrSnRKs named according to nomenclature proposed for CDPK-SnRK genes (Supplementary Table 2).

To better understand the expansion and evolutionary history of CDPK-SnRK genes in B. rapa, genes were also identified in 16 other species representing the major clades of plants. The evolutionary relationships of the species and the number of CDPK-SnRK genes are shown in Figure 1A. The data show that Glycine max contained the highest number of CDPK-SnRK genes (200), followed by Z. mays (193) and M. truncatula (188) (Figure 1A). However, A. trichopoda, a basal angiosperm species that was the single living representative of the sister lineage to all other extant flowering plants, contained the lowest number of CDPK-SnRK genes (28) in Angiospermae. The reason is that it originated prior to the split of eudicots and monocots and has not experienced any whole genome duplication (WGD), while the other 12 angiosperms had several rounds of WGDs/triplications after their split from A. trichopoda. Furthermore, the number of CDPK-SnRK genes in algae, Bryophyta and Pteridophyta was less than that in Angiospermae. This phenomenon was also caused by several WGD events that occurred during angiosperm evolution (Figure 1B). These results indicated that the expansion of the CDPK-SnRK family from angiosperms mainly relied on large-scale DNA rearrangements, namely, WGDs. The elevated duplication frequency and increased retention of CDPK-SnRK genes also contributed to neofunctionalization and caused them to gain important functions in angiosperm development.

To investigate the extent of lineage-specific expansion of the CDPK-SnRK genes in B. rapa, phylogenetic trees were constructed using the maximum likelihood method (Figure 2A). The phylogenetic tree showed that all the CDPK-SnRK genes were clustered into five distinct gene classes (CPK, CRK, PEPRK, PPCK, SnRK) (Figure 2A), while the CPK family was divided into four groups (I, II, III, and IV) and the SnRK family was classified into three groups (SnRK1, SnRK2, and SnRK3), consistent with the reports in A. thaliana (Hrabak et al., 2003). In B. rapa, the CPK, CRK, PEPRK, PPCK, and SnRK gene families contained 49 members, 14 members, 3 members, 5 members, and 56 members, respectively, whereas in A. thaliana, the CPK, CRK, PEPRK, PPCK, and SnRK families contained 34 members, 8 members, 2 members, 2 members, and 39 members, respectively. Next, the synteny of CDPK-SnRK genes between A. thaliana and three subgenomes in B. rapa was analyzed. There were 34 CPK, 8 CRK, 2 PEPRK, 2 PPCK, and 38 SnRK genes on the conserved collinear block (Supplementary Table 5). Meanwhile, 2 CPKs, 2 CRKs, 1 PEPRK, 1 PPCK, and 3 SnRKs were retained completely; conversely, 4 CPKs and 3 SnRKs from B. rapa were lost. Due to the Brassica-specific WGT event, the gene numbers of these classes in B. rapa were greater than those in A. thaliana.

Furthermore, the different domain architectures, motif compositions and gene structures of CDPK-SnRK were analyzed (Figure 2B, Figure S1). All members of the CDPK-SnRK superfamily have a kinase domain of similar length and sequence, with the kinase domains at or near the N-terminus, then the junction domains, followed by the regulatory domains (Figure 2B). Although CPK proteins have a functional kinase domain coupled with regulatory calcium-binding EF-hands, the C-terminal domains (EF-hands) of CRK proteins contain apparently degenerate calcium-binding sites with no function. Meanwhile, 10 conserved motifs were detected in BrCPK, BrCRK, BrPEPRK, BrPPCK, and BrSnRK, respectively (Figure S3). All BrCDPK-SnRK proteins had highly conserved kinase domains, which corresponded to motifs 1-4,7,9 in BrCDPK and motifs 1-4,6,7,9 in BrSnRK, whereas motif 8 was found in BrSnRK3.1, corresponding to the NAF/FISL domains (Figures S1A,B). The amino acid sequence of BrCDPK-SnRK was aligned with AtCDPK-SnRK protein sequences from five gene classes. In CPK, CRK, PEPRK, PPCK, and SnRK, higher sequence similarities were identified in the N-terminus, which corresponded to the conserved kinase domain (Figure S2). In addition, variable gene structures of BrCDPK-SnRK were observed. As shown in Figures 3A,C, the intron numbers of the BrCPK genes ranged from 5 to 9 with a median of 6, while the BrSnRK genes ranged from 0 to 15 with a median of 7. Interestingly, we found that 22 BrSnRK3 genes have no introns. The theoretical pI of the BrCPK gene family ranged from 4 to 9 with a median of 6, but BrSnRK proteins showed a pI range from 2 to 11 and with a median of 8 (Figure 3B). Other classes had complex theoretical pI ranging in value from 4 to 9 (Figure 3B).

To investigate different retention in CPK, CRK, PEPRK, PPCK, and SnRK during B. rapa WGT events, 44/49, 14/14, 3/3, 5/5, and 50/56 were located in the syntenic regions, respectively (Figures 4A,B and Supplementary Table 3). The results demonstrated that 43% (44/102) of the CPK genes were retained in the syntenic regions, relative to 44% (50/114) of the SnRK genes. The retention rates of CRK, PEPRK and PPCK are 58% (14/24), 50% (3/6), and 83% (5/6), respectively (Figure 4F). Additionally, we counted gene copies and analyzed the distribution of the three subgenomes by comparing the retention of CPK, CRK, PEPRK, PPCK, and SnRK (Figure 4E). The result showed that all PEPRK genes had more than two copies retained, which is more than the retention of the other subfamilies (42%). However, 3% of the CPK and SnRK genes were completely lost. Next, the proportions of CPK and SnRK genes retained were higher in the least fractionated (LF) subgenome than in the medium fractionated (MF1) and most fractionated (MF2) subgenome, consistent with a previous report showing that the degree of retained genes in these three subgenomes (LF, MF1, and MF2) was decreased (Wang et al., 2011). In contrast, the PEPRK and PPCK families were retained more in the MF1 subgenome than in the LF subgenome (Figure 4F). In summary, the results confirmed that PEPRK genes were more preferentially retained than other subfamilies and that CPK genes were retained similarly to SnRK genes during diploidization following WGT in B. rapa.

Furthermore, the retention rates of four CPK (CPKI, II, III, IV) and three SnRK (SnRK1, 2, 3) groups were observed. As show in Figure 4H, 67% of CPKIIIs and SnRK1s had more than two copies retained, which is greater than the retention of the other groups (Figures 4H,D). In addition, the proportion of CPKIIIs and SnRK1s retained was higher in the LF subgenome than in the other subgenomes, which once again confirmed that CPKIII and SnRK1 genes were more preferentially retained than other groups (Figures 4C,G).

All BrCDPK-SnRK genes could be mapped onto 10 chromosomes of Chinese cabbage with a non-random distribution, except BrPPCK3, which is located in Scaffold000191 (Figure 5). On every chromosome, the proportion of BrCPK genes was similar to that of BrSnRK genes. However, Chromosome 09 contained more BrCPK genes (11 genes) than BrSnRK genes, whereas chromosomes 01 had the opposite. B. rapa shares two WGDs (WGD: α and β) and one whole-genome triplication event (WGT: γ) in its evolutionary history with Arabidopsis but has undergone an additional WGT event. Therefore, the B. rapa genome was further divided into three differentially fractionated subgenomes (LF, MF1, MF2), of which LF contained more BrCPK genes and BrSnRK genes than either of the other two subgenomes. In addition, the 24 conserved ancestral genomic blocks (labeled A–X) in the B. rapa genome were reconstructed according to previous reports (Cheng et al., 2013). The color coding of these blocks was based on their positions in a proposed ancestral karyotype (AK1-8). We also found that most of the BrCPK genes cluster together in a region of AK1 (20%), whereas BrSnRK genes belonged to AK1 (18%) and AK3 (18%).

Furthermore, the types were identified by the MCScanX program, and the divergence times of the duplicated genes were estimated by calculating the synonymous substitution rates (Ks) and non-synonymous substitution rates (Ka). In total, 70 BrCDPK-SnRK duplicated gene pairs were analyzed (Supplementary Table 4). BrCPK (74%), BrCRK (86%), BrPEPRK (100%), BrPPCK (100%), and BrSnRK (70%) duplicated gene pairs were segmental duplications (Figure 6B), and all the duplicated BrCDPK-SnRK gene pairs had Ka/Ks ratios less than 1, indicating the purifying selection of these genes (Figure 6A, Supplementary Table 4). To understand the divergence time, the Ks values of the BrCPK genes ranged from 0.3 to 0.5 and had a mean of ~0.34 (~11 Myr), while the BrSnRK genes ranged from 0.2 to 0.55 and focused on ~0.25 (~8.5 Myr; Figures 2A,B, 6A). The divergence time of BrSnRK duplicated gene pairs was 8.49 MYA, which indicates that their divergence occurred during the Brassica triplication events (5 ~ 9 MYA). The divergence times obtained for the BrCPK duplicated gene pairs ranged from 10 to 16.6 MYA, indicating that these duplications occurred during the divergence of Chinese cabbage and Arabidopsis (9.6–16.1 MYA).

To investigate the evolution of the CDPK-SnRK family in the plant kingdom, we selected 13 Angiospermae (8 eudicots, 4 monocots and one basal angiosperm), 3 Gymnospermae, 1 Pteridophyta, 1 Bryophyta, and 1 Chlorophyta species for comparative analysis (Figure S3). We constructed a phylogenetic tree of the CDPK-SnRK genes to analyze the evolutionary relationships of these species. The phylogenetic tree showed that the CDPK-SnRK genes formed five gene classes (CPK, CRK, PEPRK, PPCK, and SnRK), which is consistent with the result for B. rapa and A. thaliana. Meanwhile, we found that no CRK, PEPRK, or PPCK genes were detected in Volvox carteri. Therefore, the CRK, PEPRK, and PPCK gene families may only exist in land plants. To further determine the relationship among the five groups, the genetic distance was analyzed (Figures S4A,B). The box plot shows the genetic distance of CPK vs. CRK, PEPRK, and PPCK, which was smaller than SnRK vs. these groups (Figure S4B). Notably, the genetic distance between CPK and CRK was smaller than that between CPK and PEPRK or CPK and PPCK. These results indicate that CPK has a closer relationship with CRK, PEPRK, and PPCK, especially the CRK closest to CPK, which is consistent with previous reports that plant CPK and CRK may share a common evolutionary origin (Hrabak et al., 2003).

To further investigated the footprint of the CPK CRK, PPCK, PEPRK, and SnRK families, we selected four angiosperms (C. papaya, P. trichocarpa, A. trichopoda, and V. vinifera). The reason is that Vitis vinifera, P. trichocarpa, and C. papaya did not undergo α and β duplications and A. trichopoda, a basal angiosperm, did not undergo the γ duplication event (Jiao et al., 2011b; Albert et al., 2013; Lee et al., 2013). Phylogenetic trees in each species were constructed (Figures S5, S6). In each species, the CPK family was divided into four clades (CPKI, CPKII, CPKIII, and CPKIV), and SnRK was divided into three clades (SnRK1, SnRK2, and SnRK3). CPK, CRK, and SnRK were found to exist in Amborella trichopoda, which indicates that these three groups originated from duplication events prior to the γ event. However, PEPRK appeared between the salicoid duplication and the γ event. The PPCK family exists in P. trichocarpa, which indicates that it originated from duplication events prior to the salicoid duplication. Furthermore, due to the salicoid duplication and Brassica-specific WGT events, there were more CDPK-SnRK gene family members in P. trichocarpa and B. rapa than in other species (Tuskan et al., 2006; Wang et al., 2011). In general, during the course of evolution, CPK appeared most recently and expanded most rapidly. Above all, we inferred a possible evolutionary footprint of the CPK family (Figure S7).

The family size and the percentage of CPKs in five plants suggested that CPKs expanded rapidly during evolution and further expanded in the Brassicaceae (Figure 7). WGD is known to have important impacts on the expansion and evolution of gene families in plant genomes. However, along with the gradual increase in the CPK percentage, the genes of group III were completely lost in V. vinifera (Figure 7C). Compared with other groups, the expansion of group III was more unstable. To further investigate the relationship among the four groups, the genetic distance was analyzed. The results indicated that the genetic distance between CPKI and CPKIII was shorter than that between CPKI and CPKIV or that between CPKI and CPKII, and the genetic distance between CPKII and CPKIV was shorter than that between CPKII and CPKIII (Figure 7B). In summary, we inferred a possible evolutionary footprint of the CPK family. Before the γ event, all groups in the family (CPKI, CPKII, CPKII, CPKIV) had already appeared. The gene family further expanded within Brassicaceae. Thus, the CPK family almost doubled in size in the B. rapa genome compared with that of A. trichopoda through three duplications, one triplication, and fractionation. The expansion of groups I, II, and IV played a major role in the expansion of the CPK gene family.

To investigate the divergence expression profiles of CDPK-SnRK genes in different organs between A. thaliana and B. rapa, including roots, stems, leaves, flowers and siliques, the expression patterns of all genes were investigated (Figure S8, Supplementary Tables 5, 6). BrCDPK-SnRK genes were found to be expressed in roots (99 BrCDPK-SnRKs; 77.95%), stems (104; 81.89%), siliques (106; 83.46%), leaves (101; 79.53%), and flowers (119; 93.7%) (Figures 8A–D). A total of 75 (88%) AtCDPK-SnRKs showed high expression (mean-normalized value > 1) in at least one of the five tissues (Figure S8), including roots (34 AtCDPK-SnRKs; 40.00%), stems (36; 42.35%), siliques (22; 25.88%), leaves (19; 22.35%), and flowers (16; 18.23%). Among the 116 CPKs (including 46 AtCPKs and 70 BrCDPKs), 4 (CPK28, CPK29, CPK30, and CPK33) were not expressed and 2 (BrCPK4 and BrCPK25) were only expressed in flower tissue (Figures 8A,B). In addition, a total of 34 BrCPKs, 8 BrCRKs, 4 BrPEPRKs, and 2 BrPPCKs had high expression levels (FPKM value > 10) in at least one tissue; 13 BrCPKs, 2 BrCRKs, and 1 BrPEPRK were highly expressed in all 5 tissues. However, only 2 SnRK genes were expressed in one tissue; the remains 93 SnRK genes were expressed in five tissues.

Subsequently, we selected the expression patterns of genes in five tissues on the phylogenetic tree of all CDPK-SnRK genes to investigate whether the functions of homologous genes were divergent (Figure S8). All CPKI, PEPRK, and PPCK and most CPKIII and CRK genes had high expression levels, suggesting significant roles of these genes in plant development. Most BrSnRK3s exhibited little or no expression. However, the AtSnRK3s were all expressed in five tissues, indicating that BrSnRK3 genes may have lost some functions after the duplication events (Figure S8).

CDPK-SnRKs are suggested to play an important role in the regulation of gene expression in response to abiotic stresses (Figure 9, Supplementary Table 7). To investigate the divergence information of the BrCDPK-SnRK gene family, the expression patterns following different treatments, including ABA, GA, NaCl, heat, cold, and PEG treatments, were analyzed (Figures 9, 10). The qRT-PCR results demonstrate that those BrCDPK-SnRK genes respond differentially to particular abiotic stresses. Sixteen percent of investigated BrCPK genes show increased expression levels upon GA at 6 h, while the other 84% of genes display downregulation or no significant changes. Meanwhile, 20% of genes were upregulated at 1 h and 6 h under GA treatment (Figures 9A,C). Two genes (BrCPK4 and BrCPK10) were induced under both ABA and GA treatment (Figure 9B). During four abiotic stress treatments, excluding PEG, 80% of BrCPK genes were highly responsive to cold, NaCl, and heat. We found that, with the exceptions of BrCPK2, 29, 23 and 26, all BrCPK genes were significantly induced in response to NaCl treatment (Figures 9D–F,H). Under heat treatment, we found that only two genes (BrCPK4 and BrCPK10) had no expression; other BrCPK genes were highly expressed in NaCl treatment. SnRK2 genes recognized as coding for enzymes involved in abiotic stress signal transduction in plants (Kulik et al., 2011). Therefore, we selected six BrSnRK2 genes; all of these genes had the highest expression under heat treatment (Figure 10). Furthermore, BrSnRK2.12 was induced under all stress treatments (Figure 10).

To investigate the connections between these genes, coregulatory networks were established based on the PCCs of stress-inducible BrCDPK-SnRK gene pairs (Supplementary Table 8, Figure 11). All the BrCDPK-SnRKs appeared to have different degrees of positive correlation. Next, 25 BrCPKs with PCC values that were significant at the 0.05 significance level and were greater than 0.8 were collected and visualized to construct hormone and abiotic stress coregulatory networks (Figure 11A). Six BrSnRKs and seven BrCPKs had positive significant correlations (Figure 11C). Most correlations occurred among members belonging to the same group, suggesting that the gene duplication not only led to functional divergence but also enhanced the cooperative interaction of homologs to help plants to adapt to their complex environment. In addition, based on predicted Chinese cabbage BrCDPK-SnRK superfamily protein interactions, their interactions form a complex network in both Arabidopsis thaliana and B. rapa. Moreover, to study the protein interactions of stress response genes between B. rapa and Arabidopsis, STRING 9.1 was utilized (Figure S9). Figure S9 shows four complex interaction networks of CDPK-SnRKs, providing an overall view of the relationship between B. rapa and Arabidopsis.

The B. rapa sequences were downloaded from BRAD (http://brassicadb.org/brad/) (Wang et al., 2011). The Arabidopsis sequences were retrieved from TAIR (http://www.arabidopsis.org/), and the sequences of rice were retrieved from RGAP (http://rice.plantbiology.msu.edu/; Kawahara et al., 2013). The gene information on A. trichopoda genes was retrieved from the Amborella Genome Database (http://www.amborella.org/). The others 14 species' gene information was downloaded from Phytozome v9.1 (http://www.phytozome.net/; Goodstein et al., 2012). The homologous CDPK-SnRK genes in other species were identified through comparison with Arabidopsis. First, BLASTP searches were performed against the rice protein sequences using an E-value threshold of 1 × 10⁻¹⁰. The top-ranked rice hit was used for BLASTP searches of the Arabidopsis proteins to confirm homologies. Starting with both Arabidopsis and rice homologs, BLASTP searches were performed against the proteins of other species (e < 1 × 10⁻¹⁰, identity >70%). These potential sequences were analyzed using the tool SMART (http://smart.embl-heidelberg.de/) and the NCBI database (http://www.ncbi.nlm.nih.gov/).

BLAST and the Multiple Collinearity Scan toolkit (MCScanX) were used for gene synteny analysis according to previous reports (Wang et al., 2012). An all-against-all BLASTP comparison provided the pairwise gene information and the P-value for a primary clustering. Paired segments were extended by identifying the clustered genes using dynamic programming. The potentially duplicated genes were also identified using MCScanX. The positions of B. rapa CDPK-SnRK genes on the blocks were verified by searching for homologous genes between A. thaliana and three subgenomes of B. rapa (LF, MF1, and MF2) in BRAD (http://brassicadb.org/brad/searchSynteny.php; Wang et al., 2011). The syntenic diagram was drawn using Circos software (Krzywinski et al., 2009).

Phylogenetic analyses were conducted using MEGA v5.073; maximum-likelihood (ML) trees were constructed with a bootstrap value of 1,000 replications to assess the reliability of the resulting trees. The genetic distance used in this study was also calculated with MEGA v5.0. To identify conserved motifs in the complete amino acid sequences of CDPK-SnRK proteins, we used MEME software (http://meme.sdsc.edu/meme/; Bailey et al., 2009). We researched the gene structure using Gene Structure Display (GSDS, http://gsds.cbi.pku.edu.cn/). The interaction network of CDPK-SnRK proteins in Chinese cabbage was constructed with STRING software (Search Tool for the Retrieval of Interacting Genes/Proteins, http://string-db.org/; Yamada et al., 2011).

The duplicated CDPK-SnRK genes from B. rapa were aligned using MUSCLE (Edgar, 2004). The protein alignments were translated into coding sequence alignments using an in-house Perl script. Ks (synonymous substitution rate) and Ka (non-synonymous substitution rate) values were calculated based on the coding sequence alignments using the method of Nei and Gojobori as implemented in KaKs_calculator (Zhang et al., 2006). The Ks values of all duplicated BrCDPK-SnRK genes were then plotted as the density and a boxplot using the R program (Ihaka and Gentleman, 1996). The divergence time was calculated with the formula T = Ks/2r, with Ks being the synonymous substitutions per site and r being the rate of divergence for nuclear genes from plants. The value of r was taken to be 1.5 × 10⁻⁸ synonymous substitutions per site per year for dicotyledonous plants (Koch et al., 2000).

For expression profiling of the CDPK-SnRK genes in B. rapa, we utilized the Illumina RNA-seq data that were previously generated and analyzed by Tong et al. (2013). Six tissues of B. rapa accession Chiifu-401-42, including callus, root, stem, leaf, flower, and silique, were analyzed. The transcript abundance was expressed as fragments per kilobase of exon model per million mapped reads (FPKM). The gene expression patterns of each tissue were analyzed using Cluster 3.0, and the expression values were log2 transformed. Finally, heat maps of hierarchical clustering were visualized using Tree View (http://jtreeview.sourceforge.net/). A. thaliana development expression profiling was analyzed using the AtGenExpress Visualization Tool (AVT; http://jsp.weigelworld.org/expviz/expviz.jsp) with mean-normalized values (Schmid et al., 2005). Venn diagrams were drawn using the R program.

Chinese cabbage (Chiifu-401–42) was used for the experiments. The germinated seeds were grown in plastic pots in a 3:1 soil–vermiculite mixture, and the artificial growth conditions were set at 24/16°C, with a photoperiod of 16/8 h for day/night and a relative humidity of 65–70%. Four weeks later, seedlings were subjected to various treatments. For heat and cold treatment, the pots were exposed to 38°C or 4°C; the other growth conditions were same as described earlier. Meanwhile, plants were cultured with the following four treatments: (1) control; (2) 100 μM ABA; (3) 100 μM GA; (4) (w/v) polyethylene glycol (PEG) 6,000, and (5) 250 mM NaCl. All these treatments were performed over a continuous time course (0, 1, 6, 12 h). Each treatment consisted of three biological replicates. All of the samples were frozen in liquid nitrogen and stored at −70°C for RNA preparation.

Total RNA was isolated from treated frozen leaves using Trizol (Invitrogen, San Diego, CA, USA) according to the manufacturer's instructions. Specific primers used for qRT-PCR were designed by Beacon Designer 7 and are shown in Supplementary Table 9. To check the specificity of the primers, we used BLAST against the Brassica genome. The qRT-PCR assays were performed with three biological and technical replicates. The reactions were performed using a StepOnePlus Real-Time PCR System (Applied Biosystems, Carlsbad, CA, USA). qRT-PCR was performed according to a previous report (Song et al., 2013). Relative fold expression changes were calculated using the comparative Ct-value method.

Pearson correlation coefficients (PCCs) of stress-inducible CDPK-SnRK gene pairs were calculated using an in-house Perl script based on log2-transformed quantitative real-time (qRT)-PCR data (Tang et al., 2013). All gene pairs whose PCC was more than 0.8 and was significant at the 0.05 significance level (P-value) were collected for a gene coregulatory network analysis. The coexpression networks were graphically visualized using Cytoscape based on the PCCs of these gene pairs (Shannon et al., 2003).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.