Illumina-based whole genome sequencing and assembly of strain PBC has been previously described (Gan et al., 2012b). To generate sequencing data using the MinION device, the gDNA of strain PBC was extracted from a 5-day old culture on nutrient agar using a slightly modified SDS-lysis method (Sokolov, 2000). Five micrograms of gDNA was subsequently used for library construction using the Nanopore sequencing kit SQK-NSK007 (Oxford Nanopore Technologies, Oxford, United Kingdom) according to the manufacturer’s instructions. The library was subsequently loaded onto an R9 MINION flowcell (Oxford Nanopore Technologies, Oxford, United Kingdom) and run for 48 h. 2D-basecalling was performed on the cloud-based Metrichor. Fast5 to fasta conversion used Nanocall (David et al., 2016).

SPAdes version 3.8.1 (Antipov et al., 2016) was used for an initial hybrid genome assembly incorporating both Illumina and Nanopore reads. Contigs longer than 1,000 bp were then selected for in silico scaffolding and gap-closing using npScarf (Cao et al., 2016). The assembled genome was annotated using the NCBI Prokaryotic annotation pipeline (Tatusova et al., 2016) and manually curated to include annotations for the pcaH2, pcaG2, pcaB2 and 4SLH genes coding for protocatechuate 3,4-dioxygenase alpha subunit type II, protocatechuate 3,4-dioxygenase beta subunit type II, 3-carboxymuconate cycloisomerase type II and 4-sulfomuconolactoe hydrolase, respectively, that were missed by the default annotation setting. The identification and annotation of prophage sequences were performed using the PHASTER web-server (Arndt et al., 2016). Visualization of the complete genome of strain PBC and genomic comparisons with public available Hydrogenophaga genomes (Table 1) were performed using Blast Ring Image Generator (blastN e-value setting of 1e-10) (Alikhan et al., 2011). Mapping of Nanopore reads and Illumina contigs was performed using Minimap2¹ and subsequently visualized in Integrative Genomics Viewer version 3 (Thorvaldsdóttir et al., 2013). Generation of Illumina contigs was done by assembling Illumina-only reads using SPAdes version 3.8.1 and filtering for contigs with coverage and length of more than 5× and 300 bp, respectively (Antipov et al., 2016).

Pair-wise average nucleotide identity of strain PBC against the publicly available genome sequences of Hydrogenophaga species were performed using JSpecies V1.2 (ANIm setting) (Richter et al., 2016). Strains exhibiting pair-wise ANI of more than 95% were considered as members of the same genospecies. A heatmap was plotted in R version 3 using the pheatmap library package.

Whole proteome was predicted using Prodigal (-p meta setting) and piped into PhyloPhlan for the identification of conserved proteins (Hyatt et al., 2010; Segata et al., 2013). The concatenated alignment generated by PhyloPhlAN was subsequently trimmed using trimAl version 1.9 (-automated1 setting) and used to construct a maximum likelihood tree using IqTree version 5.15 with 1000 ultrafast bootstraps (Capella-Gutiérrez et al., 2009; Nguyen et al., 2015). Visualization and annotation of trees was performed with MEGA6 (Tamura et al., 2013).

HMMsearch version 3.1b was used for the identification of proteins of interest utilizing either the PFAM (-E 1e^-5 setting) or TIGRFAM (-cut_nc or –E 1e^-5 setting) HMM profiles (Haft et al., 2003; Eddy, 2011). Protein alignment, trimming, phylogenetic tree construction and visualization were performed using MAFFT-linsi, TrimAl, IqTree and FigTree, respectively (Capella-Gutiérrez et al., 2009; Katoh and Standley, 2013; Nguyen et al., 2015). Visualization of gene neighborhood was conducted using EasyFig (blastN setting) (Sullivan et al., 2011).

A total of 90,491 1D and 51,630 2D reads were generated from a full 48-h MinION run. After filtering for reads longer than 2,500 bp, 14,404 1D and 12, 959 2D reads were retained for genome assembly. The final data output used for assembly was 173 megabases representing approximately 33× genome coverage of strain PBC. The average read length was 6,350 bp and the longest read generated was 126,959 bp. Initial Spades hybrid assembly using both Illumina and MinION filtered reads generated 5 scaffolds (>10x coverage, >1,000 bp) which were subsequently scaffolded and gap-closed into a single contig using Npscarf. The complete genome of strain PBC is 5,232,477 bp in length with a GC content of 68.4%. Based on NCBI annotation pipeline, the complete genome consists of 4,904 protein-coding genes, 43 tRNA and single-copy 5S, 16S, and 23S rRNA genes. Comparisons of strain PBC and other Hydrogenophaga species revealed two genomic regions (Figure 1A: Regions 3 and 4) uniquely present in strain PBC, that were subsequently identified as prophage in origin by Phaster. Both phage regions were incomplete and only code for two to three major phage components. Phage components 1 and 2 on the contrary are complete, uniquely shared by both H. intermedia strains (Figure 1A) and have the highest number of homologous proteins to Enterobacter phage Arya (NCBI Reference Sequence: NC_0310480) isolated from a termite gut.

