The genetic information of selected strains within the Ureaplasma sp. (Supplementary Table S1) was gathered from the GenBank database at the National Center for Biotechnology Information (NCBI), which functions as a central repository for biomedical and genomics data. In order to maintain uniformity in genome annotations, the RAST server (Rapid Annotations using Subsystems Technology) (Overbeek et al., 2014) was utilized to annotate all of the genomes. The analysis of the core proteome among 13 genomes was conducted using reciprocal best BLAST hits of all coding sequences in the EDGAR version 3.0 software framework (Dieckmann et al., 2021). For this analysis, the genome of U. urealyticum serovar 10 str. ATCC 33699 was chosen as the reference genome, and the genomes of the remaining strains were compared to the reference.

The duplicated copies or redundant sequences from the main protein set were eliminated by implementing the CD-HIT module within the CD-HIT suite (Li and Godzik, 2006). To be considered redundant, the sequence identity had to exceed a specific cutoff of 0.6 (or 60%). For the subsequent analysis, priority was given to the protein sequences that were not duplicates. To identify the essential proteins of U. urealyticum, the core set of proteins was further assessed using GEPTOP version 2.0 (Naorem et al., 2022). An essentiality score cutoff of 0.24 was utilized. The GEPTOP platform was employed to detect essential genes in prokaryotic organisms. This was achieved by comparing the orthology and phylogeny of the query proteins with the experimentally defined datasets in the database of essential genes (DEG). After identifying non-duplicated protein sequences deemed essential, a BLASTp search against the Homo sapiens genome was conducted using the threshold E-value cutoff of 10⁻⁴. Any resulting sequences showing significant similarity to the human host were disregarded, while sequences with no homology were chosen for downstream analysis.

To identify metabolic pathways in the human host and U. urealyticum, we utilized the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Automatic Annotation Server (KAAS). Using the KEGG database, we retrieved unique identification numbers for the metabolic pathways in H. sapiens and U. urealyticum. We then conducted a manual comparison to distinguish exclusive from shared pathways. Pathways exclusive to the pathogen were labeled as unique metabolic pathways, while those found in both the human host and the pathogen were identified as common metabolic pathways. In this research, the protein sequences implicated in the pathogen’s unique metabolic pathways, as well as proteins annotated with KO numbers but not associated with neither specific nor common pathways, were screened for subsequent detailed analysis.

To prioritize proteins found in the cytoplasm for assessing potential therapeutic targets, we first utilized the PSORTb version 3.0.3 server (Yu et al., 2010) to predict the subcellular localization of these proteins. To enhance the credibility of the PSORTb version 3.0.3 findings, we subsequently validated them through the BUSCA (Savojardo et al., 2018) and CELLO version 2.5 (Yu et al., 2006) servers. The most accurate results from the localization analysis were determined by aggregating the data from these three servers. The final protein localization result was the same location predicted by two or more servers.

The ability of a target to bind with drugs and drug-like molecules at a high affinity determines its druggability. To evaluate the druggability of crucial cytoplasmic proteins, Blastp was used to search for them in the DrugBank database, accessible at https://go.drugbank.com/. DrugBank is a valuable resource for bioinformatics and cheminformatics, offering extensive information on drugs and their targets. Potential drug targets or drugs were identified by considering proteins with an E value below 10⁻⁵. In instances where no matches to known drugs or drug targets were found, these proteins were classified as novel drug targets of U. urealyticum (Khan M. T. et al., 2022). The DrugBank non-homologous proteins were then prioritized for druggability assessment based on druggability probability via the PockDrug server (Hussein et al., 2015) and their protein–protein interaction characteristics retrieved from STRING version 12.0 database (Szklarczyk et al., 2023).

