An extended genotyping framework for Salmonella enterica serovar Typhi, the cause of human typhoid

The population of Salmonella enterica serovar Typhi (S. Typhi), the causative agent of typhoid fever, exhibits limited DNA sequence variation, which complicates efforts to rationally discriminate individual isolates. Here we utilize data from whole-genome sequences (WGS) of nearly 2,000 isolates sourced from over 60 countries to generate a robust genotyping scheme that is phylogenetically informative and compatible with a range of assays. These data show that, with the exception of the rapidly disseminating H58 subclade (now designated genotype 4.3.1), the global S. Typhi population is highly structured and includes dozens of subclades that display geographical restriction. The genotyping approach presented here can be used to interrogate local S. Typhi populations and help identify recent introductions of S. Typhi into new or previously endemic locations, providing information on their likely geographical source. This approach can be used to classify clinical isolates and provides a universal framework for further experimental investigations.

T yphoid fever (typhoid), caused by Salmonella enterica serovar Typhi (S. Typhi) bacteria, is a systemic human infection that affects an estimated 20.6 million people globally each year, causing an estimated 223,000 deaths [1][2][3] . Typhoid remains endemic in populations with limited access to sanitation and safe water, and is a notifiable or reportable infection in many industrialized countries, where it is generally associated with travel to endemic areas. Public health laboratories have relied on techniques such as phage typing 4,5 or pulsed-field gel electrophoresis 6 , which are phylogenetically naive and have limited discriminatory power to support epidemiological investigations and surveillance.
A genotyping scheme based on 88 single-nucleotide polymorphisms (SNPs) identified within a limited set of genes was previously developed for S. Typhi 7 . This enabled the classification of the S. Typhi population into 85 haplotypes (haploid genotypes) based on biallelic profiles and provided the first phylogenetic framework for epidemiological studies 8 . Subsequently, wholegenome sequencing (WGS) has been used to identify many more SNPs and other phylogenetically informative markers for discriminating within S. Typhi, which has limited genetic variation [9][10][11][12][13][14][15][16] . Similar progress has been made in other monophyletic clades of bacterial pathogens, such as Mycobacterium tuberculosis 17 and Yersinia pestis 18 .
We have recently reported the WGS of almost 2,000 S. Typhi isolates sourced from 63 countries 14 . This study identified 422,000 chromosomal SNPs in the core genome, which were used to build a comprehensive phylogenetic tree. Notably, the analysis confirmed the emerging dominance of the multidrug resistance-associated H58 clade, including the recent spread of H58 S. Typhi into Africa, confirming the value of SNP-based WGS analysis of S. Typhi to understand contemporary typhoid epidemiology. Here we utilize these WGS data to define a global population framework for S. Typhi and to define a new genotyping scheme comprising 68 SNPs that provides extensive coverage of typhoid-causing bacteria circulating globally. Given the increasingly widespread adoption of WGS by public health laboratories for the tracking of bacterial pathogens 19,20 , we further aimed to explore the utility of S. Typhi WGS data, analysed via genotyping, to predict the geographical source of travel-associated S. Typhi isolated in the United Kingdom. This approach gives greater discriminatory power and improved phylogenetic information than the earlier scheme 7 , and forms a robust framework for public health surveillance, epidemiological investigations and laboratory experiments of typhoid.

