Whole genome sequencing of increased number of azithromycin-resistant Shigella flexneri 1b isolates in Ontario

Azithromycin (AZM) resistance among Shigella is a major public health concern. Here, we investigated the epidemiology of Shigella flexneri serotype 1b recovered during 2016–2018 in Ontario, to describe the prevalence and spread of AZM resistance. We found that 72.3% (47/65) of cases were AZM–resistant (AZMR), of which 95.7% (45/47) were among males (P < 0.001). Whole-genome based phylogenetic analysis showed three major clusters, and 56.9% of isolates grouped within a single closely-related cluster (0–10 ∆SNP). A single AZMR clonal population was persistent over 3 years and involved 67.9% (36/53) of all male cases, and none reported international travel. In 2018, a different AZMR cluster appeared among adult males not reporting travel. A proportion of isolates (10.7%) with reduced susceptibility to ciprofloxacin (CIP) due to S83L mutation in gyrA were AZM susceptible, and 71.4% reported international travel. Resistance to AZM was due to the acquisition of mph gene-bearing incFII plasmids having > 95% nucleotide similarity to pKSR100. Plasmid-borne resistance limiting treatment options to AZM, ceftriaxone (CRO) and CIP was noted in a single isolate. We characterized AZMR isolates circulating locally among males and found that genomic analysis can support targeted prevention and mitigation strategies against antimicrobial-resistance.

Resistance to AZM among all isolates except one was primarily due to the plasmid-borne mph (A) gene encoding a macrolide 2′-phosphotransferase that inactivates macrolides.However, one AZM R isolate with MIC of 64 mg/L and lacking the mph or erm genes had a mutation in the rplD gene that encodes 50S ribosomal protein L4 (V126I), as well as a G98S mutation in rplL gene that encodes 50S ribosomal protein L7/L22, and a A146V mutation in the rplF that encodes 50S ribosomal protein L6.Though mutations at these positions have not been previously reported to cause macrolide resistance for Escherichia coli or Shigella spp., macrolides are known to bind at the 50S ribosome and recent reports show that mutations in other positions in these genes may confer resistance to macrolides, including erythromycin and AZM, in other species 12,13 .
Among AMP R isolates, either bla TEM-1B or bla OXA-1 genes were present (n = 18).A single CRO R isolate had a bla CTX-M-15 gene.Two CIP R isolates contained the qnrS gene and one isolate also had a mutation at position 87 (D87N) in gyrA gene.Seven isolates were CIP-intermediate resistant (CIP IR ) and had a S83L mutation in the gyrA gene.All SXT R had the dfrA and sul genes (Supplementary Table S1).

Phylogenetic analysis
On performing whole genome sequencing (WGS) on all S. flexneri 1b isolates (n = 65), most isolates grouped (76.9%; 50/65) into three clusters (Fig. 2, Supplementary Table S1).Among these, two clusters were AZM R and the second cluster emerged in Ontario in 2018.For phylogenetic analysis, all serotype 1b isolates were mapped to sflex002 (GENBANK accession XXX).The majority AZM R isolates (78.7%; 37/47) grouped together (cluster I) and differed by 0-9 SNP from each other.All isolates in cluster I were AZM R , harbouring mph (A) as well as aadA5, sul1 and dfrA17 resistance genes.Cluster I isolates were recovered throughout 2016-2018.All except one cluster I case, were from males aged > 18 years and the majority were from the GTA (75.0%; 27/36) (Fig. 2).None of the cluster I cases reported travel outside of Canada.
A second group of seven isolates, bearing gyrA substitutions (S83L), clustered together (cluster II, ∆SNP 4-31).Isolates were recovered in 2016 (n = 4) and 2018 (n = 3).Isolates in this cluster were recovered from both males (n = 3) and females (n = 4).All isolates in this cluster, except one, were susceptible to AZM.Five of the seven cases (71.4%) had reported travel to Cuba.A small sub-cluster (n = 4) within cluster II with 4-11 ∆SNPs was observed and single isolate was AZM R (Fig. 2).Finally a third cluster consisting of AZM R isolates recovered from males with no travel history appeared in 2018 (cluster III, ∆SNP 2-8, n = 7).These isolates had mph, bla TEM-1B and bla OXA-1 genes.www.nature.com/scientificreports/Some S. flexneri 1b isolates did not cluster together, 61.5% (8/13) of which were recovered from females.These isolates differed by 49 to 115 SNPs from the reference strain (sflex002).The majority of these isolates (76.9%; 10/13) were susceptible to AZM and recovered from cases reporting international travel (53.8%; 7/13).Among these isolates, a single non-travel related S. flexneri 1b isolate (sflex069) carried bla CTX-M-15 and mph genes and was confirmed to be phenotypically resistant to CRO and AZM.