In the field of host cell biology, proteins that act as anti-targets are of great significance due to their ability to bind with potential therapeutic compounds designed to combat corresponding pathogens. Among the human population, a comprehensive search of the existing literature discovered 210 such proteins (Fatoba et al., 2021). Noteworthy examples of these proteins include P-glycoprotein (referred to as P-gly), adrenergic receptor, dopaminergic receptor, and ether-a-go-go-related protein. In order to minimize any negative consequences resulting from the interaction between these anti-target proteins and the proposed drug targets, a thorough analysis was performed. The novel drug targets were subjected to a search using NCBI BLASTp against these 210 anti-target proteins, applying the following criteria: an E-value of less than 0.005, a query coverage greater than 30%, and an identity below 30%.

Accidental blockage or unintentional inhibition of proteins found within the gastrointestinal microflora may result in dysbiosis, significantly affecting the microenvironment and possibly causing toxicity with adverse effects on the human host (Sarker et al., 2023). Essential non-homologous proteins, identified as potential novel drug targets, underwent screening via BLASTp, applying an E-value threshold of less than 0.005, a query coverage exceeding 30%, and an identity percentage below 30%. This analysis utilized the dataset available on the Human Microbiome Project server¹ under “43,021 (BioProject)” (Khan K. et al., 2022; Peterson et al., 2009). Proteins exhibiting less than 30% similarity were classified as novel therapeutic targets and proceeded to the subsequent phase.

The target proteins analyzed in this study underwent structure prediction using the SWISS-MODEL program (Waterhouse et al., 2018) accessed through the Expasy web server, known for its comprehensive homology modeling service. The accurate evaluation of the 3D model is of paramount importance in the field of computational structure prediction (Senior et al., 2020). Recent advancements in sequencing techniques have led to groundbreaking discoveries in computational structural biology (Mendez et al., 2023; Tripathi et al., 2019). The introduction of widely accepted and efficient techniques for structure evaluation has enabled the qualitative estimation of protein structures. In this research, three distinct tools, namely PROCHECK (Laskowski et al., 1996), ERRAT (Colovos and Yeates, 1993) and ProsA-web (Wiederstein and Sippl, 2007), were utilized to assess the stereochemical quality and accuracy of the predicted models (Hooda et al., 2012).

To identify a suitable target for drug development, it is essential to investigate the interactions between the drug target and potential drug or drug-like compounds. Consequently, it is necessary to search for suitable ligand molecules that are capable of binding to the target proteins. In preparing the ligands library, we selected the FDA-approved drug library, which contains 2,470 drugs, along with an experimental drug library comprising 6,041 drugs (Gan et al., 2023). All drug molecules were obtained from DrugBank in the structured data file (SDF) format. The structures of these drug molecules were subjected to energy minimization and subsequently compiled into a ligands library designed for virtual screening.

The recently developed DrugRep server² facilitated the execution of virtual high-throughput screening as an online tool for computer-aided drug discovery. In the process of virtual screening with DrugRep, the structures of well-prepared drug target proteins underwent receptor-based virtual screening against a library of ligands that included 8,511 distinct drug molecules. For the docking procedure, the DrugRep server employed Autodock Vina (Trott and Olson, 2010) to identify potential inhibitors against the drug target proteins and returned 100 best molecules based on estimated binding affinities (kcal/mol).

To elucidate the interaction dynamics between the receptors (shortlisted drug target proteins) and the ligands (the top-ranked drug molecules), the tools AutoDock Vina and PoseEdit were employed. The PoseEdit tool was sourced from the ProteinsPlus Server.³ An examination was performed on both receptors alongside their respective ligands to evaluate interactions, binding conformations, and affinities. Furthermore, the Top-10 screened drugs underwent filtering based on Lipinski’s rule of five (RO5). Ultimately, the ligands that demonstrated the most advantageous binding affinities were selected according to their optimal docking scores and adherence to RO5.

