The protein sequences of Arabidopsis and rice were obtained from the NCBI database (https://www.ncbi.nlm.nih.gov/) database and a Hidden Markov Model file was created using the method BLAST (Altschul et al., 1997). The protein sequence of StCDPK genes was downloaded from the published potato genome database in “DM v4.04.” All candidate genes containing a kinase domain and EF-hands were identified and further verified in the Pfam (http://pfam.xfam.org/) and SMART database (http://smart.embl-heidelberg.de/) (Letunic et al., 2021). The molecular weight, theoretical pI, grand average hydropathicity, and the instability coefficient of the protein encoded by the StCDPKs gene were predicted using the ProtParam tool of ExPaSy (http://web.expasy.org/protparam/) (Gasteiger et al., 2005). The signal peptide of the StCDPK protein was analyzed using SignalP-5.0 (http://www.cbs.dtu.dk/services/SignalP/) to determine if it was a secretory protein. Interpro (http://www.ebi.ac.uk/interpro/) was used to identify StCDPK domains and predict the number of EF-hands present. The Palmitoylation and myristoylation sites of StCDPK proteins were predicted using GPS-Palm (http://csspalm.biocuckoo.org/index.php) and GPS-Lipid (http://lipid.biocuckoo.org/index.php), respectively (Ren et al., 2008; Xie et al., 2016).

A phylogenetic tree was generated by using the Neighbor-Joining method with 1,000 rapid bootstrap repeats using MEGA 7.0 (Kumar et al., 2016) using protein sequences from potato, Arabidopsis, and rice. The generated phylogenetic tree was visualized using Evolview v2 (http://www.evolgenius.info/evolview/).

The chromosomal distribution of StCDPKs was identified using the MCScanX (http://chibba.pgml.uga.edu/mcscan2/) program, and the genomic positions of the StCDPK genes were mapped using the S. tuberosum database Spud DB (http://solanaceae.plantbiology.msu.edu/pgsc%20download.shtml). The map was drafted using the Mapchart program (https://www.wageningenur.nl/en.htm). The synteny blocks were used for constructing a synteny analysis map that included the genomes of potato (S. tuberosum), A. thaliana, O. sativa, and S. lycopersicum. Figures were generated using the Circos program (version 0.69) (https://circos.ca/). Gene duplication events of CDPK genes in potato were investigated. Three criteria were used in defining gene duplication: 1) the alignment length covered >80% of the longer gene, 2) the aligned region had an identity >80%, and 3) only one duplication event was counted for tightly linked genes. All of the relevant genes identified in the potato genomes were calculated using MSCanX and visualized in Circos. Based on the phylogenetic tree results, a molecular evolutionary analysis of the StCDPK genes was performed by calculating the nonsynonymous (Ka) to synonymous (Ks) substitution ratio of the duplicated gene pairs in S. tuberosum using the KaKs_Calculator in TBtools (https://github.com/CJ-Chen/TBtools) (Chen et al., 2020).

To identify and draw the gene structure of the CDPKs genes, the genomic sequences were aligned with the corresponding coding sequences in the GSDS 2.0 server (http://gsds.cbi.pku.edu.cn). Proteins that share motifs within the StCDPK family were identified using the Multiple Expression motifs for Motif Elicitation (MEME) motif search tool (http://memesuite.org/tools/meme) (Bailey et al., 2009). They were then visualized using TBtools software.

The 1,500 bp genomic sequence upstream of the transcription start site of the CDPKs gene family members was obtained from the potato genome database (http://solanaceae.plantbiology.msu.edu/index.shtml) and the cis-acting elements in the promoter region of StCDPKs gene family members were predicted by the PlantCARE online tool (http://bioinfomatics.psb.ugent.be/webtools/plantcare/html).

A predicted protein-protein interaction network of the StCDPKs and their interacting proteins was constructed using the online program STRING V11.5 (https://string-db.org/) with the following terms; textmining, databases, co-expression, neighborhood, co-occurrence, and experimental evidence. The co-expression network was visualized using Cytoscape (Shannon et al., 2003).

Two potato cultivars QingShu 9 (QS9) and Atlantic (Atl), known to be drought-tolerant and drought-sensitive, respectively, were used in this experiment. They were provided by the Gansu Provincial Key Laboratory of Arid Habitat Crop Science/Gansu Provincial Key Laboratory of Crop Genetic Improvement and Germplasm Innovation. Uniform plantlets from two nodal cuttings were cultured on Murashige and Skoog potato growth medium (Murashige and Skoog, 1962), which contained 5 g L^-1 of agar, 30 g L^-1 of sucrose, and 0.1 g L⁻¹ of inositol. The pH of the medium was adjusted to 5.8 and autoclaved at 121°C at 15 1b psi for 25 min. In vitro propagated plants were maintained in a growth chamber at 25 ± 1°C, with a 16 h photoperiod, with active photosynthetic radiation of 45 µmol photons m⁻²s⁻¹, and with a relative humidity of 55–66% for a 30-day growth period. Healthy and good-looking plantlets were then selected for stress treatment in MS medium supplemented with 150 mM mannitol for simulated drought treatment and control conditions in sterilized glass bottles (120 × 50 mm). The experiment was replicated three times, with each glass bottle containing five cuttings of potato plantlets and subjected to growth conditions 25 ± 1°C, with a 16 h photoperiod, with active photosynthetic radiation of 45 µmol photons m⁻²s⁻¹, and with a relative humidity of 55–66%. Sampling was done on the 15th, 20th, and 25th days, respectively. Collected samples were immediately frozen in liquid nitrogen and stored at–80°C prior to subsequent RNA extraction.

Some physical and physiological parameters were measured after drought induction in the two potato cultivars. Plant height was measured using a meter ruler, i.e., above the root portion of the plant and total root length was measured with a root scanner (EPSON Scan 2 1200XL 2.2) and a root morphology and structure analysis device (Pornaro et al., 2017).

Oxidation stress-related traits such as malondialdehyde (MDA) content, catalase (CAT), and peroxidase (POD) activities, which have been reported as important traits related to drought tolerance in potato (Li et al., 2019), were measured after drought induction. Osmolytes such as proline was also measured after the plant was exposed to drought. Malondialdehyde (MDA content) was determined by the thiobarbituric acid (TBA) method (Heath and Packer, 1968); proline accumulation (Pro) after drought treatment was determined by the sulfosalicylic acid indandione method (Bates et al., 1973). Peroxidase activity (POD) was determined according to the Upadhyaya et al. (1985) method. In brief, POD activity was evaluated in a reaction solution containing 0.2 g of ground fresh plant tissue with 50 mM phosphate buffer, 25 mM guaiacol, and 20 mM H₂O_2. The amount of enzyme activity was expressed as the average change in absorbance (at a wavelength of 470 nm) per minute with readings every 1 min for a total of 3 times. Catalase activity Catalase (CAT) was measured according to the Hamurcu et al. (2013) method. In summary, an amount of 0.8 ml of stock solution was put in a 0.2 g ground sample and centrifuged for 10 min at 12,000 rpm at 4°C. The supernatant was pipetted (40 µl), and mixed with 560 µl of 0.067 M H₂O₂ for CAT determination. The absorbance value was measured at a wavelength of 240 nm within 1 min. All physiological analyses were repeated 3 times.

Total RNA was isolated from in vitro drought-stressed plants using an RNA kit (Tiangen) according to the manufacturer’s instructions. Approximately, 0.1–1 g of the plant material was collected and crushed in liquid nitrogen using a mortar and pestle. The integrity of the extracted RNA was checked on a 1% agarose gel and the concentration and purity were determined using the nucleic acid analyzer. Then, the DNA-free total RNA was used for first-strand cDNA synthesis with the kit (TOYOBO) for gene expression, and the QuantStudio 5 fluorescent real-time quantitative PCR system was used for qRT-PCR detection and analysis. The cDNA was amplified in a total reaction volume of 20 μl; i.e, SYBR Premix Ex Taq TM 10 μl, forward primer 0.8 μl, reverse primer 0.8 μl, cDNA 2 μl, ROX Reverence Dye (2X) 0.4 μl, ddH₂O 6 μl. PCR amplification conditions were: 95°C pre-denaturation for 30 s; 95°C denaturation for 5 s, 60°C annealing for 35 s, 40 cycles; 95°C denaturation for 15 s, 58°C annealing for 60 s, and 95°C for 15 s. After the reaction, the dissolution curve was analyzed and the specificity of product amplification was checked. Actin was used as an internal reference gene and the 2^−ΔΔCt method (Livak and Schmittgen, 2001) was used to calculate the relative expression of StCDPK genes. Each gene amplification result represents data comprising three technical replications for each of three biological replicates. Supplementary Table S1 shows the primers used for gene expression analysis.

Excel 2016 version was used to sort the data, and the IBM SPSS software version 24.0 (International Business Machine Corporation, United States) was used for statistical analysis. Statistically significant differences (p < 0.05) are reported in the text and shown in the figures. GraphPad Prism 7 and Excel 2016 were used for data mapping and the R software was used to draw the heat map.

A total of 25 gene loci were identified as StCDPKs and their coding genes were designated as StCDPK1∼ StCDPK25 according to their location on potato chromosomes (Table 1). In silico analysis of the chromosomal locations of the CDPK loci indicated that the 25 CDPKs were distributed among 11 chromosomes in potato (Figure 1); only chromosome 9 lacked a CDPK locus. There were 5 StCDPK genes distributed on chromosome 10 (ST4.03ch10), which had the highest number of StCDPK loci. Two chromosomes ST4.03ch01 and ST4.03ch11 had 4 loci each; 2 loci were present on chromosomes 3, 4, 6, and 12; there was only 1 locus on chromosomes 2, 5, 7, and 8.

To explore the physicochemical properties of CDPK proteins, the molecular weight (MW), isoelectric point (pI), grand hydropathy score (GRAVY), and instability index were predicted for the 25 StCDPK members. Predicted protein sequence analysis of the 25 StCDPKs revealed that the amino acid lengths ranged from 474 aa (StCDPK9) ∼ 638 aa (StCDPK17), with coding region lengths ranging from 1,425 bp (StCDPK9) ∼ 1917 bp (StCDPK17). The predicted molecular weight of the proteins ranges from 53,356.88 to 70,171.73 kDa. The differences in molecular weight can be largely attributed to the different numbers of domains. The pH at which a given molecule carries no net charge is known as the isoelectric point; this can be informative for identifying proteins with pH-dependent properties. The theoretical isoelectric points of the StCDPKs range from 5.02 to 9.34 (Table 1). With the exceptions of StCDPK5, StCDPK6, and StCDPK9, the predicted isoelectric points of the other StCDPKs proteins were below 7, a trend indicating that most of the StCDPKs are rich in acidic amino acids (Garcia-Moreno, 2009; Talley and Alexov, 2010).

Our analysis of the 25 StCDPK protein domains using the Interpro online tool (Table 1) indicated that most of the CDPK proteins were predicted to contain 4 EF-hand domains with the notable exceptions of StCDPK4/9, which were predicted to contain 1 and 2 EF-hand domains, respectively indicating that the arrangement of StCDPK protein domains is largely conserved. At the same time, the results of protein hydrophilicity analysis and instability coefficient analysis showed that the average hydrophilic GRAVY values of all StCDPKs were negative. The negative GRAVY values of all 25 proteins indicate a soluble or hydrophilic nature of the CDPK proteins (Kyte and Doolittle, 1982). Notably, 12 of the StCDPKs had protein instability coefficients predicted above 40 (Table 1), suggesting that these 12 proteins may be unstable.

We also assessed the predicted signal peptides, myristoylation sites, and palmitoylation sites for the StCDPKs. A signal peptide analysis showed that there were no obvious signal peptides. Prediction of potential myristoylation and palmitoylation sites indicated that 3 of the 25 CDPK proteins have no myristoylation or palmitoylation sites: whereas StCDPK4, StCDPK10, and StCDPK21 were predicted to have only palmitoylation sites, the remaining 21 proteins were predicted to contain both myristoylation and palmitoylation sites (Table 1).

A phylogenetic tree was constructed using 90 full-length protein sequences to compare the CDPKs of potato with those of rice and Arabidopsis (Figure 2). Using the NCBI database search tool, 34 Arabidopsis and 31 rice CDPKs protein sequences were found, and a neighbor-joining tree of CDPK protein sequences from potato, Arabidopsis, and rice was constructed using MEGA 7.0 (Figure 2). Similar to rice and Arabidopsis, the potato CDPKs were clustered into four main groups (I, II, III, and IV). The largest was group I which comprised 11 StCDPKs, 11 OsCPKs, and 10 AtCPKs. The second group (II) comprises 6 StCDPKs, 8 OsCPKs, and 12 AtCPKs; the group with the fewest CDPK members was group IV, with 3 StCDPKs, 3 OsCPKs, and 3 AtCPKs.

Gene structure analysis of the CDPK genes in potato was conducted to assess gene family expansion and divergence. The number of introns varied from 5 to 12 (Figure 3B). Our analysis showed that most of the 25 StCDPK genes contained 7-9 introns; only the group I gene StCDPK9 had 5 introns, while the group IV genes StCDPK5 and StCDPK6 each contained 12 introns. There were large differences in gene length, exon/intron structure distribution, and fragment length (Figure 3B). Using the MEME website to predict and analyze the conserved motifs of potato CDPKs, a total of 10 conserved motifs were identified (Figure 3A). Of the 10 identified motifs, motifs 1, 2, and 3 occurred in all 25 proteins, followed by motifs 5, 6, and 7, which occurred in all proteins except StCDPK9 and StCDPK4; motif 10 was present in only 5 StCDPK proteins. Groups I and II contain similar motifs, except StCDPK9, which lacks motif 5; all members of group III contain motif 10 but lack motif 9. The proteins of group IV have the lowest number of motifs (Figure 3A). More than 80% of the StCDPKs contain motifs 1–8. Moreover, one protein StCDPK4 contained only 4 of the 10 motifs.

To examine the effect of duplications on the CDPK gene family, we assessed putative tandem duplication and segmental duplication gene pairs using PGDD (Plant Genome Duplication Database), and visualized the results using Circos. We identified 4 pairs of paralogous segmental duplicated genes distributed on different chromosomes (Supplementary Figure S2). The identified CDPK gene pairs were StCDPK1/StCDPK16 found on chromosomes 1 and 10, StCDPK5/StCDPK6 located on chromosomes 2 and 3, StCDPK20/StCDPK8 located on chromosomes 11 and 4, and StCDPK13/StCDPK24 located on chromosomes 7 and 12. No tandem duplication was found among the potato CDPK genes.

To explore evolutionary relationships among the CDPK genes in potato, a synteny analysis that included Arabidopsis, rice, and tomato was performed using MCScanX. The synteny analyses help in gaining knowledge of the evolutionary and functional relationship between orthologs. This analysis found a synteny of 8 potato chromosomes with Arabidopsis, 11 potato chromosomes with tomato, and 1 chromosome with rice. Potato showed maximum synteny with tomato and Arabidopsis and least with rice because of the divergence along with the evolutionary time frame (Supplementary Figure S2). The lowest synteny found in potato and rice is due to the distant evolutionary relationship between dicots and monocots plants.

To estimate the divergence time of potato CDPKs, synonymous (Ks) and nonsynonymous (Ka) substitutions between paralogous gene pairs were calculated using the KaKs_Calculator in TBtool. Ka/Ks ratio <1, Ka/Ks = 1, and Ka/Ks > 1.00 respectively indicate purifying, neutral, and accelerated evolution with positive formation (Li, et al., 2009). In total 4 pairs of paralogous genes were found in potato CDPKs. The calculated Ka/Ks ratio of the gene pairs varied from 0.11—0.82 (Supplementary Table S2). All the paralogous StCDPK gene pairs had a Ka/Ks ratio <1, indicating a purifying formation (Yang et al., 2006). The estimated divergence time ranges from 20.32 to 39.32 million years ago (MYA).

Promoter cis-elements influence the initiation of gene transcription. We performed a bioinformatics analysis to identify possible cis-elements in the promoter sequences of StCDPKs. PlantCARE was used to identify putative cis-acting elements in the 1,500 bp upstream sequence of each StCDPK gene promoter. The results indicated that key components of the StCDPKs gene family include core promoter elements (Figure 4A) and plant-inducible promoter elements. In addition to nuclear promoter elements, we detected putative functional elements including light response elements, hormone response elements, stress response elements, and some growth and development-related regulatory elements.

Among them, the types and numbers of light-responsive elements are the largest, and these are present in the promoter regions of each CDPK gene, including box 4, G-box, GT1 motifs, TCT motifs, and GATA motifs (Figure 4B). Some StCDPKs have unique light-responsive elements, such as AAAC motifs, AT1 motifs, chs-CMA2a, GA motifs, Sp1, Box II, and Gap-box elements, which were found only in StCDPK1, StCDPK8, StCDPK9, StCDPK11, StCDPK17, StCDPK14, and StCDPK21. Interestingly, our results revealed five hormone-dependent elements in the promoters of the StCDPKs (Figure 5), including an abscisic acid-responsive element (ABRE), a methyl jasmonate responsive element (CGTCA-motif, TGACG-motif), a salicylic acid-responsive element (TCA-element), a gibberellin responsive element (TATC-box, GARE-motif, P-box), and an auxin-responsive element (TGA-element, TGA-box, AuxRR-core, AuxRE). Among these five types of hormone-responsive elements, 64% of the StCDPKs contain abscisic acid response elements, making it the most frequent hormone-responsive element; the number of methyl jasmonate response elements is second at 40%

In addition, 68% of the StCDPKs were predicted to contain anaerobically inducible elements (AREs): 52% of the members contained the WUN motif; 40% of the members contained defense and stress response elements (TC-rich repeats); 32% and 28% of the family members contained low-temperature response elements (LTR) and drought response elements (MBS), respectively. The StCDPKs with these elements appears likely to exert functions when potatoes are exposed to low temperature or drought stress (Figure 5).

Network interaction analysis can be an effective method for studying gene function (Zhao et al., 2018). We used the STRING 11.5 (https://string-db.org/cgi/) database to predict an interaction network for the StCDPK proteins (Szklarczyk et al., 2019), and Cytoscape was used to visualize the network. The network can help in connecting proteins of interest to other pathways. The network indicated that potato CDPKs apparently interact with several respiratory burst oxidase homolog proteins (RBOH), which are known to promote ROS scavenging (Foreman et al., 2003) (Figure 6). In particular, the network indicated interactions between potato CDPKs and RBOHA, RBOHB, and RBOHC, proteins with known functions in defense responses (Zhang et al., 2018) and plant development (Torres and Dangl, 2005). StCDPK5/6 was predicted to interact with StWIPK, a mitogen-activated protein kinase that may be involved in the catalytic activity as well as other cellular processes (Kamiyoshihara et al., 2010). Additionally, CDPKs interact with activated disease resistance 1 (ADR1), a protein that mediates resistance against pathogens in a salicylic acid-dependent manner (Saile et al., 2021). StCDPK12/19 was predicted to interact with glycogen phosphorylase (GlgP), an allosteric enzyme involved in carbohydrate metabolism.

Plants are regularly threatened by abiotic and biotic stresses, and CDPKs have been described as essential factors in regulating plant tolerance to biotic and abiotic stresses (Huang et al., 2018). To gain insight into the potential functions of StCDPKs in abiotic stress responses, two cultivars (QS9 and Atl) with distinct drought resistance capacities (Li et al., 2019; Bi et al., 2021) were propagated in vitro with mannitol-induced drought treatment and their expression levels measured on 15, 20 and 25 days. The expression levels of the 25 StCDPKs were assessed by qPCR, among which 20 StCDPK transcripts were detected (Figures 7A-E). Moreover, 5 StCDPK genes (StCDPK4, StCDPK10, StCDPK11, StCDPK17, and StCDPK19) were not detected in either of the examined potato cultivars grown with mannitol-induced drought treatment. The expression of 20 StCDPKs were altered in response to stress, and the expression levels of seven of these genes increased significantly in the drought-resistant cultivar QS9 (StCDPK7/9/20/21/22/24/25) as compared to the drought-sensitive cultivar Atl (p < 0.05). Expression of the genes StCDPK2/5/8/18/23 was significantly higher in QS9 than Atl at 15 and 20 days (p < 0.05) (Figures 7A-E).

We assessed changes in physiological indicators in the Atl and QS9 cultivars after 15, 20, and 25 days of drought. The plant height and total root lengths of the Atl plants were each significantly lower than those of QS9 under all the examined drought durations (p < 0.05) (Figures 8A,B). We also examined the levels of peroxidase (POD) and catalase (CAT) in Atl and QS9 under drought and normal conditions. In Atl and QS9, a high accumulation of these two antioxidants occurred after drought treatment. The POD activity of Atl and QS9 increased by 5.46%, 136.73%, 101.93%, and 171.41% 53.10%, and 75.00%, respectively. The POD activity in QS9 increased significantly (Figures 8E). On the 20th day of treatment, the CAT activities of Atl and QS9 were significantly higher than those of control (p < 0.05) (Figures 8F). The effect of cultivar x mannitol interaction on MDA activities was also significant (p < 0.05) in the two potato cultivars (Figure 8C). Compared to their corresponding unstressed control samples, the percentage increases in MDA content at 15, 20, and 25 days were respectively 98.21%, 56.95%, and 225.61% for the Atl plants and 85.89%, 115.09%, and 60.45% for the QS9 plants; note that the increases in MDA content in the Atl plants were significantly higher than for QS9 for all of the sampling days (p < 0.05). Proline is a proteogenic amino acid that is often used as a stress marker in plants where unstressed plants were reported to contained low content of proline per gram of plant tissue, whereas stressed plants contained high amount of proline (Hossain et al., 2016). We measured the proline content in the two potato cultivars after drought treatment and compared the stressed plants to their corresponding unstressed controls, the Alt plants proline increases by 118.14%, 376.64%, and 32.5% at the 15, 20, and 25-day drought stress time points, respectively, while the QS9 plants had increases of 258.91%, 302.87%, and 364.29% (Figures 8D).

We also conducted a correlation analysis for the expression of StCDPK genes in Atl and QS9 with the detected changes in physiological indicators under drought stress (Figure 9). We detected negative correlations between the expression levels of the StCDPK9/18/20/21/22/23 genes in Atl with the plant height values, whereas positive correlations for plant height were detected for QS9 plants. The StCDPK9/21/22 levels in Atl plants were positively correlated with MDA levels, whereas negative correlations were detected for these genes in the QS9 plants. In drought-stress-treated Atl plants, the expression levels of the gene pairs StCDPK13/23 and StCDPK6/20 were, respectively, positively and negatively correlated with the detected increases in MDA content (p < 0.05). Also the detected increases in MDA content (p < 0.05) were negatively correlated with the expression levels of StCDPK2/21/22 in QS9 (Figure 9). Proline accumulation in QS9 correlated positively with the expression of the gene StCDPK3/7/13/14. There was a negative correlation between the StCDPK1/3 expression in Atl and the detected increase in proline content (p < 0.05), but a positive correlation was detected between StCDPK2 expression and the change in root length (p < 0.05) in Atl. POD activity in QS9 and the expression of StCDPK24 showed a strong positive correlation (p < 0.01).