F3031 (GenBank accession no. FQ670178) is a BPF clone strain that is indistinguishable from other isolates by various typing systems, including multilocus sequence typing. F3047 (GenBank accession no. FQ670204) is a conjunctivitis isolate from Brazil that was established by typing to be unrelated to the BPF clone. F3031 and F3047 are described in more detail elsewhere (21).

Bacterial genomes were sequenced at the Wellcome Trust Sanger Institute, Cambridge, UK. The first drafts of the F3031 and F3047 genomes were assembled from sequence to ≈7-fold coverage, from pOTWI2 and pMAQ1Sac_BstXI genomic shotgun libraries, by using BigDye Terminator chemistry on an Applied Biosystems 3730 DNA Analyzer (Applied Biosystems, Foster City, CA, USA). End sequences from large insert fosmid libraries in pCC1FOS (insert size 38–42 kb) were used as a scaffold for each strain. Further sequencing was performed on the Illumina Genome Analyzer (Illumina, Inc., San Diego, CA, USA). Assemblies were created and gaps and repeat regions were bridged by read pairs and end-sequenced PCR products.

Coding sequences were predicted by using Glimmer 3 (www.cbcb.umd.edu/software/glimmer). Automated annotation by similarity was done by searching the Glimmer 3 coding sequence set against the National Center for Biotechnology Information Clusters of Orthologous Groups database and the SwissProt dataset (www.uniprot.org). Annotation by similarity was done by importing the NTHI strain 86–028NP annotation and comparing it with the F3031 coding sequence set by using reciprocal FASTA (www.ebi.ac.uk/Tools/sss/fasta). Automated annotation was confirmed by manual curation with the Artemis genome visualization tool (22). Gene definitions and functional classes were added manually by using FASTA analyses of the primary automated comparisons. tRNA genes were predicted by using tRNAScan-SE version 1.2 (23). Identification of the rRNA operons was based on similarity to homologs in the NTHI strain 86–028NP genome.

Generation of pairwise comparisons of complete genome sequences was based on alignment of basepairs in MAUVE (24), which enabled alignment of whole genome sequences despite rearrangements. For each pairwise comparison of whole-genome sequences, the length of the alignment between the 2 strains was calculated and a distance matrix was created. The distance matrix, based on the lengths of the sequence alignments, was used to create a heat map showing the clustering of strains.

Evolutionary relationships between protein-coding sequences from different strains were inferred by using MEGA version 5.02 (25). Phylogenetic trees were constructed by using sequence alignments, and a neighbor-joining tree was built under a Poisson correction substitution model assuming uniform rates of substitutions among sites.

Genome sizes ranged from 1.83 to 2.0 Mb (Table 1). The F3031 genome comprises 1,985,832 bp, is 8% larger than Rd KW20, and encodes 1,892 genes. The F3047 genome is larger (2,007,018 bp) and encodes 1,896 genes. All strains have a genome G+C content of 38%, typical of H. influenzae. HaeBPF strain F3031 contains an ≈24-MDa plasmid, previously sequenced and annotated (26), with average G+C content of 36.7%. This plasmid has been excluded from analysis.

Whole-genome alignment of Hae strains F3031 and F3047 revealed substantial colinearity with 1 major rearrangement and 3 small inversions (Figure 1). Pairwise nucleotide alignments of the 7 sequences indicated a closer relatedness of the 2 Hae strains to each other than to the 5 other H. influenzae genomes (Figure 2). A core genome of 77% was shared across all 7 strains.

F3031 shares 10.6% of its genomic sequence with 1 other strain and 88% of this shared sequence (9.3% of total) with F3047, emphasizing the closer relatedness of these 2 strains to each other than to the other H. influenzae strains. A total of 163 predicted coding sequences lie within this Hae-specific DNA. A total of 99 (61%) coding sequences lie within regions of previously characterized Haemophilus bacteriophages, encoding proteins inferred by the similarity of their deduced sequences to be phage components associated with coexpressed genes transported by the phage (phage cargo). These proteins are are either homologs of conserved hypothetical proteins in other organisms or previously unidentified proteins of unknown function. Of all Hae-specific genes, >22% encoded homologs of products identified elsewhere as being involved in host–pathogen interactions; prominent members were putative adhesins and invasins not previously found in strains of H. influenzae (Figure 3). Description of the Hae accessory genome will focus on these putative adhesins.