All-atom MD simulations were conducted using small molecule-protein complexes derived from docking as the initial structures, employing Gromacs 2023.3 software (Pronk et al., 2013; Rodríguez-Martínez et al., 2024). Both the small molecule and the protein were represented using the AMBER protein force field (Wang et al., 2004; Maier et al., 2015). The pdb2gmx tool was utilized to incorporate hydrogen atoms into the system, followed by the addition of a truncated cubic TIP3P water box at a distance of 10 Å from the system. Sodium (Na+) and chloride (Cl−) ions were included to neutralize the system’s charge (Bejagam et al., 2018). Subsequently, the topological and parameter files for the simulations were generated. MD simulations were executed with Gromacs 2023.3 software for a total simulation duration of 100 ns (Case et al., 2005). Prior to the simulations, energy minimization was performed using the mdrun command and the steepest descent method within the NVT ensemble, with an initial step size of 0.01 nm and a maximum force tolerance of 1,000 kJ/mol•nm. Following energy minimization, a 100 ps NVT (isothermal-isochoric) ensemble simulation was conducted at constant volume, with gradual heating from 0 to 310.15 K to facilitate the even distribution of solvent molecules within the solvent box. This was succeeded by a 100 ps NPT (isothermal-isobaric) ensemble simulation, employing the Berendsen barostat to equilibrate the system pressure at 1 bar. During the MD simulations, all hydrogen bonds were constrained using the LINCS algorithm, with an integration step size of 2 fs. Long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method, with a cutoff of 1.2 nm. The cutoff for nonbonded interactions was set to 10 Å, and these interactions were updated every 10 steps. The resulting trajectories underwent periodic boundary condition removal, after which root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), and radius of gyration (Rg) analyses were performed.

In our study, a total of 396 coding DNA sequences (CDs) were discovered to be present in all 13 strains. These CDs were converted into protein sequences and defined as the core proteome (Supplementary Sequence file 1). In order to verify the uniqueness of the identified protein sequences, we employed the CD-HIT web server with a 60% identity cutoff. Our analysis conclusively showed that all 396 protein sequences were distinct from any paralogous counterparts.

The 396 unique protein sequences underwent analysis using the GEPTOP server, identifying 170 sequences as essential proteins (Supplementary Sequence file 2). These vital proteins, crucial for the pathogen’s survival within the host, serve as potential targets for drug development. However, these proteins must not bear resemblance to human proteins to avoid potential drug-related complications. To address this requirement, the 170 essential proteins were compared against human (H. sapiens) proteins using BLASTp, resulting in 94 sequences identified as essential and non-homologous proteins (Supplementary Sequence file 3).

The KEGG database currently catalogs 44 metabolic pathways for U. urealyticum and 358 for human. By manually comparing these metabolic pathways between the pathogenic U. urealyticum and its human host, we identified 9 unique pathways specific to the pathogen and 35 pathways that are common to both the pathogen and the host, as detailed in Table 1. Furthermore, an analysis conducted using the KAAS server demonstrated that out of 94 proteins screened from the previous step, 92 had assigned KEGG orthology (KO) identifiers. Notably, among these 92 proteins, four were only linked with unique metabolic pathways exclusive to the pathogen, as illustrated in Table 2, and 21 were involved in neither unique nor common metabolic pathways (Supplementary Table S2). All of these 25 proteins were selected for further investigation. The remaining 67 proteins were excluded from further analysis because they were found to participate in the common metabolic pathways.

The analysis conducted with the PSORTb, BUSCA, and CELLO servers identified that among the 25 proteins screened, a significant majority, specifically 24 proteins, were categorized as cytoplasmic proteins, as detailed in Supplementary Table S3. This classification indicates that these cytoplasmic proteins might play pivotal roles within cellular processes, enhancing their relevance in drug development. Consequently, given their potential as effective drug target candidates, these proteins were selected for subsequent more in-depth investigation.

