The National Center for Biotechnology Information (NCBI) Genome database (https://www.ncbi.nlm.nih.gov/home/genomes/) was browsed by organism to collect protein genomes of 71 Enterobacter strains (17 species, complete draft). The 71 Enterobacter strains which were investigated included; 27 Enterobacter hormaechei strains, 1 Enterobacter cancerogenus strain, 6 Enterobacter ludwigii strains, 6 Enterobacter roggenkampii strains, 2 Enterobacter bugandensis strains, 13 E. cloacae strains, 4 Enterobacter asburiae strains, 3 Enterobacter kobei strains, 1 Enterobacter chengduensis strain, 1 Enterobacter sichuanensis strain, 1 Enterobacter oligotrophicus strain, 1 E. N18-03635 strain, 1 E. JBIWA003 strain, 1 E. JBIWA005 strain, 1 E. JBIWA008 strain, 1 E. DNB-S2 strain and 1 E. SGAir0187 strain.

Clp ATPases were mined from different Enterobacter strains by first obtaining the protein file for each strain from NCBI. Each file was then converted from.gz to FASTA format using an online server (https://www.convertfiles.com/file_type_description/GZ_Archive_FIle.html). The FASTA files were then individually searched for the presence of Clp ATPases. The Clp ATPase sequences were separated from the main file and used for further analyses.

The Clp ATPase protein sequences were aligned using MAFFT v6.864 embedded on the T-rex web server (http://www.trex.uqam.ca/) (Katoh et al., 2005; Boc et al., 2012). The alignments were deduced by the T-rex server, following this the file for the best tree was downloaded, visualized and colored using the interactive Tree of Life (iTOL) v6 server (https://itol.embl.de/) (Letunic and Bork, 2021). The phylogenetic tree (Neighbor-joining, BLOSUM62) of ClpV in comparison with ClpVI from Pseudomonas aeruginosa (UniProt: Q9I742) and Vibrio cholera ClpV1 (UniProt: A0A395TZY1) was constructed using JalView v 2.11.2.2 (Waterhouse et al., 2009a).

Clustal omega (https://www.ebi.ac.uk/jdispatcher/msa/clustalo) was used in four ways during the investigation. Firstly, the percentage identity between the various Clp ATPase classes was analyzed. Secondly, the percentage identity between the ClpV mined from the E. cloacae strains, P. aeruginosa ClpV1 (UniProt: Q9I742), and V. cholera ClpV1 (UniProt: A0A395TZY1) was compared. Thirdly, the percentage identity between the E. cloacae E3442 ClpV homologs were analyzed. Lastly, the percentage identity between the three E3442 ClpV proteins and V. cholera ClpV1 (UniProt: A0A395TZY1) was investigated (Sievers et al., 2011). The domains were assigned to the ClpV proteins using the annotations provided through InterPro (Blum et al., 2020).

The three ClpV proteins mined from E. cloacae E3342 were modelled using the template 1qvrA through Swiss Model (https://swissmodel.expasy.org/) (Waterhouse et al., 2018). The structures were refined using the Protein REFinement via Molecular Dynamics (PREFMD) web server (https://feig.bch.msu.edu/prefmd) (Heo and Feig, 2017). Structural assessment of the refined structures was then performed using the Swiss-Model assessment tool (https://swissmodel.expasy.org/assess) and MolProbity (http://molprobity.biochem.duke.edu/) (Waterhouse et al., 2018; Williams et al., 2018). The refined structures were superimposed with the modelling template (PDB 1qvrA) using PyMol (Schrödinger, 2020).

The ProtParam tool on the Expasy server (https://web.expasy.org/protparam/) was used to compare the physicochemical properties of the three ClpV proteins using the respective protein sequences (Gasteiger et al., 2003). ClpV protein disorder was predicted using the IUPred3 server (https://iupred3.elte.hu/) (Erdős et al., 2021). ClpV sequences were analyzed using Mafft on the Jalview software and the InterPro server (http://www.ebi.ac.uk/interpro/) was used to investigate the protein families, domains and function of the ClpV sequences (Waterhouse et al., 2009b; Apweiler et al., 2001). Normal Mode Analysis (NMA) from the DynaMut server (https://biosig.lab.uq.edu.au/dynamut/analysis) was performed to compare the conformation of the ClpV proteins (Rodrigues et al., 2018). The ClpV sequences were annotated as per the InterPro server and Lee et al. (2003). The 3D structures of the ClpV proteins were then compared using the Dali server (http://ekhidna2.biocenter.helsinki.fi/dali/) (Holm, 2022). Each ClpV protein sequence was analyzed using the STRING v11.5 server (https://string-db.org/) for protein-protein association analysis (Szklarczyk et al., 2019). The minimum required interaction score was set to medium (0.400). DogSite3 was used through the ProteinsPlus server (https://proteins.plus/) to predict binding sites on the ClpV proteins (Graef et al., 2023; Schöning-Stierand et al., 2022; Schöning-Stierand et al., 2020; Fährrolfes et al., 2017).

Molecular dynamic (MD) simulations were performed using Maestro v12.2 through the implemented GPU-enabled Desmond molecular dynamics simulation engine (Schrodinger, 2006). The modelled ClpV structures were saved as Protein Data Bank (PDB) files and submitted to the Linux (Ubuntu) desktop server for molecular simulations studies. The system builder module was used to incorporate the TIP3P solvent model with the OPLS forcefield. The proteins were placed in an orthorhombic box (distance from the box face to the outermost protein atom was set to 10 Å, the box angle was α = β = γ = 90°). The volume box containing the proteins was minimized, and counter ions were added to neutralize the system. In addition, 0.15 M NaCl was added into the solvent box for physiological conditioning. The system was then submitted for MD simulations for 250 ns.

Post dynamic analyses of the trajectories derived from the MD simulation studies were performed using Schrodinger Maestro v12.2. Firstly, Simulation Quality Analysis was used to analyze the quality of simulations through the average energy, pressure, temperature, and volume analysis. Secondly, Simulation Interaction Diagram algorithm was used to analyze the RMSD of the alpha carbon atoms (Cα), the root mean square fluctuation (RMSF) of the residues, and secondary structure element analysis. Lastly, the Simulation Events Analysis algorithm (implemented in Maestro v12.2) was used to calculate the radius of gyration (Rg).

Advancement into sequencing technologies enables researchers to study the presence and distribution of genes in microorganisms to shed light into relationships amongst species (Pareek et al., 2011). Clp ATPases are essential for bacterial survival and pathogenicity and are emerging as potential drug targets (Culp and Wright, 2017). However, it was observed that the presence and distribution of Clp ATPases in the Enterobacter species has not been fully investigated hence, genome wide mining was performed to gain an insight into the diversity of Clp ATPases in the Enterobacter species.

The data mined and tallied for the 71 Enterobacter species shows that all the studied species have both Class I (ClpA and ClpB) and Class II (ClpX and ClpV) Clp ATPase members (Figure 1). It was also observed that fewer Enterobacter species (20%) contained ClpK, a class II member of Clp ATPases that has been recently linked to thermotolerance properties of another ESKAPE pathogen member, Klebsiella pneumonia (Bojer et al., 2010; Motiwala et al., 2021). It was very interesting to observe from genome data mining that all investigated strains contained one gene copy for ClpA, ClpB and ClpX. Contrastingly, certain Enterobacter strains contained more than one copy of the ClpV gene.

It then became necessary to ascertain the evolutionary relationship of the Clp ATPase species in Enterobacter species. Subsequently, we constructed a phylogenetic tree of the mined Clp ATPases. Figure 2 showed that the ClpA, ClpB, ClpK and ClpV proteins diverge from one point. This is expected as all these proteins contain two nucleotide domains and belong to Clp ATPase Class I. Therefore, one would assume that they would be more closely related to each other than to ClpX which belongs to Clp ATPase Class II (Figure 2A). Analysis of the phylogenetic tree shows that there is not a great amount of genetic variation among the ClpA, ClpB, ClpX and ClpK proteins, however, there appears to be significant genetic variation among the ClpV proteins (Figure 2A). This correlates with Clustal omega analysis of each protein group (Figure 2B).

The role of ClpV proteins has been widely studied in Gram-negative opportunistic pathogens; V. cholera and P. aeruginosa (Smith-Roberge et al., 2018; Sun et al., 2014; Lesic et al., 2009; Li et al., 2020). The ClpV proteins mined from various E. cloacae strains were compared to ClpV, from V. cholera and P. aeruginosa through Clustal omega and phylogenetic analysis. This was done to elucidate whether the ClpV proteins from the E. cloacae strains may have a similar function to those of ClpV1 from V. cholera and P. aeruginosa. Clustal omega showed that ClpV from the various E. cloacae strains had between 43.56% and 100% similarity with ClpV1 from P. aeruginosa, and between 43.31% and 100% similarity with ClpV1 from V. cholera. ClpV from E. cloacae had a percentage identity of greater than 30% in comparison to ClpV1 from both V. cholera and P. aeruginosa, therefore indicating that the sequences are homologous (Pearson, 2013). Additionally, analysis of the phylogenetic tree comparing the evolution of ClpV from E. cloacae and ClpV1 from both V. cholera and P. aeruginosa showed that these proteins all shared a common ancestor, indicated by the number 1 (Figure 3). Figure 3 shows that a group of ClpV from E. cloacae (orange box) are more closely related to ClpVI from V. cholera and P. aeruginosa compared to the other group of ClpV from E. cloacae (blue box). Clustal omega analysis showed that three mined ClpV proteins from the E3442 strain had between 46.82% and 100% similarity, this correlated with the trend observed in Figure 3 which shows that these proteins lie on different branches of the phylogenetic tree.

To further investigate the uniqueness of the ClpV proteins identified within E. cloacae strain E3442, the three ClpV proteins were compared through homology modelling. The ClpV proteins were modelled using the 1qvrA template and the model quality was assessed (Figure 4; Table 1). The alignment of ClpV1, ClpV2, and ClpV3 with the template (1QVR-A) revealed a root mean square deviation (RMSD) score of 0.270, 0.283 and 0.347, respectively. These scores are close to 0 and indicate that there is a good alignment between the coordinates of the model and template. Therefore, the proteins are structurally similar to the template (Carugo and Pongor, 2001). Secondary and tertiary analysis of the proteins show that they are mainly α-helical (Figure 4).

The Rama Z-score provides an overall insight into the quality of a modelled protein by comparing the model to a reference set of high-resolution structures (Sobolev et al., 2020). The Rama Z-scores of the modelled proteins were less than two and therefore the modelled structures are confirmed to be of adequate quality (Table 1). The Ramachandran favored scores of the three modelled proteins are greater than 80%, this indicates that the generated models can be trusted and used for further analysis (Maxwell and Popelier, 2017). The QMeanisCo uses statistical potentials of mean force to quantify model quality and modelling errors (Studer et al., 2019). The QMeanisCo values for the three modelled proteins are close to 1, thus once again confirming that the models are of adequate quality (Table 1).

Table 2 shows the physicochemical properties of the three ClpV proteins. The instability index of all three proteins is greater than 40, the proteins are therefore considered to be unstable. However, it is essential to remember that the stability of the protein is not only determined by the intrinsic nature of the protein but also by the protein environment (Gamage et al., 2019). The high aliphatic indexes indicate that these proteins are thermally stable. Additionally, the negative Grand average of Hydropathicity (GRAVY) score indicates that the proteins are all hydrophilic (Panda and Chandra, 2012). Figure 5 shows the protein disorder prediction, the disorder score for the residues of all three proteins lie below one and the three proteins are not considered to be disordered. Therefore, these proteins can be expressed and purified (Deng et al., 2012).

Following analysis of the physicochemical characteristics and disorder prediction of the proteins, the difference in the protein sequence and 3D structure were further investigated. Figure 6A shows significant sequence similarity and conservation, especially between ClpV2 and ClpV3. The difference observed between the protein sequences indicate that these proteins may have variable enzyme activity and binding affinities. The 3D structure analysis corresponds with the primary sequence analysis in that ClpV2 and ClpV3 fall on the same branch of the phylogenetic tree (Figures 6B,C). Additionally, the Dali Z-score of 44.3 is greater than 20 and confirms that the three ClpV structures are homologous to each other (Holm et al., 2008).

The InterPro server analysis identified protein families, domains, and functional sites in ClpV proteins (Apweiler et al., 2001). ClpV1 and ClpV2 contain the repeat domain, which plays a key role in substrate specificity, a characteristic feature of the Clp A/B family (Lo et al., 2001). The nucleotide binding domains (NDB I and II) which are crucial for ATP binding and are well conserved among Clp A/B family members (Gottesman et al., 1990). They exhibit high sequence identity, with only minor variations across the three ClpV proteins. A highly conserved site (highlighted in orange) is situated in NBDI and remains identical across the three ClpV proteins (Figure 7). Similar to ClpB and ClpK, ClpV proteins contain a middle linker region that separates the NBDs and is essential for chaperone activity (Lee et al., 2003; Motiwala et al., 2021). Additionally, the AAA + lid domain, located in the middle domain of ClpV proteins, acts as a regulatory “gatekeeper” for substrate access through the central pore (Miller and Enemark, 2016).

Further conserved regions include the Walker A (highlighted in pink) and Walker B (highlighted in red) motifs, located within NBD I and NBD II. The Walker A motif, which is responsible for ATP binding, is largely conserved, though a single amino acid variation suggests potential differences in ATP interactions among the three proteins. The Walker B motif, which facilitates metal ion binding, remains identical and homologous across ClpV proteins, indicating a shared metal ion interaction mechanism (Lee et al., 2003; Motiwala et al., 2021). Conserved residues (highlighted in yellow) further support the structural and functional similarity of the ClpV proteins (Figure 7).

Potential binding pockets on the three ClpV proteins were analyzed using the DogSite3 module on the ProteinsPlus server. This was done following confirmation on the difference in the primary and 3D structure of the proteins. Interestingly, the best binding pocket was observed to be around the same area for all three proteins with ClpV1 having the biggest binding pocket (Figure 8). The difference in the size of the binding pockets indicates that ClpV1 may bind bigger proteins or ligands, while ClpV2 and ClpV3 bind to ligands or proteins of similar size. The difference in size between the binding pockets could also be attributed to the presence of the repeat domain (Figure 7).

Normal Mode Analysis (NMA) was performed to investigate the conformational changes between the ClpV proteins. The color transformations from blue to red in Figure 9 represents the conformational changes and flexibility of the three protein structures (Chen et al., 2025). The most flexible region of the protein appears to be the middle domain with ClpV1 having the most flexible middle region compared to ClpV2 and ClpV3.

Finally, the three ClpV proteins were analyzed using the STRING server to examine their protein-protein interactions (Figure 9). This analysis aimed to explore the similarities and differences in the interacting partners of each ClpV protein. Interestingly, Figure 9 reveals that all three ClpV proteins interact with ClpP, which is unexpected given that ClpV lacks the tripeptide sequence required for ClpP interaction (Schlieker et al., 2005). The mechanism of interaction requires further investigation through in vitro studies or additional in silico protein-protein binding analyses. However, the observed interaction may result from indirect binding or alternative interaction motifs (Förster et al., 2014).

Our findings indicate that the three ClpV proteins are not structurally identical. Among them, ClpV1 and ClpV2 interact with ClpX, suggesting that the presence of a ClpX-binding site is another distinguishing feature of these proteins. As noted earlier, each strain contains a single copy of ClpX. However, since cells often produce multiple copies of the same protein, ClpX could potentially interact with two different types of ClpV proteins (Alberts et al., 2002). ClpX along with ClpA and ClpB are assumed to play a role in the functioning of the T6SS of Klebsiella pneumoniae, however the function of these Clp ATPases in the functioning of this system is yet to be elucidated (Barbosa and Lery, 2019).

MD simulations were performed for the three modelled ClpV proteins. Potential energy, radius of gyration, root mean square deviation (RMSD) and root mean square fluctuation (RMSF) were analyzed to explore the dynamic nature and stability of the proteins (Figure 11). The potential energy profiles of the three proteins remained stable throughout the simulation time of 250 ns (Figure 11A). The radius of gyration (rg) gives an indication of the compactness of the protein, the rg values of all three proteins did not exceed 1 Å therefore indicating that the proteins are compact and stable (Figure 11B). The RMSF shows the residues of the proteins which are the most flexible during the period of simulation (Mohan et al., 2022). The flexibility of all three proteins is significant (>1Å), thus indicating that the proteins are highly dynamic (Figure 11C). The RMSD further confirms the stability of the three proteins, ClpV2 and ClpV3 appear to reach equilibration around 10 ns, while ClpV1 seems to stabilize around 10 ns and continues fluctuating until 100 ns (Figure 11D). The MD run was extended to 400 ns and the dynamic nature of the protein remained consistent with the 250 ns run (data not shown). This along with the RMSF, Rg, and RMSD analysis at 250 ns confirmed that the proteins are stable, and compact. Radial distribution function (RDF) analysis also showed that the proteins had reached full conformational sampling as the systems equilibrated over time (data not shown).