Comprehensive Analysis of Phaseolus vulgaris SnRK Gene Family and Their Expression during Rhizobial and Mycorrhizal Symbiosis

Sucrose non-fermentation-related protein kinase 1 (SnRK1) a Ser/Thr protein kinase, is known to play a crucial role in plants during biotic and abiotic stress responses by activating protein phosphorylation pathways. SnRK1 and some members of the plant-specific SnRK2 and SnRK3 sub-families have been studied in different plant species. However, a comprehensive study of the SnRK gene family in Phaseolus vulgaris is not available. Symbiotic associations of P. vulgaris with Rhizobium and/or mycorrhizae are crucial for the growth and productivity of the crop. In the present study, we identified PvSnRK genes and analysed their expression in response to the presence of the symbiont. A total of 42 PvSnRK genes were identified in P. vulgaris and annotated by comparing their sequence homology to Arabidopsis SnRK genes. Phylogenetic analysis classified the three sub-families into individual clades, and PvSnRK3 was subdivided into two groups. Chromosome localization analysis showed an uneven distribution of PvSnRK genes on 10 of the 11 chromosomes. Gene structural analysis revealed great variation in intron number in the PvSnRK3 sub-family, and motif composition is specific and highly conserved in each sub-family of PvSnRKs. Analysis of cis-acting elements suggested that PvSnRK genes respond to hormones, symbiosis and other abiotic stresses. Furthermore, expression data from databases and transcriptomic analyses revealed differential expression patterns for PvSnRK genes under symbiotic conditions. Finally, an in situ gene interaction network of the PvSnRK gene family with symbiosis-related genes showed direct and indirect interactions. Taken together, the present study contributes fundamental information for a better understanding of the role of the PvSnRK gene family not only in symbiosis but also in other biotic and abiotic interactions in P. vulgaris.

A better understanding of the signalling pathways that affect the productivity and sustainability of such crop systems is particularly important. A gene family such as SnRK with diverse regulatory mechanisms might play a pivotal role in biotic and abiotic interactions of legumes, and greater knowledge of these aspects may help in crop improvement. In the present study, we sought to identify and understand the diversity of the SnRK gene family in P. vulgaris and elucidate their differential expression patterns under rhizobial and mycorrhizal symbiotic conditions.

