Molecular Evolutionary Analyses of the Pseudomonas-Derived Cephalosporinase Gene

Despite the increasing evidence of the clinical impact of Pseudomonas-derived cephalosporinase (PDC) sequence polymorphisms, the molecular evolution of its encoding gene, blaPDC, remains elusive. To elucidate this, we performed a comprehensive evolutionary analysis of blaPDC. A Bayesian Markov Chain Monte Carlo phylogenetic tree revealed that a common ancestor of blaPDC diverged approximately 4660 years ago, leading to the formation of eight clonal variants (clusters A–H). The phylogenetic distances within clusters A to G were short, whereas those within cluster H were relatively long. Two positive selection sites and many negative selection sites were estimated. Two PDC active sites overlapped with negative selection sites. In docking simulation models based on samples selected from clusters A and H, piperacillin was bound to the serine and the threonine residues of the PDC active sites, with the same binding mode for both models. These results suggest that, in P. aeruginosa, blaPDC is highly conserved, and PDC exhibits similar antibiotic resistance functionality regardless of its genotype.

PDC contains more than 500 variants according to the Beta-Lactamase DataBase (http://bldb.eu/BLDB.php?prot=C#PDC, accessed on 26 December 2022) that can lead to its structural modifications [12,13]. Moreover, the PDC sequence polymorphism can alter resistance to antibiotics, thereby resulting in a major obstacle in P. aeruginosa treatment in clinical settings [12,[14][15][16]. However, the evolutionary history of bla PDC and the molecular interactions between PDC and the relevant antibiotics for each bla PDC genotype remain elusive.
Recent advances in in silico techniques and computing have resulted in impressive progress in molecular-evolutionary and docking-simulation analyses. These analyses provide opportunities for deciphering the evolution of PDCs and their molecular mechanisms of antibiotic resistance. We, therefore, analyzed the molecular evolution of bla PDC collected from different parts of the world. For the bla PDC genotypes selected from specific clonal variants (clusters), we clarified its molecular interactions with piperacillin, a major antibiotic against P. aeruginosa, via a docking simulation between PDC and piperacillin.

Sample Selection and Alignment
To investigate the molecular evolution of bla PDC , we obtained all available full-length nucleotide sequences for it from GenBank (https://www.ncbi.nlm.nih.gov/genbank/, accessed on 19 May 2021). Samples with missing data regarding the years or regions of isolation or detection or exhibiting ambiguous sequences or genetic recombination were omitted, leaving 511 samples. MEGA7 [17] was used to generate amino acid sequences for these samples, and multiple alignment was conducted using MAFFT 7 [18]. The genome data with homologous sequences (100% identity) were excluded using Clustal Omega [19], and the genome that was first isolated was also left. In total, 215 samples, isolated or detected in 40 countries between 1963 and 2021, were used. The present genome sample information is shown in Supplemental Table S1.

Time-Scaled Phylogenetic Tree Analysis
A time-scaled phylogenetic tree was constructed using the Bayesian Markov Chain Monte Carlo (MCMC) method in BEAST 2.4.8 [20], with a chain length of 40,000,000 steps and sampling every 1000 steps. First, using jModelTest 2.1.10 [21], we selected the best substitution model (TPM3uf + I + G). Subsequently, we applied the path-sampling/steppingstone sampling method [22,23] to search for the optimal combination of four clock models (Strict Clock, Relaxed Clock Exponential, Relaxed Clock Log Normal, and Random Local Clock) and two prior tree models (Coalescent Constant Population and Coalescent Exponential Population). For the analyses, we selected the Strict Clock and Coalescent Exponential Populations models. We then used Tracer 1.7.1 [24] to assess the effective sample size with respect to convergence, accepting values >200. After discarding the first 15% of iterations as burn-in, phylogenetic trees were generated using Tree Annotator 2.4.8 implemented in BEAST. Subsequently, the trees were illustrated using FigTree 1.40 (http://tree.bio.ed.ac.uk/software/figtree/, accessed on 12 October 2021). Moreover, we selected the representative bla PDC samples from each cluster determined by MCMC tree (GenBank: Cluster A, NZ_JAHBSB010000003; cluster B, NZ_RRBX01000028; cluster C, NZ_PHTA01000011; cluster D, NZ_CP053917; cluster E, NG_055276; cluster F, NZ_JAGFEQ010000011; cluster G, NZ_JAGMWR010000017; cluster H, CP020560). Then, their PDC sequences were compared. The evolutionary rate of the selected bla PDC samples was calculated using Tracer 1.7.1 by selecting a suitable model, as described.

Phylogenetic Distance Analysis
To calculate the phylogenetic distances among the samples, a phylogenetic tree was created using the maximum likelihood (ML) method in MEGA7. First, the best substitution model was selected using jModelTest 2.1.10. The phylogenetic distances of the ML tree, and for each cluster, were estimated using Patristic [25].

Selective Pressure Analyses
To reveal the positive and negative selection sites of PDC amino acid sequences, non-synonymous (dN) and synonymous (dS) substitution rates were estimated using the Datamonkey web server (http://www.datamonkey.org/, accessed on 6 December 2021). We used single-likelihood ancestor counting (SLAC), fixed-effects likelihood (FEL), internal fixed-effects likelihood (IFEL), fast unconstrained Bayesian approximation (FUBAR), and mixed-effects model of evolution (MEME) to estimate positive selection sites, and SLAC, FEL, IFEL, and FUBAR to estimate negative selection sites. The determination of positive (dN/dS > 1) and negative (dN/dS < 1) selection was based on p < 0.05 for SLAC, FEL, IFEL, and MEME, and posterior probabilities > 0.9 for FUBAR. These selection sites were mapped onto the structural models of PDCs.

Docking Simulation
To elucidate the molecular interactions between the structural models of PDC and piperacillin, docking simulations were performed using AutoDock Vina 1.1.2 [33]. Prior to docking simulation, polar hydrogen atoms and Gasteiger charges were added to PDC models using AutoDockTools 1.5.6. The grid box size was set to include all proteins. After the 20 docking models were generated using AutoDock Vina, 3D molecular interactions and conformations were visualized using PyMOL 2.3.4. Docking models with root mean square deviation ≥2, relative to the values obtained prior to docking simulation, were omitted. The optimal models in each cluster were determined from the docking models using AutoDock Vina based on the binding energy. The molecular interactions were analyzed in a 2D diagram using the BIOVIA Discovery Studio Visualiser.

Time-Scaled Phylogenetic and Evolutionary Analysis
The MCMC time-scaled phylogenetic tree of bla PDC (Figure 1 A 173 (128-268). The bla PDC belonging to cluster A were predominant at the time of analysis. Moreover, the active sites (motifs 1 to 3) of representative PDC sequences from each cluster were found to be conserved (Supplementary Figure S1). The evolutionary rate (in substitutions/site/year) of the entire bla PDC gene was 4.08 × 10 −5 (95% HPD, 2.52 × 10 −5 to 5.68 × 10 −5 ).

Time-Scaled Phylogenetic and Evolutionary Analysis
The phylogenetic distance among all of the blaPDC, estimated using an ML-based phylogenetic tree (Figure 2a), was 0.013 ± 0.013 (mean ± SD). For clusters A to G (Figure 2b), the mean distances were short, ranging from 0.005 to 0.086, with cluster A having the shortest phylogenetic distance (0.005 ± 0.000) and H having a relatively long distance (0.44 ± 0.046).

Time-Scaled Phylogenetic and Evolutionary Analysis
The phylogenetic distance among all of the bla PDC , estimated using an ML-based phylogenetic tree (Figure 2a), was 0.013 ± 0.013 (mean ± SD). For clusters A to G (Figure 2b), the mean distances were short, ranging from 0.005 to 0.086, with cluster A having the shortest phylogenetic distance (0.005 ± 0.000) and H having a relatively long distance (0.44 ± 0.046).

Selective Pressure Analysis
We estimated the positive and negative selection sites for the present PDC using the Datamonkey web server. Two amino acid substitutions (amino acids 79 and 239), common to all five methods, were identified as positive selection sites. The number of negative selection sites was 18, 54, 48, and 56 for SLAC, FEL, IFEL, and FUBAR, respectively. Of these, 18 sites were common to all four methods. No positive or negative selection sites were located in motif 1. In contrast, negative selection sites Tyr149, Thr315, and Gly316 were located in motif 2 and motif 3, respectively (Figure 3).

Molecular Interactions between PDC and Piperacillin
Docking simulation between the PDC structural models and piperacillin to clarify their molecular interactions was performed using AutoDock Vina. Based on the phylogenetic analysis results, PDC structural models were generated from the representative samples selected from clusters A and H.

Selective Pressure Analysis
We estimated the positive and negative selection sites for the present PDC using the Datamonkey web server. Two amino acid substitutions (amino acids 79 and 239), common to all five methods, were identified as positive selection sites. The number of negative selection sites was 18, 54, 48, and 56 for SLAC, FEL, IFEL, and FUBAR, respectively. Of these, 18 sites were common to all four methods. No positive or negative selection sites were located in motif 1. In contrast, negative selection sites Tyr149, Thr315, and Gly316 were located in motif 2 and motif 3, respectively (Figure 3).

Molecular Interactions between PDC and Piperacillin
Docking simulation between the PDC structural models and piperacillin to clarify their molecular interactions was performed using AutoDock Vina. Based on the phylogenetic analysis results, PDC structural models were generated from the representative samples selected from clusters A and H.
In the PDC structural models of cluster A, piperacillin interacted with Ser62 (motif 1) and Thr315 (motif 3) at the active sites via conventional hydrogen bonds (Figure 4a). An unfavorable interaction (donor-donor) was formed with Asn345 in PDC. The forces that interacted with sites other than the active sites were carbon-hydrogen bonds and hydrophobic interactions (π-π stacked and alkyl). The interacting residues other than those in the active sites were alanine, asparagine, serine, and tyrosine. The computational binding affinity was estimated to be −9.0 kcal/mol.
In the PDC structural models of cluster H, piperacillin formed conventional hydrogen bonds with Ser62 (motif 1) and Thr315 (motif 3) at the active sites (Figure 4b). There was an unfavorable interaction (donor-donor) between Arg348 in PDC and piperacillin. Carbon-hydrogen bonds and hydrophobic interactions (π-π stacked and alkyl) were formed as attractive forces at sites other than the active sites. The interacting residues other than those in the active sites were alanine, arginine, asparagine, serine, and tyrosine. The binding affinity was computed as −9.0 kcal/mol. In the PDC structural models of cluster A, piperacillin interacted with Ser62 (motif 1) and Thr315 (motif 3) at the active sites via conventional hydrogen bonds (Figure 4a). An unfavorable interaction (donor-donor) was formed with Asn345 in PDC. The forces that interacted with sites other than the active sites were carbon-hydrogen bonds and hydrophobic interactions (π-π stacked and alkyl). The interacting residues other than those in the active sites were alanine, asparagine, serine, and tyrosine. The computational binding affinity was estimated to be −9.0 kcal/mol.

Discussion
Sequence polymorphism of PDC has generated increasing concern [34,35]. To elucidate the molecular evolution of blaPDC, we investigated this using blaPDC gene samples from various parts of the world. Based on the blaPDC gene clusters revealed by this analysis, we conducted a docking simulation between PDC and piperacillin to elucidate their molecular interactions. This revealed that PDC evolved constantly and formed eight clusters, whereas the phylogenetic distances of blaPDC were short. These results, suggesting low levels of genetic divergence among the current blaPDC, are consistent with a prior study of clinically isolated samples, in which the sequences were relatively conserved despite variations in gene size [36]. At the same time, the number of PDC with altered resistance to In the PDC structural models of cluster H, piperacillin formed conventional hydrogen bonds with Ser62 (motif 1) and Thr315 (motif 3) at the active sites (Figure 4b). There was an unfavorable interaction (donor-donor) between Arg348 in PDC and piperacillin. Carbonhydrogen bonds and hydrophobic interactions (π-π stacked and alkyl) were formed as attractive forces at sites other than the active sites. The interacting residues other than those in the active sites were alanine, arginine, asparagine, serine, and tyrosine. The binding affinity was computed as −9.0 kcal/mol.

Discussion
Sequence polymorphism of PDC has generated increasing concern [34,35]. To elucidate the molecular evolution of bla PDC , we investigated this using bla PDC gene samples from various parts of the world. Based on the bla PDC gene clusters revealed by this analysis, we conducted a docking simulation between PDC and piperacillin to elucidate their molecular interactions. This revealed that PDC evolved constantly and formed eight clusters, whereas the phylogenetic distances of bla PDC were short. These results, suggesting low levels of genetic divergence among the current bla PDC , are consistent with a prior study of clinically isolated samples, in which the sequences were relatively conserved despite variations in gene size [36]. At the same time, the number of PDC with altered resistance to antibiotics is increasing [34,35,37]. These findings indicate that even minor amino acid substitutions may affect the functions of PDC.
Our selection-pressure analysis of PDC revealed two positive selection sites and many negative selection sites. Moreover, negative selection sites overlapped with motifs 2 (Tyr149) and 3 (Thr315 and Gly316) in PDC active sites. Negative selection acts as a sieve, which leads to the elimination of harmful mutations [38][39][40]. This suggests that P. aeruginosa with mutations in PDC active sites do not adapt to their environment, leading to the presence of highly conserved active sites. Indeed, the active sites of representative PDC sequences from each cluster were found to be conserved.
Molecular evolutionary analysis of P. aeruginosa gyrA revealed that key residue substitutions in GyrA, with positive selection due to quinolone drugs, altered the binding mode between GyrA and ciprofloxacin [41]. Although PDC is an antibiotic resistance-conferring enzyme rather than a target, its molecular mechanisms of action may evolve under selective pressure. In our PDC-piperacillin docking simulation, piperacillin interacted with serine and threonine residues at the active sites for clusters A and H. These results suggest that PDC shows similar antibiotic resistance mechanisms regardless of its genotype. This is partially due to the conservation of active sites, which play an essential role in β-lactam antibiotic hydrolysis, among the clusters [42]. Furthermore, deeper insight into the molecular mechanisms of PDC in evolutionarily selected variants will pave the way towards developing promising β-lactam/β-lactamase inhibitors against P. aeruginosa.
This study has computing capability-related limitations that need to be addressed in the future. The first is the reliability of structural models. The homology modeling generates a good stereochemical protein model but is limited in predicting the native structure accurately from its sequence and template structure alone [43]. Quite recently, a new powerful tool of structural modeling using deep learning techniques, Alphafold2, has become available [44]. AlphaFold2 may exhibit a smaller structural gap between native and target proteins than homology modeling [44]. However, we could not perform AlphaFold2 because it requires high computing capability. The second limitation concerns the lack of a dynamic process through which PDC recognizes piperacillin in our simulation. Subtle changes in proteins may cause drug resistance by maintaining substrate recognition and catalysis while preventing the antibiotic from exerting its full effect [15,45]. In PDC, various amino acid substitutions at sites other than the active sites, such as those in the short peptide strand traversing the active sites (Ω loop) and the C-terminal domain, can affect resistance to β-lactam antibiotics [29,35]. We were unable to include these catalytic residues in our docking simulation because our method applies a "snap-shot" approach that does not present the molecular mechanisms of target molecules as a dynamic image. Although molecular dynamic simulations using extremely high-performance computer systems provide opportunities to explore dynamic molecular mechanisms [46][47][48], relatively few laboratories have the facilities to conduct them. Although the methodology of the present study is subject to these limitations, our molecular evolution and docking simulation analyses will guide the designing of effective therapies targeting PDC.

Conclusions
To the best of our knowledge, this is the first comprehensive evolutionary analysis of PDC gene. These findings will improve our understanding of time-scaled genetic and functional changes in PDC. Based on the MCMC phylogenetic tree, a common ancestor of bla PDC diverged ca. 4660 years ago, giving rise to eight clusters. The phylogenetic distances within the samples were short. Our analysis revealed two positive selection sites, and many negative selection sites, in PDC. These were not located in the active sites motif 1 but in motifs 2 (Tyr149) and 3 (Thr315 and Gly316). Docking simulations between piperacillin and the PDC structural models based on clusters A and H revealed similar binding modes. These findings imply that bla PDC is highly conserved, particularly at the active sites, resulting in similar antibiotic-related functions for PDC, regardless of the genotype.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest in association with the present study.