A highly related plasmid bearing mph (A) gene conferring AZM resistance
Most of AZM R isolates had incFII plasmids bearing mph gene with high similarity to pKSR100.Plasmids bearing mph (A) gene from cluster I (isolates sflex002 and A-6726) had an incFII replicon and were highly similar to each other (99.9% identity) though recovered from patients in different years (Fig. 3A).The mph (A) gene was flanked by IS26 and IS6100 belonging to the IS6 family of insertion sequences.Additional antibiotic resistance genes (aadA5, sul1 and dfrA17) that are part of a class I integron were present (Fig. 3A).This plasmid was highly similar to the pKSR100 plasmids (99.7% nucleotide similarity, 85% coverage).All cluster I isolates had plasmidborne mph (incFII) with high similarity to psflex002 and pKSR100 4 .
The plasmids bearing the mph gene from cluster III isolates (pA-6817 and pA-10383) had 93% homology and 60% length coverage, with a similar backbone as pKSR100 (Fig. 3B).Cluster III plasmids bearing mph had two replicon types, incFIB and incFII and all harboured multiple resistance genes (dfrA14, bla TEM-IB , aph (6)-Id, and aph (3″)-Ib and sul1) (Fig. 3B).Units (PHUs), year of receipt, gender (M = male, F = female, U = other/unknown), travel-status (colored purple for international travel, light grey = if international travel was not indicated in iPHIS) and travel destination and azithromycin susceptibility (AZM; R = Resistant, S = Susceptible).SNPs were called against reference strain sflex002 belonging to cluster I (labelled in red font).Maximum Likelihood Tree (ML-Tree) was generated using IQ-TREE with the GTR + ASC method and ultrafast bootstrap over 5000 replicates.The antimicrobial resistant genes present are indicated as dark grey were detected by using SRST2 with ResFinder as the choice of database.Each cluster is labelled and colored with a different set of colors and SNP differences within each cluster are indicated.Bootstrap values of major nodes are shown.
Vol A plasmid from a single isolate (sflex069) that did not belong to any cluster also had the mph (A) gene flanked by IS26 and IS6100 (Fig. 3C).Interestingly, this pKSR100-like plasmid also contained bla CTX-M-15 and qnrS genes conferring resistance to CRO and CIP, respectively, included in a DNA fragment flanked by different transposons that could be involved in its acquisition.Comparison of this plasmid with mph-carrying plasmids from clusters I as well as to the pKSR100 showed high similarity (> 95% nucleotide homology) as they all shared a similar plasmid backbone and also harbored dfrA17, aadA5 and sul1 genes.