Results
Defining phylogenetically informative genotypes for S. Typhi. In order to develop a comprehensive genotyping system, we used WGS data from 41,800 globally representative S. Typhi 14 to identify phylogenetically informative clades and subclades based on SNP architecture 21 . A summary of the isolates is shown in Table 1 and full details are provided in Supplementary Data 1 and Supplementary Table 1. Using a combination of phylogenetic tree topology and population genetic methods (using BAPS; Bayesian Analysis of Population Structure 21 ), we defined 16 S. Typhi clades that could be further divided into 49 subclades (Fig. 1, see Methods). Most of the clades could be grouped into four nested clusters (1-4, which we refer to as 'primary clusters'), each with 100% bootstrap support and defined by 420 SNPs (coloured branches in Fig. 1a). The median pairwise distances between isolates were as follows: 25 SNPs within subclades, 109 SNPs within clades and 243 SNPs between clades. We labelled these primary clusters, clades and subclades using a structured hierarchical nomenclature system similar to that used for M. tuberculosis 17 , whereby cluster 1 is subdivided into clades 1.1 and 1.  Fig. 1a). However, mapping the Roumagnac haplotypes to the wholegenome phylogeny showed that the older scheme provides highly uneven resolution across the S. Typhi phylogeny ( Supplementary  Fig. 1b), with a lack of resolution in some cases (11 Roumagnac haplotypes span two or more distinct subclades each; for example, H52 comprises clades 3.4, 3.5, 4.1 and 4.2) and excessive resolution in others (24 subclades are further divided into two or more haplotypes in the Roumagnac scheme).
A new SNP-based genotyping framework for S. Typhi. We identified a minimum set of 68 SNPs that can be used to genotype S. Typhi into the four primary clusters, 16 clades and 49 subclades. For each of these groups, we identified all SNPs that were unique to members of the group, and selected one such SNP to be used for genotyping. We prioritized the inclusion of synonymous intragenic SNPs (that is, located within a protein-coding sequence, but with no change to the encoded amino acid), within genes that showed evidence of genetic stability within the S. Typhi population (that is, nucleotide diversity o1% and dN/dS o0.7 across the global data set, with no inactivating mutations identified). Details of the genotyping SNPs are given in Supplementary Table 2. This genotyping scheme has greater discriminatory power than the original Roumagnac haplotyping scheme (D ¼ 0.96 versus 0.78), is phylogenetically informative by design and the hierarchical nomenclature of genotypes is intrinsically informative with respect to phylogenetic relationships between clades and subclades.
Geographical distribution of S. Typhi clades and subclades. Next, we examined the geographical distribution of S. Typhi genotypes. For these analyses, isolates of the same subclade, country and year were collapsed to a single representative to reduce the impact of localized outbreaks on our collection; this resulted in 541 unique isolates for analysis. Primary clusters 2, 3 and 4 were broadly distributed across continents (greens, blues and reds, respectively, in Fig. 1c), likely reflecting the relatively ancient spread of S. Typhi across the globe. Isolates outside these clusters, which result from deep branching closer to the root of the S. Typhi whole-genome tree, were rare in our collection (n ¼ 24 unique isolates) and mostly found in Africa (n ¼ 16). While the three common clusters (2)(3)(4) were present in most regions we analysed, cluster 2 predominated among American isolates (n ¼ 18/23 unique isolates, 78%). Most clades were detected on multiple continents (n ¼ 11/16) and included isolates from Asia (n ¼ 13/16) and/or Africa (n ¼ 10/16), which together  made up 78% of our isolate collection (Table 1). However, there were differences in the geographic distributions of clades, with most clades being dominated by unique isolates from a single continent (Asia, Africa or Oceania; see Supplementary Fig. 2).
In contrast, at the subclade level, only 22% of subclades (n ¼ 11) were found on more than one continent, and most were dominated by unique isolates from a single country or region: 40 subclades (82%) had Z50% of non-outbreak isolates from a single country (Fig. 2) and 44 subclades (90%) had Z50% of non-outbreak isolates from a single region (Fig. 2). A total of 28 subclades comprised five or more non-outbreak isolates each, and of these common subclades, 12 (43%) were detected in a single region only (six in Oceania, five in Southeast Asia and one in South Asia; Fig. 2). In total, 16 common subclades (57%) were highly restricted to a region (490% of isolates drawn from a single region) and 20 (71%) were generally associated with one region (470% of isolates drawn from a single region; Fig. 2). These data suggest that most S. Typhi subclades represent localized bacterial subpopulations with barriers to geographical dispersion, and that transfers to new locations rarely result in long-term establishment of local populations. In contrast with this general pattern, subclade 4.3.1 (previously H58) was found in nine different regions across Africa, Asia and Oceania. Only 10 other subclades (20%) were found on more than one continent, and the majority of these were dominated either by Asian, African or Oceanian isolates (Fig. 2). Thus, the recent global dissemination of subclade 4.3.1, which spread out of South Asia B30 years ago and has established successful local clonal 0   50  83  92  99  100  100  100  100  100  100 0    Genomic prediction of the geographical origins of S. Typhi by comparison with the global framework. Since most S. Typhi subclades were associated with a narrow geographical source, we hypothesized that genotyping of S. Typhi isolates could be used to predict the likely geographical origins of typhoid cases. As this is clearly challenging for the more widely distributed subclades, we also sought to examine whether specific SNPs could be used to predict origins down to the country level. For 1,501 out of 1,831 (82%) isolates in our global collection, the genetically closest isolate was from the same country. Where the closest isolate was 0-1 SNPs away, this frequency was 95% and for o10 SNPs, 90% ( Supplementary Fig. 3). Since our current global genome collection includes groups of isolates that were frequently collected from the same time and place, this should not be taken as a reliable measure of the general predictive power of SNP distance for S. Typhi. In order to further explore the power of our global genomic framework to predict geographic origins of travel-associated typhoid, we sequenced and genotyped 99 novel S. Typhi that were isolated from patients attending a hospital in East London, United Kingdom between 2005 and 2010 ( Table 2). A total of 13 genotypes were identified. Epidemiological interviews were able to link 81 of these cases with travel to a specific country; the remaining 18 cases were not associated with travel. The median SNP distance between these novel isolates and genomes in our global collection was 21 SNPs (interquartile range, 18-25 SNPs), posing a challenge for prediction of their geographical origin. Among the 81 travel-associated UK isolates, 53 were genotyped as 4.3.1; these were all linked to travel to countries within South Asia (Table 3), and clustered along with South Asian isolates from the global collection ( Fig. 3 and Supplementary Fig. 4n). For the 28 non-4.3.1 travel-associated UK isolates, the location of travel generally matched the geographical origin of the closest isolate (in terms of number of SNPs) in the global collection: travel location and closest global isolate source matched at the region level in all cases, and at the country level in most cases (n ¼ 20, 71%). That is, prediction of geographical origin based on the  Closest genome country match, number of travel-associated isolates whose country of travel matched that of the closest genome in the global collection (based on lowest number of SNPs); N (global), number of isolates in the global collection that were assigned to this subclade; N (travel), number of travel-associated isolates that were assigned to the subclade; Region frequencies, frequency of each geographic region among isolates of this subclade from the global collection (note groups of isolates from the same subclade, country and year were classified as outbreaks and represented only once per group in the frequency calculations); SNP, single-nucleotide polymorphism; S. Typhi, Salmonella enterica serovar Typhi. *Highlights the most frequent region for this subclade among the global collection, where this matches the region of travel.

