Genome-wide identification of calcium-dependent protein kinase (CDPK) family members in Phaseolus vulgaris L. and expression analysis during abiotic stresses

Identification of PvCDPK gene family and analysis of basic parameters

Amino acid (aa) sequences of the P. vulgaris L. CDPK gene family were retrieved from the Phytozome v12.1 database (https://phytozome-next.jgi.doe.gov/) under accession number PF03492 (http://pfam.xfam.org). The genomes of A. thaliana and Glycine max were examined in the same database to identify potential CDPK proteins (Lamesch et al., 2012; Valliyodan et al., 2019). The default configurations of the Hidden Markov Model (HMM) validated the CDPK protein sequences. Table S1 enumerates the CDPK protein sequences of various plants. The HMMER database (http://www.ebi.ac.uk) was utilized to examine the CDPK domains within the sequences. The aa count, molecular weight (kDa), and other properties of the CDPK proteins were assessed utilizing the “ProtParam tool” (https://web.expasy.org/protparam/). The phylogenetic studies employed the Neighbor-Joining (NJ) technique with a bootstrap value of 1,000 replicates. The ClustalW algorithm was employed to align the PvCDPK protein sequences (Thompson et al., 1997). Evolutionary diagrams were produced via MEGA v7 (Tamura et al., 2011). The iTOL database was utilized to construct the phylogenetic tree (Letunic & Bork, 2011).

Structure, chromosomal localization, duplication and comparative analysis of PvCDPK gene family members

The coding and non-coding sections of the PvCDPK genes were retrieved utilizing the Gene Structure Display v2.0 web tool using the genomic and coding sequence (CDS) sequences (http://gsds.gao-lab.org/) (Hu et al., 2015). The positions of PvCDPK genes on the chromosome were derived from the Phytozome v12.1 database (https://phytozome-next.jgi.doe.gov/). PvCDPK genes were delineated on each chromosome of P. vulgaris L. utilizing MapChart (Voorrips, 2002). MCScanX (The Multiple Collinearity Scan Toolkit) (Wang et al., 2012) was used with default settings to determine the orthologous relationship between CDPK genes of P. vulgaris L. and G. max.

The substitution ratios (Ka, Ks, and Ka/Ks) between duplicate pairs of PvCDPK genes were estimated using PAL2NAL (http://www.bork.embl.de/pal2nal/#Ref) (Suyama, Torrents & Bork, 2006) and the AML interface tool (https://github.com/abacus-gene/paml) (Yang, 2007). To estimate the duplication and divergence time (MYA) of CDPK genes, synteny maps were generated using TBtools (Chen et al., 2020), and the formula T = Ks/2λ (λ = 6.56E−9) was applied (Yang & Nielsen, 2000; Lynch & Conery, 2003).

To uncover more conserved PvCDPK protein motifs, the Multiple Em for Motif Elicitation (MEME Tool) (https://nbcr.net/meme/meme/meme-download.html) was used (Bailey et al., 2006). The parameters 2, 50, and 10 were used for minimum and maximum width and maximum number of motifs, respectively. The motif regions were set in the range of 2–300. Area distribution repetitions might be any number. Motifs were analyzed with the InterPro database as outlined by Quevillon et al. (2005). The WEBLOGO online web tool (http://weblogo.berkeley.edu/logo.cgi) produced CDPK domain sequence logos for conserved area sequence analysis (Crooks et al., 2004).

Subcellular localization and analysis of cis-acting elements of PvCDPK gene family

The upstream sections (Table S1) containing two kilobase pair (kbp) DNA segments of each PvCDPK gene family member were analyzed using the PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) database for cis-acting element analysis (Lescot, 2002). A figure showing cis-acting elements was generated using TBTools (Chen et al., 2020). WoLF PSORT (https://wolfpsort.hgc.jp) predictor was used to predict the subcellular localization of PvCDPK proteins (Horton et al., 2007).

Common bean homology modeling for CDPK proteins

The Phyre2 database (http://www.sbg.bio.ic.ac.uk/ phyre2/html/page.cgi?id=index) was used to acquire the 3D structures, and protein homology modeling was obtained using previously identified CDPK protein sequences (Kelley et al., 2015). The best 3D image was obtained by comparing the protein models’ reliability rates.

An analysis of the ontology of genes and the links between CDPK proteins in P. vulgaris L.

Protein–protein interactions were examined to ascertain their functional and physical relationships using the STRING database (https://string-db.org). The obtained information was categorized and integrated with the confidence level for every interaction between proteins. Cytoscape software was used to visualize the interactions of CDPK proteins, both among themselves and with other proteins (Shannon et al., 2003). The implementation of functional genomics techniques is a crucial prerequisite for the functional annotation of novel sequence data in plant biotechnology research. Ontology data for PvCDPK genes were acquired using the Blast2GO program, and this information was utilized to assess the functional characteristics of PvCDPK proteins (Conesa et al., 2005).

In silico gene expression analysis

RNA-seq data of P. vulgaris L. under salt and drought stress was taken from the National Center for Biotechnology Information (NCBI)’s SRA collection. The used accession numbers were SRR957668 (leaf subjected to salt stress), SRR958469 (leaf salt control) (Hiz et al., 2014), SRR8284481 (leaf subjected to drought stress), and SRR8284480 (leaf drought control). Gene expression data were normalized utilizing reading per kilobase of transcript per million mapped reads (RPKM) (Mortazavi et al., 2008). The Orange software (Demsar et al., 2013) was employed to transform the RPKM data to log2 and generate a heatmap.

Experimental plant materials and treatments

The P. vulgaris L. cultivars “Elkoca-05” and “Serra” employed in this study were obtained from the Molecular Biology and Genetics Department, Erzurum Technical University. The genotype-specific seeds underwent surface sterilization for 5–7 m with a 1% (v/v) NaOCl solution. Subsequently, perlite was utilized for the germination process. The seedlings were relocated to a hydroponic medium comprising 0.2 L of modified 1/10 Hoagland solution upon attaining the developmental stage specified by Büyük et al. (2019). P. vulgaris L. seedlings were grown at 25 °C and 70% relative humidity in a controlled cultivation room with light and a photosynthetic photon flux of 250 mmol m⁻² s⁻¹. After common bean seedlings reached the first trifoliate stage in the growth chamber, the control group was treated with 0 mM NaCl, and the stress treatment group was subjected to salt stress for nine days using Hoagland’s solution and 150 mM NaCl (for medium salinity stress). Concurrently, drought-stressed common bean plants grown in the same conditions were kept for 24 h in a Hoagland solution that was treated with either 0 (control) or 20% PEG6000 (Aygören et al., 2023). Root and leaf tissues from two distinct common bean cultivars were collected following the ninth day of stress treatment. Following the specified duration, the leaf tissue of the common bean genotypes was stored in liquid nitrogen and preserved at −80 °C until the analysis was performed. Three biological replicates of the common bean genotypes utilized in the study were cultivated, and these replicates were employed for qRT-PCR analysis. The root and leaf tissues were subjected to distinct qPCR analyses.

In vitro qRT-PCR analysis

Trizol Reagent (Invitrogen Life Technologies, Waltham, MA, USA) was used to extract total RNAs. The Multiskan Go spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) was employed to quantify RNA, while a 1.5% agarose gel was utilized to evaluate the quality of the sample. To perform complementary DNA synthesis, the SensiFAST cDNA Synthesis Kit (Cat No: Bio-65053, UK) was utilized, following the instructions provided by the manufacturer. The qRT-PCR study was focused on five PvCDPK genes that were selected from the RNA-seq data. The qRT-PCR reactions were conducted using the RotorGene Q Real-Time PCR System (Corbett Research, Qiagen GmbH, Hilden, Germany) and ABT SYBR Green Mix (Cat. No.: Q03-02-01, Ankara, Turkey). A total of 20 µL of qRT-PCR mix was used, including 10 µL of ABT SYBR Green Mix (2x), 0.4 µL of each primer (one µM forward and reverse), and 200 ng of cDNA. The reaction was carried out as follows: 10 min at 95 °C to be 1 cycle; 15 s at 94 °C, 30 s at 60 °C, and 30 s at 72 °C to be 40 cycles.

Statistical analyses

The housekeeping gene used was the β-actin gene from P. vulgaris L. The 2^−ΔΔCT technique for relative quantification was used to standardize the qRT-PCR data (Livak & Schmittgen, 2001). Information on the primers used in this study is presented in Table S2. A two-way analysis of variance (ANOVA) with Dunnett’s test at the 0.05 significant level was utilized to conduct statistical studies in GraphPad Prism 7.

Identification and physicochemical characteristics of CDPK gene family in P. vulgaris L.

Here, 25 CDPK s were identified in the P. vulgaris L. genome through bioinformatics tools and renamed PvCDPK1.1 to PvCDPK32 according to their locations on chromosomes. The amino acid number, protein molecular weight, and theoretical pI (isoelectric point) of the proteins were identified (Table 1). The number of amino acid sequence lengths of PvCDPK s ranged from 298 (PvCDPK4.2) to 582 (PvCDPK2). The molecular weights ranged from 33.43 kDa (PvCDPK4.2) to 65.13 kDa (PvCDPK2), and the pI values ranged from 4.82 (PvCDPK4.2) to 9.21 (PvCDPK16).

Chromosomal location and duplication events

The 25 PvCDPK genes were distributed unevenly on eight chromosomes: Chr1 (four genes), Chr2 (four genes), Chr3 (three genes), Chr6 (one gene), Chr7 (six genes), Chr8 (four genes), Chr9 (two genes), and Chr11 (one gene) but not on the other chromosomes (4, 5, 6, and 10) of the common bean (Table 2).

The investigation of gene duplication revealed that the PvCDPK1.1/PvCDPK2 and PvCDPK3.1/PvCDPK3.2 genes had a synonymous substitution (Ks) value of 0.64, the PvCDPK8/PvCDPK32 genes had a Ks value of 0.81, and the PvCDPK10.1/PvCDPK10.2 genes with 0.86 Ks, PvCDPK11.1/PvCDPK11.2 genes with 0.74 Ks, and PvCDPK17.1/ PvCDPK17.2 genes with 0.69 Ks. PvCDPK4.1/PvCDPK4.2 and PvCDPK29.1/PvCDPK29. 2 genes were found to be tandemly duplicated with Ks values of 0.42 and 1.11, respectively (Table 3). These genes had Ka/Ks ratios ranging from 0.08 to 0.28. Natural selection during duplication events is represented by values equal to 1, purifying selection is indicated by values less than 1, and positive selection in the evolutionary process is shown by Ka/Ks values greater than 1 (Juretic et al., 2005). Accordingly, it can be inferred that all PvCDPK genes have undergone purifying selection (Fig. 1). In addition, the differentiation time of six pairs of segmentally duplicated genes ranged from 48.94 million years ago (MYA) to 65.57, while tandem pairs ranged from 32.09 to 84.95 MYA.

Phylogenetic relationships and synteny analysis of PvCDPKs in P. vulgaris L. and different plants

Using the protein sequences from the common bean, A. thaliana, and G. max, a phylogenetic tree was created to ascertain the phylogenetic relationships for the common bean’s CDPK gene family. In the phylogenetic analysis, researchers used a total of 98 protein sequences: 25 from P. vulgaris L., 34 from Arabidopsis, and 39 from G. max. As illustrated in Fig. 1, the genes were clustered into three major subfamilies: A, B, and C. Group A is the largest with 69 genes, while Group C is the smallest with nine genes. PvCDPK4.2 was noticed relatively independently of other PvCDPKs. Within the three groups formed in this phylogenetic tree, which played a crucial role in elucidating the molecular evolutionary process, PvCDPK genes were observed to be homologously distributed, especially with AtCDPK genes (Fig. 2).

Synteny analysis was performed to examine shared structural changes in the genome, including chromosomal fission and fusion. The analyses revealed 57 syntenic relationships between P. vulgaris L. and G. max and 23 syntenic relationships between P. vulgaris L. and A. thaliana. While a syntenic relationship was found between all PvCDPK genes and G. max genes, no syntenic relationship was found between AtCDPK genes and PvCDPK genes only in PvChr-6. This indicates a strong evolutionary similarity between G. max and P. vulgaris L. It can also be said that there is a strong syntenic relationship between A. thaliana and P. vulgaris L. in terms of chromosomal significance. In addition, CDPK genes are equally distributed in these genomes, indicating that these gene pairs are widely distributed within the genomes (Fig. 3).

Gene structure and motif analysis of PvCDPK gene family

The analysis of gene structure revealed that the 25 PvCDPK gene family members together harbor 169 introns and 195 exons. PvCDPK16 and PvCDPK28 have 11 introns; PvCDPK3.1, PvCDPK3.2, PvCDPK8, PvCDPK9, PvCDPK17.1, PvCDPK17.2, PvCDPK21, PvCDPK24, PvCDPK29.1, PvCDPK29.2, and PvCDPK32 have seven introns; PvCDPK1.1, PvCDPK1.2, PvCDPK2, PvCDPK4.1, PvCDPK6, PvCDPK10.1, PvCDPK10.2, PvCDPK11.1, PvCDPK11.2, PvCDPK13, and PvCDPK20 have six introns; and PvCDPK4.2 has five introns (Fig. 4).

A total of 10 conserved motifs were identified with lengths ranging in length from eight to 50 amino acids using MEME (Fig. 5; Table S1). All proteins were found to contain motifs 2, 5, and 7, while additional motifs were detected solely in particular subgroups (Fig. 4). As an example, motif 8 is only found in one subgroup (Fig. 4). Except for PvCDPK16 and PvCDPK28, all PvCDPKs have four motif 2.

Analysis of PvCDPKs promoter cis-elements

The regulatory mechanisms of PvCDPK genes were explored by analyzing the two kb upstream sequences from their start codons for cis-regulatory element composition.

As a direct result, a total of 341 cis-acting regulatory elements were identified in the promoters of PvCDPK genes (Fig. 6; Table S1). The promoters of PvCDPK genes also contained a total of 15 cis- acting regulatory elements. The cis-acting elements were divided into four main categories: abiotic/biotic stress-responsive, including 12 elements (ABRE, MYB, MBS, LTR, etc.); development, including elements (CCAAT-box), core elements, and binding sites (W box); and hormonal-responsive (as-1), including only one element.

The greatest number of cis-acting elements were determined in MYC and MYB, and these elements, together with the MBS element, are associated with drought stress defense. While MYC and MYB were found in all genes, MBS was found to be associated with PvCDPK-1.2, -3.1, -4.1, -4.2, -9, -10, -11.1, -13, -17.1, -20, -29.1, and -29.2 genes. Biotic and abiotic stress is represented by the highest number of elements, indicating PvCDPK genes have important roles in response to biotic and abiotic stress (Fig. 6; Table S1). Considering this viewpoint, it is likely that these PvCDPKs participate in multiple biological functions.

Homology modeling, 3D structural prediction, protein interaction network construction, and gene ontology annotation of PvCDPK proteins

Proteins interact to fulfill their functions; therefore, comprehending the connections and mechanisms of complex biological processes is essential. PvCDPK protein sequences were utilized to ascertain their protein-protein interactions via the STRING interface. Here, it was found that 25 proteins interacted with five different common bean proteins. These proteins were V7BRL4_PHAVU-Phvul.006G142500 and V7BS13_PHAVU-Phvul.006G157600 (calcium binding protein 39), V7BH71_PHAVU-Phvul.007G22340 (heat stress transcription factor A-3), V7CME5_PHAVU-Phvul.002G160700 (respiratory burst oxidase homolog protein F-related), and V7CPP5_PHAVU-Phvul.002G293700 (PTHR11972//PTHR11972:SF81 - NADPH oxidase) (Fig. 7; Table S1). All PvCDPK proteins interacted with V7CME5_PHAVU-Phvul.002G160700. PvCDPK17.1 and PvCDPK17.2 interacted with PvCDPK24. The remaining proteins showed no interaction with each other.

Gene ontology facilitates the comprehension of genes through the analysis of biological processes, molecular functions, and cellular components. In the biological process, PvCDPK genes are enriched in the peptidyl-serine phosphorylation, protein autophosphorylation, and intracellular signal transduction. The cellular component category included the nucleus, cytoplasm, and intracellular anatomical structure. Protein serine kinase activity, calcium ion binding activity, calmodulin-dependent protein kinase activity, calmodulin binding activity, and calcium-dependent protein serine/threonine kinase activity were categorized in molecular function (Fig. 8; Table S1). CDPK proteins were identified using the Phyre2 database, and homology modeling was visualized via the 3D modeling technique. Figure 9 illustrates the three-dimensional homology models of the proteins identified in the study.

In silico expression profiles of PvCDPK gene family drought and salt stress

Throughout their period of development and growth, plants are greatly impacted by a wide variety of environmental conditions, including low temperatures, high salt, and drought. Expression pattern analysis can help understand the biological functions of PvCDPK in tissue-specific or abiotic stresses such as salt and drought. To comprehensively analyze the mRNA expressions of PvCDPK genes, RNA-seq data from four normal and treatment samples from the NCBI SRA database were obtained, and FPKM values of 25 PvCDPK genes were evaluated. Five different tissues were taken for evaluation in this study. All the PvCDPK genes were expressed in at least one tissue. Different PvCDPK genes revealed differential expression patterns. For example, only PvCDPK17.1 and PvCDPK24 displayed expression in flowers but no other tissues. PvCDPK11.2 was the most expressed PvCDPK in leaves and stem, whereas PvCDPK6 was significantly expressed in flowers, nodules, root, and stem (Fig. 10A; Table S1). With these findings, it can be said that CDPK genes actively contribute to common bean organ development.

In this study, in silico gene expression analysis under drought and salt stresses was determined. PvCDPK16 and PvCDPK6 expressions were higher in control plants compared to drought- and salt-treated plants. PvCDPK11.2, PvCDPK10.2, PvCDPK32, PvCDPK21, PvCDPK13, and PvCDPK10.1 genes expressed higher than control under drought stress. However, these gene expressions showed lower expression under salt stress (Fig. 10B; Table S1). PvCDPK28 and PvCDPK3.2 were induced after salt treatments, but their expressions reduced after drought treatment. No important change was determined in the expression patterns of other genes. On the other hand, among these genes, only PvCDPK10.2 expression displayed the same expression level between control and treated plants. These findings suggest that PvCDPKs may be involved in the response to a range of abiotic stresses, with different genes displaying distinct responses to stress.

qRT-PCR analyses

In this study, common bean seedlings were treated with drought and salt to examine the function of PvCDPK gene members. Here, two cultivars, Elkoca-05 and Serra, were used. qRT-PCR analyses were performed for five specific primers (PvCDPK1, PvCDPK4, PvCDPK10, PvCDPK20, and PvCDPK29) designed using RNA-seq data, and the results are shown in Fig. 11. Firstly, it was determined that no non-specific results were obtained in the negative control analyses performed in qPCR. Under drought stress, the PvCDPK1, PvCDPK4, PvCDPK10, and PvCDPK29 gene expressions were increased in Elkoca-05. However, there was no significant change in the expressions of PvCDPK analyzed, nor of the PvCDPK genes in Serra, neither in leaf nor in root.

Under salt stress, different expression patterns were observed compared to drought stress. PvCDPK1 and PvCDPK29 were induced in leaf, while PvCDPK4 was induced in root in the Serra cultivar. Besides, PvCDPK4 expression also increased in Elkoca-05 root. No significant change was observed for PvCDPK4 expression in leaf under salt stress in both Serra and Elkoca-05.

Consequently, although gene expression levels of the PvCDPK gene family varied among cultivars, the same genes analyzed under drought stress exhibited no variation in the Elkoca-05 cultivar across different tissues. Although the expression levels of PvCDPK1 and PvCDPK29 in the Elkoca-05 cultivar increased in drought stress treatment, their expression decreased in salt treatment in all tissues. PvCDPK4 was induced under both drought and salt stress in both two cultivars. These findings are in agreement with in silico analyses.