The 5 CU strains used in this study were the only strains available at the time the study was initiated (Table 1); their associated clinical features are listed in S1 Table. The class I and class II strains used in this study were chosen because these strains had been previously analyzed by multilocus sequencing (Table 1) [10]. 35000HP, whose genome has been sequenced (GenBank accession no. NC_002940.2), was used as the reference strain in this study; 35000HP was isolated from a volunteer who was experimentally infected on the arm with strain 35000 and has been extensively characterized in human inoculation experiments [12, 13]. The H. ducreyi strains were grown on Columbia agar plates or in Columbia broth supplemented with 1% bovine hemoglobin (Sigma-Aldrich), 1% IsoVitaleX, and 5% fetal bovine serum (Hyclone) at 33°C with 5% CO₂.

Genomic DNA was extracted from H. ducreyi strains using the DNeasy Blood & Tissue kit (Qiagen) and quantified using the Quant-It High Sensitivity dsDNA Assay kit (Life Technologies).

The sequencing libraries were prepared using the NexteraXT DNA Library Preparation kit (Illumina, Inc.) following the manufacturer’s instructions. Samples were multiplexed using the NexteraXT Dual Index Primer kit. Equimolar concentrations of indexed libraries were combined into a single pool and were sequenced at the Tufts University Genomics Core Facility. Paired-end 250-bp sequencing was performed on the Illumina MiSeq platform using the MiSeq V2 500 cycles chemistry. The de novo assembly was performed using Edena, with a customized bash script that optimizes the assembly process by optimizing three key Edena parameters [14]. The assembled contigs were annotated using the RAST online annotation tool [15].

A flow chart of comparative genome analysis of CU and GU strains is depicted in S1 Fig. For all comparative genome analyses in this study, the genome sequence of 35000HP was used as the reference. The de novo assembled contigs were ordered into Locally Collinear Blocks (LCBs) by Mauve Contig Mover (MCM) [16]. The breakpoints between LCBs were resolved by using BLAST analysis of the unaligned contigs produced by MCM, the breakpoint regions in the de novo assembled contigs, and by alignment of raw reads against 35000HP. After resolving the breakpoints, the ordered contigs were concatenated into draft genomes using Emboss 6.3.1. Pairwise genome conservation distances, which represent both gene content and sequence similarity, were estimated from draft genomes using ProgressiveMauve and plotted as heat map using CIMminer [17, 18].

The draft genome sequences for the 14 H. ducreyi strains NZS1, NZS2, NZS3, NZS4, 82–029362, 6644, HD183, HMC46, HMC56, NZV1, 33921, CIP542, DMC64, and DMC111 were deposited in GenBank under the accession numbers CP011218, CP011219, CP011220, CP011221, CP011222, CP011223, CP011224, CP011225, CP011226, CP011227, CP011228, CP011229, CP011230, and CP011231, respectively.

Genome rearrangements were identified from multiple alignments of the draft genomes generated by ProgressiveMauve, BLAST Ring Image Generator, and nucleotide BLAST and from the assembly of raw reads against 35000HP by SeqMan NGen [17, 19, 20]. Because reference-based alignment can miss additional genes that might be absent in the 35000HP genome, the de novo assembled contigs that did not align to the 35000HP genome by ProgressiveMauve were aligned against other microbial genomes using translated nucleotide BLAST.

SNPs and small insertions and deletions (indels, <10 bp) were detected using DNASTAR Lasergene (DNASTAR, Inc., Madison, WI). Briefly, the sequenced reads were assembled by SeqMan NGen against 35000HP. SNPs and indels were discovered by Seqman Pro using default parameters except that a minimum frequency of 90% reads and a minimum coverage of 50 reads were used for the analysis. SNPs were grouped as non-coding, synonymous, or nonsynonymous. Nonsynonymous SNPs were further categorized as substitutions, no-start, no-stop, nonsense, or frameshifts. All SNPs in the genomes of the CU strains were manually verified for accuracy.

Diversity analyses of whole-genome nucleotide sequences and translated concatenated coding sequences were performed using Mega 6.0 [21]. The reliability of the diversity analyses was tested using 1000 bootstrap replicates.

Recombination analysis was performed using the Phi test implemented in PhiPack and the likelihood ratio test implemented in TOPALi v2 [22, 23]. For both tests, a threshold P < 0.05 was used to define a recombination event.