Discussion
Our study describes the epidemiology of S. flexneri serotype 1b isolates in Ontario from 2016 to 2018.Though there is limited literature on distribution of S. flexneri serotypes, geotemporal differences in resistance and serotype distributions exist.When we compare the prevalence of serotype 1b in comparison to other S. flexneri serotypes, serotype 1b constituted 20.6% of all S. flexneri in Ontario during our study while percentage of serotype 1b in the UK from 2004 to 2015 was 10% 17 .AZM R was observed in 72.3% of the S. flexneri 1b isolates included in this study.There is limited data available where resistance levels is reported for S. flexneri for each serotype.During the same period, percentage of AZM R isolates varied from 32.9 to 53.1% in other high-income countries for S. flexneri, however these don't report resistance for each serotype 5,10,18 .
We report dissemination of resistance to AZM among S. flexneri 1b isolates that are not travel related in Ontario.The AZM R cases in Ontario (clusters I and III) were male-dominated and none reported foreign travel.We also showed that plasmids among S. flexneri 1b in this study had genetic similarity to the pKSR100 plasmid that caused several MSM outbreaks globally during 2008-2018 4,5,7,14 .Since isolates in cluster I had a SNP distance of 0-9, this suggests that most cases of cluster I were part of a cluster that was determined to be domestically acquired.
Plasmids are important vehicles that carry various mobile genetic elements (MGEs) and antimicrobial resistance genes.Bacteria that acquire these plasmids may selectively persist and adapt under selective pressure.MGEs such as IS26/IS6100 (as observed in mph gene insertion) have driven plasmid acquired antibiotic resistance in incFII plasmids 19,20 .Cluster I isolates also possessed dfrA17-aadA5 and sul1 genes, which were present in the pKSR100 plasmids circulating in certain sublineages of S. flexneri 2a, 3a and S. sonnei that were predominantly in males or known MSM clusters during 2010-2018 1,4,5,7,14 .
Cluster III plasmids had a different gene arrangement including antibiotic resistance genes.A single nonclustering isolate had a pKSR100 plasmid that was also carrying additional genes (bla-CTX-M15 and qnrS) rendering this isolate resistant to AZM, CIP and CRO.Recovery of plasmid-borne resistance to AZM, CIP and CRO is alarming as these are three primary antibiotics used to treat shigellosis.It is more concerning as this multidrug resistance was found on a pKSR100 like plasmid which has been reported as widely disseminated in multiple MSM outbreaks.
As a high percentage of AZM R cases were phylogenetically closely related males without travel history, one may speculate them as MSM, though our study did not capture sexual behaviour and route of disease exposure or acquisition.Resistance to AZM was observed to be highest in 2017 and incidentally none of cases in 2017 report international travel and 94.1% cases in 2017 were males.Multivariate regression analysis showed that compounded effect of travel, gender and year of recovery together were important to AZM susceptibility outcome (multinom Wald's test P value < 0.05).This suggests that S. flexneri 1b isolates circulating locally in Ontario were AZM R and travel-related cases and cases among females in years 2016 and 2018 introduced AZM S isolates.We acknowledge that our study is limited to a time period of 3 years as trend related to travel and importation of cases change over time.However, 71.4% (5/7) of isolates in cluster II were travel-related suggesting that this lineage may be imported.Contact history between patients was also unknown.We found that 13.8% (9/65) of the isolates were non-susceptible to CIP, and 66% (6/9) were recovered from returning travellers from India or Cuba.While the importation of non-susceptible CIP Shigella spp.from South Asia has been reported many times 1,3 , and to our knowledge, there are limited reports from other high-income countries of importation of non-susceptible CIP shigellosis cases from the Caribbean region.Of concern, we observe importation of CIP intermediate-resistant cases from the Caribbean which is historically not associated to be high-risk region to acquire CIP resistant Shigella, indicating a possible shift in the epidemiology of global lineages.
Using WGS, we established that a dissemination of AZM R S. flexneri 1b bearing plasmids with a high percentage of homology to pKSR100 was prevalent in Ontario in the years 2016 to 2018.Multidrug-resistant shigellosis among men and the MSM population has been recorded world-wide 1 .Our data concurs with previous data showcasing antimicrobial resistance is increasing among male shigellosis cases [1][2][3][4][5][6][7][8][9][10] .There was limited data highlighting increased resistance from Ontario among Shigella isolates and therefore, our study is important in imparting awareness to high percentage antimicrobial resistance for improving treatment strategies.By bringing awareness to disease and high antimicrobial resistance rates and risks-groups, transmission events and outbreaks in future may decrease.Susceptibility testing for each shigellosis case is urged.As our study suggests that AZM R serotype 1b isolates are circulating locally in Ontario, future use of AZM to treat shigellosis needs re-evaluation.More comprehensive data showcasing antimicrobial-susceptibility trends for different Shigella spp.and their serotypes and sub-serotypes is needed.
Our study highlights the importance of continued surveillance of Shigella spp.The widespread prevalence of drug-resistant Shigella spp.necessitates public health action for increasing awareness, especially among high risk communities such as MSM and travellers.Spreading awareness of high percentage of azithromycin resistance among males will help clinicians in treatment and management of patients.

Shigellosis cases
Shigellosis is a notifiable disease in Ontario.All positive isolates are reported to 34 local Public Health Units (PHUs) that apply the provincial surveillance case definition to identify cases for provincial reporting using the web-based iPHIS.Data on case demographics (full name, health number, age at onset, gender, and PHU of residence) and travel status were extracted from iPHIS for all laboratory confirmed shigellosis cases reported from January 1, 2016 to December 31, 2018.Travel associated cases were defined as those reporting any international travel history up to 4 days prior to onset of symptoms.For our study, we included five PHUs covering regions of Toronto, Peel, Halton, York and Durham as GTA.

Bacterial isolates
As the provincial reference laboratory, PHO's laboratory confirms species identification including serotype for all Shigella isolates using an in-house protocol that ultilises culture, biochemical assays and serotype assays based on previously published protocols [21][22][23] .Demographic data associated with each isolate including age, date of birth, full name, gender, reporting public health unit (PHU) and unique health number and laboratory values including serotype and MIC were extracted from the laboratory information system.

Serotyping and susceptibility testing
Shigella isolates at the PHO laboratory are serotyped using anti-sera raised to specific surface antigen (O) as per an in-house protocol similar to previously published documents [21][22][23] .In addition, all isolates undergo susceptibility testing using agar dilution method following the Clinical and Laboratory Standards Institute (CLSI) guidelines 24 .As recommended by CLSI, E. coli ATCC25922 and Pseudomonas aeruginosa ATCC 27853 served as quality control strains in each test panel.To determine exact MIC, susceptibility testing for AZM and CRO was repeated using agar dilution method.The MIC results for all susceptibility testing were interpreted as per CLSI clinical breakpoints published in M100-S28 25 .For the purpose of this study, we present the results for antibiotics that are commonly used to treat shigellosis including AZM, CRO, CIP, SXT and AMP.For AZM, MIC values of ≤ 8 mg/L were interpreted as susceptible and ≥ 16 mg/L were interpreted as resistant in accordance to 28th edition of CLSI 25 .

WGS and bioinformatics analysis
Short-read sequencing All S. flexneri 1b (n = 65) were subjected to WGS.Genomic DNA for Illumina sequencing was extracted using QIAamp DNA mini kit (Qiagen Inc., Germany).Extracted nucleic acids were quantified using the Qubit fluorometer and quality of DNA was assessed by agarose gel electrophoresis and Nanodrop.Paired end libraries were prepared from the extracted DNA using the Nextera XT DNA Library Prep kit (Illumina Inc., San Diego, CA) per the manufacturer's instructions.Libraries were analysed using Agilent 2100 bioanalyser and pooled and sequenced on the Miseq platform using v3 chemistry (2 × 300 bp) in different runs.Post-completion of illumina run, bcl2fastq2 conversion software onboard the sequencer was used to demultiplex and remove barcodes and convert the base call files to fastq files (Illumina Inc., San Diego, CA).

Short-read mapping and phylogenetic analysis
Illumina reads from each isolate were subjected to quality check using FASTQC (http:// www.bioin forma tics.babra ham.ac.uk/ proje cts/ fastqc/) and Confindr 26 .Reads were trimmed using Trim-galore.For S. flexneri 1b, there was no complete reference.A randomly chosen S. flexneri serotype 1b isolate 002 was used as reference to improve the percentage of reads aligned and the de novo assembly was put into order by mauve 27 .The genome and assembly were further completed by combining Nanopore MinION long-reads and Illumina short-reads using unicycler (version 0.4.7) 28.Following assembly, a single round of error correction was performed by aligning the hybrid assembly to the de novo assembly and any variants on the hybrid were corrected to match the variant observed with the short-reads.Using sflex002 (GENBANK accession number XXXX) as reference, we obtained high-quality single nucleotide variant alignment using a custom approach that utilised an in-house pipeline and other open sourced bioinformatics tools.This in-house pipeline uses SMALT (http:// www.sanger.ac.uk/ scien ce/ tools/ smalt-0) to map reads against a reference genome and converts the smalt generated SAM files to BAM format using SAMtools and calls SNPs using Freebayes with a minimum coverage of 15, minimum base quality of 30 and minimum mapping quality of 30 to contain the variant and finally excludes repetitive regions.Additional variant confirmation using SAMtools mpileup with -BQ0 -d100000000 and bcftools view -cg is performed.Positions were filtered to be consistently called by both FreeBayes and SAMtools.MGEs (genomic islands, phages, transposons, insertion sequences, and other mobile elements) were screened using open-sourced tools like IslandViewer, Phaster, IS-finder, Prokka, Integron-Finder [29][30][31][32][33] .SNPs present in MGEs and areas of recombination identified using Gubbins were removed 34 .This resulted in multiple sequence alignments of 895 variant sites for all the isolates.Phylogenetic tree analysis was performed using IQ-TREE (version 1.6.2) with the GTR + ASC model with 5000 ultrafast bootstraps pseudo replicates to yield a maximum likelihood tree (ML Tree) 35 .SNP distance between isolates was computed using P-distance in MEGA X 36 .
Short-Reads were evaluated for antimicrobial resistance for genes presence and absence using SRST2 (version 0.1.8)with the curated ResFinder database (version 3.0) 39,40 .We initially screened SNPs for mutations that may cause resistance to AZM and CIP manually.Genetic determinants were further validated by using RGI (version 4.2) 41 .The ML-tree was visualized along with meta-data using Phandango 42 .

Long-read sequencing on MinION
In order to characterise AZM resistance carrying plasmids, select AZM R isolates representing the AZM R clusters were selected and subjected to long-read MinION sequencing to get complete plasmid sequence (two from clusters I [A-6726, sflex002] and III [A-10383, A-6817] and a single non-clustering isolate-sflex069  43 .Hybrid assemblies were derived by using the corresponding short-reads and long reads to yield complete plasmid sequences using Unicycler 28 (Supplementary Table S2).The completed plasmid assemblies were annotated using Prokka (version 1.3.3) 31.The homology of Ontario plasmids to pKSR100 plasmid (LN624486) was compared using BLAST and visualized using BRIG 44 .The presence of insertion sequences, integrons were searched using bioinformatics tools as described previously [31][32][33] .Whole genome assemblies from short-reads were generated using SPAdes (version 3.9) 45 .A combination of different approaches was performed to ascertain if Ontario short-reads had pKSR100.These included mapping short-reads against pKSR100 as reference as previously described 4,46 .We further evaluated if pKSR100 was present in short-reads by performing BLAST on de novo assemblies created using SPAdes with pKSR100 as reference 47 .The SPAdes assemblies were interrogated further using MOB-suite (version 1.4) which outputs plasmids sequences, replicon factor 48 .The output of MOB-suite also predicts inc Factor, relaxase typing and conjugation potential for each of the predicted plasmid contig aggregates.Pseudoplasmids were created by combining the BLASTn results and the aggregated contigs generated by MOB-suite 47,48 .All the pseudoplasmids were evaluated for antibiotic resistance genes as described earlier and compared with the sequence of pKSR100 (> 90% nucleotide match) to confirm presence of pKSR100 42,47,48 .

Data management and statistical analysis
Specimen duplicates and isolates received from the same case within 90 days were removed from the S. flexneri serotype 1b isolates dataset.S. flexneri cases with episode dates within 90 days of another case were also removed from the iPHIS case dataset.After applying these criteria, the two datasets were linked using deterministic matching based on health number, followed by supplemental probabilistic matching based on first and last names and date of birth.Sixty-five S. flexneri serotype 1b isolates received at PHO from January 1, 2016 to December 31st, 2018 and 63 shigellosis cases with episode dates from January 1, 2016 to December 31st, 2018 were successfully linked to iPHIS.The linked laboratory and iPHIS data was then de-identified to remove any personal information prior to further analysis.Statistical significance was assessed by Fisher's exact test, where a p value < 0.01 was deemed significant.We used the function 'multinom' in R (version 4.1) to perform multi-variate regression analysis and to calculate the Wald's p value.

Ethical review
Data linkage from iPHIS to PHO's laboratory database was done as part of regular surveillance.All personal identifiers linked to laboratory specimen numbers were removed.Whole genome sequencing was done on deidentified specimens as part of surveillance investigation and therefore did not require review by ethical review board.

Figure 2 .
Figure 2. Phylogenetic analysis of S. flexneri 1b and distribution of antibiotic resistant genes, Public health Units (PHUs), year of receipt, gender (M = male, F = female, U = other/unknown), travel-status (colored purple for international travel, light grey = if international travel was not indicated in iPHIS) and travel destination and azithromycin susceptibility (AZM; R = Resistant, S = Susceptible).SNPs were called against reference strain sflex002 belonging to cluster I (labelled in red font).Maximum Likelihood Tree (ML-Tree) was generated using IQ-TREE with the GTR + ASC method and ultrafast bootstrap over 5000 replicates.The antimicrobial resistant genes present are indicated as dark grey were detected by using SRST2 with ResFinder as the choice of database.Each cluster is labelled and colored with a different set of colors and SNP differences within each cluster are indicated.Bootstrap values of major nodes are shown.

Figure 3 .
Figure 3.Comparison of mph bearing representative plasmids from each of azithromycin-resistant (AZM R ) clusters.Complete plasmid sequences obtained by hybrid assembly (Unicycler) after combining long reads and short reads were annotated using Prokka.Insertion sequences were verified using the Prokka annotated file and the sites indicated on ISFinder website (https:// isfin der.bioto ul.fr/).Arrangement of the genes is shown as obtained in the hybrid assembly.(A) Comparison of mph bearing cluster I plasmid with psflex002 (recovered in 2017 in Ontario) as reference and showing synteny to pA-6726 (another cluster I isolate recovered in Ontario in 2018) and with plasmid pKSR100 (recovered from S. flexneri 3a isolates in Ontario in 2012).(B) Characteristics of cluster III plasmids recovered in 2018 and their synteny to pKSR100.The inner most ring represents the plasmid pA-6817 which serves as the backbone, followed by plasmid pA-10383 and plasmid pKSR100.(C) Figure showing the genetic arrangement of plasmid from isolate sflex069 (a non-clustering azithromycin resistant isolate that was multidrug-resistant to ceftriaxone, ciprofloxacin, trimethoprim-sulfamethoxazole and ampicillin).The figure also compares its relative homology to similar plasmids recovered from other 1b isolates, pKSR100 and pl183660 (a pKSR100 plasmid with mph (A) and bla CTX circulating in the UK among MSM during 2015-2016).The last outer ring for each plasmid in (A-C) shows the gene annotations.The mph (A) gene locus and the bla CTX-M-15 gene is colored green, while the class 1 integron int1-dfrA17-aadA5-sul1 gene loci is colored yellow in (A).The bla TEM-1B -aph (6′)-Id, aph (3′)-sul2 gene loci for plasmid p6817 in (B) is colored brown.Annotations in yellow color indicate the corresponding genes identified using integronFinder software for (A) and (C). https://doi.org/10.1038/s41598-023-36733-wwww.nature.com/scientificreports/