Figure 1

Figure 1. Origin and year of collection for 437 clinical isolates collected and sequenced from different countries and patient communities as part of study of travel- and community-based transmission of multidrug-resistant Shigella sonnei...

We performed whole-genome sequencing on 437 S. sonnei isolates as part of this study. These isolates were from patients associated with OJCs outside of Israel (n = 171), from 221 patients in Israel (200 OJC, 21 of unknown ethnicity), or from patients in the United Kingdom not associated with OJCs (n = 45) (Technical Appendix 1). The isolates were collected from 6 countries (Israel, the United Kingdom, France, Belgium, the United States, and Canada) during 1992–2014 (Figure 1). The collection included isolates from most previously reported OJC-associated outbreaks of S. sonnei shigellosis; we defined cases as being OJC-associated separately for each public health agency (Table 1).

DNA was extracted at multiple sites by using the Wizard genomic DNA extraction kit (ProMega, Madison, WI, USA) according to manufacturer’s instructions. DNA was sequenced by using the MiSeq and HiSeq 2000 platforms (Illumina, San Diego, CA, USA) at multiple institutes according to in-house protocols (Technical Appendix 1) (25–27). Sequencing data for all isolates described passed internal quality control and were assembled into <687 contiguous sequences with a total length of <5.0 MB. Sequence data are available in the European Nucleotide Archive (http://www.ebi.ac.uk/ena; accession numbers in Technical Appendix 2).

Analysis of sequencing data was similar to that previously described (28). Multiple sequence alignment for phylogenetic analysis was generated by mapping to reference isolate S. sonnei Ss046 (GenBank accession no. CP000038), then masking mobile and repetitive elements (16) and stripping sites of recombination (29). Analysis of remaining variable sites was performed by using maximum-likelihood analysis in RAxML version 7.8.6 to create phylogenetic trees (30).

For isolates in the OJC-associated lineage for which the sampling date was known (n = 333; Technical Appendix 1), BEAST version 1.8 software (http://beast.bio.ed.ac.uk/) was used to estimate the emergence date of the lineage (31). Root-to-tip distances were generated by using Path-O-Gen version 1.4 (32). BEAST results shown are from 4 chains of 100 million Markov chain Monte Carlo generations run according to a general time reversible plus gamma substitution model, with a relaxed normal clock and Bayesian Skyline Population model, previously used for this pathogen (16,17). Chains were sampled every 1,000 generations with a 10% initial burn-in for root-height (time to most recent common ancestor) analysis. The maximum clade credibility tree was generated with a 10% burn-in and sampling every 100,000 generations. These results were consistent with those generated similarly by using a constant population growth model (Technical Appendix 2 Table).

De novo assembly, annotation, and antimicrobial resistance gene detection in the isolates was done as previously described (28) (accession numbers for annotated draft genome assemblies in Technical Appendix 1). Contiguous sequences containing antimicrobial resistance genes were extracted from assemblies and the presence of plasmid incompatibility groups on these contiguous sequences was determined by using PlasmidFinder (33). The presence of the Tn7/Int2 cassette was confirmed by mapping, and synteny detected by using ACT (34).

Figure 2

Figure 2. The OJC-associated lineage of multidrug-resistant Shigella sonnei in context with other global lineages. These background (non-OJC) isolates were contemporaneously collected and selected on the basis of phage typing; that is, including...

To determine the relationships among S. sonnei from OJCs inside and outside of Israel, we constructed a phylogenetic tree including whole-genome sequence data from 437 isolates of S. sonnei alongside the 118 isolates from the global context dataset (Table 1; Figure 2). This analysis showed the existence of a large, unique monophyletic sublineage (n = 396 isolates) of the global III lineage that was almost exclusively (388/396, 98%) composed of isolates from OJC-associated outbreaks and samples from persons in Israel (OJC-associated lineage; Figure 2; Technical Appendix 2 Figure 2). This lineage contained nearly all isolates (217/221, 98%) sequenced from Israel and collected during1992–2014 (Figure 2); 170/171 (99%) of those identified from samples collected during the same time frame from OJCs in the United States, Canada, France, Belgium, and the United Kingdom; and the 1 isolate from the global context dataset that originated in Israel (Figure 2; Technical Appendix 1). The clustering of most (387/392, 99%) of the strains from Israel and the other OJC-associated strains in the OJC-associated lineage is remarkable considering that the lineage represented approximately 10% of the diversity of the S. sonnei: the largest intra-lineage pairwise distance was 8.8-fold less in the OJC-associated lineage relative to the remainder of the tree (Figure 2).

The inclusion of isolates from patients in the United Kingdom that were not associated with OJCs provided further illumination of the association of this lineage with OJCs. These non-OJC samples were almost entirely (37/45, 82%) located outside the OJC-associated lineage, distributed elsewhere in the global III and II lineages (Figure 2). Ultimately, this resulted in a statistically significant association of lineage with sample designation (i.e., OJC or non-OJC) among UK isolates (p<0.0001 by Fisher exact test). This association correlated better with phylogenetic position than phage type, which was a comparatively poor indicator of genome level phylogeny (Technical Appendix 2 Figure 3). The relative phylogenetic positions of non-OJC and OJC isolates from the United Kingdom when viewed in an international context, i.e., including the other strains (Figure 2) showed that strains from UK OJCs were more likely to be related to strains from Israel than to strains circulating in non-OJCs. For example, strains from OJCs sampled in the United Kingdom in 2014 were phylogenetically adjacent to strains from Israel sampled in 2014, rather than to non-OJC strains sampled in the United Kingdom in 2014 (Figure 2.)

Figure 3

Figure 3. The OJC-associated lineage of multidrug-resistant Shigella sonnei across time and associated antimicrobial drug resistance. The Bayesian-inferred phylogenetic tree shows the evolutionary relationships of 333 isolates (those for which a fixed date...

Consistent with this finding, phylogenetic relationships within the OJC-associated lineage were defined more by time than geography (Figure 2; Technical Appendix 2 Figures 2, 4) Bayesian phylogenetic analysis showed that the OJC-associated lineage emerged in 1988 (95% highest posterior distribution 1985–1990) (Figure 3). Since that time, contemporaneously collected isolates were phylogenetically proximate with subsequent evolution, resulting in strain replacement rather than coexistence over time (Figure 3; Technical Appendix 2 Figure 2). Contrasting with the clear time signature in the lineage, geographic admixing of isolates occurred within the lineage (Figure 3). For example, >5 and 7 monophyletic clusters of isolates from OJC-associated outbreaks in France and the United Kingdom, respectively, were encompassed within the diversity of strains characterized in Israel (Figure 3). Samples from a single epidemiologic OJC-associated outbreak in Belgium were separated from each other by 44 single-nucleotide polymorphisms, making them phylogenetically distinct (Figure 2; Technical Appendix 2 Figure 2). This finding is consistent with contemporaneous OJC-associated outbreaks in different geographic areas representing real-time transmission events.

Because of the potential consequences of this intercontinental transmission to the transfer of AMR to S. sonnei, we determined the AMR characteristics of the OJC-associated lineage. We found that the lineage had a unique AMR profile relative to isolates outside of the OJC-associated lineage (Table 2). The OJC-associated lineage belonged to the global III lineage of S. sonnei, and every isolate in the OJC-associated lineage contained a Tn7/Int2 cassette that encoded the MDR thought to have facilitated the global dispersal of global III (16). This cassette contains the aadA1, sat2, and dfrA1 genes that confer resistance to aminoglycosides, streptothricin, and trimethoprim and was chromosomally integrated adjacent to the glmS gene in these isolates (indicating a single acquisition event). This region is identical to a Tn7-like island, also adjacent to the glmS gene, in the newly emerging Xv serotype of S. flexneri (reference strain 2002017 [35]; Technical Appendix 2 Figure 5). No mutations known to confer quinolone resistance were found in gyrA or parC sequences of isolates in the OJC-associated lineage, and no plasmid-encoded quinolone resistance genes were detected (Technical Appendix 1).

To consider AMR mechanisms that had potential to mobilize among bacteria, we further examined antimicrobial resistance genes that were inconsistently present (in 5%–95% of isolates) (Figure 3; Technical Appendix 1; Technical Appendix 2 Figure 2). The 7 genes found to be inconsistently present across the lineage were second copies of aadA1 present in some isolates; the aminoglycoside resistance–conferring strA and strB genes; the sulphonamide resistance genes sul1 and sul2; the tetracycline resistance gene tetA; and the ampicillin resistance gene bla_TEM (Figure 3; Technical Appendix 2 Figure 2). Attempts were made to determine the coinheritance and genetic carriage elements of these genes (within the limitations of genome assembly). In isolates that had additional copies of aadA1, the gene was typically co-inherited with the sul1 gene (Figure 3; Technical Appendix 1; Technical Appendix 2 Figure 2); this combination was found on plasmids of 2 different incompatibility groups, I1 and P, as well as on contiguous sequences where no plasmid incompatibility groups were identified, shown as unknown (Figure 3; Technical Appendix 2 Figure 2). Similarly, strA, strB, and sul2 were frequently co-inherited and found on plasmids of 4 different incompatibility groups. Isolates collected earlier in the lineage’s evolution tended to carry strA/strB/sul2 on B/O/K/Z, Q1, and P incompatibility group plasmids, whereas later isolates carried the genes on I1 plasmids (Figure 3; Technical Appendix 2 Figure 2). Similarly, the tetA gene appeared to have had 2 major introductions into the lineage, earlier on a P group plasmid and later on an I1 plasmid (Figure 3; Technical Appendix 2 Figure 2). Last, bla_TEM genes were found in 86% of isolates in the lineage, compared with 14% outside of the lineage (Table 2; Technical Appendix 1); these genes were carried on plasmids of 5 different incompatibility groups, with sporadic coinheritance patterns with other resistance genes (Figure 3; Technical Appendix 2 Figure 2).

We used whole-genome sequencing to develop a high-resolution picture of the international transmission of S. sonnei and its AMR determinants among OJCs over several decades. These analyses offer insight for the epidemiology of shigellosis inside and outside of Israel as well as for the broader transmission of AMR enteric pathogens. We showed that, in countries outside of Israel, outbreak strains in OJCs were distinct from strains circulating in the general population and that OJC-associated strains were more closely affiliated with outbreaks associated with OJCs in other countries (irrespective of geographic distance) and strains circulating in Israel. Strains from Israel and strains from nearly all previously reported OJC outbreaks elsewhere formed a distinct OJC-associated sublineage that emerged ≈30 years ago. Unlike other described emergent Shigella sublineages (17,28), the OJC-associated lineage lacked a defining association with AMR.

Isolates collected during outbreaks of S. sonnei in OJCs outside Israel were phylogenetically linked to contemporaneous isolates from Israel. This finding was suspected from previous studies that used pulsed-field gel electrophoresis, the results of which supported that samples from outbreaks among OJCs in the United States and Belgium were distinct from samples of S. sonnei circulating locally in non-OJCs (22,24) and were related to strains from Israel (22). We confirmed this link in an analysis of specimens collected in the United Kingdom that showed that strains from OJC-associated outbreaks were distinct from other circulating strains in that country but related to contemporaneous strains from Israel (Figure 2). The broader analysis, expanded in time and geography, showed that local epidemics in OJCs in France, Belgium, and North America (21,22,24) were also linked with contemporaneous isolates from Israel (Figure 3; Technical Appendix 2 Figure 2). This pattern of phylogenetic clustering by community affiliation was also recently demonstrated for S. flexneri 3a strain transmission among a global epidemiologic community of men who have sex with men, through which a unique MDR sublineage spread during ≈20 years (28). These studies demonstrate the speed with which AMR Enterobacteriaceae can be transmitted among persons in an internationally linked community rather than by contiguous geographic spread.

These findings also have implications for the epidemiology of shigellosis within Israel. The isolates from Israel in this study derive primarily from OJCs, which drive cyclic S. sonnei epidemics in Israel (19); here, they were shown to belong to a single, low-diversity sublineage. The lineage was monophyletic and had a strong time signature, consistent with a point introduction and subsequent epidemic emergence. This pattern was similar to that observed in Vietnam after the introduction of another global III S. sonnei sublineage (15,17). The date of the emergence of the OJC-associated lineage (1988 [95% highest posterior distribution 1985–1990]) is consistent with that estimated in a previous study where a 2-isolate lineage emerged in the Middle East during 1983 (16). Considering the timing and context of the emergence, it is possible that the OJC-associated lineage emerged from the large waterborne epidemic that occurred in Israel in the mid-1980s (19) or was potentially introduced from the first reported OJC-associated outbreak of shigellosis outside of Israel, which was an outbreak of >13,000 cases across the United States that also occurred in the 1980s (18). Samples from this period were not available to explore these origins of these outbreaks, but it is clear that the OJC-associated lineage is now endemic to OJCs in Israel and is causative of OJC-associated outbreaks elsewhere.

The AMR profile of the OJC-associated lineage is consistent with the phenotypic information and is likely influenced by antimicrobial resistance selection pressures in the 0–4-year age group, which is primarily affected by shigellosis in OJCs. Sulfonamide and tetracycline resistance determinants were in flux across the lineage (Figure 3), and these antimicrobial classes have been reported as being phenotypically dynamic over time among S. sonnei isolated in Israel (19). Furthermore, trimethoprim resistance was chromosomally encoded in all isolates, and plasmid-mediated ampicillin resistance was a common finding (86% of isolates) (Figure 3). These antimicrobial classes are key for the treatment of children with shigellosis. Similarly, resistance to tetracyclines, which can be used in children >7 years of age (8), and macrolides (mphA gene [21]) (Table 2) were also found. Despite being reported in other global III S. sonnei strains (16,17), resistance to quinolones was not found in the OJC-associated lineage, possibly because the use of quinolones is contraindicated in children.

The acquisition of the Tn7/Int2-encoded MDR in this lineage may have facilitated its epidemic emergence, as has been hypothesized for the broader global III lineage (16) and as is possible for S. flexneri Xv (35), although these possibilities cannot be explored fully by using these data. The presence of this gene cassette throughout the lineage and other phylogenetic clusters of antimicrobial resistance genes demonstrates that AMR can spread through this geographically dispersed, but closely associated, community. However, in contrast to the emergence of S. flexneri 3a among men who have sex with men or of S. sonnei in Vietnam, further acquisition of AMR does not appear to have shaped the subsequent evolution of the lineage. The presence of additional resistance genes (Figure 3) was not correlated with later time points, and the same genes were carried on distinct mobile genetic elements, consistent with sporadic reintroduction, rather than maintenance of the additional resistance genes in the population. This absence of a defining association with AMR suggests that it is probably primarily the epidemiologic suitability of OJCs to the transmission of S. sonnei that supports its maintenance in these communities. This likelihood is consistent with previous studies of OJC-associated shigellosis, which suggest that transmission is largely driven by communal childcare arrangements in a host population that is young and densely structured (19,20,22).

This study documents the speed with which MDR enteric bacteria can transmit intercontinentally through travel within a geographically dispersed, but closely linked, community. Awareness of this mode of sustained, geographically noncontiguous transmission must inform public health practice, including targeted control and reduction of consequences of the pathogen through developing effective relationships with affected communities; identifying specific risk factors; and designing, piloting, and eventually implementing specific culturally appropriate interventions with the participation and support of the community while avoiding stigmatization. Effective ongoing surveillance is also vital, and we demonstrated the resolution required for that purpose and the need for effective data sharing to track these otherwise silent transmission phenomena. Furthermore, the repeated application of high-resolution tools supports identification of parallels and contrasts with previous genomic epidemiology studies, supporting construction of a complete picture of the global transmission of AMR enteric pathogens, advancing our position for tackling this global public health issue.