Bacterial isolates of S. flexneri from 330 cases were submitted to the Gastrointestinal Bacterial Reference Unit between August 2015 and January 2016, from local and regional hospital laboratories in England and Wales. This strain set comprised the following phenotypic serotypes (numbers of isolates belonging to each serotype in parenthesis): 1a (3), 1b (16), 1c (25), 2a (185), 2b (16), 3a (41), 3b (7), 4av/E1037 (9), 6 (19) X (4), Y (4), X and Y (1). Epidemiological data on age, sex, and region of residence were available from laboratory report forms. Travel history was available for 176 of cases. All isolates were serotyped and WGS.

Phenotypic identification of S. flexneri isolates were confirmed using the Omnilog GenIII MicroPlate (Biolog, Hayward, CA, United States). Serotyping was carried out using standard methods by slide agglutination using both commercially available monovalent antisera (Denka Seiken, Japan) and monoclonal antibody reagents (Reagensia AB, Sweden) and in-house antisera raised in rabbits (Gross and Rowe, 1985). Molecular PCR was carried out on the discrepancies between phenotypic and WGS typing results as previously described (Gentle et al., 2016) and a serotype assigned according to the gene combination detected (Supplementary Table S1).

Genome sequencing and sequencing analysis were carried out as previously described (Dallman et al., 2015). Genomic DNA extracted using the QiaSymphony DNA extraction platform (Qiagen) from 330 S. flexneri was fragmented and tagged for multiplexing with Nextera XT DNA Sample Preparation Kits (Illumina) and sequenced using the Illumina HiSeq 2500 at PHE. A reference database containing the gene sequences encoding the 12 O-antigen synthesis or modification genes described by Sun et al. (2011, 2012a,b), including wzxc1-5, wzxe1-5, wzx6, gtrI, gtrII, gtrIV, gtrV, gtrX, gtr1c, oac oac1b and opt, was constructed. Using the GeneFinder tool (Doumith, unpublished), FASTQ reads were mapped to the S. flexneri O-antigen synthesis or modification genes using Bowtie 2 (Langmead and Salzberg, 2012) and the best match to each target was reported with metrics including coverage, depth, mixture and nucleotide similarity in XML format for quality assessment. Only in silico predictions of serotype that matched to a gene determinant at >80% nucleotide identity over >80% target gene length were accepted. FASTQ sequences were deposited in the National Center for Biotechnology Information Short Read Archive under the bioproject PRJNA315192 (see Supplementary Table S2 for SRA identifiers).

Short reads were quality trimmed and mapped to the reference S. flexneri serotype 2a strain 2457T (AE014073.1) (Wei et al., 2003) or the reference strain NC_007613 if S. flexneri serotype 6, using BWA v0.75 (Li et al., 2009; Bolger et al., 2014). The Sequence Alignment Map output from BWA was sorted and indexed to produce a Binary Alignment Map (BAM) using Samtools (Li and Durbin, 2010). GATK v2.6.5 was used to create a Variant Call Format (VCF) file from each of the BAMs, which were further parsed to extract only SNP positions which were of high quality (MQ > 30, DP > 10, GQ > 30, Variant Ratio > 0.9) (McKenna et al., 2010). Gubbins v2.0.0 (Croucher et al., 2015) was used to identify recombinant regions of the genome which were subsequently masked for phylogenetic analysis. Pseudosequences of polymorphic positions were used to create maximum likelihood trees using RAxML v8.1.17 (Stamatakis, 2014). De novo assembly was carried out using Spades 3.5.0 using ‘–careful’ and ‘ -k 21,33,55,65,77,83,91’ options (Bankevich et al., 2012).

To proactively detect outbreaks from WGS data, SNP typing was carried out on S. flexneri isolates belonging to clonal complex (CC) 245 and CC145. CC145 mostly comprises S. boydii serotypes but includes S. flexneri serotype 6, as this S. flexneri serotype was misidentified historically (Wirth et al., 2006; Chattaway et al., 2017). Hierarchical single linkage clustering was performed at seven descending thresholds of SNP distance (Δ250, Δ100, Δ50, Δ25, Δ10, Δ5, Δ0) as previously described (Dallman et al., 2016). This clustering results in a discrete seven digit code where each number represents the cluster membership at each descending SNP distance threshold. The resultant SNP addresses describes an isolates position in the S. flexneri population structure where two isolates with the same SNP addresses have 0 SNP differences between them.