Phylogenetic Analysis of PvSnRK Family Genes
Multiple sequence alignment of PvSnRKs and AtSnRKs was performed using ClustalW. Based on the alignments, phylogenetic analysis of the aligned sequences was carried out using Molecular Evolutionary Genetics Analysis (MEGA XI) with the neighbour-joining (NJ) method and the JTT + I + G substitution model with 1000 bootstrap replicates and default parameters [52].
Conserved motifs in the PvSnRK gene family in P. vulgaris were identified using the Multiple Expectation Maximization for Motif Elicitation (MEME) online program (http://meme.sdsc.edu/meme/itro.html, accessed on 2 October 2021) with the following parameters: number of repetitions = any, maximum number of motifs = 20; and optimum motif length = 6 to 100 residues. The gene structure of the PvSnRK gene family was analysed using the Gene Structure Display Server online program (GSDS: http://gsds.cbi.pku.edu.ch, accessed on 10 October 2021) [53].
Sequences 2 kb upstream of PvSnRKs were downloaded from the Phytozome database. The plant transcriptional regulatory map (http://plantregmap.gao-lab.org/, accessed on 25 October 2021) was used to analyse the promoter sequences.

Calculation of Ka/Ks and Dating of Duplication Events
To identify putative orthologues between two different species (A and B), each sequence from species A was searched against all sequences from species B using BLASTN; each sequence from species B was also searched against all sequences from species A. The two sequences were defined as orthologues when reciprocal best hits were each within ≥300 bp of the two aligned sequences. The duplication period (Million Years ago, MYA) and divergence of each PvSnRK gene were calculated using the following formula: T = Ks/2λ (λ = 6.56 × 10 −9 ) [54]. The calculation of Ka and Ks was performed using TBtools. Graphs were developed with R studio with the R commander package [55].

Transcriptome Profiling and RT-qPCR Analysis
Data on differential expression of SnRK genes in P. vulgaris tissues under nitrogen treatments and after inoculation with Rhizobium tropici (CIAT899) were obtained from the PvGEA website (https://plantgrn.noble.org/PvGEA/, accessed on 12 January 2022). Previously, we performed global transcriptome profiling in P. vulgaris L. cv. Negro Jamapa roots colonized by Rhizophagus irregularis spores or Rhizobium tropici strain CIAT899 [57]. The present study uses the same transcriptomic data to obtain expression profiles of PvSnRK family genes under both types of symbiotic conditions. Heatmaps were constructed with fold-change values applying the R package (https://www.r-project.org/, accessed on 14 January 2022).
To validate the RNA-seq data, we surface-sterilized P. vulgaris L. cv. Negro Jamapa seeds and germinated them as described by Nanjareddy et al. [57]. Two-day-old germinated seedlings were transplanted into sterile vermiculite and inoculated with R. irregularis or R. tropici according to Nanjareddy et al. [58]. Total RNA preparation and RT-qPCR analysis were carried out according to Quezada et al., 2019 [59].

Identification of PvSnRK Protein Orthologues in P. vulgaris
A BLAST search was carried out using Arabidopsis SnRKs as a reference, and a total of 42 genes were identified based on conserved domains in each sub-family (Table S1). The sub-family PvSnRK1 was identified by the domain KA1 (PF02149), the PvSnRK2 sub-family by the OST domain [60], and the PvSnRK3 sub-family by the NAF/FISL domain (PF02149). The genes were named PvSnRK1.1, PvSnRK1.2, PvSnRK2.1-PvSnRK2.11 and PvSnRK3.1-PvSnRK3.29 based on sequence homology with the Arabidopsis SnRK family. The amino acid length of the 42 PvSnRK gene family members ranged from 310 aa (PvSnRK2.11) to 528 aa (PvSnRK1.2), corresponding to molecular weights of 60.54 to 35.67 kDa (Table 1). The theoretical isoelectric point of PvSnRKs (PI) ranged from 4.7 to 9.24, with PvSnRK1 sub-family members showing a basic PI, the PvSnRK2 sub-family being mostly acidic (4.7-6.65) and the PvSnRK3 sub-family being slightly acidic to highly basic (6.4-9.37). Subcellular localization analysis was carried out using ProtComp v.9.0, and the results showed PvSnRK1s to localize to the extracellular space; PvSnRK2s mostly localized to the nucleus and plasma membrane and PvSnRK3 sub-family members to the plasma membrane (Table 1). However, the subcellular localization analysis through different software showed some variation, as depicted in Table S3.

Phylogenetic and Structural Analyses
To determine the evolutionary relationship among Arabidopsis and Phaseolus SnRK superfamily genes, a phylogenetic tree was constructed using the protein sequences of 42 PvSnRK genes and 34 AtSnRK genes using the neighbour-joining (NJ) method with 1000 bootstrap replications (Figures 1 and S2). The accession numbers or locus IDs of the SnRK genes are listed in Table 1 and Supplementary Table S2. The resulting tree categorized the PvSnRKs and AtSnRKs into four clades, indicating that the ancestral genes of these two clades diverged before Brassicaceae and Fabaceae separated. In the P. vulgaris phylogeny, clade I contains PvSnRK2 sub-family members represented by the OST domain, clade II contains the PvSnRK1 sub-family containing the KA1 domain, and clades III and IV contain the NAF/FISL domain and 3 and 26 members of the PvSnRK3 sub-family, respectively ( Figure S1). The combined phylogenetic analysis of Phaseolus and Arabidopsis SnRK proteins also showed a similar distribution, whereby the proteins from two species appear scattered across the branches of the evolutionary tree, suggesting that they experienced duplications after the lineages diverged. The theoretical isoelectric point of PvSnRKs (PI) ranged from 4.7 to 9.24, with PvSnRK1 sub-family members showing a basic PI, the PvSnRK2 sub-family being mostly acidic (4.7-6.65) and the PvSnRK3 sub-family being slightly acidic to highly basic (6.4-9.37). Subcellular localization analysis was carried out using ProtComp v.9.0, and the results showed PvSnRK1s to localize to the extracellular space; PvSnRK2s mostly localized to the nucleus and plasma membrane and PvSnRK3 sub-family members to the plasma membrane (Table 1). However, the subcellular localization analysis through different software showed some variation, as depicted in Table S3.

Phylogenetic and Structural Analyses
To determine the evolutionary relationship among Arabidopsis and Phaseolus SnRK superfamily genes, a phylogenetic tree was constructed using the protein sequences of 42 PvSnRK genes and 34 AtSnRK genes using the neighbour-joining (NJ) method with 1000 bootstrap replications ( Figure 1, Figure S2). The accession numbers or locus IDs of the SnRK genes are listed in Table 1 and Supplementary Table S2. The resulting tree categorized the PvSnRKs and AtSnRKs into four clades, indicating that the ancestral genes of these two clades diverged before Brassicaceae and Fabaceae separated. In the P. vulgaris phylogeny, clade I contains PvSnRK2 sub-family members represented by the OST domain, clade II contains the PvSnRK1 sub-family containing the KA1 domain, and clades III and IV contain the NAF/FISL domain and 3 and 26 members of the PvSnRK3 sub-family, respectively ( Figure S1). The combined phylogenetic analysis of Phaseolus and Arabidopsis SnRK proteins also showed a similar distribution, whereby the proteins from two species appear scattered across the branches of the evolutionary tree, suggesting that they experienced duplications after the lineages diverged.  To examine the evolution of Phaseolus SnRK genes, their chromosomal distribution was determined. PvSnRK superfamily genes are distributed across 10 pairs of homologous chromosomes among 11 pairs in the P. vulgaris genome. PvSnRK1 genes are located on chromosomes 4 and 8 and PvSnRK2 genes on chromosomes 2, 3, 6 and 8. PvSnRK3s are located on nine chromosomes, where chromosomes 1 and 11 are the exceptions, and chromosome 3 contains the highest number of PvSnRK genes (nine PvSnRKs), followed by chromosome 8 (eight PvSnRKs). Chromosomes 1, 4, 5 and 7 each have one gene each, as shown in Figure 2 and Table 1.
protein sequences of SnRK family genes of two species. The phylogenetic tree was constructed using MEGA XI software with the neighbour-joining tree method with 1000 bootstrap values.
To examine the evolution of Phaseolus SnRK genes, their chromosomal distribution was determined. PvSnRK superfamily genes are distributed across 10 pairs of homologous chromosomes among 11 pairs in the P. vulgaris genome. PvSnRK1 genes are located on chromosomes 4 and 8 and PvSnRK2 genes on chromosomes 2, 3, 6 and 8. PvSnRK3s are located on nine chromosomes, where chromosomes 1 and 11 are the exceptions, and chromosome 3 contains the highest number of PvSnRK genes (nine PvSnRKs), followed by chromosome 8 (eight PvSnRKs). Chromosomes 1, 4, 5 and 7 each have one gene each, as shown in Figure 2 and Table 1. Intron-exon analysis was carried out to obtain better insight into the structure of PvSnRK genes. PvSnRK gene family members exhibited a great variation, from 1 to 14 introns, as shown in Figure 3 and Table 1. In the PvSnRK1 sub-family, PvSnRK1.1 and PvSnRK1.2 have 9 and 10 introns, respectively, and all members of the PvSnRK2 sub-family have 8 introns each. A great variety in introns numbers was found in the sub-family PvSnRK3, with 16 PvSnRK3 members having no introns, PvSnRK3.19 having 1 intron, PvSnRK3.21 and 3.28 having 11 introns, and PvSnRK3.4 having 14 introns; the remaining eight PvSnRK3 sub-family members have 13 introns. This divergence in introns numbers indicates that exon gain, and loss occurred during evolution of the PvSnRK gene family. These findings are corroborated by the clades in the phylogenetic analysis, where clade I contains all PvSnRK2 members, clade II has PvSnRK1 individuals, and PvSNRK3 members are divided into clade III, with genes comprising 11 and 14 introns. Finally, clade IV includes all the remaining members of PvSNRK3 ( Figure S1). Intron-exon analysis was carried out to obtain better insight into the structure of PvSnRK genes. PvSnRK gene family members exhibited a great variation, from 1 to 14 introns, as shown in Figure 3 and Table 1. In the PvSnRK1 sub-family, PvSnRK1.1 and PvSnRK1.2 have 9 and 10 introns, respectively, and all members of the PvSnRK2 sub-family have 8 introns each. A great variety in introns numbers was found in the sub-family PvSnRK3, with 16 PvSnRK3 members having no introns, PvSnRK3.19 having 1 intron, PvSnRK3.21 and 3.28 having 11 introns, and PvSnRK3.4 having 14 introns; the remaining eight PvSnRK3 sub-family members have 13 introns. This divergence in introns numbers indicates that exon gain, and loss occurred during evolution of the PvSnRK gene family. These findings are corroborated by the clades in the phylogenetic analysis, where clade I contains all PvSnRK2 members, clade II has PvSnRK1 individuals, and PvSNRK3 members are divided into clade III, with genes comprising 11 and 14 introns. Finally, clade IV includes all the remaining members of PvSNRK3 ( Figure S1). A search for conserved motifs in all 42 PvSnRK proteins using the MEME program revealed a total of 25 conserved motifs, named from 1 to 25. The identified motifs were annotated in Pfam, and the details of the putative motifs are shown in Table S4. Motifs 1-4 are designated protein kinase domains and are found in all three sub-families. Motif 8, a kinase domain associated with kinase1 and KA1, and motifs 10, 11, and 20, which encode an NAF domain, are only present in PvSnRK3. Ubiquitin-associated domain 20 was found in PvSnRK1 and PvSnRK3. Furthermore, motifs 5 and 13, designated protein superfamily kinase domains, were found only in the PvSnRK1 and PvSnRK3 sub-families ( Figure 4). The remaining motifs were not annotated functionally. A search for conserved motifs in all 42 PvSnRK proteins using the MEME program revealed a total of 25 conserved motifs, named from 1 to 25. The identified motifs were annotated in Pfam, and the details of the putative motifs are shown in Table S4. Motifs 1-4 are designated protein kinase domains and are found in all three sub-families. Motif 8, a kinase domain associated with kinase1 and KA1, and motifs 10, 11, and 20, which encode an NAF domain, are only present in PvSnRK3. Ubiquitin-associated domain 20 was found in PvSnRK1 and PvSnRK3. Furthermore, motifs 5 and 13, designated protein superfamily kinase domains, were found only in the PvSnRK1 and PvSnRK3 sub-families ( Figure 4). The remaining motifs were not annotated functionally.

Ka/Ks and Gene Duplication
To further explore evolutionary constraints on Phaseolus PvSnRK genes, synonymous (Ks) and nonsynonymous (Ka) substitutions per site and their ratio (Ka/Ks) and divergence time of paralogous and orthologous SnRK family genes were calculated for AtSnRK orthologues of PvSnRKs (Tables 2 and S5). The Ka/Ks ratio among all SnRK sequences was lower than 1, indicating purifying selection. These Ka/Ks ratios suggest the conservation of SnRK homologues in terms of both sequence and biological function [61].

Ka/Ks and Gene Duplication
To further explore evolutionary constraints on Phaseolus PvSnRK genes, synonymous (Ks) and nonsynonymous (Ka) substitutions per site and their ratio (Ka/Ks) and divergence time of paralogous and orthologous SnRK family genes were calculated for AtSnRK orthologues of PvSnRKs (Table 2 and Table S5). The Ka/Ks ratio among all SnRK sequences was lower than 1, indicating purifying selection. These Ka/Ks ratios suggest the conservation of SnRK homologues in terms of both sequence and biological function [61].

Cis-Elements in Promoter Regions of PvSnRKs
To determine the gene expression pattern of PvSnRKs, the 2 kb region upstream of the CDS was analysed using the PlantRegMap database. Among all transcription factors recorded, ERF and MYB were found to be the most abundant. PvSnRK3.7 contained the greatest number of cis-elements (607) in the examined regulatory region, with 286 ERF binding sites and 49 MYB and 44 C2H2 sites. PvSnRK 3.20 has 338 TF sites; this was followed by PvSnRK 3.19 with 318 TFs, with 111 TFs being ERF TFs and 56 being bHLH TFs, and PvSnRK 1.2, with 312 TFs, with 100 being NAC TFs, 45 being MYB TFs and 40 being ERF TFs ( Figure 5, Table S6, Figure S3). The most abundant TFs, ERF, C2H2, bHLH, NAC and MYB identified in PvSnRKs were also found by symbiosis related studies in other species such as M. truncatula and L. japonicus (Table 3).

Cis-Elements in Promoter Regions of PvSnRKs
To determine the gene expression pattern of PvSnRKs, the 2 kb region upstream of the CDS was analysed using the PlantRegMap database. Among all transcription factors recorded, ERF and MYB were found to be the most abundant. PvSnRK3.7 contained the greatest number of cis-elements (607) in the examined regulatory region, with 286 ERF binding sites and 49 MYB and 44 C2H2 sites. PvSnRK 3.20 has 338 TF sites; this was followed by PvSnRK 3.19 with 318 TFs, with 111 TFs being ERF TFs and 56 being bHLH TFs, and PvSnRK 1.2, with 312 TFs, with 100 being NAC TFs, 45 being MYB TFs and 40 being ERF TFs ( Figure 5, Table S6, Figure S3). The most abundant TFs, ERF, C2H2, bHLH, NAC and MYB identified in PvSnRKs were also found by symbiosis related studies in other species such as M. truncatula and L. japonicus (Table 3).  bHLH Controls a diverse processes from cell proliferation to cell lineage establishment Nodule organogenesis [64] bZIP Plant bZIP proteins preferentially bind to DNA sequences with an ACGT core Nodule organogenesis [65] C2H2 The majority of such proteins characterized to date are DNA-binding transcription factors, and many have been shown to play crucial roles in the development of plants, animals and fungi Symbiosome development [66] Dof Regulation of gene expression in processes such as seed storage protein synthesis in developing endosperm, light regulation of genes involved in carbohydrate metabolism, plant defense mechanisms, seed germination, gibberellin response in post-germinating aleurone , auxin response and stomata guard cell specific gene regulation No report  bHLH Controls a diverse processes from cell proliferation to cell lineage establishment Nodule organogenesis [64] bZIP Plant bZIP proteins preferentially bind to DNA sequences with an ACGT core Nodule organogenesis [65] C2H2 The majority of such proteins characterized to date are DNA-binding transcription factors, and many have been shown to play crucial roles in the development of plants, animals and fungi Symbiosome development [66] Dof Regulation of gene expression in processes such as seed storage protein synthesis in developing endosperm, light regulation of genes involved in carbohydrate metabolism, plant defense mechanisms, seed germination, gibberellin response in post-germinating aleurone, auxin response and stomata guard cell specific gene regulation No report  No report MYB transcription factors promote expression of genes involved in cell proliferation and differentiation. ERFs are involved in regulation of developmental processes in response to stimuli, and NAC, C2H2 and bHLH are involved in the pathogen response, cell proliferation and development.

GO Analysis
Gene Ontology analysis of all PvSnRK genes showed them to be involved in biological processes and molecular functions, but not cellular components. The biological processes involved related to PvSnRKs are cellular signalling, phosphorylation, and metabolic processes such as protein, macromolecule, and phosphorous metabolism. Among molecular functions, the majority are associated with binding and catalytic activities, followed by transferase and nucleotide-binding functions (Figure 6a). RT-qPCR data are the averages of three biological replicates (n > 9). The statistical significance of differences between mycorrhized and nodulated roots was determined using an unpaired two-tailed Student's t-test (** p < 0.01). Error bars represent means ± Standard error mean (SEM).

Expression Profiles of PvSnRK Genes in Different Tissues
Differential expression data for PvSnRK genes were obtained from PvGEA: Common Bean Gene Expression Atlas and Network Analysis (https://plantgrn.noble.org/PvGEA/, accessed on 16/11/2021). The expression patterns of all 42 PvSnRK members were analysed in 25 different tissues, including leaves, stems, shoots, pods, seeds, roots (inoculated and uninoculated with Rhizobium) and nodules, at different developmental stages and treatments (Figure 6b, Table S7). The lowest expression of all PvSnRK genes was found in young flower, seed and shoot tissues; at least one of the 42 genes was expressed in the remaining tissues. Among nodulation-related tissues, most of the PvSnRK genes showed high expression, except for the PvSnRK1 sub-family, in whole roots separated from fix+ nodules collected at 21 days after inoculation (PvRE). Pre-fixing (effective) nodules collected at 5 days after inoculation (PvN5) showed very high expression of only PvSnRK3. 5 and PvSnRK3.29 (Figure 6b).

Expression Patterns of PvSnRK Genes under Symbiotic Conditions
Additionally, differential expression patterns of PvSnRKs under rhizobial and mycorrhizal symbiotic conditions were analysed by transcriptomics at the early infection stages 5 dpi and 7 dpi, respectively. Most of the genes encoding PvSnRK1s and PvSnRK3s Candidate genes were selected and corresponding transcript accumulation under mycorrhized and nodulated conditions was quantified by RT-qPCR. RT-qPCR data are the averages of three biological replicates (n > 9). The statistical significance of differences between mycorrhized and nodulated roots was determined using an unpaired two-tailed Student's t-test (** p < 0.01). Error bars represent means ± Standard error mean (SEM).

Expression Profiles of PvSnRK Genes in Different Tissues
Differential expression data for PvSnRK genes were obtained from PvGEA: Common Bean Gene Expression Atlas and Network Analysis (https://plantgrn.noble.org/PvGEA/, accessed on 16 November 2021). The expression patterns of all 42 PvSnRK members were analysed in 25 different tissues, including leaves, stems, shoots, pods, seeds, roots (inoculated and uninoculated with Rhizobium) and nodules, at different developmental stages and treatments (Figure 6b, Table S7). The lowest expression of all PvSnRK genes was found in young flower, seed and shoot tissues; at least one of the 42 genes was expressed in the remaining tissues. Among nodulation-related tissues, most of the PvSnRK genes showed high expression, except for the PvSnRK1 sub-family, in whole roots separated from fix+ nodules collected at 21 days after inoculation (PvRE). Pre-fixing (effective) nodules collected at 5 days after inoculation (PvN5) showed very high expression of only PvSnRK3.5 and PvSnRK3.29 (Figure 6b).

PvSnRK Protein-Protein Interaction Network Prediction
A protein interaction network was constructed for PvSnRK proteins based on orthologues in Arabidopsis using the STRING database with the highest bit score (0.9). While predicting interactions among the PvSnRK proteins, we used experimentally proven, coexpressed and co-occurring protein interactions, revealing a total of 40 interacting proteins. The majority of them are phosphoprotein phosphatase (PPP) family proteins, PP2A (protein phosphatase PP1-α catalytic subunit), PP2CA/PP2CB (protein phosphatase 2A catalytic subunit α/β isoform), PPP2RA/PPP2RB (protein phosphatase 2A 65 kDa regulatory subunit A β isoform) and PPP1C (protein phosphatase PP1-α catalytic subunit)-type proteins. Phosphatases are proteins involved in substrate recognition, plant signalling pathways such as stress regulation, light, pathogen defence and hormonal signalling, the cell cycle, differentiation, metabolism, and plant growth.
The next important interactions were with cell cycle proteins and cyclin B proteins. Cyclin B proteins have been implicated in plant growth and development. All PvSnRK protein sub-families were found to interact with PPP family members and cyclin B proteins through PKGII, indicating their involvement in various cellular and developmental processes. SnRK2 sub-family members were found to specifically interact with PP2C proteins (Figure 7). were found to be upregulated under both symbiotic conditions, though PvSnRK2 was less induced during symbiosis with rhizobia. It was particularly interesting to find highly induced PvSnRK1.2 in comparison to PvSnRK1.1, and among the PvSnRK3 sub-families, PvSnRK3.10 was highly induced, followed by PvSnRK3.7, PvSnRK3.20, and PvSnRK3. 22. Under mycorrhizal inoculation conditions, PvSnRK1s and PvSnRK3s were all induced, and some PvSnRK2 genes were also induced. In the SnRK1 sub-family, PvSnRK1.2 was more induced; among PvSnRK2 genes, PvSnRK2.6, PvSnRK2.7 and PvSnRK2.11 were least induced. In the PvSnRK3 sub-family, PvSnRK3.1, PvSnRK3.2, PvSnRK3.3, PvSnRK3.4, PvSnRK3.5 and PvSnRK3.6 were less induced than other members (Figure 6c & 6d).

PvSnRK Protein-Protein Interaction Network Prediction
A protein interaction network was constructed for PvSnRK proteins based on orthologues in Arabidopsis using the STRING database with the highest bit score (0.9). While predicting interactions among the PvSnRK proteins, we used experimentally proven, co-expressed and co-occurring protein interactions, revealing a total of 40 interacting proteins. The majority of them are phosphoprotein phosphatase (PPP) family proteins, PP2A (protein phosphatase PP1-α catalytic subunit), PP2CA/PP2CB (protein phosphatase 2A catalytic subunit α/β isoform), PPP2RA/PPP2RB (protein phosphatase 2A 65 kDa regulatory subunit A β isoform) and PPP1C (protein phosphatase PP1-α catalytic subunit)-type proteins. Phosphatases are proteins involved in substrate recognition, plant signalling pathways such as stress regulation, light, pathogen defence and hormonal signalling, the cell cycle, differentiation, metabolism, and plant growth.
The next important interactions were with cell cycle proteins and cyclin B proteins. Cyclin B proteins have been implicated in plant growth and development. All PvSnRK protein sub-families were found to interact with PPP family members and cyclin B proteins through PKGII, indicating their involvement in various cellular and developmental processes. SnRK2 sub-family members were found to specifically interact with PP2C proteins (Figure 7).

Figure 7.
In silico prediction of protein-protein interactions among SnRK1s, SnRK2s and SnRK3s using the Cytoscape tool based on the Pearson correlation coefficients of the relative expression of the gene. Each node represents a protein, and each edge refers an interaction. The red-coloured box represents SnRK1s, blue represents SnRK2s, and green represents SnRK3s.

Prediction of Coregulatory and Interaction Networks of PvSnRK and Symbiotic Genes
The symbiotic interaction between legumes and rhizobia is unique and involves a set of common symbiotic genes that regulate root nodule symbiosis and mycorrhizal Figure 7. In silico prediction of protein-protein interactions among SnRK1s, SnRK2s and SnRK3s using the Cytoscape tool based on the Pearson correlation coefficients of the relative expression of the gene. Each node represents a protein, and each edge refers an interaction. The red-coloured box represents SnRK1s, blue represents SnRK2s, and green represents SnRK3s.

Prediction of Coregulatory and Interaction Networks of PvSnRK and Symbiotic Genes
The symbiotic interaction between legumes and rhizobia is unique and involves a set of common symbiotic genes that regulate root nodule symbiosis and mycorrhizal symbiosis. Another 190 genes have been implicated in regulating symbiotic interactions. A total of 200 genes known to be involved in symbiosis were chosen to predict an interaction network with each of the PvSnRK sub-families in the STRING database. Interaction between PvSnRK1 sub-family members and symbiosis-related proteins shows that PvSnRK1.1 and PvSnRK1.2 do not interact directly with any of the symbiosis-related genes. However, they do interact with TOR, PI3K and ATG-RP, which interact with other symbiotic genes represented by 54 nodes and 119 edges. The network mostly involves nucleoporins (NUPs), RAB proteins, auxin-responsive ARP proteins and many proteins involved in vesicle transport (Figure 8a). The largest PvSnRK3 sub-family, with 29 members, interacts with TOR, PI3K and ATG-RP at the first level, similar to PvSnRK1 members. The predictions showed 102 nodes and 319 edges, indicating a larger network. CCS interacts with nucleoporins (NUPs), and through GNBI, they interact with several RPSs (ribosomal proteins) and a variety of symbiosis-related proteins. Network prediction showed that PvSnRK3 sub-family members may be involved in regulating symbiosis in Phaseolus via different signalling pathways ( Figure 9). The largest PvSnRK3 sub-family, with 29 members, interacts with TOR, PI3K and ATG-RP at the first level, similar to PvSnRK1 members. The predictions showed 102 nodes and 319 edges, indicating a larger network. CCS interacts with nucleoporins (NUPs), and through GNBI, they interact with several RPSs (ribosomal proteins) and a variety of symbiosis-related proteins. Network prediction showed that PvSnRK3 sub-family members may be involved in regulating symbiosis in Phaseolus via different signalling pathways ( Figure 9).

Figure 9.
In silico prediction of protein-protein interactions among SnRK3 subfamily proteins with symbiosis-related proteins using the Cytoscape tool based on the Pearson correlation coefficients of the relative expression of the gene. Each node represents a protein, and each edge refers an interaction. The different colours represent the different gene families.

Discussion
The SnRK gene family, serine/threonine kinases, and its orthologues in animals and yeast are highly conserved. In plants, they have been identified as regulators of abiotic and biotic stress [80][81][82]. In recent years, genome-wide analysis of this gene family in an array of both model and economically important plant species has focused primarily on abiotic stress. Phaseolus vulgaris is the most important legume consumed by humans worldwide, as it is an affordable source of proteins, vitamins, minerals, antioxidants, and bioactive compounds [83]. Climate change has impacted the world's crop yield and quality, leading to socioeconomic and food insecurity [84]. The yield quantity and quality of N2-fixing crops can be improved by several agronomic practices, such as irrigation, sowing density and Rhizobium application. Since common bean does not need exogenous N fertilizer application, productivity is cost effective. Any efforts towards the betterment of rhizobial association to improve N fixation in crops such as Phaseolus should be undertaken. The present investigation encompasses the identification, classification, and analysis of the expression patterns of the SnRK gene family under symbiotic conditions. Genome-wide identification studies have been carried out, with early reports documenting the presence of various numbers of members in both monocots and dicots. A total of 39, 114, 30, 34, 26, 149, 46, 48 and 44 SnRK genes have been identified in A. thaliana [1], Brassica napus [30], C. sativus [31], E. grandis [32], F. ananassa [33], T. aestivum [34], H. vulgare [35], O. sativa [17] and B. distachyon [36], respectively. In Phaseolus, we identified 42 SnRK genes with 2 SnRK1s, 11 SnRK2s and 29 SnRK3s with the characteristic domains of each sub-family. As in any other species, the SnRK2 and SnRK3 sub-families in Phaseolus are larger than SnRK1, supporting the view that these two sub-families originated from duplication of SnRK1 [35]. The nonmotile nature of plants exposes them to biotic and abiotic factors, and plants adopt such changes by expanding their genes and gene families. Gene and genome duplications are important events that contribute to polyploidy and genome evolution. One or multiple polyploidies are prevalent in angiosperms [85,86] and explain the large number of SnRK2 and SnRK3 sub-family members in all reported plant species.
In Phaseolus, the SnRK gene family is distributed on 10 of 11 chromosomes. This distribution is unlike in other species, in which all chromosomes in the genome contain at

Discussion
The SnRK gene family, serine/threonine kinases, and its orthologues in animals and yeast are highly conserved. In plants, they have been identified as regulators of abiotic and biotic stress [80][81][82]. In recent years, genome-wide analysis of this gene family in an array of both model and economically important plant species has focused primarily on abiotic stress. Phaseolus vulgaris is the most important legume consumed by humans worldwide, as it is an affordable source of proteins, vitamins, minerals, antioxidants, and bioactive compounds [83]. Climate change has impacted the world's crop yield and quality, leading to socioeconomic and food insecurity [84]. The yield quantity and quality of N2-fixing crops can be improved by several agronomic practices, such as irrigation, sowing density and Rhizobium application. Since common bean does not need exogenous N fertilizer application, productivity is cost effective. Any efforts towards the betterment of rhizobial association to improve N fixation in crops such as Phaseolus should be undertaken. The present investigation encompasses the identification, classification, and analysis of the expression patterns of the SnRK gene family under symbiotic conditions. Genome-wide identification studies have been carried out, with early reports documenting the presence of various numbers of members in both monocots and dicots. A total of 39, 114, 30, 34, 26, 149, 46, 48 and 44 SnRK genes have been identified in A. thaliana [1], Brassica napus [30], C. sativus [31], E. grandis [32], F. ananassa [33], T. aestivum [34], H. vulgare [35], O. sativa [17] and B. distachyon [36], respectively. In Phaseolus, we identified 42 SnRK genes with 2 SnRK1s, 11 SnRK2s and 29 SnRK3s with the characteristic domains of each sub-family. As in any other species, the SnRK2 and SnRK3 sub-families in Phaseolus are larger than SnRK1, supporting the view that these two sub-families originated from duplication of SnRK1 [35]. The nonmotile nature of plants exposes them to biotic and abiotic factors, and plants adopt such changes by expanding their genes and gene families. Gene and genome duplications are important events that contribute to polyploidy and genome evolution. One or multiple polyploidies are prevalent in angiosperms [85,86] and explain the large number of SnRK2 and SnRK3 sub-family members in all reported plant species.
Each of the SnRK sub-families has a characteristic domain; however, the N-terminal kinase domain is highly conserved across species and sub-families. Exon-intron structural diversification and motif composition play an important role in the evolution and function of many gene families [87]. PvSnRK1 sub-family members have 10-11 exons, such as BdSnRK1s, EgrSnRK1s and CsSnRK1s. All PvSnRK2s have nine exons, similar to most HvSnRKs, HcSnRK2s, EgrSnRK2s, VvSnRK2s, AtSnRK2s and BdSnRK2s, indicating the conserved nature of these sub-families. The sub-family PvSnRK3 can also be subdivided into exon-rich and exon-poor types, as can all other species reported thus far. Reports suggest the origin of the SnRK3 sub-family from green algae, and the intron-poor group first appeared in seed plants [88]. These results are consistent with phylogenetic tree studies showing that SnRK exon-intron numbers are highly conserved during the evolution of each sub-family. The phylogenetic tree and exon-intron structure showed that most paralogous gene pairs contain the same exon number, though some gene pairs have different exon numbers. Motif analysis revealed a close evolutionary relationship within sub-groups due to the conserved nature of motif composition among sub-families. The motif structure of each sub-family might define the biological function in which they are involved. The gene structure and sequence conservation were similar to most of the previously studied species [1,17,[30][31][32][33][34][35][36].
Gene expression is regulated by external factors that are perceived through signalling mechanisms. Such signals activate specific transcription factors that, when combined with cis-acting elements in the promoter regions of genes, alter gene expression. In most genomewide analysis studies of SnRK gene families, detection of cis-acting elements has focused on abiotic stress conditions in which the presence of hormone-, salt-, temperature-, and drought-specific cis-elements [1,17,[30][31][32][33][34][35][36]. As the aim of our investigation was to predict the possible role of the PvSnRK gene family in symbiosis, we analysed symbiosis-related cis-elements as demonstrated in Medicago and L. japonicus. All PvSnRK gene sub-families contain symbiosis-related cis-elements. Among all cis-elements, ERF and MYB are the most abundant, followed by C2H2 and bHLH.
Ka/Ks analysis showed that PvSnRK gene family duplications either occurred slowly or are highly conserved [89]. Gene Ontology studies revealed that the PvSnRK genes function mostly in molecular and biological processes, specifically in cellular and metabolic processes followed by nucleotide binding activity. When we analysed expression of PvSnRK gene family members in different tissues, the lowest expression of any PvSnRK was found in aerial tissues such as young pods, flowers, seeds, and shoots. Most of the SnRKs were found to be expressed in root or root nodules at some stage of their development. Furthermore, transcriptomic data under early symbiotic conditions showed elevated expression levels of PvSnRK1 s and PvSnRK3s, with some members being more highly induced than others. For some genes, these expression patterns were found under both root nodule and mycorrhizal symbiotic conditions, suggesting a possible role for PvSnRKs in the establishment of symbiotic relationships.
To predict possible interaction networks among the PvSnRKs and PvSnRKs with symbiosis-related genes, we carried out in situ interaction network building using the STRING database and Cytoscape. The results were interesting, as most of the PvSnRKs interact among themselves, and interaction is mediated by master regulators of cellular processes such as TOR and PKGII. Through these proteins, PvSnRKs interact with several protein phosphatases and cell cycle-regulating cyclins. We chose symbiosis-related genes based on a previously published article and identified a total of 200 genes in the Phaseolus genome. PvSnRK1s mostly interact with other major proteins, such as TOR and PI3Ks, which are connected to the NUPs, RAB proteins, ARPs and proteins involved in vesicle transport. On the other hand, PvSnRK2s interact with PKGII, which interact with the cell cycle regulatory proteins cyclins, CDC, and APC. In contrast, PvSNRK2s shows very few symbiotic gene interactions. Of the largest PvSnRK sub-families, PvSnRK3 interacts with NUPs, RPSs and many symbiotic genes.
Taken together, our extensive analysis of the PvSnRK gene family revealed structural conservation of genes across species and possible functional conservation as well. Furthermore, the principle aim of this study was to understand the putative role of PvSnRKs in regulating symbiotic interactions between Phaseolus and Rhizobium or Phaseolus; mycorrhiza are promising, as some genes contain specific cis-elements and showed transcript upregulation in response to symbionts. Finally, the identified in situ protein-protein interactions may help in predicting candidate genes for functional characterization to obtain a clear picture of the regulatory mechanisms involved.