Phylogenetic analyses were performed using Mega 6.0 and Realphy [21, 24]. Briefly, whole-genome alignments were imported into Mega 6.0 and subjected to model testing to identify the best-fit models of nucleotide substitution. Model testing identified Hasegawa-Kishino-Yano plus invariant sites plus gamma-distributed model as the best-fit nucleotide substitution model for our data. Using the best-fit model, phylogenetic analyses were performed with both whole-genome alignments and alignments of translated amino acid sequences from concatenated protein-coding regions using different methods of phylogeny reconstruction, including Maximum Likelihood, Maximum Parsimony, Minimum Evolution, and Neighbor Joining with different gap treatment approaches. We also inferred phylogenies using Realphy, which generates phylogenetic trees by merging alignments obtained by mapping to multiple reference genomes. A rooted Maximum Likelihood tree was reconstructed by including other Pasteurellaceae members (Actinobacillus pleuropneumoniae, Mannheimia haemolytica, Pasteurella multocida, Aggregatibacter actinomycetemcomitans, and Haemophilus influenzae) as outgroups. The reliability of all the trees generated was verified by 1000 bootstrap replicates.

The times to the most recent common ancestor (MRCA) were estimated by Bayesian molecular clock method using Beast v1.8.1 [25]. Hasegawa-Kishino-Yano plus invariant sites plus gamma-distributed model and a relaxed clock model were used to account for variation in substitution rates. The results from the Beast analysis were visualized using Tracer v1.6. A best-fit tree was identified from the tree data generated by Beast using TreeAnnotator and visualized using FigTree v1.4.2. As described previously, we used a substitution rate of 4.5 × 10⁻⁹ per site per year to calibrate the tree [26].

Selection analyses were performed using Mega 6.0 and Hyphy 2.1 [21, 27]. Briefly, protein-coding regions were extracted from the annotated genomes, ordered against 35000HP using MCM, concatenated using Emboss 6.3.1, and aligned using ProgressiveMauve [16, 17, 28]. The alignments were manually edited for accuracy to obtain a codon-delimited alignment, which was used for all the selection analyses. Rates of nonsynonymous (d_N) and synonymous (d_S) substitutions are widely used as a sensitive measure of selection occurring in a protein with d_N = d_S, d_N > d_S, and d_N < d_S indicating neutral, positive, and negative selection, respectively. Alignment-wide evidence for selection was tested using the codon-based Z test. For the codon-based Z test, we first calculated d_N and d_S and their variances using 1000 bootstrap replicates. We then used this information to test the null hypothesis of neutrality (d_N = d_S) versus alternative hypothesis of positive (d_N > d_S) or negative (d_N < d_S) selection using a Z-test. A branch-site random effects likelihood test was used to test whether any of the branches in the tree are evolving under positive selection. A branchTestDNDS test was performed to test whether a prespecified branch of the tree is evolving under different selection strength than the rest of the tree [27]. Individual sites under positive or negative selection were identified using the single likelihood ancestral counting and fixed effects likelihood methods [27].

The draft genomes were interrogated for the presence of genes that are required for the virulence of strain 35000HP in the human inoculation experiments, using nucleotide BLAST [13]. For identifying sequence variation, the nucleotide sequences of virulence genes were translated into amino acids and the translated sequences were aligned using Clustal Omega [29].

The draft genomes were searched for the presence of known antimicrobial resistance genes using ResFinder with default parameters [30]. AST was performed using the agar dilution method as described previously with some modifications [31–33]. Briefly, H. ducreyi strains were grown on Columbia agar (Difco) containing 1% hemoglobin (BBL), 0.2% activated charcoal (Sigma-Aldrich), 5% fetal bovine serum (Atlanta Biologicals), and 1% IsoVitaleX (BBL) for 48 h at 33°C under microaerophilic conditions. The colonies were suspended into Mueller-Hinton (BBL) broth containing 1% IsoVitaleX and 0.002% Tween-80 (Sigma-Aldrich), passed through a 22-gauge needle and left at room temperature for 15 min. The optical density of the culture was adjusted to that of a 0.5 McFarland standard using a Spectronic 20 Plus spectrophotometer (Milton Roy). AST was performed on Mueller-Hinton II medium (BBL) containing 33% lysed horse blood (Remel), 5% fetal bovine serum, and 1% IsoVitaleX. The following antibiotics were tested: amoxicillin (AMX), amoxicillin/clavulanic acid (AMC; 2:1), azithromycin (AZT), ciprofloxacin (CIP), ceftriaxone (CRO), doxycycline (DOX), erythromycin (ERY), and penicillin (PEN; all from Sigma-Aldrich). The H. ducreyi strains CIP542, 35000HP, and the H. influenzae strain 49247 were used as controls. A 10⁴/ml suspension of each strain was delivered onto each plate with a Steer’s Replicator (CMI-Promex, Inc.), and the plates were dried for 15 min at room temperature. The minimal inhibitory concentrations were recorded after incubating the plates for 48 h at 33°C under microaerophilic conditions. The presence of three or fewer colonies was recorded as no growth.