There were 330 isolates reported between 1st August 2015 and 18th January 2016, 236 (71%) were from males, 85 (26%) from females, and for nine (3%) cases the sex was not stated. Travel history was not provided for 154 (47%) cases, 85 cases reported travel 7 days prior to onset of symptoms and 91 (27%) reported that they did not travel during the 7 days prior to onset of symptoms. The most frequently reported destinations were India (n = 16, 5%) and Pakistan (n = 11, 3%) with other counties accounting for less < 1% each.

Of the 330 cultures tested prospectively by WGS, 306 (93%) had concordant results with phenotypic serotyping. Three of the mismatched results between the phenotypic and genotypic tests were attributed to mutations identified in O-antigen synthesis or modification genes (Table 1), including nonsense mutations resulting in early stop codons (n = 2) and a frameshift mutation. Repeat testing of sample 4 revealed the discrepancy was due to an auto-agglutination reaction of the strain with the sera (Table 1).

The remaining mismatches were novel serotype gene profiles detected by WGS, both associated with outbreaks (Cluster 2, n = 13; Cluster 4, n = 5) (Tables 1, 2). Cluster 2 comprised 13 isolates associated with a local community outbreak. The isolates failed to agglutinate any of the serotype specific S. flexneri antisera and were positive for wzx1-5 and gtrIc gene targets in the PCR but were negative for gtrI. Typical strains of serotype 1c are positive for wzx1, gtrI, and gtrIc (Supplementary Table S1). These strains were designated 1c variants (1cv). Cluster 4 comprised five isolates of S. flexneri associated with an outbreak of gastrointestinal symptoms in five captive chimpanzees phenotypically identified as serotype 3a. The isolates had wzx1-5, oac, and grtX but were also positive for gtrII and were designated 3a variant (3av).

Analysis of the WGS data organized 161 of the 330 isolates into 36 five SNP single-linkage clusters of two or more isolates associated with CC245 and one cluster in CC145. The median number of cases in these clusters was 2 and ranged from 2 to 27; 32 (89%) clusters investigated comprised less than 5 isolates (Table 2). It was not possible to identify epidemiological links associated with these small clusters using only the limited epidemiological data available from laboratory report forms. However, 9/32 (28%) comprised at least one case reporting recent travel abroad prior to onset of symptoms.

During the period of the study, four outbreaks were identified following routine surveillance of local hospital reports of gastrointestinal symptoms caused by S. flexneri. SNP typing confirmed that all the isolates belonging to each outbreak cluster were closely related and monophyletic. Cluster 1 was the largest cluster and had a high male to female ratio with the 97% of cases being adult males (for two cases the gender was not stated) (Table 2). Isolates from this cluster were observed throughout the study period. The minimum SNP distance between isolates in this cluster was zero with the median distance 17. The maximum SNP distance between any two isolates was 64, however the majority of these SNPs were caused by a transposase mediated recombination of pic, encoding a serine protease. This demographic of adult males has previously been shown to be characteristic of clusters linked to sexual transmission among the MSM community (Borg et al., 2012; Gilbart et al., 2015; Baker et al., 2015; Simms et al., 2015). Cluster 1 was part of larger outbreak of S. flexneri serotype 2a previously described by Simms et al. (2015).

Clusters 2 and 3 were community outbreaks and cases were geographically linked. They were temporally restricted and genetically homogenous, exhibiting 0–1 SNPs difference in the core genome between isolates. Despite a thorough epidemiological investigation, the source and route of transmission associated with Cluster 2 could not be determined. Cluster 3 was linked to consumption of contaminated food at a restaurant, most likely due to an infected food handler. Cluster 4 was associated with an outbreak of gastrointestinal symptoms in a group of captive chimpanzees with a minimum SNP distance between isolates of one and a maximum SNP distance of four.

The S. flexneri population structure clusters into seven distinct phylogenetic groups (PGs) (Connor et al., 2015). Phylogenetic analysis of the diversity of CC245 in the PHE collection is shown in Figure 1. Four out of seven clades are represented by samples received by PHE through routine surveillance. With respect to the clusters described in this study, Clusters 1, 3, and 4 fall within PG3 and Cluster 2 falls within PG1.