These new Hae-specific adhesins include 4 novel fimbrial operons, unique high-molecular-weight (HMW) proteins, and a 10-member family of trimeric autotransporter adhesins (TAAs). Many of these coding sequences are associated with simple sequence repeats (SSRs), indicating that phase variation may confer the potential for antigenic variation and immune response evasion during infection.

The presence of duplicated hafABCDE operons (27) was confirmed in the F3031 and F3047 genomes. We also identified 4 more Hae-specific fimbrial gene clusters, aef1–aef4 (Figure 4). Clusters aef1–aef3 were present in both strains, although not identical (55%–100% similarity on gene-by-gene comparison), but aef4 was not present in HaeBPF F3031. Each aef operon encodes 4–6 proteins and has modest sequence identity to products of corresponding haf genes (38%–57%) and to F17 fimbrial adhesins (25%–64%) produced by pathogenic Escherichia coli associated with septicemic diarrheal diseases. Three clusters (aef1, aef3, aef4) are associated with mononucleotide SSRs of 10–17 nt located in the putative promoter region upstream of the aefA gene (Figure 4), conferring capacity for phase-variable expression through expansion and contraction of the SSR, altering efficiency of promoter binding.

The Hae genomes each encode a much richer repertoire of autotransporter adhesins than is found in other sequenced Haemophilus spp. Monomeric (classical) and novel trimeric autotransporter adhesins are present (MAA and TAA, respectively). Of the established Haemophilus autotransporter adhesins, the MAA Hap (Haemophilus adhesion and penetration protein), widely distributed in H. influenzae and proposed as a candidate NTHI vaccine antigen, is present as a pseudogene in F3031 and F3047, as previously reported by Kilian et al. (28). IgA1 protease, previously identified in the BPF clone, is also present in conjunctivitis strain F3047. Sequence alignment to other H. influenzae demonstrated that IgA1 from F3047 is more closely related to IgA1 from Rd KW20 (88% aa identity) than from F3031 (65% aa identity). Homologs of the HMW adhesins HMW1 and HMW2 (MAAs) and of the TAA H. influenzae adhesin Hia are present in F3031 and F3047. In contrast to the many NTHI strains for which substantial sequence information is available, where HMW1/HMW2 or Hia have almost always been alternatives, both are found in these Hae strains. HMW1 and HMW2, encoded at loci each consisting of 3 genes (hmwABC), were first identified in NTHI strain R2846 as HMW surface-exposed proteins, mediating attachment to human epithelial cells (29). More than 75% of NTHI encode HMW1 and HMW2, present at the same chromosomal locations in almost all HMW-containing NTHI isolates examined. hmw1A and hmw2A encode adhesins with different receptor binding specificities resulting from domains in variable regions comprising amino acid residues 114–237 of mature Hmw1A and 112–236 of mature Hmw2A (30). Despite conservation in binding specificity, hmw1A or hmw2A alleles from different isolates are highly polymorphic in the receptor binding domains (30). Hae Hmw1A– and Hmw2A–binding domain sequences (deduced from comparison with R2846 sequence) were aligned by using ClustalW (www.ebi.ac.uk/Tools/msa/clustalw2) with those from homologs in other NTHIs, regardless whether they were Hmw1A or Hmw2A, and the alignment was used to construct a phylogenetic tree (Figure 5). The Hae HmwA–binding domains are distinct from those in other NTHIs, suggesting that in Hae these proteins have diverged separately from other NTHIs.

In Hae F3031 and F3047, the HMW clusters are not at the homologous NTHI chromosomal position; they are elsewhere, with a 22-kb bacteriophage insertion directly downstream of hmw2ABC. The hmwA alleles are further differentiated from those found in other NTHI strains by their associated SSRs. The putative promoter region upstream of the hmw1A and hmw2A homologs contains the octanucleotide repeat unit 5′-GCATCATC-3′; there are 14 and 15 copies, respectively, in F3031 and 13 and 12 copies, respectively, in F3047. This repeat pattern contrasts with all hmwA genes so far sequenced in different NTHI strains, in which 7 basepair SSRs of either 5′-ATCTTTC-3′ or 5′-TGAAAGA-3′ in varying copy numbers are located upstream of the genes (31,32).