A u s t r a li a O c e a n ia E a s t A s ia S o u t h e a s t A s ia S o u t h A s ia W e s t e r n A s ia N o r t h A f r ic a W e s t A f r ic a C e n t r a l A f r ic a S o u t h A f r ic a S o u t h e r n A f r ic a E a s t A f r ic a
closest strain of known location in the current global framework would have yielded the correct region of origin in all cases, and the correct country of origin in 71% of cases (95% confidence interval (CI), 66-76%). Furthermore, for non-4.3.1 subclades, genotyping alone was predictive of geographical origin at the regional level for the same proportion of isolates (71%). It is likely that power to predict the geographical sources of UK isolates would be improved by wider geographical coverage in the reference genome collection. Two of these isolates were genotyped as subclade 3.1.1 and linked with travel to Ghana and Nigeria; the closest isolates in our global collection were 16-17 SNPs away and were not from these precise locations, but likely originated from bordering countries in West Africa (Supplementary Fig. 4e). It is likely that a deeper coverage of West African isolates in our global framework would provide greater power to resolve geographic associations within this region, which comprised less than 2% of our current global collection (n ¼ 30 isolates). Similarly, for the other travel-associated isolates for which the recorded country of travel did not match the closest genome in the global tree, the closest genome was also from a neighbouring country (for example, Pakistan, India, Bangladesh; see Supplementary Fig. 4).
Genomic predictions of the geographical origins of 18 non-travel-associated UK isolates are shown in Table 4 and Supplementary Fig. 4. Thirteen isolates were 4.3.1 and clustered together with travel-associated isolates from South Asia, within a broader group of South Asian 4.3.1 isolates (Fig. 3). This suggests that S. Typhi imported into the United Kingdom from these regions have likely been transmitted onwards within the United Kingdom to individuals with no recent travel history ( Supplementary Fig. 4n). Two additional isolates were from subclades that were dominated by a single region in our global collection-3.1.1 (68% West Africa) and 3.3.0 (83% South Asia). Notably, while the 4.3.1 isolates were closely related to travelassociated isolates recently obtained in London, they were Z17 SNPs away from any isolates in the global collection. Thus, the diversity captured by the global collection does not provide the resolution to precisely identify the origin of these isolates 23 .