Genes responsible for the catabolism of 4-aminobenzenesulfonate in strain PBC are clustered in two distinct genomic regions. Interestingly, the sad operon appears to be fairly conserved in all currently reported Hydrogenophaga genomes, suggesting that the ability to deaminate 4-aminobenzenesulfoate may be a defining metabolic trait of Hydrogenophaga strains. However, given the diverse isolation sources of Hydrogenophaga strains (Table 1) and their presumably limited exposure to xenobiotic compounds, it is likely that the original substrate of the 4-aminobenzenesulfonate dioxygenase is a naturally occurring compound with close structural homology to 4-aminobenzenesulfonate.

In contradistinction, the genomic region containing pcaH2G2, pcaB2 and 4-SLH involved in the conversion of 4-sulfocatechol to maleylacetate is conserved only in the Hydrogenophaga intermedia strains (Figure 1B). In addition, it is also worth noting that this region was identified as a phage region and flanked by genes coding for IS6 transposable elements (Figure 1C), an indication of the mobility of the pcaH2G2, pcaB2 and 4-SLH gene cluster presumably through phage transduction (Hickman et al., 2010). As expected, the presence of multiple transposase genes in this genomic region led to a fragmented initial Illumina read assembly (Figure 1C) due to the inability of short reads to resolve the repetitive regions (Phillippy, 2017). This shortcoming was overcome by integrating Nanopore long reads into the assembly as evidenced by their ability to span repeats in the genome (Figure 1C).

It is tempting to speculate that the acquisition of this gene cluster may be a key factor leading to complete 4-aminobenzenesulfonate biodegradation in Hydrogenophaga. Recently, the complete genome of Novosphingobium resinovorum SA1, another well-studied 4-aminobenzenesulfonate degrader, has been reported and comparison of this genomic region and that of the homologous region in Novosphingobium resinovorum SA1, revealed a comparable gene organization, e.g., pcaH2-pcaG2-tnp-pcaB2-4SLH-MAR in H. intermedia PBC (Figure 1C) and pcaB2 (scaA)-4SLH (scaB)-MAR(scaC)-oxidoreductase-pcaH2(scaE)-pcaG2(scaF) in Novosphingobium resinovorum SA1 (Hegedüs et al., 2017), suggesting recently shared ancestry of the currently known 4-sulfocatechol biodegradation pathway and possibly inter-class (alphaproteobacteria-betaproteobacteria) horizontal gene transfer of this gene cluster.

Hydrogenophaga flava was initially designated as the type species of the genus Hydrogenophaga but due to its slow and unreliable growth, H. pseudoflava has been proposed for use as an alternative reference taxon for the genus as it exhibits similar genotypic and protein electrophoretic profiles to that of H. flava (Willems et al., 1989). In silico genome-genome hybridization indicates that although H. flava and H. pseudoflava are closely related to one another, they are clearly two distinct genospecies with a pairwise ANIm of 93% (Figure 2). A more in-depth analysis of the H. flava genome may be useful to identify its additional nutritional requirements for optimal and consistent growth on culture medium. Such a genomic approach has been successfully demonstrated in Clostridium tyrobutyricum and other Clostridia associated with butyric acid fermentation leading to the identification of key vitamins and several amino acids essential for growth (Storari et al., 2016). Furthermore, in silico genome-genome hybridization analysis also indicates that H. sp. H7, which was previously isolated from a coal mine in China (Table 1), could be reclassified to the species H. pseudoflava given its high pairwise ANIm (97%) to H. pseudoflava NBRC 102511^T (Figure 2).

Maximum likelihood tree based on 400 universal conserved proteins provides maximal support for the monophyletic clustering of all currently sequenced Hydrogenophaga strains with H. intermedia strains and H. sp. PML 113 being basal to the rest of Hydrogenophaga (Figure 3). The genus Hydrogenophaga clade shares a sister relationship with a group of highly alkaliphilic hydrogen-utilizing Comamonadacea bacteria isolated from an alkaline serpentinizing springs at The Cedars, California (Suzuki et al., 2014). Beyond the genus Hydrogenophaga, however, considerable taxonomic incongruence was observed in the genera Acidovorax and Comamonas. Some members of the Acidovorax, such as Acidovorax sp. strains 121606, 202149 and 12322-1, could be reclassified as members of the genus Comamonas based on their phylogenetic affinities. Within the genus Comamonas, Comamonas badia DSM 17552^T (Suzuki et al., 2014) is peculiar as it does not cluster with a majority of other Comamonas, but instead forms a sister group with members of the genus Alicycliphilus, albeit with moderate bootstrap support. In addition, Comamonas granuli NBRC 101663^T appears to be considerably divergent from other Comamonas species given its basal position in the clade containing all currently sequenced Comamonas strains (Suzuki et al., 2014). A more in-depth taxonomic investigation based on the percentage of conserved proteins (POCP) to specifically re-define the genus boundary between Comamonas and Acidovorax should be undertaken (Qin et al., 2014).