The Hae accessory genome includes a 10-member gene family that encodes proteins with the sequence characteristics of TAAs (Table 2). These TAAs are distinct from the Haemophilus TAAs Hsf (33), Hia (34), or the recently described Cha (35). In strains F3031 and F3047, a total of 8 genes (1–6,8,9) are present as homologs, termed tabA (for the HaeBPF trimeric autotransporter [bpf] alleles) or tahA (for the regular hae [conjunctivitis] alleles). tahA7 has no homlog in HaeBPF. The tenth gene, tabA10, is the recently described adhesin/invasin gene hadA (36). This gene is found only in HaeBPF F3031; F3047 has no corresponding gene. Each gene appears to be locus specific, sharing the same flanking regions, but sequences differ substantially between homologs 1 and 2 in particular. tabA4/tahA4, tabA5 and tabA9/tahA9 seem to be pseudogenes, carrying frameshift mutations within the coding sequence. All TAAs except tabA8/tahA8 and tabA10 (hadA) are associated with SSRs located either within the coding sequence or upstream in the putative promoter region, indicating that expression may be modulated by phase variation.

All these TAAs share the characteristic 3-domain structure of N-terminal signal peptide and C-terminal outer membrane translocator domain, separated by an internal passenger domain. However, comparison of orthologous TAAs revealed striking differences between their passenger domains for TabA1/TahA1 and TabA2/TahA2, suggesting different functions of these proteins in the 2 strains (Figure 6). The passenger domains of these proteins vary in the number of binding domains (hemagglutinin and Hep_Hag domains) and in possession of different-sized, low-complexity spacer regions consisting of approximate heptapeptide repeats. TabA1 from F3031 contains 90 copies of tandemly duplicated AASSSAS with occasional T, N, or other substitutions in many copies; TahA1 from F3047 contains 48 copies of tandemly duplicated AETAKAG with occasional R, V, or other substitutions in many copies. In the prototypic TAA YadA, a series of 15-residue repeats appears to have such a spacer function between the protein head and its anchor in the outer membrane (37), holding any receptor-binding domains away from the bacterial cell surface.

In the context of the unusual virulence of the HaeBPF clone, the tabA1 locus is particularly intriguing. Comparison with the tahA1 locus indicates not only the substantial difference between the genes themselves, in the sequence encoding the putative stalk domain, but also (in F3031) an additional gene, HIBPF06250, encoding a conserved hypothetical protein, homologous to an uncharacterized gene product in the Haemophilus cryptic genospecies strain 1595 (35). In this strain, the gene (tandem duplicated) lies downstream of the TAA Cha. In F3031, HIBPF06250 is interposed between tabA1 and IS1016 (Figure 7), and the gene (like the insertion sequence) is absent in F3047. Association of IS1016, first described as the Haemophilus capsulation locus–associated insertion sequence, with unusual and invasive virulence of NTHI strains has been suggested elsewhere (17,18), although no specific gene association has been identified.

The part of the Hae accessory genome unique to HaeBPF amounted to 102,304 bases (5.2% of its genome). Ten HaeBPF-specific loci ranged in size from 370 to 20,002 bases and in G+C content from 27.9% to 44.5%. Deviation from the Haemophilus average of 38% suggests that these are more recently acquired regions. Much of this DNA is located within 5 bacteriophage domains, containing all 219 coding sequences (12 Hae specific, 11 HaeBPF specific) (Table 3) and including 1 (phage region 1) now termed HP3, similar in size and gene content to Haemophilus bacteriophage HP2, found in NTHI strains associated with unusual virulence (38). The HaeBPF-specific accessory genome comprises these and another 10 coding sequences (Table 4), which remained apparently BPF specific after BLASTP analysis of their deduced amino acid sequences against the nondegenerate public databases (October 2011), which include many more Haemophilus sequences from incomplete genome sequencing projects (19). The nearest matches to these sequences were mainly homologs in other pathogenic bacterial species that occupy the same ecologic niche. One gene (hadA at BPF-specific locus 10) has recently been characterized as encoding an epithelial adhesin/invasin plausibly contributing to HaeBPF virulence (36), but the function of the others, and any part their products may play in the serum resistance of the HaeBPF clone that endows it with pathogenic potential, remains to be established. Eleven genes appear to be phage cargo (Table 3); these are either homologs of conserved hypothetical proteins identified in other organisms or entirely unknown and might represent novel virulence factors. Four F3031-specific gene products do not have homologs in any other bacterial species and cannot be assigned a putative function. Novel genes have generally formed a much larger part of newly sequenced bacterial genomes, and identification of so few unknown genes in HaeBPF strain F3031 reflects the current availability of a large amount of Haemophilus sequence data, in particular from strains of NTHI.