A total of five proteins exhibited similarity to the available drug targets documented in the DrugBank database (Table 3). Conversely, the other 19 proteins displayed no matches with any entries currently identified in DrugBank. This observation indicates their potential as novel drug targets for U. urealyticum. Among those 19 proteins, 18 possess suitable binding pockets that have a druggability probability exceeding 0.5 (Table 4). Additionally, protein–protein interaction (PPI) analysis recognized 15 of these proteins as hub proteins, fulfilling the requisite threshold of a node degree of at least 5 (Aslam et al., 2021), as indicated in the STRING version 12.0 database (Table 4). Following the elimination of the five proteins that did not meet the specified parameters, each of the remaining 14 hub proteins was found to interact with several proteins (Figure 2). These hub proteins are integral to a variety of functions, and the inhibition of their activities may affect the functions of other interacting proteins. Moreover, these 14 unique therapeutic proteins lacked similarity to 210 identified human “anti-targets.” Furthermore, a BLASTp analysis conducted on all microbial strains from the Human Microbiome Project (HMP) utilizing the NCBI blast server demonstrated that only two out of the 14 proteins showed a similarity of less than 30%. Given their unique properties and the absence of overlaps with common host-pathogen pathways and human “anti-targets,” it is advisable to assess these two proteins (B5ZC96 and B5ZAH8) as promising candidates for drug targeting.

Due to the unavailability of the desired protein 3D structures in the Protein Data Bank (PDB), the homology modeling for each protein was conducted by utilizing the FASTA sequences sourced from the UniProt database, specifically for the proteins with accession numbers B5ZC96 and B5ZAH8. The AlphaFold DB model of C0AGN1 was selected as the template for the homology modeling of B5ZC96 (Figure 3A) based on optimal GMQE values and sequence similarities, while the AlphaFold DB model of C0AGQ9 served as the template for modeling B5ZAH8 (Figure 3B). With these chosen template proteins, the 3D structures of both proteins were successfully modeled.

The validation of the developed models was performed using the PROCHECK, ERRAT, and ProSA-Web tools. For the protein B5ZC96, the Ramachandran plot revealed that approximately 88.1% of the residues are situated in the most favored regions, with 11.1% located in the additionally allowed regions, and both 0.4% of the residues found in the generously allowed and disallowed regions, as illustrated in Figure 3C. Conversely, concerning B5ZAH8, the Ramachandran plot indicated that nearly 92.9% of the residues resided in favorable areas, 7.1% in additionally allowed regions, and neither the generously allowed nor disallowed regions contained any residues (Figure 3D). Furthermore, the ERRAT plot provided an overall quality factor of 87.2038 for B5ZC96 and a perfect score of 100 for B5ZAH8, suggesting that the proposed models exhibit excellent quality. Additionally, the ProSA tool’s cross-validation of the modeled structures yielded a Z-score of −6.11 for B5ZC96 and −8.36 for B5ZAH8, indicating that these models are consistent with structures derived from NMR/X-ray crystallography.

The DrugRep server conducted a virtual high-throughput screening involving a set of 8,511 unique drug molecules. Our docking analysis indicated that the docking scores for the top 10 candidates ranged from −8.7 to −7.8 kcal/mol at the B5ZC96 active site, implying a significant affinity for this target protein (Table 5). Regarding RO5, the compounds 2-[3-({Methyl[1-(2-Naphthoyl)Piperidin-4-Yl]Amino}Carbonyl)-2-Naphthyl]-1-(1-Naphthyl)-2-Oxoethylphosphonic Acid, Irinotecan, and 3-(1 h-Indol-3-Yl)-2-[4-(4-Phenyl-Piperidin-1-Yl)-Benzenesulfonylamino]-Propionic Acid showed violations associated with molecular weight (MW), exceeding the ideal upper limit of 500 g/mol. Additionally, both 2-[3-({Methyl[1-(2-Naphthoyl)Piperidin-4-Yl]Amino}Carbonyl)-2-Naphthyl]-1-(1-Naphthyl)-2-Oxoethylphosphonic Acid and Pimozide went beyond the LogP constraint, exceeding the suggested threshold of 5. Furthermore, 2-[3-({Methyl[1-(2-Naphthoyl)Piperidin-4-Yl]Amino}Carbonyl)-2-Naphthyl]-1-(1-Naphthyl)-2-Oxoethylphosphonic Acid also violated the criteria for rotatable bonds, surpassing the preferred number of less than 10. Nevertheless, since Irinotecan and Pimozide have previously received FDA approval, it is unlikely that these two inhibitors will adversely impact pharmacokinetic properties. The remaining drugs complied closely with the RO5 criteria.

