Comparative genomic analysis of emerging non-typeable Haemophilus influenzae (NTHi) causing emerging septic arthritis in Atlanta

Oxford nanopore sequencing and hybrid assembly with Illumina data

One Haemophilus influenzae isolate was randomly chosen from each of the C1 and C2 clusters for hybrid assembly of Oxford Nanopore minIOn and Illumina data to produce complete reference sequences (referred to as C1-1 and C2-1). Genomic DNA was extracted using the Promega Wizard Genomic DNA Purification Kit. Sequencing libraries prepared using the SQK-LSK109 1D ligation sequencing kit and sequenced on a FLO-FLG001 Flongle flow cell, yielding 496.9 Mb and 552.6 Mb of raw reads (~267x and ~297x coverage) for C1-1 and C2-1, respectively. The H. influenzae C1-1 and C2-1 genomes were then assembled from Nanopore and Illumina paired-end reads using Unicycler (Wick et al., 2017). Assemblies were deposited in the NCBI assembly database under GCF_044096575 (C1-1) and GCF_044096445 (C2-1). Raw data was deposited to NCBI as part of accession PRJNA544724.

H. influenzae virulence protein database

We created a non-redundant BLAST database of all of the H. influenzae virulence factors defined in the Virulence Factor Database (http://www.mgc.ac.cn/VFs/main.htm) and the Victors database (https://phidias.us/victors/) (Sayers et al., 2019; Chen, 2004). We also manually curated the recent literature to identify genes and sequences associated with virulence that may not have been included in the VFDB or Victors databases. In April 2024 we performed a search in pubmed for articles containing the term ‘“Haemophilus influenzae” [All Fields] AND (“pathogenicity” [MeSH Subheading] OR “pathogenicity” [All Fields] OR “virulence” [All Fields] OR “virulence” [MeSH Terms] OR “virulences” [All Fields] OR “virulent” [All Fields]) ‘published between 2004–2024. Gene names associated with virulence were identified from reading the articles and were matched to uniprot accessions from H. influenzae or related gram-negative bacteria. From this search, we added gene and protein sequences of 105 H. influenzae reference genes described in Pinto et al. (2019).

We created a BLAST database of all the protein sequences in our table. To match a query protein against our database we used blastp v2.10.1+ with default parameters and used a threshold of 80% identity to define a hit.

Processing whole genome data

We used “bactopia search” from the Bactopia software package (Petit & Read, 2020) to identify 4,842 publicly available H. influenzae genomes from the Sequence Read Archive database. Whole genome shotgun Illumina paired-end fastq data files were downloaded from NCBI SRA database on 05-April-2021 and processed using Bactopia (v1.6.0) (Petit & Read, 2020). The multi-locus sequence type (MLST) schema for H. influenzae from PubMLST.org (Jolley, Bray & Maiden, 2018) was included and completed genomes for C1-1 and C2-1 were added as optional datasets.

Within the Bactopia pipeline the reads were cleaned and error-corrected using BBDuk (v38.86, (Bushnell, 2014)) and Lighter (v1.1.2, (Song, Florea & Langmead, 2014)). Processed reads were assembled with SKESA (v2.4.0, (Souvorov, Agarwala & Lipman, 2018)) using Shovill (v1.1.0, (Seemann, 2014)). Assembly quality metrics were determined with assembly-scan (v0.3.0 (Petit, 2017) and CheckM (v1.1.3, (Parks et al., 2015)). The species composition of the assembly was determined by screening against a minmer sketch of GenBank (Clark et al., 2016) using sourmash (v3.5.0, (Titus Brown & Irber, 2016)). The MLST was determined using BLAST+ (V2.10.1, (Camacho et al., 2009)) and Ariba (v2.14.6, (Hunt et al., 2017)).

Generating a representative dataset

Due to the large number of samples collected and diverse sequence types among them, we decided to create a smaller representative dataset for downstream analyses. This set would include at least one high quality sample from each sequence type, as well as all high-quality Georgia-specific samples. Both capsulated and unencapsulated samples were included in this set to enable a broader pan-genome comparison when assessing unique traits of C1 and C2. First, we generated an intermediate group of high-quality samples that would include candidates for the final representative dataset. We used the “summary” tool from Bactopia to aggregate the quality control results for all genomes into a single table. This included assignment of a rank, either gold, silver, or bronze, based on coverage, mean per-read quality, mean read length, and total number of contigs. The cutoffs for each category were:

Gold:

Coverage >= 100x

Quality >= Q30

Read Length >= 95 bp

Total Contigs < 100

Silver:

Coverage >= 50x

Quality >= Q20

Read Length >= 75 bp

Total Contigs < 200

Bronze:

Coverage >= 20x

Quality >= Q12

Read Length >= 49 bp

Total Contigs < 500

Any samples lower than a bronze quality rank were dropped. The average nucleotide identity (ANI) between C1 and all other genomes was calculated with FastANI (Jain et al., 2018) (v1.32 (Souvorov, Agarwala & Lipman, 2018)) to determine relatedness between samples.

To generate the intermediate data from the post-QC samples, a custom Python script (‘generate-representative-set.py’) was created. The script tags representative samples by prioritizing Georgia samples first, then samples with ANI greater than a specified cutoff. For any remaining non-Georgia non-C1/C2 samples, it selects representatives based on the sample rank (gold, silver, bronze), CheckM completeness, and the number of contigs, prioritizing higher completeness and fewer contigs. Samples are excluded if they do not meet the criteria for paired-end reads, fail quality checks, or have undesirable sequence types, and the exclusion reasons are documented in an output file. The script updates exclude and representative files accordingly, utilizing Python library argparse for command-line argument parsing, os for file operations, and collections for ordered dictionary management.

For the final representative dataset, we further reduced the intermediate data to the best single representatives of each sequence type. Any genomes with more than 500 assembled contigs, those that did not have a sequence type determined, or were not screened as “Haemophilus influenzae” by sourmash (Titus Brown & Irber, 2016), were excluded. The representative set included all high-quality C1, C2, and non-cluster isolates from the GA EIP investigation from this study and all genomes that had an ANI greater than ANI of the most distant C2 genome from C1-1 (a member of C1). This was done to ensure that no highly divergent samples were included. If an ST was identified, but not yet included, a genome was picked prioritizing its quality rank, assembly completeness, and total number of contigs. The final representative set was used to determine the pan-genome using PIRATE (v1.0.4, (Bayliss et al., 2019)). A recombination masked core genome alignment was produced using ClonalFrameML (Didelot & Wilson, 2015) using a maximum-likelihood calculation. A phylogenetic tree based on this alignment was created with IQ-TREE (v2.0.3, (Nguyen et al., 2015; Hoang et al., 2018)) with default parameters.

Intra-clade comparisons

SNPs were called using the bactopia-tools workflow for snippy (Seemann, 2014). To understand if there was any variability across the entire genomes of C1 and C2 samples, reads from each sample were individually mapped to the completed reference genome of each and potentially deleted regions were assessed by lower-than-expected coverage in regions. In 1,000 bp windows, sliding every 250 bp, reads were counted, adjusted for reads per kilobase per million (RPKM), and selected if one of the samples had 0 reads over that bin. Adjacent regions of zero coverage were combined into a single region of interest and coverage was converted into a binary presence or absence.

Genome-wide association study (GWAS)

We used the pangenome GWAS tool Scoary (Brynildsrud et al., 2016) to investigate whether blood or sputum-associated genes were enriched in C1s and C2s. We randomly choose one representative sample from each ST collected from either blood or sputum to reduce bias introduced by oversampling a small number of STs with many genomes. This resulted in a set of 146 blood and 87 sputum-associated genomes that were used in the GWAS.

Recombination analysis

The core-genome alignment of 50 randomly chosen STs (including C1 and C2) with the sample() function in R with set.seed = 9,696 and recombination events were predicted by ClonalFrameML. The output includes a recombination-free alignment and a set of predicted recombination events, which were subsequently analyzed to understand the patterns and impact of horizontal gene transfer on the core-genome evolution of the selected STs

AMR analysis

We used the AMRFinder tool (Feldgarden et al., 2021) with default settings to identify antimicrobial resistance (AMR) genes in the analyzed sequences. The presence of these AMR genes was then plotted onto the phylogenetic tree to visualize their distribution across different strains, highlighting potential evolutionary relationships and patterns of resistance.

Inter-species comparative genomic analyses

The MAUVE (Darling et al., 2004) tool was used for comparative analyses between NTHi and other Haemophilus reference species to identify shared or unique ancestral sequences. The reference species H. aegyptius (GCA_900475885), H. parainfluenzae (GCA_000191405), H. haemolyticus (GCA_900477945) and H. seminalis (GCA_008605885) were compared to our NTHi samples using MAUVE.

Data availability

Commands, code, and data used in this study (including the “generate-representative-set.py” script, PIRATE gene family results, and curated blastp virulence gene database) are available at https://doi.org/10.5281/zenodo.14861587 and https://doi.org/10.5281/zenodo.14861572.

Genetic variation within the C1 and C2 clusters

There were 26 C1 and 23 C2 isolates with whole genome shotgun (WGS) Illumina data from the original Atlanta investigation from 2008–2018. One isolate was picked at random from each group (“C1-1” and “C2-1”) for closing using Oxford nanopore sequencing techolony. The final C1-1 hybrid assembly included one circular 1.875 Mb contig, while C2-1 included one circular 1.885 Mb and one circular 37.6 Kb plasmid. Table S1 shows hits of the virulence gene database to coordinates on the assembled C1-1 and C2-1 genomes.

The C1 isolates were found to be highly related, with a maximum distance of 132 SNPs in the 1.55 Mb core genome alignment (Fig. 1A). All C1 isolates were of the same sequence type, ST164. There were 29 unique regions of low coverage found in at least one sample (hereafter “deletion regions”) in the C1 isolates compared to the reference genome, ranging from 1 to 17 kb and encompassing 48 genes (Fig. 2A). Of these deletion regions, 15 genes had a significant match to H. influenzae virulence genes (ghfB, ghfC, ghfD, ghfE, hafB, hafC, hafD, hafE, hgpB, hgpD, hifB, hifC, hifD, hifE, and tbpB). The C2 isolates were also found to be highly related, all representing the same sequence type, ST1714, and having a maximum distance of 149 SNPs in the 1.60 Mb core genome alignment (Fig. 1B). The C2 cluster was separated into 2 subclades, with sample SNP distances as low as 1 and as high as 35 SNPs within these sub-groups. In a similar manner to C1 analysis, there were 13 large deletion regions (~1 to 32 kb), spanning 24 annotated genes (Fig. 2B), 6 of these genes matched suspected virulence genes: (hgpB, hgpD, licA, licB, licC, mrsA/glmM). Of the genes that had annotated gene names, five were deleted in at least one sample in both C1 and C2 (eamA, ninG, sRNA-Xcc1, sRNA-Xcc1, and tolB). None of these genes are found in the virulence factor database or the literature-curated database.

Few accessory genes distinguish the C1 and C2 clusters from other H. influenzae strains

Our public data search yielded 536 distinct MLST sequence types from a total of 4,842 available H. influenzae genomes. From this we randomly selected 536 ST representatives (one for each ST), in addition to 26 C1 isolates, and 23 C2 isolates for the pangenome analysis. This final group of 583 samples is referred to as the “representative dataset”. The pangenome consisted of 6,560 gene families, of which 1,368 were core (>= 95% genomes), 1,107 intermediate accessory (5% < x < 95%) and 4,085 rare accessory (x = <5%) (Fig. S1). Overall, we found C1 and C2 were part of a closely related subclade of NTHi strains (Fig. 3).

We found that there were no accessory genes unique to either C1 or C2, or both C1/C2 (i.e., present in one or both clades and not found in other MLSTs). There were also no core genes missing in only C1 or C2, or both. While there were no genes absent in C1 and/or C2 that were present in 100% of the rest of the H. influenzae pangenome, there was one gene family, identified as pxpB, that was lost in 100% of the C1 and C2 isolates and present in more the 90% of rest of the population of strains (Fig. 3). There were 20 gene families that were ‘rare’ in the context of the H. influenzae pangenome (i.e., found in < 10% STs) present in C1 but not C2 isolates (Table 1). Further, there were seven ‘rare’ gene families unique to C2, that were not identified in C1 strains (Table 2). There were no rare genes shared by both C1 and C2.

We used blastp to investigate the presence or absence of virulence factors in the predicted proteins from each of the 583 representative samples. We found the number of sequences present was comparable across most STs, with a total of 167 hits across all samples out of 207 total tested sequences (Fig. S2A). When comparing C1 and C2, there were two sequences unique to C1, three unique to C2, and 129 common between them (Fig. S2B). None of these rare unique genes were in our set of known virulence genes (Fig. S3). There were no known antimicrobial resistance (AMR) genes in the C1, C2 or closely related STs (Fig. S4).

C1 and C2 isolates possess a unique ancestral mobile cassette

While examining the H. influenzae pangenome data, we observed an interesting pattern of gains and losses that correspond to a previously undescribed mobile cassette inserted in the same ancestral region of C1, C2, and a small subset of other sequence types. There were eight genes identified as present in both C1 and C2 but in less than 90% of the rest of the pangenome (Fig. 3). All eight gene families had significant homoplasy, with a consistency index less than or equal to 0.2 on the species core genome tree (Table S1). Within this subset of rare genes, we identified five genes that were likely acquired as a cassette in the ancestor and inserted at the ancestral site of the one lost gene family. The gene that was the apparent site for cassette integration in C1 and C2 encoded a protein annotated as 5-oxoprolinase subunit PxpB (Niehaus et al., 2017). Using comparative genomic analysis with the MAUVE tool, we found that the five inserted genes were on a cassette of 9,444 bp in C1 and C2. In the same region of the outgroup that lacked the cassette (GCF_014701215.1_ASM1470121v1), there was an intact pxpABC gene cluster. When that region is compared to H. aegyptius, H. parainfluenzae, H. haemolyticus and H. seminalis using MAUVE, this region was not like any sequence in these reference pathogens.

Four of the five genes within the boundary of the cassette have been assigned a gene name and function by bakta annotations (Schwengers et al., 2021). One of those gene families is the gene tnpA, which is crucial for IS200/IS605 family transposition (Barabas et al., 2008; Morero et al., 2018). The remaining three were potentially part of a sugar metabolic operon based on their annotations. These were (ptsEII)-a sugar transmembrane transporter; malQ, 4-α-glucanotransferase important in maltose metabolism; and treR the repressor of the trehalose metabolic pathway. The presence of intact pxpB was highly homoplasic in H. influenzae (consistency Index of 0.14), which suggested a history of frequent gain and loss in the species, commensurate with cassette insertion. There was a pattern of lost and gained gene families found in other sequence types on the tree. We observed that there were 42 STs that have the exact same pattern of the lost pxpB gene family and concurrently gained the five gene families within the cassette. The pattern of gains and losses could be explained by multiple independent insertions of a cassette at the same location that disrupted pxpB. The metadata available for these 42 STs included isolates from both blood and respiratory infections. Many isolates were from infections of populations at risk of pulmonary infection such as patients with chronic obstructive pulmonary disease (COPD) and cystic fibrosis.

C1 and C2 have unique accessory gene profiles closer to H. influenzae isolates from blood than sputum

Instances of NTHi infection are usually respiratory in nature, but C1 and C2 infections were found to be largely systemic. This prompted us to do a GWAS of C1 and C2 to the larger NTHi pangenome to determine whether gene C1 and C2 profiles more closely matched samples collected from sputum or blood. Out of 4,842 H. influenzae isolates with public genomes, 1,624 had metadata indicating the isolates were collected from a blood or system infection and 1,441 were labeled as being isolated from sputum (the other genomes did not have identifiable metadata). When this metadata is plotted on the phylogenetic tree, we do not observe notable representation of blood or sputum samples in the clade that includes both C1 and C2 samples (Fig. S5). Using Scoary, we identified 24 accessory genes that have a Bonferroni p-value less than 0.05 and odds ratio greater than 2 associated with blood vs. sputum (Table 3). Compared to the H. influenzae 4,842 pangenome set, representative C1 and C2 genomes were enriched in the presence of these accessory genes associated with blood infections, having 16 and 14 genes (of the 24 identified), respectively (Fig. 4).

After identifying this enrichment of potential bloodstream infection genes, we wanted to determine whether these enrichment patterns were unique to NTHi or were reflected in the larger Haemophilus genus. Of the 24 genes that were identified by Scoary as potentially discriminatory for bloodstream infection, only two were known virulence factors in the Haemophilus genus according to the VFDB. The lsgB gene, found in both C1 and C2 samples, is associated with H. parasuis virulence through its involvement in lipooligosaccharide biosynthesis sialylation (Wang et al., 2018). It is one of several virulence genes in H. parasuis, contributing to the bacterium’s pathogenicity by influencing sialylated lipooligosaccharide production. Another of the 24, an igaA1 gene, found in C1 but not C2 genomes, is a homolog of Salmonella membrane protein IgaA (Tierrez & García-del Portillo, 2004). IgA regulates bacterial regulons like RcsC-YojN-RcsB and PhoP-PhoQ, with the igaA1 allele (due to an R188H mutation) altering the expression of PhoP-PhoQ-activated (pag) genes, such as ugd, which is linked to lipopolysaccharide modification and colanic acid capsule synthesis (Tierrez & García-del Portillo, 2004).

Rates of recombination were similar in C1 and C2 to the rest of the H. influenzae species

Previous work in Neisseria meningitidis has shown that homologous recombination events can contribute to a novel pathogenicity phenotype (Tzeng et al., 2017). We hypothesized that allelic variation introduced by homologous recombination may have played a role in the pathogen evolution in the C1 and C2 clusters. We created a core genome alignment of 50 randomly chosen H. influenzae STs and the C1-1 and C2-1 genomes to predict potential regions of recombination using ClonalFrameML (Kwong, 2017) (Fig. 5A). As found in previous studies (Carrera-Salinas et al., 2021; Gonzalez-Diaz et al., 2022) recombination was common across H. influenzae genomes. There were 44 recombination events in the C1-1 genome with a rho/theta of 1.6, (rho/theta is the frequency of recombination events relative to the mutation rate). There were no recombination events identified in the C2-1 genome (rho/theta = 0.58). When we looked at the rho/theta values of the other H. influenzae samples we found that both C1-1 and C2-1 were within the ranges of other clades (Fig. 5B), suggesting that patterns of recombination were not atypical for the species.