CSNK2A1 and CSNK2B variants were compiled from the OCNDS literature (Iossifov et al., 2014; Okur et al., 2016; Trinh et al., 2017; Yuen et al., 2017; Akahira-Azuma et al., 2018; Chiu et al., 2018; Owen et al., 2018; Duan et al., 2019; Nakashima et al., 2019; Martinez-Monseny et al., 2020; van der Werf et al., 2020; Wang et al., 2020; Wu R et al., 2020; Xu et al., 2020; Dominguez et al., 2021; Wu et al., 2021), and the POBINDS literature (Poirier et al., 2017; Sakaguchi et al., 2017; Fernández-Marmiesse et al., 2019; Li et al., 2019; Nakashima et al., 2019; Ernst et al., 2021; Rahman and Fatema 2021; Selvam et al., 2021; Valentino et al., 2021; Yang et al., 2021). We also compiled CSNK2A1 and CSNK2B variants from the DECIPHER database v11.6 (https://www.deciphergenomics.org; last accessed 1/5/22) (Firth et al., 2009), Simons Searchlight Single Gene Dataset v8.0, AutDB Autism Informatics Portal (http://autism.mindspec.org/autDB/; last accessed 1/5/22) (Basu et al., 2009), and website/url to ClinVar, “(https://www.ncbi.nlm.nih.gov/clinvar/) (last accessed 1/5/22) (Landrum et al., 2018).

To review the literature, we used the queries: OCNDS, POBINDS, CSNK2A1, CSNK2B. These queries identified papers focused on the two disorders, large-scale sequencing efforts, and other articles that did not contain OCNDS and POBINDS variants. In addition, we used the AutDB Autism Informatics Portal to find additional research articles that contain variants in CSNK2A1 and CSNK2B. For the gene variants described in the literature, we recorded each patients’ variant on CSNK2A1 or CSNK2B. For the ClinVar database, we compiled the gene variants described as OCNDS, POBINDS and also variants described as “neurodevelopmental syndrome/inborn genetic disease,” provided that the phenotypic description closely-related to OCNDS or POBINDS. Some of the patients in the DECIPHER database were published in research articles. These patients were identified by comparing clinical manifestations and/or sequencing facilities, to ensure that each patient data was only counted once.

We applied for and obtained the Simons Searchlight Single Gene Dataset v8.0 which includes CSNK2A1 gene variants of patients. The collection of gene variants was reviewed and deemed non-human subjects research by the Institutional Review Board (IRB) of Boston University Medical Campus and Boston Medical Center (IRB # H-41033). To avoid duplicates, we selected variants from the Simons database not already included in ClinVar, as it was highly possible that the variants were already submitted to ClinVar by the original sequencing teams. These variant data were combined and presented in Tables.

CK2α and CK2β gene and protein data were obtained from UniProt (Accession P68400; Accession P67870). Protein domains were determined based on the most consequent domain definition (according to the CATH database) given in (Niefind and Battistutta 2013), the crystal structure of protein kinase CK2α from Zea mays at 2.1 A resolution (Niefind et al., 1998), and the crystal structure of CK2β (Chantalat et al., 1999). For schemes and histograms, we utilized similar color palettes to facilitate reading and understanding of the data.

We utilized R (v4.0.3) and the Ggplot package (v3.3.5) to plot the number of unique patients per residue in the primary protein sequences, including data for missense, nonsense, and frameshift mutations (excluding splice-site variants). To determine the average number of point mutations per residue in each domain we calculated the number of unique patients per domain (End Residue − Start Residue + 1), and divided the total number of unique patients in each domain by the domain length in Microsoft Excel (v16.16.27).

To determine the probability of finding the number of patients with a particular mutation(s) in each residue, we calculated an exact chi-square test statistic comparing the observed mutation probabilities to those expected by chance and present the corresponding p-value. For this, we assumed that each of the 9 variants per codon are equally likely to occur, thus the counts across the 9 mutations will follow a multinomial distribution with the probability of observing each variant = 1/9 (=0.11).

The nonrandom mutation clustering (NMC) algorithm was utilized to identify nonrandom clusters of mutated residues on each protein (Ye et al., 2010). This probability model assumes that mutations along a protein follow a uniform distribution, that mutations are independent of each other in and within samples, and that the number of samples exceeds the number of mutations. NMC utilizes the differences between pair-wise order statistics to derive probabilities. In the original study, NMC was performed on missense mutations to identify activating mutations in cancers that could be targeted for pharmacological intervention. We have expanded the statistical method to identify residues that are found mutated non-randomly among the point mutations that we collected. We utilized all mutations, except splice, as our study differs from cancer studies where missense mutations are key for disease development and progression. All the mutation analyzed in CK2α and CK2β in this study appear in the patient populations that we are studying and therefore, all may have a role in disease (note that it has not been yet demonstrated that any of these mutations causes disease in vivo). For each analysis we obtained clusters ranging from 1 residue to almost full-length protein. Out of these data, significant (p < 0.05) clusters of only one residue were depicted in a table. R code used to generate histograms and NMC analysis is provided in R markdown upon request.

To determine the conservation of the changed residues across eukaryotic species through evolutionary alignment, the primary sequences of CK2α and CK2β from diverse eukaryotic organisms were downloaded. We utilized the species selected in Homologene, and utilized their curated sequences, except for Xenopus laevis and Danio rerio, for which we found sequences with high similarity to the human sequence using BLASTP. Multiple sequence alignment of the sequences utilized was built with MUSCLE (https://www.ebi.ac.uk/) (Edgar 2004) and displayed in HTML format. In this format, conserved residues appeared in blue (conservation across most species), gray (conservation across some species) or white (no conservation).

In silico tools were used to obtain predictions of the impacts of CK2α and CK2β missense mutations based on their evolutionary conservation (PANTHER, MutationTaster2), functional impact (SIFT, PROVEAN, Polyphen-2, I-Mutant 3.0 Disease, MutationAssessor, SNAP2), and importance in protein stability (PremPS, Kinact, I-Mutant ΔΔG). Two consensus programs (PredictSNP, REVEL) were also used for functional predictions. Our classification of these programs into categories of evolutionary conservation, functional impact, and effect on protein stability are not strict classifications, as many of these programs generate their predictions based on many of these factors. Rather, these classifications are used to more easily compare programs that consider similar protein characteristics when evaluating mutations. The following steps (exemplified for CSNK2A1) were followed to obtain predictions from each of these prediction programs:

1) PANTHER is a protein classification system that predicts the functional consequence of a single-nucleotide polymorphism (SNP) based on the preservation time of the mutated amino acid (Tang and Thomas 2016). Click on the cSNP Scoring tab. Input the FASTA sequence for CSNK2A1 and a list of all missense substitutions. Select “Homo sapiens” for the organism and submit. Prediction outputs are “Benign,” “Possibly Damaging,” or “Probably Damaging.”

The McNemar’s test was used to compare the categorical outputs of the functional programs (SIFT, Polyhphen-2, PROVEAN, Mutation Assessor, SNAP2, I-Mutant3.0 Disease, REVEL, and PredictSNP) and the stability programs (PremPS, Kinact mCSM, Kinact SDM, Kinact DUET, and I-Mutant3.0 ΔΔG). The Kinact suite of tools is only available for kinases such as CK2α therefore for CK2β we only compared PremPS and I-Mutant3.0 ΔΔG. We calculated two values to compare each pair of tests that were performed on the mutations. First, the Kappa Coefficient (which calculates the degree of agreement beyond what would be expected by chance.) Second, the p-value from McNemar’s test, an assessment of whether there is a significant difference between tests in the ratings or “benign” vs. “effect”. Third, the p-value from McNemar’s test, an assessment of whether there is a significant difference between tests in the ratings or “benign” vs. “effect.”

For prediction of changes in protein-protein binding affinity on mutations we utilized BeAtMuSiC (http://babylone.ulb.ac.be/beatmusic). BeatMusic predicts changes in binding free energy (ΔΔG) induced by point mutations (Dehouck et al., 2013) based on known protein structures, the strength of interactions at the interface and the overall stability of the complex. As input for this prediction we used the CK2 holoenzyme structure (PDB ID: 4DGL; (Lolli et al., 2012)).

The CK2α^1-335 structure (PDB ID: 2PVR (Niefind et al., 2007)) or CK2β^1-193 structure (PDB ID: 3EED (Raaf et al., 2008)) were used as reference structures. The protein structures were downloaded from the RCSB Protein Data Bank (Berman et al., 2000). To visualize the conserved residues of CK2α and CK2β the ConSurf server (Glaser et al., 2003; Landau et al., 2005; Ashkenazy et al., 2010; Celniker et al., 2013; Ashkenazy et al., 2016) was used. Multiple Sequence Alignments were built using MAFFT (Katoh et al., 2018). The Homologues were collected from UNIREF90 using the homolog search algorithm HMMER (Eddy 2011) (E-value: 0.0001; No. of HMMER Iterations: 1). As maximal identity between sequences 95% and as minimal identity for homologs 40% were chosen. The calculation was performed on a sample of 150 sequences that represent the list of homologues to the query (Supplementary Table S3). Conservation scores were calculated using the Bayesian method (Mayrose et al., 2004). For structural analysis the programs COOT (Emsley et al., 2010) and PyMOL were used. Figures 1, 4A,C, 6, 7 were drawn using PyMOL (Schrödinger 2013).

The first publications for OCNDS and POBINDS associated these syndromes with variants in CSNK2A1 and CSNK2B, respectively, and provided clinical manifestations that defined the syndromes. Since then, other variants in these two genes have been published and classified as OCNDS and POBINDS. In addition, a number of unpublished CSNK2A1 and CSNK2B variants have been placed in data repositories which include clinical manifestations similar to the published variants. In view of the spread of the data, we collected all of the CSNK2A1 and CSNK2B variants from publications and data repositories. Chromosomal variants were not collected in this study. The results of this data collection are found in Tables 1, 2, where the data is divided into non-conservative amino acid mutations and splice mutations. Table 1 includes the 36 published CSNK2A1 variants (some of which were also in data repositories), and 32 unique variants from data repositories (Simons, DECIPHER, ClinVar, AutDB). Table 2 includes 40 published CSNK2B variants (some of which were also in data repositories) and 27 unique variants from data repositories. By pooling these diverse variant data resources, we have increased the number of variants analyzed by 89% in CK2α and by 68% in CK2β. Interestingly, no genetic diseases are known so far to be associated with mutants of CSNK2A2, a CSNK2A1 paralog, located on a different chromosome (Yang-Feng et al., 1991).

For CSNK2A1, we collected 7 splice variants and 61 single nucleotide variants (SNVs) (Table 1). The SNVs were: 1 start site, 47 missense variants (S51R was coded by 2 different variants), 9 nonsense, and 4 frameshift variants. For CSNK2B, we collected 14 splice variants and 49 SNVs: 4 different start site variants, 20 missense, 13 nonsense (Y101* was coded by 2 different variants), 11 frameshift, and 1 codon duplication variant (Table 2). Most of the SNVs in CSNK2A1 were substitutions leading to missense mutations while in CSNK2B there was a high number of duplications and deletions, which resulted in a larger number of nonsense and frameshift mutants compared to CSNK2A1. Figure 2 depicts the location of the CK2α and CK2β mutants along their primary protein structures. Indicated in the Figure are distinct structural domains, highly conserved functional domains and putative functional residues from the literature and from PhosphositePlus (Hornbeck et al., 2012) that we will analyze and discuss below.

We collected the number of patients with each variant to identify hotspots (excluding splice sites). First, we represented the total number of unique patients per mutation in each residue (Figure 3). Then, we utilized a non-random mutation clustering (NMC) approach to determine clusters of mutated residues that cannot be explained by random mutation in the nucleotide sequence. We found a number of significant clusters of residues of varying lengths (Supplementary Table S1), of which we selected single-residue clusters (Table 3). For CK2α, these were (in order of significance): K198, Y50, R47, R80, R191, R312, H160, D156, S51, Y39, R195 and D175. For CK2β, these were D32, R47, M1, C137 and R86. Since mutations in these residues were determined to be non-random in the NMC analysis, these residues may be potential mutational hotspots.

For some residues, all the patient mutations were the same (e.g., K198R in CK2α, and R47* in CK2β) while for others, there were multiple different mutations in the same residue (Supplementary Table S2). For CK2α, 24% of residues had 2 or more distinct mutations (11 of 45 residues), and 18% for CK2β (7 of 39 residues) (Supplementary Table S2). Therefore, we calculated the probability of finding the number of patients with mutation(s) in each residue. If we assume that each of the 9 variants per codon are equally likely, then the counts across the 9 mutations will follow a multinomial distribution with the probability of observing each variant once in nine (1/9 = 0.11). We calculated an exact chi-square test statistic comparing the observed mutation probabilities to those expected by chance and present the corresponding p-value (Supplementary Table S2). For CK2α, significant results were found for R47G/Q, Y50C/F/S, R80C/H, H160R, R191Q/*, R195*, K198R and R312Q/W; and for CK2β for D32A/N, R47*, R86C, R111P, and H165R.

Understanding the genetic mechanisms (population, gene characteristics, etc.) underlying the higher mutability of these sites may lead to important insights that could help us understand the etiology of these diseases.

Studies on neurodevelopmental diseases (NDD) find that missense mutations cluster in or near the functional domains of NDD-associated proteins, in contrast with rare missense mutations found in the 1000 Genomes project (Geisheker et al., 2017). CK2α and CK2β missense mutants were found in both functional domains and interdomain sequences (Figures 2, 3).

We compiled the number of patients in the different functional domains (excluding splice sites) and calculated the ratio of mutations to residues (number of mutations in a protein domain divided by the number of residues in the protein domain) (Table 4; Supplementary Table S1). All functional domains had mutations, albeit at different ratios. None of the Ser and Thr phosphorylated sites in CK2α or CK2β were mutated, but as we will discuss in the structural section, phosphorylation may be affected. For CK2α, the highest ratio was found in the Gly-rich-loop (4.3) followed by the P+1 loop (3.2) while the lowest ratio was found in the activation loop (0.08, with a similar ratio to interdomain sequences). For CK2β, the highest ratios were in the Cys bound to Zinc (1), putative destruction box (0.75) and interface one to CK2α (0.71) while the lowest was in the acidic loop (0.13). Intriguingly, the ratio of the N-terminal segment was among the four highest, suggesting that this region contains functional domains. The interdomain sequences in CK2β showed a higher ratio (0.17) than one known functional domain, suggesting that there may be unidentified functional domains in CK2β located in these sequences. Both these hypotheses are supported by the fact that most of the clusters of one residue in the NMC analysis for CK2β fall in known functional domains, except for D32 and R47 in the N-terminal segment, and R86 in interdomain sequences.

Therefore, as the regulation of CK2α and CK2β is still not fully understood, we examined the NMC data to identify short clusters of mutations outside the known functional domains that could be novel functional domains (Supplementary Table S1). For CK2α, we found clusters including the known functional domains (or shorter versions), except for the hinge region that appears in clusters including the activation segment. Residues in the N-terminal segment were found at minimum in clusters that included at least the Gly-rich loop (e.g.: cluster 21-53), suggesting a functional interaction between the N-terminal and the Gly-rich loop. We also found clusters from residue 107 until the second CK2β binding domain, and from residue 147 to the catalytic loop, the activation segment and residues 312/326. These residues (107, 147 and 312/326) are outside of the known functional domains. This suggests that these residues have functional interactions with known functional domains and/or form part of novel functional domains, particularly the Ct residues that are essential for CK2α’s tertiary structure. This is supported by the fact that most of the clusters of one residue in the NMC analysis fall in known functional domains, except for R312 in CK2α (discussed in our structural analysis).

For CK2β, the significance found in the NMC analysis was lower than in CK2α, most likely due to the limited sample number. We did not find short clusters that included the residues of the known functional domains. A number of clusters started at residues 1, 20, 21 and 27 of different lengths (highest significance ending in 32–47). We found significant short clusters around residue 32 (e.g., 31–35), and clusters including residues 47 (34/42-47) and 86 (86–88). These data suggest that these may be one to three novel functional domains, which will be further discussed in our structural analysis.

Overall, these data suggest that protein function will be affected at least by some of these mutations. We will need experimental data to show that the identified potential domains in CK2α and CK2β are bona fide functional domains. If so, NMC analysis could be used to identify novel potential functional domains for NDD-associated proteins. Our next analyses will assess, via diverse methods, the potential impact of these mutations in protein structure and function, and will indicate potential biochemical mechanisms.

Disease-causing mutations are more likely to occur in evolutionary conserved positions in proteins (Michaelson et al., 2012; Yuen et al., 2016; Geisheker et al., 2017; Lelieveld et al., 2017; Wang et al., 2020); therefore, we assessed whether the residues mutated in CK2α and CK2β in OCNDS and POBINDS were conserved across species including vertebrates, invertebrates, plants, yeast and fungi. Supplementary Figures S1, S2 show an alignment of a subset of species using MUSCLE (highest conserved residues highlighted in blue), with the residues mutated in OCNDS and POBINDS highlighted in red. Supplementary Figures S3, S4 show an alignment of 150 species (highest conserved residues highlighted in raspberry red) using MUSCLE. The accession numbers of the 150 sequences can be found in Supplementary Table S3. Figure 4 shows a 3D representation of alignment of 150 species using a 3D surface representation of the conserved residues calculated using the ConSurf server (Glaser et al., 2003; Landau et al., 2005; Ashkenazy et al., 2010; Celniker et al., 2013; Ashkenazy et al., 2016), and Supplementary Figure S5 shows the 3D surface representation using the subset of sequences in Supplementary Figures S1, S2. The alignments showed that the highest conservation across species was found in the first 324 CK2α residues. The CK2β sequence is less conserved than the CK2α sequence and there was a higher degree of conservation between vertebrate and invertebrate sequences and between plant and unicellular organisms (Supplementary Figure S5). Chantalat et al. (Chantalat et al., 1999) found 40 identical residues in the seven species they investigated, which are highlighted in Supplementary Figure S2. Among these conserved residues, 32, 34, 35 and 86 have missense mutations in POBINDS (defined by Chantalat as surface cluster 1, and discussed in the structural section). As discussed above, based on the NMC analysis, these residues could form part of novel functional domain(s).

Prediction programs based on evolutionary conservation are frequently used to determine the possible consequences associated with mutations. We used two well-known prediction programs, MutationTaster2 and PANTHER, to analyze CK2α and CK2β mutations (Supplementary Table S4). In MutationTaster2, all CK2α mutants were classified as disease causing and in PANTHER, all mutants were classified as probably damaging except for S356T (possible damaging). For CK2β mutants, both MutationTaster2 and PANTHER classified the mutants as disease causing/probably damaging. There was a concordance between the two prediction programs even for residues not fully conserved across species.

We assessed the functional effects of CK2α and CK2β mutations with six functional prediction programs (SIFT, PROVEAN, I-Mutant 3.0 Disease, MutationAssessor, Polyphen-2, and SNAP2) and two consensus prediction programs (REVEL and PredictSNP). The majority of the programs only analyze missense mutations; therefore, the computational analyses below are restricted to these mutations. Figure 5 represents the residues and number of patients with missense mutations.

For CK2α, the mutations predicted to be functionally impacted across all programs were D156H, D156E, D156Y, K158E, N161S, N161D, P231R, R312Q, and R312W (Table 5A). Mutations Y50C, Y50S, S51I, S51R, V53L, C147Y, K158R, G177S, S194F, and F197I may also be highly affected as none of the programs predicted these mutations to be neutral, tolerated, or benign. In addition to these mutations, M1?, G46V, M153R, D175G, D175E, and L178W were predicted to be damaging by both consensus programs. Mutations S356T, P363H and P382L were predicted neutral/benign/non-disease causing in most of the programs. For the rest of the mutations there were discrepancies between the programs. For CK2β, the mutations predicted to be functionally impacted across all programs were D32A, N35K, R86C, R111P, C137R, C137G, C137F, and H165R (Table 5B). The mutations D32N, F34S, F106V, L187P, and L187R may also be highly affected because none of the programs predicted these mutations to be neutral/tolerated/benign. Six prediction programs also provided numerical values associated with the predictions (Supplementary Table S5).

We performed McNemar’s tests to determine whether there were statistically significant differences in the categorial values between the programs. For this, each categorical value was coded as either as “benign” (neutral, benign, tolerated, non-disease causing) or “effect” (deleterious, effect, probably damaging, possibly damaging, not tolerated, disease, low/medium/high impact, disease causing). We then calculated percent agreement, the Kappa coefficient (to address agreement beyond what would be expected by chance) and the p-value in McNemar’s test (Table 6A,B). For CK2α functional predictions, REVEL had the least agreement with most other programs followed by Mutation Assessor Functional Impact, while SIFT and PredictSNP had agreement with the most number of programs. For CK2β there were fewer significant agreements between programs, probably due to limited sample size.

Since disease-causing mutations are more likely to have common structural features, these computational analyses could identify the best mutant candidates for experimental testing of the functional impact. A complexity in using these programs is that they each have their own thresholds for determining what is damaging and the extent to which a mutation is damaging. For example, some programs predict nearly all mutations to be damaging, while other programs categorize mutations based on predicted severity. In addition, not all programs consider other factors, such as position of the mutation, which have a significant impact on the consequence of the mutation. For instance, in the case of M1?, some of the programs predicted this mutation to be neutral/benign/tolerated and did not take into consideration that this is a mutation at the start codon, and therefore cannot not be neutral. Therefore,

Our next aim was to compare the results from these computation analyses with experimental data to determine which prediction program(s) can best help guide us in understanding the potential consequences of CK2α and CK2β mutations. We acknowledge that we have limited experimental data to date, particularly for CK2β, nonetheless, we correlated experimental data, in particular, kinase activity and subcellular localization with the functional predictions. The in vitro kinase activity of 18 CK2α recombinant mutant proteins purified from bacteria ranges from 0% to 50% of wild-type protein (Dominguez et al., 2021). Table 5C displays the CK2α mutants sorted from low (0% for D175G, S51R and F197I) to medium activity (50% for V73E and R191Q) compared with the wild type protein. Given that all of these mutants had less than 50% activity, PROVEAN and SNAP2 could be used to find mutants that have impacted activity since these two programs properly categorized all of these mutants as having an effect (except M1?). For the rest of the programs, we assessed whether their categorical values may have discrimination power between mutants with the less activity (0%–25%) compared to those with higher activity (25%–50%) (bold line in the table). Based on this cut-off, PredictSNP appears to most accurately predict most of the 0%–25% range to be deleterious (eight out of nine mutations) and it predicts less of the 25%–50% to be deleterious (only 2 out of five mutations). This suggests that predictSNP may be able to pinpoint the mutants most affected and may be utilized to predict large activity changes. We also investigated the numerical predictions to determine whether we can set a range of values for mutants more or less affected (Supplementary Table S5). For CK2α, SNAP2 and REVEL showed somewhat overlapping ranges. For SNAP2: mutants with 0%–25% activity were in the range 96–6 and those with 25%–50% activity in the range: 79–31, suggesting that the most affected mutants could be found in the 96–79 range. For REVEL: mutants with 0%–25% activity were in the range 0.902–0.361 and mutants with 25%–50% activity were in the range: 0.552–0.357, suggesting that the most affected mutants could be found in the 0.902–0.552 range.

We also compared these prediction results with experimental data on nuclear/cytoplasmic distribution of CK2α mutant proteins (Dominguez et al., 2021) as the subcellular distribution could be an indicator of the similarity to the wild-type protein structure (e.g., folding). CK2α wild-type is predominantly nuclear in 80% of cells. To assess the predictions, we categorized the mutants into three groups: high, medium and low similarity to wild-type depending on the percentage of cells with predominantly nuclear distribution. High similarity mutants to wild-type had approximately 80% nuclear distribution (E27K, R80H, I174M, D175G, R191Q, K198R, P231R), medium similarity mutants had 50%–70% nuclear distribution (R47G, R47Q, Y50S, S51R, V73E, D156H and F197I), and low similarity mutants had 20%–30% nuclear distribution (R321Q/W). Table 5D displays the CK2α mutants from low to high similarity to wild-type nuclear distribution. There does not seem to be a clear correlation between the categorial predictions and CK2α subcellular distribution for any of the programs. We also investigated the numerical predictions to determine whether we can find a range of values to distinguish mutants more or less affected (Supplementary Table S5). SNAP2 and MutationAssessor showed somehow overlapping ranges, as follows. For SNAP2, these were the range 96–90 (20%–30% nuclear distribution); range: 91–31 (50%–70% nuclear distribution), and range 89–6 (80% nuclear distribution) suggesting that the most affected mutants could be found in the 96–91 range. For MutationAssessor, these were range >4.61 (20%–30% nuclear distribution), range 4.21–0.525 (50%–70% nuclear distribution), and range 3.865–0.095 (80% nuclear distribution), suggesting that the most affected mutants would be ranked >4.61 and the less affected mutants ranked <0.525.

More data is needed to determine the predictive power of these programs as it relates to CK2α and CK2β, including computing of more complex characteristics that can be affected by changes in the primary structure of these proteins, such as interaction and oligomerization.

We utilized PremPS, Kinact, and I-Mutant 3.0 ΔΔG to predict protein stability changes of missense mutations. All three of these programs were utilized for CK2α mutations. Only PremPS and I-Mutant 3.0 ΔΔG were used for CK2β mutations since Kinact is specifically for kinase proteins. The prediction programs provide the Gibbs Free energy (kcal/mol) for change in protein stability. As in the above analyses only missense mutations were analyzed (Table 7).

For CK2α, mutations Y39S, G199D, and R312Q were predicted to be highly destabilizing across all programs (Table 7A; Supplementary Table S6). Mutations Y39C, Y50S, R80C, R80H, K158E, L178W, and P231R may also be destabilizing as several programs predicted them to be highly destabilizing and none of the programs gave them a stabilizing prediction. Predictions were somehow different between the programs. Kinact mCSM and I-Mutant 3.0 ΔΔG predicted all mutants to be destabilizing at some level, and PremPS predicted all mutants to be destabilizing at some level except for E27K and D156Y. Both Kinact SDM and Kinact DUET predicted all mutations to be destabilizing except for a few, including G46V, S51I, T127M and D175G that both programs predicted stabilizing. For CK2β, both PremPS and I-Mutant 3.0 ΔΔG predicted all CK2β mutations to be destabilizing at some level (Table 7B; Supplementary Table S6). Mutations F34S, F34C, F106V, C137G, L187P, and L187R were predicted to be highly destabilizing by both programs while D32A and E77K were predicted to have a lower destabilizing effect. We performed McNemar’s tests for the differences in the categorial values between the programs (Tables 6C,D). For CK2α stability, Kinact SDM ΔΔG and Kinact DUET ΔΔG had the least agreement with other programs (significant p-value).

A complexity of the use of these programs is that these prediction programs required a structure input, except for I-Mutant 3.0 ΔΔG that used protein sequence. Overall, there are 287 CK2α structures in the PDB; 228 of them are from the two human paralogs CK2α and CK2α. Selecting the structure to use may be difficult for the newcomer to the field. Indeed, we obtained varied results when using three different structures for human CK2α to run Kinact (PDB ID’s 2PVR, 5OMY, 2ZJW) (not shown). For the data in the table, we used PDB ID 2PVR as it was the structure with the highest resolution.

There is no experimental study of the stability of CK2α mutants to be compared with the results from these program analyses. However, there is data on expression levels from a subset of CK2α mutant proteins that were expressed in cell lines (Dominguez et al., 2021). Three mutants appear less expressed than others in cell lines: R312Q, R312W and R47Q. As we will discuss below, R312 is a key residue to maintain the 3D structure of CK2α as it forms an ion pair with Q201 in α-EF helix, coupling the a-GHI-subdomain with the α-EF-helix; therefore, it is predicted to be less stable. The mutant with the lowest expression, R312Q, was predicted to be destabilizing by all programs (Table 7C). We also investigated the numerical predictions to determine whether we can find a range of values to distinguish mutants more or less affected (Supplementary Table S6). For CK2α, the three less expressed mutants had I-Mutant 3.0ΔΔG scores from −1.27 to 0.73. For CK2β, T37Yfs*5 had low expression levels (Nakashima et al., 2019). However, we could not assess this mutant, as the prediction programs only analyze missense mutations. More experimental data is needed to determine specific programs that are better at predicting stability properties of the CK2α and CK2β mutants.