Whole-genome sequencing generated between 0.5 and 4 million reads for each of the 14 strains (Table 1). The estimated genome sizes ranged from 1.52 Mb to 1.74 Mb, with an average GC content of 37.8% to 38.6% (Table 1). The total number of contigs for each strain ranged from 40 to 129 (Table 1). The estimated average genome coverage ranged from 92 to 596 fold (Table 1). Contig ordering generated 2 to 6 LCBs for the CU strains, 5–15 LCBs for the class I strains, and 12–16 LCBs for the class II strains. Inspection of the LCBs revealed that the majority of the putative breakpoints between LCBs occurred in genes that share high homology with other genes in the H. ducreyi genome such as lspA1 and lspA2, genes encoding rRNAs, and bacteriophage-related genes and that the majority of the breakpoints did not contain any rearrangements.

Analysis of pairwise genome conservation distance of the draft genomes showed that CU strains form a subcluster within class I strains and that class II strains form a separate cluster from CU and class I strains (Fig 1).

Compared to 35000HP, all the CU strains consistently contained ~20-kb deletion (HD1528 to HD1565) in a bacteriophage locus that is homologous to Pseudomonas aeruginosa bacteriophage B3 and five small deletions that ranged in size from 30–767 bp (Fig 2 and S2 Table). The class I strain HMC56 contained a 50-kb deletion (HD0897 to tRNA-Lys-1) in a region homologous to the H. influenzae ICEHin1056 integrative conjugative element (Fig 2 and S3 Table). All the class II strains contained 3 major deletions of 37 kb (HD0087 to HD0161), 35 kb (HD0478 to HD0495) and 50 kb (HD0897 to tRNA-Lys-1), which are homologous to Escherichia coli bacteriophage D108, Haemophilus bacteriophage SuMu, and H. influenzae ICEHin1056, respectively (Fig 2 and S4 Table). The class II strains also contained several deletions (between HD1528 and HD1618) in a region that is homologous to P. aeruginosa bacteriophage B3 (Fig 2 and S4 Table). All the GU strains also contained several other small deletions as listed in S3 and S4 Tables.

Compared to 35000HP, we did not find any inversions in CU strains with the exception of NZV1, which contained an inversion of ~428 kb that spanned from HD0054 (tuf) to HD0659 (S2 Fig). Among the class I strains, HMC56 contained an inversion of ~300 kb that spanned from glpA (HD1157) to lspA1 (HD1505) (S2 Fig). HD183 contained a ~161 kb inversion that spanned from hhdA (HD1327) to lspA1 (HD1505) (S2 Fig). All the class II strains contained an inversion of ~17 kb that spanned from HD1532 to HD1565 (S2 Fig). However, BLAST analysis of the inversion breakpoints showed no major changes in their genetic content.

Compared to 35000HP, the CU strains, the class I strain 82–029362, and the class II strain CIP542 contained no additional genes in their genomes. All the remaining class I and class II strains contained several additional genes as listed in S5 Table.

To get a deeper understanding of the relationship of CU strains to GU strains, we next performed whole-genome SNP analysis using 35000HP as a reference. CU strains differed from 35000HP by ~400 SNPs (Table 2). The class I strain HD183 differed from 35000HP by ~160 SNPs, while all other class I strains differed by ~2,000 SNPs (Table 2). The class II strains differed from 35000HP by ~30,000 SNPs (Table 2).

Analysis of within lineage genetic diversity showed that CU strains had the least nucleotide and amino acid divergence followed by class I and class II strains (Table 3). Analysis of interlineage diversity showed that there was little divergence between CU and class I strains; however, a greater amount of divergence was observed between CU and class II strains (Table 3). Interlineage diversity analysis showed that there was high divergence between class I and class II strains (Table 3).

