Bacterial genomes and gene reannotation

Whole genome sequences for S. pneumoniae TIGR4 uid57857, S. pneumoniae R6 uid57859, H. influenzae Rd KW20 uid57771, H. influenzae 86 028NP uid58093 and M. catarrhalis BBH18 uid48809 were obtained from the National Centre for Biotechnology Information (NCBI) Genbank File Transfer Protocol (FTP) website (ftp://ftp.ncbi.nih.gov/genbank/). All open reading frame (ORF) annotations were updated using Rapid Annotation using Subsystem Technology (RAST) [16]. In this analysis, all locus coordinates in original Genbank genomes release were retained without adjustments for frame-shifts.

Orthology and gene essentiality predictions

We clustered the reannotated protein sequences into putative orthologous groups using the OrthoMCL standalone software Version 2.0.2 [17]. Most studies have consistently deciphered essential genes under ideal conditions, that is, in the richness of all necessary nutrients and without environmental stress. For the purpose of this study, we define the “essentiality” of a gene as its indispensability under rich media conditions, unless stated otherwise. The caveat with this approach is that essential genes required for metabolism within the host may be missed.

Transposon mutant libraries used were either created in-house for this study, or obtained from literature and reanalysed. The M. catarrhalis BBH18 marinerT7 transposon mutant libraries consisting of 28,000 and 7,000 independent transformants were previously described [18, 19], and the 12,500 transformants library was generated using the previously described protocol [18]. The 40,000 transformants S. pneumoniae R6 and the 11,000 transformants H. influenzae 86 028NP library were previously described [20, 21]. Libraries for the 15,000 transformants S. pneumoniae R6 and H. influenzae Rd KW20 were also respectively constructed as previously described [20, 21]. The Tn-seq technology was used to profile the relative abundance of each mutant in all libraries after growth as described previously [22], except for S. pneumoniae TIGR4. Tn-seq data for S. pneumoniae TIGR4 were obtained from literature [23]. We then performed essentiality predictions for individual genes using the in-house developed web-tool, ESSENTIALS [24], which enabled us calculate a statistical essentiality metric for each ORF, and precisely delineate the optimal boundary between essential and nonessential ORFs in each of the 5 strains. Analysis data can be found at http://bamics2.cmbi.ru.nl/websoftware/essentials/links.html.

Overrepresented metabolic pathways and subsystems

Pathways and subsystems for the strains under study were obtained from the Kyoto Encyclopedia of Genes and Genomes orthology, and the SEED databases respectively [25, 26]. Using a Fisher’s exact test, we performed functional categories enrichment for the pathways and subsystems, while incorporating the statistical essentiality value (the fold-change value predicted by ESSENTIALS) for each ORF. We corrected for multiple testing using Bonferroni correction and obtained q-values for corresponding p-values [27].

Proteins subcellular localization (SCL)

The subcellular localizations (SCL) of all proteins in this study were determined using publicly available SCL prediction tools. First, we analysed all Gram-positive and Gram-negative strains using pSORTdb version 2.0 [10] and CELLO version 2.5 [28]. Further complementation SCL predictions were performed using LocateP and GnegPloc for Gram-positive and Gram-negative strains respectively [29, 30]. Additionally, the presence of integral Gram-negative outer membrane proteins (OMP) was determined using β-barrel outer membrane protein predictor (BOMP) [31]. Proteins that showed different SCL predictions in the different predictors used were denoted “Unknown”, together with those predicted to be of unknown SCL by majority of the predictors used.

Selecting potential drug targets

To identify and eliminate essential genes with close undesirable orthologs, we performed separate unidirectional protein-protein BLAST (BlastP) searches, using an E-value cut-off of 1e-10, and minimum 70% sequence identity over 75% sequence coverage; against the human genome, and the metagenomics catalogue of non-redundant human gut microbiome genes by Qin et al.[7].

Determination of antimicrobial activity

Selection of potential drug targets for in vivo validation was mainly based on their novelty, that is, they have not been described as targets to existing antimicrobials. Commercial availability of inhibitory compounds without resorting to customized chemical synthesis was also key; all inhibitory compounds used were supplied by Sigma Aldrich. 1-Methyluric acid, 5, 5′-Dithio-bis-(2-nitrobenzoic Acid), and 5′-deoxyadenosine were dissolved in water at 5 mg/ml. When necessary, the pH was neutralized (to pH7) using 10 M NaOH solution or 1 M HCl. Antimicrobial activity of the compounds was tested by Kirby-Bauer/disk diffusion assay [32], by applying 10 μg of the inhibitory compounds to 6 mm filter paper discs at concentration ranging from 10000 to 0.05 μg/ml in 10-fold stepwise dilutions. As for (R)-6-fluoromevalonate diphosphate 2 μl of (R)-6-fluoromevalonate diphosphate was diluted in 1 ml of Milli-Q (MQ). 10 μl and 100 μl of the dilution was used in separate disk diffusion assays. Columbia III agar with 5% sheep blood medium was used for S. pneumoniae. Brain heart infusion (BHI) agar medium, and a combination medium of BHI, hemin, and NAD were used for M. catarrhalis and H. influenzae respectively. MIC calculations were performed as described by Wiegand and colleagues [33]. Experiments were performed in quadruplicate, and outliers were removed using the Grubbs test [34].

Toxicity assays on epithelial cell lines

Cellular toxicity of (R)-6-fluoromevalonate diphosphate was tested using the CellTox Green Cytotoxicity Assay (Promega, WI) on Detroit 562 (ATCC CCL-138) and A549 (ATCC CCL-185) cell lines according to the manufacturer’s instructions. The two cell lines were exposed to (R)-6-fluoromevalonate diphosphate at its effective MIC concentration of 26.6 μg/ml for 24 hours at 37°C with 5% CO₂. Fluorescence was measured on a Perkin Elmer 1420 Victor 3 V multi-label plate reader.

Genome reannotation and gene clustering

We sought to determine potential drug targets in S. pneumoniae, H. influenzae, and M. catarrhalis following the selection criteria outlined in Figure 1. For these species, five strains with the required Tn-seq data were available; S. pneumoniae strains R6 and TIGR4, H. influenzae strains Rd KW20 and 86 028NP, and M. catarrhalis strain BBH18. Altogether, genomes of these strains in their initial annotations constituted of 10,072 open reading frames (ORFs). These annotations were updated using RAST to ensure consistency and comparability among strains in subsequent analyses. This analysis resulted in putative annotations for about 50% of all ORFs originally annotated with a hypothetical function (Table 1; Additional file 1). Next, we clustered the updated protein sequences using OrthoMCL, producing 1,798 orthologous groups/clusters (OGs) with, and 2,729 without singletons respectively (Additional file 1). This clustering of orthologous proteins allowed for the determination of species and/or strain specific proteins, as well as determining the metabolic potential of the strains. For example, the “Gram-negative specific” periplasmic chaperones (SurA) were clustered in OG_756 (cluster 756), while the “Streptococci-specific” transcriptional regulators (LytR) were clustered in OG_2554. On the other hand, 300 OGs, including OG_184, OG_186, OG_216, and OG_224, among others, contained genes conserved in all the five strains. All protein in individual OGs constituted of similar or identical functional annotations. This consistency in grouping and annotation was observed across all OGs, suggesting a reliable clustering. Confirmatory clusters and respective annotations derived from the clusters of orthologous genes (COGs) database were consistent with our OrthoMCL clusters. Additionally, using the OG’s, we were able to curate annotations for the HI1586 locus in Haemophilus influenzae Rd KW20, which was possibly misannotated in the initial release, as an isoleucyl-tRNA synthetase instead of a NaP⁺P/HP⁺P antiporter.

Essential and conserved protein-coding genes

Loss of mutant readouts from a transposon library after in vitro transposition and genetic transformation of the wild-type isolate is a strong indicator of gene essentiality [35]. Although some essential genes tolerate disruptive insertions in the 3′ regions, generally, insertions in essential genes lead to lethal phenotypes [36]. For our analysis, mutant libraries and/or Tn-seq data were constructed in in-house experiments or obtained from literature (Table 1). We separately analysed the Tn-seq datasets using ESSENTIALS [24]. This analysis resulted in the identification of 532 essential genes in H. influenzae 86-028NP, representing 28% of the genome, a higher number as compared to the other Gram-negative strains; H. influenzae Rd KW20 and M. catarrhalis BBH18, in which we identified 431 and 445 essential genes respectively. In S. pneumonia, we identified 325 and 414 essential genes for the R6 and TIGR4 strains respectively (Table 1; Additional file 1). These values showed that on average, about 20% of all genes in the five strains are essential. This is consistent with earlier studies which have reported 15-25% of all genes in a genome being essential [23, 36, 37].

Differences in the number of essential genes could be explained by various factors that hamper precision in transposon mutagenesis experiments, including short gene lengths and unsaturated transposon libraries; “saturation” being the presence of at least one insertion in every gene. In practice, short genes are less susceptible to disruptive transposon insertions, hence, more likely to be misclassified as essential. In unsaturated transposon mutant libraries, dispensable genes are also more likely to be devoid of transposon insertions, and therefore misclassified as being essential genes. The low-density transposon mutant library (approximately 11,000 colony forming units; CFU) used for H. influenzae 86-028NP, and a substantial number of short genes in its genome could, therefore, explain the apparently overestimated (532) essential genes. Relatively saturated libraries of approximately 20,000 CFU and 40,000 CFU were used for H. influenzae Rd KW20 and M. catarrhalis BBH18 respectively (Table 1). A rarefaction analysis on our data confirmed that the S. pneumoniae, M. catarrhalis, and H. influenzae Rd KW20 transposon libraries approached saturation (Additional file 2). Additionally, based on derivations of Poisson’s law, there is a 99.6% probability that genes with a size of 1 kb are hit in the 1.9 Mb H. influenzae 86-028NP genome and an 11,000 CFU mutant library. Similar statistics on the 1.79 Mb H. influenzae Rd KW20 genome with a 20,000 CFU mutant library shows a 99.99% probability. Therefore, H. influenzae 86-028NP could have suffered slightly more false positive predictions due to its less saturated mutant libraries.

We selected 705 OGs containing at least one essential gene from any of the five strains for further analysis. These essential OGs mainly consist of proteins with annotated functions, participating in diverse core cellular processes, such as DNA replication, DNA transcription, protein translation, cell wall biosynthesis, signal transduction, and metabolism. Eighteen OGs, however, contained conserved proteins of uncharacterized function (Additional file 1). Functional characterization of these genes will aid in achieving the optimal set of targets that can be used to develop antimicrobials against RTI causing bacteria. The distribution and overlap of the essential genes within the three species is outlined in Figure 2. From the 705 OGs, we identified 128 OGs that constituted of genes conserved and essential in all five strains, representing targets particularly attractive for developing broad-spectrum antimicrobials to treat RTI, since they encode components of basal cellular functions in respiratory pathogens. Importantly, collective analysis of the five strains revealed species-specific and/or “Gram-category” specific essential genes, best suited for narrow-spectrum or specialized inhibition.

Essential metabolic pathways and subsystems

Protein subcellular localization

Out of the 705 OGs selected, the majority (526) consists of cytoplasmic proteins. Cellular localization of the other OGs were predicted to be: 96 in the inner membrane, 11 in the outer membrane, 12 in the periplasm, and 4 in the extracellular space. In addition, 21 OGs are non-categorically predicted to contain membrane proteins, whereas 35 are of unknown localizations. Of the 11 outer membrane OGs, 7 contained β-barrels (Additional file 1).

Orthologs in human and human gut microflora

The human gut is home to microbiota whose proper composition and functioning collectively influence human nutrition, protection against pathogens and development of disease [7]. Perturbing this microbiota with antibiotics could cause adverse side effects. Furthermore, interference with human cell physiology by antibiotics as a consequence of non-specific targeting can cause severe cellular cytotoxicity [3], which may result in organ failure or even death. We used blastP analyses against the human genome (Genome Reference Consortium) and the human gut microbial gene catalogue [7], to identify targets that would likely have off-target effects. It is noteworthy that targets with as few as 10 matches in the non-redundant gut microbial gene catalogue were allowed in the final selection, as we hypothesised that these would have no effects on the gut microbiome preventing disruption of gut health. This decision was motivated by the observation from our analysis that well known targets for both clinically approved antimicrobials and experimental small molecule inhibitors collated in DrugBank (Additional file 1; column 9) maintained on average fewer than 10 blast hits against the human gut microbial gene catalogue (Additional file 1; column 20). On the other hand, the majority of the targets with numerous blast hits were aminoacyl-tRNA synthetases (aaRSs) and ribosomal protein, including rpsL, a well-known target that had 249 hits for pneumococci, 156 for H. influenzae, and 151 for M. catarrhalis. One shortcoming of using such filtering criteria is that novel targets that have more than 10 blast hits are not effectively retained in the final selection. Nevertheless, we identified 96 OGs with orthologs in human, and 127 OGs with orthologs in human gut microflora, that is, with >10 blast hits (Additional file 1). All 20 aminoacyl-tRNA synthetases (aaRSs), essential for protein synthesis, were particularly conserved in both human and human gut microflora. Studies have shown that aaRSs can be selectively targeted as most bacterial aaRSs recognize and aminoacylate only cognate tRNA [38]. However, possible side effects are expected from drugs targeting aaRSs. RNA molecules and ribosomal proteins were also highly conserved in gut microbiota and humans. Additionally, the relatively short lengths and the presence of highly repetitive DNA in RNA sequences also rendered their essentiality predictions unreliable. All these molecules were therefore not included in the final selection of drug targets. Moreover, blast comparison between finally selected targets and their human orthologs showed minimal sequence identities (<35%) over short sequence coverage.

Drug targets selection and validation

We identified 249 potential drug targets in the five strains (Additional file 5), including key enzymes in pathways such as fatty acid biosynthesis [39–41], vitamin biosynthesis [42–45], and isoprenoid biosynthesis pathways [46–48], which have gained interest in drug discovery research, as well as 67 known targets inhibited by 75 FDA-approved antimicrobial drugs and 35 other researched small molecule inhibitors collated in the DrugBank database [49]. To validate our target prediction, we selected four novel targets with commercially available novel inhibitors of their predicted essential functions, that is, inhibitors not yet approved as clinical drugs and don’t require to be custom synthesized: We tested whether exposure to these compounds inhibited growth of the target organisms.

As an essential amino acid, methionine is not synthesized de novo in humans, who must rely on dietary intake. Enzymes involved in microbial methionine biosynthesis therefore offer highly specific and selective drug targets. We used 1-methyluric acid (CAS: 708-79-2) to target S-adenosylmethionine synthetase (EC: 2.5.1.6); a key enzyme in methionine biosynthesis, whose drug target potential has been explored in various pathogens [53, 54]. Contrary to expectations, no growth inhibition was observed in Gram-negative strains (Table 3; Figure 3): growth inhibition was only observed in S. pneumoniae. Since 1-methyluric acid formed a precipitate in concentrations above 312 μg /ml, no MIC values could be calculated. This lack of growth inhibition in Gram-negative strains may possibly be due to their double layered cell walls which are less penetrable [55], or the bacteria have expanded their resistance mechanisms to evade killing by antimicrobials [55, 56]. It is also possible that the two Gram-negative species have alternative mechanisms for methionine biosynthesis, further complicating screening for effective drugs.

The microbial fatty acid synthesis (FAS) pathway is an attractive target for drug discovery [41, 57]. This pathway is subdivided into type I and II, whereby human FAS proteins predominantly belong to type I FAS, and the bacterial ones are predominantly type II FAS. Proteins from the two FAS types generally possess distinctive molecular organization of the active site allowing for selective targeting [39, 40]. Although Gram-positive pathogens could compensate FASII inhibition by assimilating environmental fatty acids; particularly unsaturated fatty acids [58, 59], several clinical and household antimicrobials targeting key FAS enzymes, e.g. Platensimycin and Platencin have been successfully developed [41, 60]. In our analysis, we identified various genes conserved in all five strains, for example genes in OGs 085, 143, and 653, whose products play key roles in the FAS pathway. With 5′-Deoxyadenosine (CAS: 4754-39-6), we targeted lipoate synthase (LipA; EC: 2.8.1.8), a key enzyme in the lipoic acid metabolism [61], using product-level inhibition. Surprisingly, we observed growth inhibition in all three species (Figure 3; Table 3), despite the target cluster (OG_653) comprising of orthologs from only Gram-negative strains (Additional file 1). This observations are also reflected in the MIC, which ranged from 29.3 to 205.1 μg/ml (Table 3). A blastP comparison showed that the closest ortholog of the Gram-negative LipA in S. pneumoniae is the non-lipoic pathway enzyme fructose-6-phosphate aldolase I, sharing about 32% sequence identity. Moreover, a comparison between LipA and lipoate-protein ligase (LplA), the key lipoylation enzyme in S. pneumoniae[61], revealed that the two proteins are non-orthologous, as they share very low sequence identity (<25%). They however have conserved domain which may explain the observed growth inhibition.

Isoprenoids are natural products involved in many biochemical functions, such as supplying quinones for the electron transport chains, components of membranes, and subcellular targeting and regulation [47]. Humans employ the mevalonate pathway, whereas most microbes follow a non-mevalonate (1-deoxy-d-xylulose 5-phosphate/2-C-methyl-d-erythritol 4-phosphate) pathway. Functional roles of key enzymes in the isoprenoid biosynthesis pathway are well characterized, opening prospects for the discovery of novel drug targets [46, 48]. Fosmidomycin is a promising isoprenoid-based anti-malarial drug which is currently in clinical trials [48]. Using 6-fluoromevalonate (CAS: 2822-77-7) to target diphosphomevalonate decarboxylase (EC: 4.1.1.33), we observed selective growth inhibition only in S. pneumoniae as expected (Figure 3; Additional file 1; Table 3). Additionally, no effects on growth were observed in the Gram-negative strains, which was also as expected. We determined an average MIC of value 26.6 μg/ml for the S. pneumoniae growth inhibition (Table 3). At 26.6 μg/ml, no toxicity was observed in cell toxicity assays on epithelial cell lines (data not shown). Moreover, in patent WO 1995013058 A1, no cytotoxic effects of 6-fluoromevalonate were observed on T-lymphocytes. Previous literature also shown that 6-fluoromevalonate could potentially function the same as statins, as they inhibit the same pathway [62]. Diphosphomevalonate decarboxylase could therefore be a promising target for developing novel antibiotics against S. pneumoniae[63].

Availability of supporting data

The data sets supporting the results of this article are included within the article and its additional files. Tn-Seq data sets are available in the European Nucleotide Archive repository, [http://www.ebi.ac.uk/ena/data/view/PRJEB7553].