Characterization of ribonucleotide reductases of emerging pathogens Elizabethkingia anophelis and Elizabethkingia meningoseptica and streptonigrin as their inhibitor: a computational study

Abstract Antibiotic resistance is a global concern. Two members of the bacterial genus Elizabethkingia, namely, E. anophelis and E. meningoseptica have raised much concern in recent years because of their resistance to multiple commonly used antibiotics. Identification of multidrug resistant and pan-drug resistant bacteria has propelled the search for new antibiotics that can act on unconventional targets. Researches are going on to find out the possibility of using bacterial ribonucleotide reductases as a novel target for antibiotic development. Through in silico evaluations, this study aims for characterization and functional annotation of ribonucleotide reductase enzymes of E. anophelis and E. meningoseptica. Binding affinities with these enzymes of the compounds that have shown promising results in inhibiting Pseudomonas aeruginosa growth by acting on its ribonucleotide reductase were also assessed by molecular docking and dynamics simulations. Insights from this study will help in battling these infections in the near future. Communicated by Ramaswamy H. Sarma


Introduction
Antibiotic resistance is a slowly progressing pandemic of which most people are unaware. Indiscriminatory uses of antibiotics, unnecessary uses of broad-spectrum antibiotics, reluctance of patients to adhere to prescribed doses of antibiotics, over-the-counter antibiotic availability in many countries, extensive uses of antibiotics in agriculture, poultry and livestock feeds -all these factors are hastening the development of antibiotic-resistant bacteria. According to the US Center for Disease Control (CDC), at least 2.8 million American people fall victim to antibiotic-resistant bacteria or fungi, and the death toll is more than 35,000 (https://www. cdc.gov/drugresistance/about.html). The world statistics are far scarier. According to the World Health Organization (WHO), each year at least 700,000 people die of antibioticresistant bacteria, of which 230,000 people die only from multidrug-resistant tuberculosis. The projected deaths due to antibiotic-resistant bacteria by the year 2050 are 10 million (https://www.who.int/news-room/detail/29-04-2019-new-reportcalls-for-urgent-action-to-avert-antimicrobial-resistance-crisis).
Relentless endeavors for finding new antibiotics to combat drug resistance are ongoing and the search for finding new targets for antimicrobial drug development is now more pronounced than ever before. Ribonucleotide reductase (RNR) is considered such a novel target (Tholander & Sjoberg, 2012). Ribonucleotide reductase catalyzes the conversion of four ribonucleotides (ADP, GDP, CDP and UDP) to DNA building blocks, the deoxyribonucleotides, an indispensable reaction for ensuring proper DNA synthesis. It is the only enzyme responsible for such conversions and is the rate-limiting enzyme for DNA synthesis and repair (Chapman, 2012). Ribonucleotide reductase enzymes act as protein freeradicals and generate a cysteine thiyl radical in their active sites. They are classified into three types (I, II and III) depending on the metallo-cofactor needed to generate the cysteine thiyl radical. Class Ia and Ib RNRs use a diferric (Fe III 2 ) or dimanganese (Mn III 2 ) coupled tyrosyl cofactor whereas class Ic and Id RNRs exploit directly an oxidized metal cluster (Mn IV Fe III or Mn III Mn IV ) as the cofactor. Another class I RNR, class Ie utilizes dihydroxyphenylalanine radical without using any metal. Class I RNRs can only continue their actions under aerobic conditions. Alpha subunits of class I RNRs are larger than the beta subunits and contain catalytic sites and binding sites for allosteric effectors. The beta subunits contain the metallocofactors or/and organic radicals. For Class II RNRs, adenosylcobalamin and for class III RNRs, S-adenosylmethionine and an iron-sulfur cluster act as cofactors. Class II and Class III RNRs act in oxygen-independent ways. Class I and Class III RNRs are found as heterotetramers (a 2 b 2 ) of two homodimers, a 2 and b 2 . Class II RNRs are mainly homodimers (a 2 ). All these three types of RNRs are found in bacteria and a bacterium may contain more than one type of RNR (Jordan & Reichard, 1998;Rose et al., 2019;Torrents et al., 2002).
Two members of the genus Elizabethkingia-Elizabethkingia anophelis and Elizabethkingia meningoseptica have recently emerged as threats because of their multidrug resistance. These bacteria are aerobic, non-motile, gramnegative rods. Immunologically vulnerable patients, like neonates, patients in the post-operative wards, and old people with comorbidities are mainly infected by them. These two bacteria have been found to cause meningitis, hospitalacquired pneumonia, sepsis, and bacteremia (Govindaswamy et al., 2018;Swain et al., 2015;Wang et al., 2019).
The current study aims to identify structural and functional characteristics as well as phylogenetic relationships of drug target RNR enzymes of these two bacteria using computational approaches. This study also endeavors to find out their possible inhibitors through molecular docking followed by molecular dynamics (MD) simulations. Tholander and Sjoberg (2012) showed that streptonigrin (NSC45383) which is a quinone-containing antitumor antibiotic, 7-Nitro-4-(1-oxidopyridin-1-ium-2-yl) sulfanyl-2,1,3-benzoxadiazole (NSC228155), and toluidine blue (NSC36758) are the most potent inhibitors of Pseudomonas aeruginosa class I ribonucleotide reductase among a group of 1,364 compounds. We selected streptonigrin and NSC228155 for the docking and dynamics simulations. Toluidine blue was excluded from the docking studies because it was proposed that toluidine blue doesn't inhibit RNR by pure competitive inhibition, rather by redox-related activities (Tholander & Sjoberg, 2012).

Analysis of primary structure
ExPASy's ProtParam Tool (https://web.ExPASy.org/protparam/) was employed to determine the physical and chemical properties of the proteins under study. The physicochemical parameters it computes are deduced from the amino acid sequence of the query protein, and include molecular weight, theoretical pI (isoelectric point), amino acid composition, atomic composition, extinction coefficient, estimated half-life, instability index, aliphatic index, and grand average of hydropathicity (GRAVY) (Gasteiger et al., 2005). Input to the ProtParam Tool can be a Swiss-Prot/ TrEMBL accession number or a sequence identifier, or an amino acid sequence in one-letter code. We used one-letter code amino acid sequences of the desired proteins as inputs.

Prediction of domains and motifs
Conserved domains in the proteins under study were searched against the CDD v3.18-55570 PSSMs database using NCBI Conserved Domain Search (https://www.ncbi.nlm. nih.gov/Structure/cdd/wrpsb.cgi) (Lu et al., 2020). Amino acid sequences of proteins in FASTA format were used as inputs. Notes: Extinction Coefficient 1 (assumed all pairs of Cysteine residues form Cystines), Extinction Coefficient 2 (assumed all Cysteine residues are reduced). Extinction coefficients are in units of M À1 cm À1 , at 280 nm measured in water.

Secondary structure prediction
Secondary structures of the alpha and beta subunits of E. anophelis and E. meningoseptica RNRs were predicted by PSIPRED 4.0 server (http://bioinf.cs.ucl.ac.uk/psipred/). PSIPRED (PSIblast-based secondary structure PREDiction) first runs PSI-BLAST searches against a custom-built non-redundant protein sequence data bank to obtain a position-specific scoring matrix. This matrix is next processed sequentially by two artificial neural networks having feed-forward back-propagation network architectures with double hidden layers to predict the final secondary structure (Buchan & Jones, 2019;Jones, 1999).

Template selection for homology modeling
HHPRED is a web tool (https://toolkit.tuebingen.mpg.de/ tools/hhpred) that uses hidden Markov models (HMMs) for protein homology detection and structure prediction. For homology detection, it searches against precomputed HMM profile databases (Zimmermann et al., 2018). In this study, protein sequences in FASTA format were provided as inputs to HHPRED. PDB_mmCIF70 database was chosen as the target database for identifying homologs with known structures of the desired proteins.

Homology modeling & model energy minimization
HHPRED outputs were sent directly to MODELLER, an efficient homology modeling tool ( Sali & Blundell,1993). The models generated by MODELLER were next energy minimized by YASARA energy minimization server (http://www. yasara.org/minimizationserver.htm). YASARA (Yet Another Scientific Artificial Reality Application) server uses its own YASARA force field for energy minimization which is an integration of AMBER (Assisted Model Building with Energy Refinement) all-atom force field equation, multi-dimensional knowledge-based torsion angle energy potentials, and a consistent set of force field parameters (Krieger et al., 2009). PyMOL software (The PyMOL Molecular Graphics System,   Open-source PyMOL, Schr€ odinger, LLC.) was used for visualization of the 3D models.
Protein subcellular localization and function prediction, and phylogenetic tree construction Subcellular localization of E. anophelis and E. meningoseptica RNRs were predicted by CELLO2GO server (http://cello.life. nctu.edu.tw/cello2go/) (Yu et al., 2014). Phylogenetic trees of the RNR proteins were generated to reveal their evolutionary history. First, BLAST searches were conducted against the non-redundant protein sequences (nr) database using NCBI Protein BLAST tool. Standard Protein BLAST (protein-protein BLAST algorithm) was used for retrieving the homologous sequences. During BLAST searches, model proteins (XM/XP), non-redundant RefSeq proteins (WP), and proteins from the same organism were excluded to avoid repetition. The top 50 sequences from the BLAST result were selected for generating multiple sequence alignment (MSA) by COBALT (Constraint-based Multiple Alignment Tool) (Papadopoulos & Agarwala, 2007) incorporated in the NCBI platform. Phylogenetic trees were later built from MSA applying 'fast minimum evolution' method, and the 'maximum allowed sequence difference' parameter was set to 0.85.
Ligand preparation for docking 3D structures of Streptonigrin and NSC228155 were downloaded in Structure Data File (SDF) format from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). They were converted to Protein Data Bank (PDB) format using PyMOL (The PyMOL Molecular Graphics System, Open-source PyMOL, Schr€ odinger, LLC.). The ligands were next optimized

Molecular docking and dynamics study
Molecular docking studies can give us an approximation of binding affinity between a ligand and its target. In our current study, molecular docking simulations were performed using AutoDock Vina (Trott & Olson, 2009  Information about the grid boxes used for docking can be found in Dassault Syst emes, 2020) was used for visualization of receptor-ligand interactions. It was found that streptonigrin has a better binding affinity to the RNR proteins than NSC228155. Next, MD simulations were conducted for 50 ns for the E. anophelis RNR-streptonigrin and E. meningoseptica RNR-streptonigrin complexes. The conditions applied during MD simulations were previously described by us (Mahfuz et al., 2020). Triplicate MD simulations were carried on for each complex.

Results and discussion
Sequence retrieval and analysis of primary structure Alpha subunit of E. anophelis RNR contains 550 amino acid residues and that of E. meningoseptica contains 551 residues. Beta subunits of both E. anophelis and E. meningoseptica RNRs have 324 residues. Aligning the sequences of RNR alpha and beta subunits of these bacteria revealed that the alpha subunits share 98% (541/551) and the beta subunits share 88% (284/324) identical amino acid residues. Knowledge of protein primary structure can help in predicting post-translational modifications and identifying antigenic sites responsible for specific antibody binding (Eisenhaber et al., 1995). ExPASy's ProtParam Tool predicted various physical and chemical properties of the proteins under study. Theoretical isoelectric points of the proteins indicate that these proteins are acidic in nature. Extinction coefficient is a measurement of the amount of light absorbed by a protein which is relative to its amino acid composition at a particular wavelength. ProtParam predicted extinction coefficients of the proteins at 280 nm wavelength. These values would help in the purification of these proteins. The instability indices indicate that beta subunits are more enduring than alpha subunits in the test tube environment. Aliphatic index of a protein is an indicator of the relative volume of its aliphatic side chains (alanine, valine, isoleucine, and leucine). Thermostability of globular proteins usually increases with the increase in the aliphatic index value. It was found that aliphatic index of the beta subunits are more than the alpha subunits. This also implies that the beta subunits are more stable. A protein's GRAVY (GRand AVerage of hYdropathicity) value indicates its hydrophobicity. The negative GRAVY values of E. anophelis and E. meningoseptica RNR proteins indicate that these proteins are hydrophilic. Table 1 summarizes ProtParam output.

Prediction of domains and motifs
Domains are important determinants of protein structure and functions. Domains identified from NCBI Conserved Domain Search are listed in Table 2. Only domains with specific hits have been considered. No hit was identified from the ScanProsite search for motifs.

Prediction of secondary and tertiary structures
PSIPRED predicted secondary structures of the RNR proteins under study. It was predicted that alpha subunits have mostly helical and coiled regions with few strands. Beta subunits were predicted to contain mostly helical with few coiled regions but no strand. Figure 1 depicts secondary structure predictions by PSIPRED.
There are currently no crystal structures available for the proteins under study. So, we decided to build their homology models. The suitable templates for homology modeling were chosen using HHPRED. HHPRED accepts a protein sequence or multiple sequence alignment in A3M, CLUSTAL, FASTA or STOCKHOLM format as input. The most suitable templates chosen by HHPRED were-C chain of Flavobacterium johnsoniae class Id RNR alpha subunit (PDB ID: 6DQW) for both E. anophelis and E. meningoseptica RNR alpha subunits, and A chain of Flavobacterium johnsoniae dimanganese (II) RNR beta subunit (PDB ID: 6CWO) for both E. anophelis and E. meningoseptica RNR beta subunits (supplementary material, Table 1). Homology modeling can provide satisfactory results if sufficient similarity is present between the query and the template sequence. MODELLER, a homology-based 3D structure prediction tool, predicted 3D models of the proteins under investigation using the abovementioned templates. An alignment of the sequence to be modeled with a homologous solved protein structure's sequence is used as the input to MODELLER. MODELLER derives a set of spatial restraints on the query sequence from the alignment and the template structure. These restraints are the probability density functions of a set of geometrical criteria for the location of every atom in the sequence. MODELLER then calculates the 3D model by satisfying the restraints with the least-possible violations ( Sali & Blundell,1993). Energy minimization of a generated model is required to ensure that chemical structures are energetically favorable and they achieve a proper molecular arrangement in space (Roy et al.,2015). Our models were energy minimized by applying YASARA all-atom forcefield. Figures 1-4  (supplementary material) show the cartoon presentations of MODELLER generated and YASARA server energy minimized models.

Verification and validation of the 3D models
The generated 3D models of E. anophelis and E. meningoseptica RNR alpha and beta subunits were assessed by PROCHECK, ERRAT and VERIFY3D web tools. PROCHECK examines residue-by-residue geometry and overall structural geometry to find out the stereochemical quality of a query protein structure. PROCHECK calculates Ramachandran plot which is a very reliable indicator of protein structure quality, whether it is obtained from X-ray crystallography or NMR or built by homology modeling. A protein structure is considered good if it has more than 90% of residues in the most favored regions on the main-chain Ramachandran plot (Laskowski et al., 2012). PROCHECK calculated Ramachandran plots of the generated 3D models where every model was found to have >90% of residues in the most favored regions, a reliable criterion fulfilled for being a good model. ERRAT can differentiate between correctly and incorrectly determined regions in a protein structure. It was developed based on the notion that atoms of different types have a nonrandom distribution in a good protein structure due to energetic and geometric factors, and random distributions are expected in incorrect regions. ERRAT analyzes the statistics of six non-covalently bonded atom-atom interactions involving C (carbon), N (nitrogen), and O (oxygen) atoms (CC, CN, CO, NN, NO, and OO). It then uses a quadratic error function to characterize the set of pairwise interactions against a database obtained from 96 high-confidence protein structures (Colovos & Yeates, 1993). ERRAT-calculated overall quality factor of each of the models was >90. VERIFY3D assesses the quality of a protein structure by comparing the 3D structure with its 3D profile. This 3D profile is derived from the protein's own amino acid sequence. VERIFY3D characterizes every amino acid in the query structure by its environment and checks which amino acid fits that environment in accordance with solved protein structures. The profile score is high for correct structures and tends to be low for faulty structures (L€ uthy et al., 1992). VERIFY3D considers a protein structure acceptable if it has at least 80% of residues scoring !0.2 in the 3D/1D profile. All of the models met this criterion. Table 3 and Figures 5-16 (supplementary material) depict results from the analyses done by PROCHECK, ERRAT and VERIFY3D. The results obtained from these three tools ascertain the good quality of the generated models.

Protein subcellular localization and function prediction
CELLO2GO server determined subcellular localization and provided functional annotations of the proteins under study. This server also provided functional gene ontology (GO) annotations of the genes encoding the proteins. It conducts BLAST search to identify already annotated sequences homologous to the query sequence for function prediction (Yu et al., 2014). Predictions were identical for both the RNRs but with different probabilities (not shown). These annotations give an idea of the complex biological activities handled by these proteins in a comprehensive way. Subcellular localization and molecular functions of the selected proteins, biological processes performed by them, and their cellular distributions are summarized in Table 4.

Phylogenetic tree construction and analysis
RNR is one of the most primitive proteins. Lateral gene transfers as well as novel mutations have played critical roles in the evolution of different RNR classes (Torrents et al., 2002). Phylogenetic trees of E. anophelis and E. meningoseptica RNR were constructed to elucidate their evolutionary relationships with RNRs of other bacteria. It was found from the constructed phylogenetic trees that E. anophelis RNR alpha subunit forms a clade with RNR alpha subunits of E. miricola and Elizabethkingia sp. YR191. But, E. anophelis RNR beta subunit doesn't form a clade. The closest neighbors (members of sister clades) of E. anophelis RNR beta Subunit are several members from the genus Chryseobacterium and Chlorobi bacterium. E. meningoseptica RNR alpha and beta subunits also do not form a clade with other bacteria. Members from sister clades of E. meningoseptica RNR alpha subunit include different strains of E. anophelis and E. miricola. But, the sister clades of E. meningoseptica RNR beta subunit are formed by different members of the genus Chryseobacterium and Chlorobi bacterium. Elizabethkingia is a recently created genus (2005) whose members were previously placed in the genus Chryseobacterium (Etymologia: Elizabethkingia. 2016). The phylogenetic analyses exhibit close ties of Elizabethkingia with Chryseobacterium and other members of the Flavobacteriaceae family. They also suggest alpha subunits of E. anophelis and E. meningoseptica RNR are more evolutionarily conserved than the beta subunits. The RNRs of E. anophelis and E. meningoseptica have not yet been assigned to a class. Recently it was identified that one member of Flavobacteriaceae, Flavobacterium johnsoniae has class Id RNR (Rose et al., 2019). It is highly possible that E. anophelis and E. meningoseptica RNRs also belong to class Id. Figures 2  and 3 show the generated phylogenetic trees.

Molecular docking study
Molecular docking studies were performed to assess the inhibitory effects of streptonigrin and NSC228155 on E. anophelis and E. meningoseptica ribonucleotide reductase enzymes. Homology models of both E. anophelis and E. meningoseptica RNRs were built upon the crystal structure of Flavobacterium johnsoniae RNR which is a class Id RNR (Rose et al., 2019). The catalytic Cys266 residue is conserved between F. johnsoniae (Rose et al., 2019) and E. anophelis. E. meningoseptica has Cys residues in 267 and 268 positions. We postulate that Cys266 in E. anophelis RNR and either one or both of the Cys residues (267 and 268 positions) in E. meningoseptica RNR act as the catalytic Cys responsible for producing cysteine thiyl radicals. Only alpha subunits of the proteins were chosen for docking. Grid boxes for docking simulations were prepared in such a way that the putative catalytic Cys residues remain within the designed grid boxes. AutoDock Vina calculated 9 binding poses for every docking simulation. For each receptor-ligand complex cluster, the most favorable conformation (that with the least binding energy and RMSD) was considered as the best binding pose. It was found that streptonigrin binds to both E. anophelis and E. meningoseptica RNR alpha subunits with more affinity (À8.1 kcal/mol and À7.8 kcal/mol, respectively) than NSC228155 (-6.3 kcal/mol and À6.1 kcal/mol, respectively). Figures 4 and 5 show streptonigrin bound to the receptors in the best-docked poses and Table 5 represents the ligand-receptor interactions.

Molecular Dynamics study
MD simulations were run in triplicate and the obtained results were found reproducible. Changes in the RMSD values of the Ca atoms of streptonigrin-bound E. anophelis and E. meningoseptica RNR a subunits during triplicate simulations can be visualized in Figure 6. It can be found from Figure 6a that the RMSD values of the Ca atoms of streptonigrin-bound E. anophelis RNR a subunit increases up to 8 ns and then decreases and becomes stable at 12 ns. The RMSD values remain 3.0 ± 0.5 Å for the rest of the simulation. The ligand, streptonigrin, becomes stable at about 9 ns and maintains its stability afterward. From Figure 6b, it can be observed that the streptonigrin-bound E. meningoseptica RNR a subunit becomes stable within the first 5 ns and the RMSD values then oscillate between 2.9 ± 0.5 Å up to the end of the simulation time. The ligand, streptonigrin, also becomes stable within the first 5 ns and its RMSD fluctuations remain within 1 Å afterward.
RMSF variations of the a subunits of RNR proteins and streptonigrin, as obtained from triplicate MD simulations, are shown in Figure 7a and b, respectively. Residues in the Nand C-terminals and loop regions show higher variations in RMSF value. Except for them, most residues of both the RNR  Figure 8a and b shows interactions during MD simulations between the residues of E. anophelis and E. meningoseptica RNR a subunit proteins and the ligand streptonigrin, respectively. E. anophelis RNR a subunit interacts mostly by hydrogen bonds, hydrophobic interactions, and water bridges with streptonigrin ( Figure 8(a)). The hydrogen bonds are generated mainly by polar amino acids such as Arg153, Asn499, Tyr527 and present over 80% of the simulation time. The hydrophobic interactions are shown mainly by Tyr528 and for a brief duration by Pro66, Arg153, Arg154, Leu250 and Tyr527. The most representative water bridges are formed by Arg153 and Tyr528 residues. Also, ionic interactions are present for more than half of the simulation period between Arg153 and the carboxy moiety of streptonigrin (Figure 9(a)).
Most of the contacts between E. meningoseptica RNR a subunit protein and the ligand are water bridges, hydrogen bonds, and ionic interactions (Figure 8(b)). Most important hydrogen bonds are generated by Ser65, Ser66, Ile107, Arg154 and Ser396 residues, remain present >50% of the simulation time. Arg154 is in the center of multiple interactions. It also has p-cation interactions with the C ring of streptonigrin and forms two water bridges (Figure 9(b)). These interactions are responsible for accommodation and stabilization of the ligand within E. meningoseptica RNR a subunit.

Conclusion
Bacterial ribonucleotide reductase is an attractive target for drug development. In this study, we have analyzed structural and functional properties of the ribonucleotide reductase enzymes of two important emerging pathogens. These properties would be valuable for purification and detailed in vitro and in vivo studies of these unexplored proteins. Results from the molecular docking and dynamics simulations indicate that streptonigrin can bind to E. anophelis and E. meningoseptica RNR alpha subunits with good affinities and can be considered a potential inhibitor of these RNRs.
This study further suggests streptonigrin is a superior inhibitor of the two selected bacterial RNRs to NSC228155. Since Pseudomonas aeruginosa and Elizabethkingia (anophelis and meningoseptica) are phylogenetically distant bacteria, this study also corroborates the importance of streptonigrin's being a lead compound for developing novel antibiotics.