Discussion
Our data show that the global S. Typhi population consists of 49 distinct subclades that are strongly geographically clustered, with many locations harbouring subpopulations of S. Typhi established over long periods of time. We show how these ARTICLE subclades can be identified through a simple genotyping scheme consisting of 68 SNPs. Importantly, while we show that this scheme is highly phylogenetically informative, it can be readily inferred from raw sequence data without the need for multiple genome comparisons, phylogenetic analysis or any other complex or computationally intensive steps. Such properties make this universal SNP-based system a valuable tool upon which researchers can develop future studies. The S. Typhi genome is highly stable and exhibits minimal genetic variation and virtually no recombination 9,14 , and we recently estimated the substitution rate to be slower than one SNP per genome per year 14 ; therefore, the genotyping framework is expected to be robust to future evolution.
Owing to the strong geographical clustering of the various subclades, whole-genome comparison of novel S. Typhi isolates to the existing global population framework is strongly predictive of geographic origin at the regional level and has the potential to accurately predict origins to the country level. This has important public health implications for typhoid surveillance and control in endemic and non-endemic areas; however, ongoing updates to the global genomic framework will be important to ensure the utility of genomic surveillance for typhoid. For example, we found that the origin of travel-associated 4.3.1 isolates could not be resolved using the prior global framework alone, but benefitted from updated information provided by other recent travelassociated isolates of known geographical origin. This illustrates the importance of expanding and updating the global genomic framework through sequencing of novel isolates and suggests that, while ongoing surveillance in endemic areas is undoubtedly important, the use of clinically well-characterized travel-associated organisms isolated in non-endemic countries may also provide a valuable source for improving the granularity of data in the framework for genome-based surveillance of S. Typhi 23 . In addition, it will be important to expand the current global framework to include more recent isolates (the most recent in our current collection was from 2013) as well as isolates from regions that are currently under-represented (including Africa, the Americas and northeast Asia).
WGS-equipped reference laboratories provide a highly accessible source to expand the global genomic framework for typhoid, with potential benefits to local but also global typhoid control. For example, in England, Wales and Northern Ireland B520 typhoid cases are reported annually to the national reference laboratory (Public Health England). These cases are investigated in order to determine whether they are associated with travel to typhoid endemic regions 23 . However, approximately one-fifth of typhoid cases in the United Kingdom cannot be traced to a country of origin. At present, Public Health England provides molecular typing, which since April 2015 includes WGS as well as antimicrobial susceptibility profiling, for S. Typhi isolated from such cases. The resulting data are considered important for local epidemiology. However, we propose that this could also serve as a proxy for informal surveillance of typhoid molecular epidemiology in endemic regions. This may prove particularly valuable when supported by our genotyping framework for simplified attribution.
Diego, CA, USA) according to the manufacturer's protocols to generate tagged 100 base pair (bp) paired-end reads.
SNP analysis. For analysis of SNPs, the paired-end reads were mapped to the reference genome of S. Typhi CT18 (ref. 25), using SMALT (version 0.7.4; http:// www.sanger.ac.uk/resources/software/smalt/). SNPs were identified as previously described 14 , using samtools mpileup 26 and filtering with a minimum mapping quality of 30 and a quality ratio cutoff of 0.75 (ref. 24). SNPs located within phage regions, repetitive sequences or recombinant regions were excluded as previously described 14 , resulting in a final set of 22,143 chromosomal SNPs in an alignment length of 4,275,037 bp for the global collection of 1,831 S. Typhi isolates. An expanded alignment comprising 22,673 SNPs from S. Typhi isolates from the global collection (1,831) plus 99 traveller-associated UK isolates was generated using the same procedures as above. Pairwise SNP distances between isolates (that is, the number of core genome SNP loci at which pairs of isolates had discordant alleles) were extracted from each alignment i using the ape package 27 for R (v3.2; function call: dist.dna(i,model ¼ "N",pairwise.deletion ¼ T)).
Phylogenetic analyses. The maximum likelihood (ML) phylogenetic tree shown in Fig. 1 was built from the 22,143-SNP alignment of all 1,831 isolates using RAxML (version 7.8.6) 28 with the generalized time-reversible model and a Gamma distribution to model site-specific rate variation (the GTR þ g substitution model; GTRGAMMA in RAxML). The tree was outgroup-rooted by including a pseudosequence comprising S. Paratyphi A alleles in the alignment. Support for the ML phylogeny was assessed via 100 bootstrap pseudo-analyses of the alignment data.
The backbone topology of the global ML tree, showing relationships between subclades (Fig. 1b), was recovered by randomly selecting one isolate from each subclade to retain, and removing all other tips from the tree (using drop.tips() in the ape package 27 for R (v3.2)). A ML phylogenetic tree was also generated separately from 22,673 SNPs of S. Typhi isolates from the global collection (1,831) plus the 99 East London traveller-associated isolates, using the same procedures as above. All ML trees were visualized and annotated using Python (https://github.com/katholt/plotTree/#python-code).
Identification of phylogenetically informative clades and subclades. In addition to the whole-genome phylogenetic analysis outlined above, we investigated the population structure of the global S. Typhi collection using a phylogeny-free population genetics approach, implemented in BAPS v.6.0 (ref. 21). Hierarchal clustering analyses were conducted on identified clusters until singlemember clusters were obtained, thus allowing the discovery of nested genetic population structures 21 . Ten nested levels of molecular variation were fitted to the data using 10 independent runs of the stochastic optimization algorithm with the a priori upper bound of the number of clusters varying over the interval 50-300 across the runs 30 .
As our goal was to identify genotypes that were both phylogenetically and epidemiologically informative, we explored the homogeneity (1-Simpson's diversity) of geographical source within BAPS clusters (as an indicator of the potential power of genotyping to identify geographical origin of travel-associated isolates) at different levels of clustering ( Supplementary Fig. 5). This showed that within-cluster homogeneity increased up to the sixth level of clustering and then reached a plateau, with deeper clustering providing no greater resolution of geographical origin ( Supplementary Fig. 5). The third level of clustering resulted in most clusters being dominated by a single continent (14/17 clusters with 480% of isolates from one continent), while sixth-level clustering resulted in most clusters containing isolates from a single country (60/89 clusters with 480% of isolates from one country; Supplementary Fig. 6). We therefore used the BAPS clusters to guide the definition of clades (BAPS level 3) and subclades (BAPS level 6).
In order to maintain compatibility with the phylogeny, some minor modifications of the raw BAPS clusters were required (this consisted of subdividing some BAPS clusters and merging others, but not reassigning members between clusters; see Supplementary Fig. 7). The modified level-3 BAPS clusters were designated 'clades' and were assigned labels of the form [ Some BAPS clusters were polyphyletic and consisted of isolates belonging to rare phylogenetic lineages whose common ancestor in the phylogenetic tree coincided with the common ancestor of an entire clade (n ¼ 9) or primary cluster (n ¼ 2). These groups contain isolates that, given increased numbers, may emerge as distinct BAPS clusters that form sister taxa within the parent clade (or primary cluster), and were thus designated [z] ¼ 0 (or [y] ¼ 0) to indicate non-equivalence with the properly differentiated sister clades (n ¼ 16) or subclades (n ¼ 49). For example, while the genotypes 2.1 and 2.2 represent distinct sister clades that are each monophyletic, isolates assigned to 2.0 are paraphyletic and include multiple lineages that could not be further subdivided by BAPS analysis ( Supplementary  Fig. 7). SNP-based genotyping. We identified a minimum set of 68 SNPs with which to rapidly genotype S. Typhi into the 16 clades and 49 subclades, as described above (Supplementary Table 2). Short read alignment (BAM) files, generated by mapping Illumina reads to the CT18 reference genome (accession AL513382), were used to assign genotypes for each novel read set using a custom Python script (available at https://github.com/katholt/genotyphi). Briefly, the script uses samtools mpileup to extract from each BAM file the consensus base calls at the SNP loci. The resulting variant call format file is then processed to identify the presence of cluster-, cladeand/or subclade-defining SNP alleles (defined in Supplementary Table 2) that pass a minimum quality threshold (default consensus base Phred score Z20) and uses these to assign the read set to a cluster, clade and subclade. Discriminatory power was calculated using the method outlined in ref. 31.
Data availability. Raw sequence data are available in the European Nucleotide Archive under accession ERP001718. Supplementary Data 1 lists accession numbers for each isolate. The software for Microreact interactive tree viewer is available at: http://microreact.org/project/styphi 22 . SMALT is available at: http:// www.sanger.ac.uk/resources/software/smalt/. Python script to visualize and annotate trees is available at https://github.com/katholt/plotTree/#python-code. Python script to call SNPs is downloadable at https://github.com/katholt/ genotyphi.