While the Phi test showed no evidence of recombination (P = 0.44), the likelihood ratio test identified five putative recombination events in CU and GU strains (S6 Table). Removal of recombination regions had no major effect on the overall topology of the phylogenetic tree described in the following section, except for minor differences in bootstrap values and positioning of individual species within class clades (S3 Fig).

In general, all methods showed that class I and class II strains formed two separate phylogenetic clusters and that CU strains formed a subcluster within the class I clade with minor differences in bootstrap values and positioning of individual species within class clades. A rooted tree generated by the Maximum Likelihood method and Pasteurellaceae members as outgroups was used as the final tree (Fig 3).

To determine the approximate time to the MRCA of the CU strains, we performed a molecular clock analysis using the Bayesian method and the mutation rates proposed by Ochman et al. for calibration [25, 26]. The divergence time of the CU strains from the MRCA of the class I strains 35000HP and HD183 was estimated as 180,000 years ago (Fig 4). The divergence time of the CU strains, 35000HP, and HD183 from the MRCA of other class I strains was estimated as 450,000 years ago (Fig 4). The divergence time of class I strains from the MRCA of class II strains was estimated as 1.95 mya (Fig 4). Molecular clock analysis also showed that the CU strains began to diversify from each other around 27,000 years ago (Fig 4). Thus, CU strains appear to have recently diverged from class I GU strains.

Pairwise analysis of rates of nonsynonymous (d_N) and synonymous (d_S) substitutions and their variances showed that the Z-test rejected the null hypothesis of neutrality (d_N = d_S) in favor of the alternative hypothesis of negative selection (d_N < d_S) (Table 4). The d_N-d_S value averaging over all sequence pairs was -75.55 (P = 0.0000000001). Utilizing the rates of nonsynonymous (d_N) and synonymous (d_S) substitutions, we also calculated the overall mean and pairwise mean d_N/d_S ratios; the overall mean d_N/d_S ratio for all genomes was 0.31 and the pairwise mean d_N/d_S ratios for most comparisons were less than 1 (Table 4). The pairwise mean d_N/d_S ratio between the CU and GU lineages was 0.35, between CU and class I lineages was 0.38, and between CU and class II lineages was 0.33. Consistent with these analyses, the single likelihood ancestral counting and the fixed effects likelihood analyses identified 141 and 132 negatively selected sites, respectively.

To determine whether CU strains evolved under different selection strength than GU strains, we performed a TestBranchDNDS analysis. This analysis showed that the strength of selection in CU strains was not significantly different than in GU strains (likelihood ratio difference = 3.8; P = 0.58).

We determined whether the genomes of CU strains contained the genes that are required for the virulence of strain 35000HP in the human challenge model of infection and whether there were variations in these virulence determinants compared to GU strains [13]. BLAST analysis showed that all the CU and GU strains contained all of the genes known to be required for virulence in the human challenge model (S7 Table). Alignment of amino acid sequences of the virulence determinants showed that the DsrA, LspA1, and LspA2 proteins of the CU strains differed by at least 1 amino acid from class I strains (S7 Table).

To determine whether CU strains were resistant to clinically relevant antimicrobials, we performed AST using the agar dilution method. The CU strains from Samoa and Vanuatu were AZT susceptible, and had similar susceptibility patterns as the type strains 35000HP and CIP542 (Table 5). With the exception of 82–029362 and 35000HP, all the class I strains were resistant to penicillin (MIC, >256 μg/ml), amoxicillin (MIC, 64–256 μg/ml), and doxycycline (MIC, 8–16 μg/ml) (Table 5). With the exception of CIP542, all class II strains were resistant to amoxicillin and penicillin (MIC, 128–256 μg/ml) (Table 5). The class II strains 33921 and DMC64 were also resistant to doxycycline (MIC, 8–16 μg/ml) (Table 5). All the strains were susceptible to ciprofloxacin, azithromycin, erythromycin and ceftriaxone (Table 5).

Consistent with their susceptibility to clinically relevant antimicrobials, the CU strains contained no horizontally acquired genes encoding antimicrobial resistance determinants in their genomes. Consistent with their resistance to penicillin/amoxicillin and doxycycline, the genomes of GU strains contained genes that confer resistance to penicillin/amoxicillin (blaTEM-1B) and doxycycline [tet(B), tet(32) or tet(M)] (Table 5).