|Home | About | Journals | Submit | Contact Us | Français|
The emergence of multidrug-resistant (MDR) typhoid is a major global health threat affecting many countries where the disease is endemic. Here whole-genome sequence analysis of 1,832 Salmonella enterica serovar Typhi (S. Typhi) identifies a single dominant MDR lineage, H58, that has emerged and spread throughout Asia and Africa over the last 30 years. Our analysis identifies numerous transmissions of H58, including multiple transfers from Asia to Africa and an ongoing, unrecognized MDR epidemic within Africa itself. Notably, our analysis indicates that H58 lineages are displacing antibiotic-sensitive isolates, transforming the global population structure of this pathogen. H58 isolates can harbor a complex MDR element residing either on transmissible IncHI1 plasmids or within multiple chromosomal integration sites. We also identify new mutations that define the H58 lineage. This phylogeographical analysis provides a framework to facilitate global management of MDR typhoid and is applicable to similar MDR lineages emerging in other bacterial species.
S. Typhi, the primary global cause of human typhoid (enteric fever), is a monophyletic serovar of S. enterica. Unlike many Salmonella, S. Typhi are highly restricted to infection of humans and are associated with systemic infection, prolonged fever and an asymptomatic carrier state1. Typhoid is still a common disease in many regions of the world with poor infrastructure and limited economic development and is also a risk for travelers who visit such regions2. It is estimated that 20–30 million cases of typhoid occur annually, although deaths are less frequently reported than before the availability of effective antimicrobials3,4.
In addition to improvements in access to clean water and sanitation, typhoid can potentially be controlled by other interventions such as vaccination5–7 and antimicrobial therapy8. Chloramphenicol, ampicillin and trimethoprim-sulfamethoxazole were traditional first-line drugs commonly used to treat acute typhoid, and these agents continue to be used in areas of the world where S. Typhi are deemed susceptible. However, since the 1970s, S. Typhi have emerged that display multidrug resistance, defined as resistance to the above antimicrobials, compromising treatment9–11. Since the 1990s, alternative treatment options have included fluoroquinolones, third-generation cephalosporins (such as ceftriaxone) and the azalide azithromycin1. The early emergence of MDR S. Typhi was driven in large part by the acquisition of IncHI1 plasmids carrying antibiotic resistance genes12 and, more recently, by chromosomal mutations associated with resistance to fluoroquinolones, and MDR strains have been reported across Asia and Africa13–16.
Phylogenetic analysis, initially based on subgenomic DNA sequences but later on whole-genome DNA sequences, showed that the global S. Typhi population is highly clonal and likely originated from a common ancestor that moved into the human population several thousand years ago17–19. It also indicated that the population is relatively small and that recombination between S. Typhi and other Salmonellae is rare12,19,20. Simple SNP-based typing schemes have been developed that stratify the S. Typhi population into haplotypes, and these schemes are now used to unequivocally map new isolates to the phylogeny17,19,21,22. Notably, this approach identified a single emerging, highly clonal MDR haplotype of S. Typhi, H58, which is being reported with increasing frequency from many countries in Africa and Asia12,17,19,23. Within the H58 lineage, IncHI1 MDR plasmids of the restricted subtype PST6 (ref. 23) and chromosomal point mutations conferring quinolone resistance are common14,24–26. However, relatively little is known about the emergence and evolutionary history of the H58 lineage or how it is moving across endemic regions. Here we have used phylogenetic analysis based on the whole-genome sequences of a global collection of S. Typhi from 63 countries to investigate the genomic architecture of this highly successful S. Typhi lineage.
Of a global collection of 1,832 sequenced S. Typhi (listed in Supplementary Tables 1 and 2), 853 (47%) belonged to haplotype H58, initially defined by the SNP glpA-C1047T (position 2,348,902 in S. Typhi CT18, BiP33 in ref. 17). The earliest H58 isolates in our collection were from 1992 (Fiji) and 1993 (Fiji and Vietnam), and H58 isolates were represented every year from 1992 to 2013, at a mean rate of 40% per year (Fig. 1a). H58 isolates formed a tight cluster within the whole-genome maximum-likelihood phylogeny (Fig. 1b), forming a unique lineage separated by 151 SNPs from the nearest neighboring non-H58 cluster, which consisted exclusively of isolates from Fiji (n = 140) and Tonga (n = 2). Individual H58 isolates differed from the most recent common ancestor (MRCA) of the H58 lineage by a median of just two SNPs, and the median distance between pairs of H58 isolates was six SNPs, strongly indicative of recent clonal expansion. Nearly all of the H58 isolates (93%; 797/853) had ≤5 isolate-specific SNPs, consistent with frequent transmission relative to substitution mutations. This finding was in contrast to that for the rest of the S. Typhi tree, which included a wide diversity of isolates (Fig. 1b), among which only 66% (642/979) had ≤5 isolate-specific SNPs (P < 0.0001, Fisher’s exact test) (Supplementary Fig. 1).
The population structure within H58 was consistent with our previous work defining two major sublineages of H58 (I and II)22,25,27 but provided much greater resolution of substructure with a strong phylogeographical signal (Fig. 2). Our H58 S. Typhi isolates were collected from 21 countries across Asia, Africa and Oceania. We observed strong phylogeographical clustering within 13 countries (Fig. 2), indicating transmission of H58 within these locations. Although our sample spans distinct time periods in different locations, in most cases, the localized subclades were isolated from the same country over ≥4 years, indicating the establishment of long-term local reservoirs (Fig. 3a,b).
These data demonstrate that H58 is now widely disseminated across distinct geographical areas, and the phylogeny provides several insights into the spatial patterns of its spread. There were numerous instances of very closely related isolates from different countries (Fig. 2), which indicate likely transfer events or regional outbreaks and identify routes for geographical dissemination (Fig. 3a,b). Maximum-likelihood analysis of inter-region transfers based on the maximum-likelihood phylogeny and locations where isolates were collected highlighted several candidate intercontinental transfers (Fig. 4). These data suggest that South Asia was an early hub for H58, from which it was propagated to many locations around the world, including countries in Southeast Asia, western Asia and East Africa, as well as Fiji (Figs. 2–4). Most of the diversity in lineage II was present among Indian isolates, with unique local subclusters detected in neighboring countries (Nepal and Pakistan) and in Africa, indicative of occasional transfers out of Asia. In contrast, lineage I was associated mainly with Southeast Asia (Vietnam, Cambodia and Laos), with evidence of transmission to Thailand, Pakistan, Fiji and Africa.
There have been sporadic reports of the emergence of S. Typhi in Africa14,28–30. Indeed, H58 isolates were predominant among the eastern and southern African S. Typhi isolates (63%) (Supplementary Fig. 2); in contrast, the H58 lineage was relatively rare in northern, western and central Africa. H58 lineages I and II were detected in Kenya, Tanzania, Malawi and South Africa, providing compelling evidence for multiple introductions of H58 S. Typhi from South Asia into the continent (Figs. 2–4). In addition, we uncovered evidence of an unreported recent wave of transmission of H58, based on 138 isolates, from Kenya to Tanzania and on to Malawi and South Africa (Figs. 2 and and3).3). These isolates differed from one another by an average of 10 (range of 0–30) SNPs, consistent with a recent clonal expansion. Therefore, this analysis demonstrates an ongoing epidemic of H58 typhoid across countries in eastern and southern Africa.
Estimating mutation rates and divergence dates within the S. Typhi population has been challenging, as S. Typhi is known to establish persistent asymptomatic carriage, during which time it likely evolves at a different rate than during acute infection, disrupting the molecular clock17. Here a temporal signal was barely detectable across the full S. Typhi maximum-likelihood tree, assessed via linear regression of root-to-tip branch lengths on the basis of year of isolation (correlation coefficient (R) = 0.09 (95% confidence interval (CI) = 0.04–0.13); P = 0.0002, Fisher’s exact test). However, a moderate signal was evident within the H58 subtree (R = 0.60 (95% CI = 0.56–0.64); P < 1 × 10−6, Fisher’s exact test; Supplementary Fig. 2). This temporal signal was entirely destroyed by randomization of isolation dates (mean R = 0.01), indicating that uneven sampling of the H58 lineage across space and time was not solely responsible for the observed association. We propose that a temporal signal was detectable within the H58 tree because these data capture epidemic spread over a relatively short time span (2–3 decades), whereas the wider S. Typhi tree represents much more variable population dynamics over thousands of years of evolution, including recent periods of endemic transmission, that differ greatly from the clonal expansion of H58.
We therefore proceeded to estimate the divergence date of the H58 lineage via Bayesian phylodemographical modeling of the H58 population, implemented in Bayesian Evolutionary Analysis Sampling Trees (BEAST)31. To limit potential bias due to highly variable sampling intensities in different geographical locations, we performed BEAST analyses on several cross-sections of 114 H58 isolates (13% of the total), each sampled from 21 countries and spanning the years 1992 to 2013. The combined estimate for the median substitution rate within the H58 population was 1.42 × 10−7 substitutions per site per year (95% highest posterior density (HPD) = 1.0 × 10−7 to 1.8 × 10−7), equivalent to the accumulation of 0.63 SNPs per genome per year (95% HPD = 0.59 to 0.67). The analyses predicted that the MRCA of all extant H58 strains existed ~25 years ago (median calendar year for divergence, 1989; 95% HPD = 1985–1992) and that the effective population size of H58 increased dramatically after 1993 (Supplementary Fig. 2). Although the BEAST analysis was limited because of the moderate strength of the temporal signal, these results are consistent with the very low numbers of SNPs within the H58 lineage, the clear evidence of clonal expansion from the maximum-likelihood tree (Figs. 1b and 3a,b) and epidemiological data reporting increasing rates of multidrug resistance in Asia in the early 1990s (refs. 13,32).
The H58 lineage is associated with high levels of multidrug resistance and reduced susceptibility to fluoroquinolones12,14,25. Acquired resistance genes were identified in 671 of the 1,832 S. Typhi isolates, including 68% of the H58 isolates in comparison to just 9% of non-H58 isolates (P < 1 × 10−16, Fisher’s exact test). In our collection, which included 15 countries with ≥2 MDR isolates, H58 was significantly associated (P < 0.01) with multidrug resistance in nearly all these locations, with the exception of central and western Africa, where multidrug resistance was detected but H58 was not (Supplementary Fig. 3). The most common resistance genes detected were blaTEM-1 (ampicillin resistance), dfrA7, sul1 and sul2 (resistance to trimethoprim and sulfonamides, respectively, and to trimethoprim-sulfamethoxazole collectively), catA1 (chloramphenicol resistance) and strAB (streptomycin resistance). These genes were each found in >540 H58 isolates, including in 525 isolates that carried all 7 genes. These genes are encoded within a Tn2670-like complex transposable element comprising transposon Tn6029, which carries blaTEM-1, strAB and sul2, inserted into transposon Tn2670, which itself comprises Tn21 carrying a class I integron (including sul1, with dfrA7 in the gene cassette) inserted into Tn9 carrying catA1 (refs. 12,33) (Fig. 5). In addition, 405 H58 isolates harbored the tetB gene located in Tn10. Other acquired resistance genes were rare, identified in <1% of the H58 isolates.
Previous reports linked multidrug resistance in H58 S. Typhi to the Tn2670-like element described above encoded on IncHI1 plasmids of the PST6 genotype12,23. Here we identified IncHI1-PST6 plasmids in 74% of the H58 isolates harboring the MDR element, including isolates from Southeast Asia, East Africa and South Asia (Supplementary Table 3), indicating intercontinental transmission of the IncHI1-PST6 MDR plasmid with its H58 S. Typhi host (Fig. 6). However, the remaining MDR H58 isolates lacked IncHI1 plasmid sequences, indicating that the resistance-conferring genes must be located elsewhere in the genomes of these 139 isolates. We screened all H58 isolates for known plasmid replicons and identified non-IncHI1 MDR plasmids in 23 isolates but not among isolates carrying the Tn2670-like MDR element (Supplementary Table 3).
Nearly all H58 isolates carried a copy of insertion sequence IS1 between the chromosomal genes STY3618 and STY3619, near cyaA (Fig. 6). We propose that this insertion was acquired early in the H58 lineage, originating from the IncHI1-PST6 plasmid. Consistent with this hypothesis is the present analysis showing that IS1 was absent from only a few isolates that were basal in the tree (Fig. 6). Long-read sequencing of 2012 Indian isolate ERL12960 (Supplementary Table 4) confirmed that the entire MDR locus, flanked by copies of IS1, was integrated at this location near cyaA (Fig. 5), in agreement with other recent reports15,16. The nucleotide sequence was identical to that in the IncHI1-PST6 plasmid of H58 isolate 10425_1_48_Viety3-193_1997, which we also sequenced for comparison. Further in silico analysis of our data identified new insertion sequences in MDR isolates that lacked the IncHI1 plasmid, including (i) 25 phylogenetically related isolates (mostly from Bangladesh, Pakistan and Iraq) with IS1 in the yidA gene, (ii) 9 related isolates from Fiji with IS1 in STY4438 and (iii) 1 isolate from India with IS1 in the fbp gene (Fig. 6). Long-read sequencing of the 2012 isolate 12148 from India confirmed integration of the MDR locus in the yidA gene (Fig. 5); we confirmed integration at the other two sites by PCR. The distribution of isolates in the H58 tree indicates single integration events at each of the yidA, STY4438 and fbp loci but numerous independent integrations at the cyaA site (Fig. 6). The latter suggests that the Tn2670-like element, which is flanked by IS1 sequences, may target existing IS1 sequences during its mobilization.
In addition, we found that four isolates harbored genes associated with azithromycin resistance. These isolates were 10593_2_14_Alg05-8683_2005 from Algeria carrying ereA and three isolates from Indonesia carrying either msrA (10349_1_90_Indo404ty_1983) or msrD (9953_5_22_IndoA340_2010 and 9953_5_48_IndoA377_2010).
The primary targets of the fluoroquinolones are the DNA gyrase subunits (gyrA and gyrB) and the topoisomerase IV components (parC and parE)34. Nonsynonymous mutations in the quinolone resistance–determining regions (QRDR) of each gene can decrease susceptibility to fluoroquinolones such as ciprofloxacin, which is commonly used in the treatment of typhoid35. Here we found that nonsynonymous changes in the QRDR of these four genes were far more common in the H58 isolates (59%) than in other S. Typhi (13%; P < 1 × 10−6, Fisher’s exact test; Supplementary Table 5). The most frequent QRDR mutations were changes in codon 83 of gyrA encoding p.Ser83Phe (45% of H58 isolates) and p.Ser83Tyr (9% of H58 isolates) substitutions (Supplementary Table 5). The distribution of gyrA substitutions within the H58 phylogeny indicates that these mutations have arisen independently on multiple occasions, consistent with our previous observations of convergent evolution17,19 and confirming that this region of gyrA is under strong positive selection (Supplementary Fig. 4). The accumulation of multiple mutations within the gyr and par genes can result in a higher minimum inhibitory concentration (MIC) for fluoroquinolones35. Additionally, we detected multiple mutations in these genes in 199 isolates: 190 H58 isolates (predominantly from Cambodia and India) and 9 non-H58 isolates (mainly from India) (Supplementary Fig. 4 and Supplementary Table 5).
Transmissible fluoroquinolone resistance can occur in Salmonella via plasmid-mediated acquisition of qnr genes36. Here we identified such genes in seven H58 isolates, which carried both gyrA mutations and the qnrS1 gene that confer high-level resistance. The qnrS1 gene was present within a mobile element also containing blaTEM-1, sul2 and catB4, in association with and possibly mobilized by IS26, on an IncFIB(K) plasmid (Fig. 6 and Supplementary Table 3).
The data show some geographical differences in patterns of antimicrobial resistance within the H58 lineage (Supplementary Fig. 5). Multidrug resistance was common among H58 isolates from Southeast Asia in the 1990s, and in recent years gyrA mutations have arisen on this background, resulting in high rates of MDR H58 with reduced susceptibility to fluoroquinolones. This observation likely reflects the therapeutic use of fluoroquinolones to treat typhoid over this period. Although we have few examples of H58 isolates from South Asia before 2000, the pattern is clearly different in this region, with the majority of isolates from recent years almost all harboring gyrA mutations but with low rates of multidrug resistance35. The situation appears to be different yet again in Africa, where the majority of recent isolates (mainly from Malawi, Kenya and South Africa) were identified as MDR but without gyrA mutations, potentially reflecting the continued use of traditional antimicrobial agents.
Because the S. Typhi H58 lineage emerged rapidly over the past 30 years, we searched for any distinctive genetic signatures, other than those for multidrug resistance, that might be facilitating its dissemination. Prophage-like elements are known hotspots for variation within S. Typhi and other S. enterica37; however, H58 genomes shared five of the seven prophage-like elements previously identified in the S. Typhi reference isolate CT18 (haplotype H1), and only rare acquisitions of new phages were found (five phages, affecting 10% of the H58 isolates; Supplementary Fig. 6 and Supplementary Table 6). The SNPs that define the H58 lineage include several nonsynonymous changes in genes associated with pathogenicity, adaptation and chaperones (Supplementary Table 7). The affected genes consist of ssaP, encoding a Salmonella pathogenicity island–2 (SPI-2)-associated protein involved in intracellular survival and persistence38, and the regulatory genes sirA and csrB, which have been implicated in Salmonella virulence39,40. Additionally, H58-associated SNPs also included changes in genes involved in central or intermediary metabolism, including trpE, rlpB, betC, rtn, metH, yfbT, pub, nuoG and iaaA, and genes involved in membranes or structures, for example, lip1, yhdA, kefA, kcsA, yegT, lsrC, yajI and SBOV18161 (hyaE)41. All S. Typhi display substantial genome degradation, via the accumulation of deletions and inactivating mutations within coding sequences, to form pseudogenes; each S. Typhi genome has >200 pseudogenes, reflecting a loss of ~4% of protein-coding capacity19,20. Here we found that all H58 isolates additionally harbored a point mutation in the sptP gene (STY3001) resulting in a premature stop codon at position 185, effectively rendering H58 strains deficient for SptP protein. SptP is an SPI-1 effector protein known to have a role in modulating the host cell actin cytoskeleton via its GAP domain that targets CDC42 and RAC-1 (ref. 42). We previously identified this nonsense mutation in seven sequenced H58 isolates and a second distinct nonsense mutation in the same gene in the H50 isolate E98-3139. Convergent loss-of-function mutations such as this are quite rare in S. Typhi and may reflect a selective advantage for inactivation of this gene20.
Here, we provide the first comprehensive global phylogeographical analysis of the emerging MDR-associated S. Typhi clade known as H58, covering many of the key geographical regions where typhoid remains endemic. This analysis indicates a major ongoing clonal replacement of resident non-H58 S. Typhi haplotypes by this clade and identifies previously unappreciated inter- and intracontinental transmission events. Smaller regional studies performed in different Asian countries and in Kenya have described the emergence of the H58 haplotype at local and country levels12,14,15,17,19,25. However, here we show the true global impact of H58, which is transforming the S. Typhi population structure across the world. Indeed, we show definitively that H58 has expanded dramatically since the early 1990s and that the MDR phenotype of H58 strains is likely influenced by different regional antibiotic usage11. Further, our analyses show an ongoing epidemic of MDR typhoid moving across Africa, potentially driven by this antimicrobial usage. The existence of this epidemic is supported by the available epidemiological data, which include increasing numbers of reports of MDR typhoid in Africa, and by observations from members of this consortium14,28–30.
The H58 lineage has previously been associated with multidrug resistance, which may be a key factor driving its current expansion14,25. Here we show that this association holds across numerous countries in Asia and Africa, such that the majority of the global burden of MDR typhoid can be attributed to the H58 lineage. Intriguingly, our data indicate that multidrug resistance in H58 is tightly linked to the presence of a single Tn2670-like element that was probably first introduced via the IncHI1-PST6 plasmid but has since transferred to the S. Typhi chromosome in numerous distinct integration events, each affecting different sublineages of H58. Such integrations have recently been noted in isolates from Bangladesh (cyaA and yidA sites)15,43 and Zambia (cyaA site)16; however, our data provide important context for these observations, showing that integrations are relatively frequent and have been occurring since the emergence of H58. Integration of the MDR locus into the chromosome may facilitate loss of the large IncHI1 plasmid, thus moderating any potential fitness burden while maintaining the MDR phenotype. A similar phenomenon has been observed in Shigella sonnei, where chromosomal integration of an MDR transposon in the late 1970s appears to be associated with global dissemination of a single successful sublineage44, and in Salmonella that have acquired the Salmonella Genomic Island45. H58 strains also harbor a higher frequency of quinolone resistance–associated mutations than other S. Typhi, potentially owing to enhanced exposure of the large, diversifying and frequently MDR H58 population to fluoroquinolones, coupled with the lack of a fitness cost associated with these mutations46.
These results provide, to our knowledge, the first global, large-scale, genome-based study of an MDR clade of S. Typhi. The global dissemination of H58 requires urgent international attention. Indeed, the arrival of S. Typhi H58 in Africa appears to be transforming the epidemiology of the disease, with MDR outbreaks of typhoid being reported where the disease was previously unappreciated or absent. It will be particularly important to control antimicrobial prescribing practices, including the use of prophylactic antibiotics such as trimethoprim-sulfamethoxazole47 in this region, as such use likely promotes multidrug resistance. This study highlights the need for longstanding routine surveillance to capture epidemics and monitor changes in bacterial populations as a means to facilitate public health measures, such as the use of effective antimicrobials and the introduction of vaccine programs, to reduce the vast and neglected morbidity and mortality caused by typhoid.
URLs. Sprai, http://zombie.cb.k.u-tokyo.ac.jp/sprai/; SMALT, http://www.sanger.ac.uk/resources/software/smalt/; Path-O-Gen, http://tree.bio.ed.ac.uk/software/pathogen/; Velvet Optimizer, http://www.ebi.ac.uk/~zerbino/velvet/; ISmapper, https://github.com/jhawkey/IS_mapper.
A total of 1,832 S. Typhi isolates were included in the study. These were isolated between 1905 and 2013 and originated from 63 countries spanning 6 continents (Asia, Africa, North and South America, Europe, and Australia and Oceania) (Supplementary Table 1). This collection included 14 of the 19 S. Typhi isolates previously sequenced by Holt et al.19, including the 2 isolates with finished reference genomes, CT18 (AL513382) and Ty2 (AE014613) (Supplementary Table 2). Seven other public S. Typhi references sequences were downloaded from public databases and included in the analysis (Supplementary Table 2). Information on the source was available for 1,531 of our S. Typhi isolates. The isolates included in our study were supplied by numerous contributing laboratories and were cultured from a wide range of clinical specimens, including blood (90.9%; 1,391/1,531), stool or rectal swab (8.3%; 127/1,531), urine (0.4%; 6/1,531), pus (0.3%; 5/1,531), gallbladder fluid (0.1%; 2/1,531), pleural fluid (0.1%; 1/1,531) and cerebrospinal fluid (0.1%; 1/1,531). There were data on the age of the patient for 826 S. Typhi isolates, and these comprised 330 (40%) isolates originating from children (0 to ≤16 years old) and 496 (60%) isolates cultured from adults (>16 years old).
Each of the collaborating laboratories used their own individual methodologies for isolation of whole-genomic DNA. For 231 isolates, DNA was prepared using the Wizard Genomic DNA kit (Promega) according to the manufacturer’s instructions. Index-tagged paired-end Illumina sequencing libraries were prepared as previously described49. These were combined into pools each containing 96 uniquely tagged libraries and sequenced on the Illumina HiSeq 2000 or HiSeq 2500 platform according to the manufacturer’s protocols to generate tagged 100-bp paired-end reads.
In addition, the genomes of five H58 isolates were sequenced on the PacBio RS II platform (Pacific Biosciences) for better resolution of the integration of the large composite transposable element into the chromosome (Supplementary Table 4). Genomic DNA (3 μg) was sheared using the HydroShear Plus (Digilab), and a library was prepared using DNA Template Prep Kit 2.0 (Pacific Biosciences), according to the manufacturer’s instructions. Sequencing was performed on SMRT cells with XL polymerase and DNA Sequencing Kit C2 (Pacific Biosciences). De novo assembly was performed with Sprai v0.9.5 (see URLs) and HGAP v2.1.0 (ref. 50) with default parameters. The contigs from Sprai were circularized with a script in the Sprai package when the script detected a significant overlap between the beginning and end of contigs.
For analysis of SNPs, paired-end Illumina reads were mapped to the CT18 reference genome of S. Typhi, including the chromosome and pHCM1 and pHCM2 plasmids41, using SMALT (version 0.7.4) (see URLs) as previously described51,52. Candidate SNPs were identified as previously described19,53, using SAMtools command mpileup -d 1000 -DSugBf ref bam > results.bcf; bcftools view -cg results.bcf54. SNP calls with quality scores above 30 were collated into a single list of variant sites, and the allele at each SNP site in each isolate was determined by reference to the consensus base call for that genome (using SAMtools and removing low-confidence alleles with consensus base quality <50, SNPs contained in <75% of reads and those with mapping quality <30, read depth <4, <2 reads per strand, strand bias P < 0.001, mapping bias P < 0.001 or tail bias P < 0.001). SNPs called in phage regions, repetitive sequences (354 kb; ~7.4% of bases in the CT18 reference chromosome, as defined previously19) or recombinant regions (~180 kb; <4% of the CT18 reference chromosome, identified using an approach described previously49) were excluded, resulting in a final set of 22,145 chromosomal SNPs identified in an alignment of length 4,275,037 bp.
The maximum-likelihood phylogenetic tree shown in Figure 1b was built from the 22,145-SNP alignment of all 1,832 isolates, plus a S. Paratyphi A strain included as an outgroup for tree rooting, using RAxML (version 7.8.6)55 with the generalized time-reversible model and a Gamma distribution to model site-specific rate variation (the GTR+Γ substitution model; GTRGAMMA in RAxML). Support for the maximum-likelihood phylogeny was assessed via 100 bootstrap pseudoanalyses of the alignment data. A maximum-likelihood phylogenetic tree was also inferred separately from the SNP alignment of 853 H58 S. Typhi isolates using the same parameters as above (Fig. 2). The H58 phylogenetic tree was rooted using an S. Typhi isolate from the nearest neighboring cluster of non-H58 isolates (isolate Fij107364). All maximum-likelihood trees were displayed and annotated using iTOL56,57. To simplify visualization of the H58 phylogeny (Fig. 3a), clades containing only isolates from a single country were collapsed manually in R. Geographical transitions were inferred from this collapsed H58 tree using discrete trait transition modeling, implemented in the make.simmap function in the phytools R package. Briefly, each genome was assigned to a geographical region on the basis of country of isolation (or presumed region of inoculation for travel-associated isolates), these regions were treated as discrete tip states on the H58 maximum-likelihood tree and a Markov model for the evolution of this state (i.e., geographical region) was fitted to the tree, as proposed in ref. 58. Entries in the resulting transition matrix were interpreted as likely geographical transfers between regions, drawn as arrows on the map in Figure 4 (directionality inferred from visual confirmation of the phylogeny).
To investigate temporal signal in the maximum-likelihood phylogeny for S. Typhi, we used Path-O-Gen (see URLs) to extract root-to-tip dates and analyzed their linear relationship with year of isolation using R. To assess the robustness of the H58 temporal signal, analysis of the H58 subtree was repeated 100 times with randomly permuted tip dates. The evolutionary dynamics of the H58 lineage were investigated via Bayesian analysis with BEAST (v1.6)31. Initial analyses were conducted on 114 isolates from across the H58 maximum-likelihood tree, covering the full temporal and geographical range of H58. Bias in isolate distribution was reduced by selecting a maximum of eight isolates from each geographical location, with as much temporal diversity as possible within that location. For countries represented by fewer than eight H58 isolates, all isolates at that location were included in the BEAST analysis. The concatenated SNP alignments of these 114 isolates were subjected to multiple BEAST analyses using both constant population size and Bayesian skyline models of changes in population size, in combination with either a strict molecular clock or a relaxed clock (uncorrelated lognormal distribution), to identify the model that best fitted the data44. For the BEAST analysis, the GTR+Γ substitution model was selected, and tip dates were defined as the year of isolation. For all model combinations, 3 independent chains of 100 million generations each were run to ensure convergence, with sampling every 1,000 iterations. The 3 runs were combined with LogCombiner31, following removal of the first 10 million steps from each as burn-in. In all cases, the relaxed, (uncorrelated lognormal) clock model, which allows evolutionary rates to vary among the branches of the tree together with the skyline demographic model, proved a much better fit for the data (Bayes factor > 200).
For the final analyses reported here, ten independent runs were conducted with different samples of isolates in the alignment but using identical substitution (GTR+Γ), clock (uncorrelated, lognormal relaxed) and demographic (Bayesian skyline) models. Each permutation contained a different set of randomly selected isolates from India, Pakistan, Nepal, Cambodia, Vietnam, Tanzania, Kenya, Laos, Malawi, Bangladesh, Iraq and Fiji. The countries of Sri Lanka, Thailand, South Africa, Afghanistan, Indonesia, Myanmar, Lebanon, Palestine and Australia had five or fewer isolates, and all these isolates were therefore included in each of the ten permutation data sets. Each set of isolates was subjected to triplicate BEAST runs. Maximum–clade credibility (MCC) trees were generated using TreeAnnotator31. Estimates reported as median values with 95% HPDs and posterior probability values (PPVs) were used as support for the identification of ancestral nodes and their associated geographical locations. The Bayesian skyline plot was calculated and visualized using Tracer (v1.6), to investigate changes in the effective population size of the H58 lineage over time31. The effective sample sizes (ESSs) of the parameters were estimated to be >200 for all 10 independent runs of the analysis.
The reads for each isolate were assembled de novo using the short-read assembler Velvet59 with parameters optimized with Velvet Optimizer (see URLs) to provide the highest N50 value. Contigs that were less than 300 bp long were excluded from further analysis. The assemblies were constructed and annotated using Prokka60 by the Pathogen Informatics team at the Wellcome Trust Sanger Institute (Cambridge, UK) using an automated pipeline.
The de novo–assembled contig sets were mapped iteratively to the pan-genome reference set (initialized as the concatenation of the S. Typhi CT18 chromosome and plasmids) using MUMmer (nucmer algorithm)61 as previously described19. At each iteration stage i, sequences not aligning to the current pan-genome Pi – 1 set were incorporated into an extended pan-genome, Pi. The final pan-genome P was annotated using both annotation transfer (for S. Typhi reference sequences) and de novo annotation with Prokka60. Paired-end reads were then aligned to the pan-genome using bwa62 with default mapping parameters. SAMtools54 was used to produce a pileup for each aligned read set, and this was used to summarize, for each annotated gene in the pan-genome P, the coverage (percentage of bases covered) and the presence of inactivating mutations (nonsense SNPs or non-triplet insertions and/or deletions (indels) resulting in frameshifts) in each genome. The annotated pan-genome was specifically examined for the acquisition in H58 of phage sequences and plasmids of defined incompatibility groups (see “Plasmid analyses”). New phages identified in the pan-genome study were confirmed using analysis of high-quality de novo–assembled contigs with the web server PHAST (Phage Search Tool)63.
Acquired antimicrobial resistance genes were detected, and their precise alleles were determined, using the mapping-based allele typer SRST2 (ref. 64) together with the ARG-Annot database65. SRST2 was also used to identify mutations in the gyrA, gyrB, parC and parE genes that have been associated with resistance in Gram-negative bacteria (including Salmonella) to quinolones34,35,66–68. Note, however, that our data were not suitable for assessing the presence of azithromycin resistance–related mutations in the 23S rRNA sequences, which have been occasionally reported in other bacterial pathogens69–71. The horizontal transfer of resistance genes associated with a transposon from the IncHI1 plasmid into the chromosome of H58 isolates was analyzed using a combination of the data from the Illumina and PacBio sequencing, inspected using a combination of tools, including Genome Browser Artemis and Artemis Comparison Tool (ACT)72, which allowed comparison of the genome assemblies against finished sequences for IncHI1 plasmids (pAKU1 (ref. 33) and pHCM1 (ref. 41)) and the CT18 chromosome41. IS1 insertion sites were investigated using ISmapper (see URLs), which uses bwa62 to map reads to the IS1 sequence and identify those reads that flank the IS (defined as those reads that do not map within the IS but whose paired reads do). These IS1-flanking reads were then mapped to the CT18 chromosome reference sequence, to identify the sites of IS1 insertion in the chromosome.
Presence of the IncHI1 plasmid was confirmed by (i) BLASTN search of contig sets with the sequence of conserved backbone genes (164.1 kb) and antimicrobial resistance genes within Tn10 and a composite transposon, Tn2067, generated from comparative analysis of the nucleotide sequences of R27 (AF250878), pHCM1 (AL513383) and pAKU1 (AM412236)33 (Supplementary Table 2) using ACT72, and (ii) detection of an IncHI1 replicon directly from reads using SRST2 and the PlasmidFinder database (v 1.2)73. IncHI1 plasmid multilocus sequence typing (MLST) types23 were also determined using SRST2 (ref. 64). H58 S. Typhi isolates that contained fewer than 38% (63/168) of the backbone genes were excluded from further analysis. Additional plasmid replicons were identified from the SRST2 analysis of plasmid replicons, and the location of resistance genes was determined by manual investigation of the assemblies using BLASTN, Artemis and ACT72.
Nonsynonymous SNPs that were identified in >99% of the H58 isolates and in no other lineage were analyzed. The SNPs and their associated genes were studied using Genome Browser Artemis72. Functional categories were as annotated in the S. Typhi CT18 genome (AL513382.1).
Simple, descriptive statistics were used to compare the geographical distributions of lineages, prevalence of plasmids, and resistance genes and mutations and to calculate 95% confidence intervals. The significance of differences between studied groups of variables was calculated using Fisher’s exact test. All statistical tests were two-sided at α = 0.05, and analyses were performed using STATA (version 12.1, StataCorp) and R.
This work was supported by the Wellcome Trust. We would like to thank the members of the Pathogen Informatics Team and the core sequencing teams at the Wellcome Trust Sanger Institute (Cambridge, UK). We are grateful to D. Harris for his superb work in managing the sequence data. We also thank L. Fabre for her excellent technical assistance.
This work was supported by a number of organizations. The authors affiliated with the Wellcome Trust Sanger Institute were funded by Wellcome Trust award 098051; N.A.F. was supported by Wellcome Trust research fellowship WT092152MA. N.A.F., R.S.H. and this work were supported by a strategic award from the Wellcome Trust for the Malawi-Liverpool Wellcome Trust Clinical Research Programme (101113/Z/13/Z). C.M.P. was funded by the Wellcome Trust Mahidol University Oxford Tropical Medicine Research Programme, supported by the Wellcome Trust (Major Overseas Programmes–Thailand Unit Core Grant), the European Society for Paediatric Infectious Diseases and the University of Oxford–Li Ka Shing Global Health Foundation. D.D., P.N. and V.D. were supported by the Wellcome Trust (core grant 089275/H/09/Z). K.E.H. was supported by the National Health and Medical Research Council of Australia (fellowship 1061409) and the Victorian Life Sciences Computation Initiative (VLSCI; grant VR0082). C.A.M. was supported by a Clinical Research Fellowship from GlaxoSmithKline, and P.J.H. was supported by a UK Medical Research Council PhD studentship. This work forms part of a European Union Framework Programme 7 Marie Curie Actions Industry Academia Partnerships and Pathways (IAPP) Consortium Programme, entitled GENDRIVAX (Genome-Driven Vaccine Development for Bacterial Infections), involving the Wellcome Trust Sanger Institute, KEMRI Nairobi and the Novartis Vaccines Institute for Global Health. The authors affiliated with the Institut Pasteur were funded by the Institut Pasteur, the Institut de Veille Sanitaire and the French government ‘Investissement d’Avenir’ program (Integrative Biology of Emerging Infectious Diseases Laboratory of Excellence, grant ANR-10-LABX-62-IBEID). C.H.W. was supported by the UK Medical Research Council (MR/J003999/1). C.O. was supported by Society in Science and the Branco Weiss Fellowship, administered by ETH Zurich. A.K.C. was supported by the UK Medical Research Council (G1100100/1). J.J. was supported by the antibiotic resistance surveillance project in the Democratic Republic of the Congo, funded by project 2.01 of the Third Framework Agreement between the Belgian Directorate General of Development Cooperation and the Institute of Tropical Medicine (Antwerp, Belgium). F.M. was supported by a research grant from the Bill and Melinda Gates Foundation. The findings and conclusions contained within this publication are those of the authors and do not necessarily reflect positions or policies of the Bill and Melinda Gates Foundation. J.A. Crump was supported by the joint US National Institutes of Health–National Science Foundation Ecology and Evolution of Infectious Disease program (R01 TW009237), the UK Biotechnology and Biological Sciences Research Council (BBSRC; BB/J010367/1) and UK BBSRC Zoonoses in Emerging Livestock Systems awards BB/L017679, BB/L018926 and BB/L018845. S.K. was supported by US National Institutes of Health grant R01 AI099525-02. S.B. is a Sir Henry Dale Fellow, jointly funded by the Wellcome Trust and the Royal Society (100087/Z/12/Z). S.O. was supported by the National Institute of Allergy and Infectious Diseases of the US National Institutes of Health (R01 AI097493). The content is solely the responsibility of the authors and does not necessarily represent the official views of the US National Institutes of Health. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. C.D. was supported by the University of Oxford–Li Ka Shing Global Health Foundation.
Accession codes. Raw sequence data have been submitted to the European Nucleotide Archive (ENA) under accession ERP001718.
Note: Any Supplementary Information and Source Data files are available in the online version of the paper.
AUTHOR CONTRIBUTIONSStudy design: V.K.W., S.B., J. Parkhill, N.R.T., K.E.H. and G.D. Sequencing data generation: A.J.P., J.A.K. and E.J.K. Data analysis: V.K.W., K.E.H., J. Parkhill, N.R.T., A.J.P., J.A.K., D.J.E., J. Hawkey, S.R.H., A.E.M., A.K.C., J. Hadfield, C.O., R.A.K., E.J.K., D.A.G. and D.J.P. Isolate acquisition and processing and clinical data collection: D.J.P., S.B., N.A.F., N.R.T., F.-X.W., P.J.H., N.T.V.T., R.F.B., C.H.W., S.K., M.A.G., R.S.H., J.J., O.L., W.J.E., C.M., J.A. Chabalgoity, M.K., K.J., S. Dutta, F.M., J.C., C.T., S.O., C.A.M., C.D., K.H.K., A.M.S., C.M.P., A.K., E.K.M., J.I.C., S. Dongol, B.B., M.D., D.B., T.T.N., S.P.S., M.H., P.N., R.S.O., L.I., D.D., V.D., G.T., L.W., J.A. Crump, E.D.P., S.N., E.J.N., D.P.T., P.T., S.S., M.V., J. Powling, K.D., G.H., J.F. and K.E.H. Manuscript writing: V.K.W., S.B., K.E.H. and G.D. All authors contributed to manuscript editing. Project oversight: S.B., K.E.H. and G.D.
COMPETING FINANCIAL INTERESTS
The authors declare no competing financial interests.