The genus Hydrogenophaga was originally proposed by Willems et al. (1989) to group “yellow-pigmented hydrogen-oxidizing [Pseudomonas] species belonging to the acidovorans rRNA complex.” However, by identifying currently known HMM profiles related to hydrogenase system, genomic potential for hydrogen oxidation was only identified in 7 Hydrogenophaga strains of which 6 strains are closely related as shown by their monophyletic clustering in the phylogenomic tree (Figure 3). The 7 Hydrogenophaga strains include the original strains that were initially used to describe the genus Hydrogenophaga, such as H. palleronii, H. taeniospralis, H. flava and H. pseudoflava, corroborating their previously demonstrated in vitro hydrogen oxidizing ability (Willems et al., 1989; Contzen et al., 2000). The reported inability of H. intermedia S1^T and H. atypica DSM 15342^T to oxidize hydrogen in vitro (Contzen et al., 2000; Kampfer et al., 2005) correlates with the absence of key genes associated with hydrogen metabolism (Figure 3). However, the absence of genomic potential for hydrogen catabolism in nearly half of the currently sequenced Hydrogenophaga strains is unexpected and indicates that the description genus of Hydrogenophaga as “hydrogen eating bacteria” is not sustainable and warrants revision.

Mapping of Hydrogenophaga genome sequences to the contig of H. palleronii NBRC102513 (BCTJ01000048, 24 kbp) containing pabB revealed that only 8 out of 14 genomes exhibit nucleotide similarity (>40% identity and E-value < 1E^-5) to the pabB gene region (Figure 4A). Of the eight genomes, only partial pabB hits were observed for strains H7, NBRC102514, Root209, NBRC102511 and A37, hinting divergence at the nucleotide level. Interestingly, only two proteins belonging to H. palleroni and H. sp. PML113 achieved an HMMscore that is above the domain noise cutoff score for bona fide PabB (TIGR00553) (Supplementary Figure S1). On the contrary, strain IBVHS2 despite exhibiting nucleotide similarity to the pabB gene region, does not code for protein with HMMscore above the TIGR00553 noise cutoff. However, it is worth noting that multiple proteins with HMMscore slightly below the TIGR00553 noise cutoff could be identified from the annotated Hydrogenophaga genomes (Supplementary Figure S1). These proteins appear to be more closely related to PabB than TrpE as evidenced by their extremely low HMMscore to TrpE (TIGR00564) in comparison with other putative TrpE proteins and their phylogenetic affiliation to the PabB clade (Supplementary Figure S2 and Data Sheet 1). TrpE is involved in the synthesis of anthranilate from chorismate and ammonia. Functional characterization of these divergent PabB homologs will be useful to explore their role in pABA synthesis as of them belong to Hydrogenophaga pseudoflava that has been shown to grow on minimal medium without pABA supplementation (Povolo et al., 2013).

The absence pabB in H. intermedia strains PBC and S1^T is consistent with their previously reported reliance on either external nutrient supplementation or helper strain for pABA (Dangmann et al., 1996; Gan et al., 2011b; Kim and Gan, 2017). A recent genomic analysis of Hydrogenophaga intermedia PBC helper strain, Ralstonia sp. PBA, showed the presence of two PabB proteins of which one of them was fused to PabC (Gan et al., 2012a; Kim and Gan, 2017), potentially enabling the overproduction of pABA to sustain the growth of both strains when they were co-cultured in vitamin-free minimal medium containing 4-aminobenzenesulfonate as the sole carbon and nitrogen source. A closer comparison of the genomic region containing genes spatially associated with pabB in H. intermedia PBC and H. palleroni NBRC102513 suggests that pabB may be lost through gene deletion (Figure 4B). Surprisingly, strains LPB0072, RAC07, Root209 and IBVHS1 appear to also lack the genomic potential for pABA synthesis. A future study investigating the ability of these strains to grow on defined medium without pABA supplementation will be necessary to verify their predicted pABA auxotrophy as most Hydrogenophaga strains were isolated on nutrient rich medium, which will contain traced amount of pABA. For example, strains Root209 and RAC07 were maintained on tryptic soy agar and minimal medium with yeast extract supplementation, respectively (Bai et al., 2015; Fixen et al., 2016).