In addition, the binding energy for the leading 10 drugs concerning the B5ZAH8 receptor was significantly higher than that for the top 10 drugs targeting B5ZC96 (Table 6). The energy values for the top 10 drugs were established to be between −9.8 and − 9.0 kcal/mol. All top 10 drugs fulfilled the required RO5 parameters, with the exceptions of Zk-806450, 3-(2-aminoquinazolin-6-yl)-1-(3,3-dimethylindolin-6-yl)-4-methylpyridin-2(1H)-one and Arotinoid acid, each of which demonstrated a LogP violation exceeding 5, and Fluazuron, which violates both the LogP threshold greater than five and the MW limit of 500 g/mol (Table 6).

Following the docking and sorting processes, the Top-10 molecules evaluated for each drug target demonstrated the lowest docking scores. According to RO5 filtering criteria, the compound 1-(5-{2-[(1-methyl-1H-pyrazolo[4,3-d]pyrimidin-7-yl)amino]ethyl}-1,3-thiazol-2-yl)-3-[3-(trifluoromethyl) phenyl]urea was identified as the best ligand for B5ZC96, while ZK-806711 was recognized as the optimal ligand for B5ZAH8. The two-dimensional structure and binding conformations of the best ligand for each target protein were illustrated in Figure 4. These binding configurations were generated utilizing the DrugRep server. Following the completion of the docking process, further analysis revealed that the compound 1-(5-{2-[(1-methyl-1H-pyrazolo[4,3-d]pyrimidin-7-yl)amino]ethyl}-1,3-thiazol-2-yl)-3-[3-(trifluoromethyl)phenyl]urea mediated two pi bonds with Phe111 and Tyr114 (Figure 5A), whereas ZK-806711 mediated three pi bonds with Phe55 and single hydrogen bond with Ser111(Figure 5B).

The receptor-ligand complexes underwent 100 ns MD simulations using Gromacs. RMSD was employed to assess stability, with lower values indicating greater stability. For B5ZC96, the RMSD curve initially increased as the small molecule adapted within the binding cavity, stabilizing around 10 ns at 0.2 nm. Oscillations were observed around 50 ns but stabilized by 75 ns. Throughout the simulation, the small molecule maintained stable interactions, with minor adjustments reflected in the RMSD fluctuations (Figure 6A). In the case of B5ZAH8, the RMSD curve raised to 0.5 nm before stabilizing around 0.45 nm, subsequently dropping to 0.3 nm near the end of the simulation. The small molecule exhibited no significant displacement relative to the protein, maintaining stable interactions. Minor adjustments and the rotation of an alkyl chain contributed to the fluctuations observed in the RMSD curve (Figure 6B).

RMSF, assessed by the average deviation for each residue, was analyzed to evaluate receptor-ligand binding and protein dynamics. For B5ZC96, overall RMSF values were low, with increased values observed at the N- and C-termini due to their less constrained positions. Elevated RMSF in other regions may result from small molecule interactions or intrinsic flexibility (Figure 6C). Similarly, for B5ZAH8, RMSF values also remained low, with higher values at the N- and C-termini for analogous reasons. Increased RMSF in other regions might stem from small molecule interactions or inherent flexibility (Figure 6D). Low RMSF values suggest minimal fluctuations, indicating stable vibrations in the solvent environment and consistent sampling and analysis throughout the simulation.

Rg was a metric that reflected the compactness of protein structures and can also be utilized to assess changes in the looseness of the protein’s polypeptide chain during simulations. For B5ZC96, the analysis of the Rg curve indicated that during the initial 60 ns of the simulation, Rg exhibited fluctuations, gradually decreasing from approximately 2.20 nm to around 2.15 nm. In the subsequent 40 ns, Rg stabilized, exhibiting minimal fluctuations around 2.15 nm until the end of the simulation (Figure 6E). In contrast, for B5ZAH8, the examination of the Rg curve demonstrated that the protein’s Rg values oscillated within a narrow range of 1.72 nm to 1.76 nm throughout the entire duration of the simulation (Figure 6F).