IS1-related large-scale deletion of chromosomal regions harbouring the oxygen-insensitive nitroreductase gene nfsB causes nitrofurantoin heteroresistance in Escherichia coli

Nitrofurantoin is a broad-spectrum first-line antimicrobial used for managing uncomplicated urinary tract infection (UTI). Loss-of-function mutations in chromosomal genes nfsA, nfsB and ribE of Escherichia coli are known to reduce nitrofurantoin susceptibility. Here, we report the discovery of nitrofurantoin heteroresistance in E. coli clinical isolates and a novel genetic mechanism associated with this phenomenon. Subpopulations with lower nitrofurantoin susceptibility than major populations (hereafter, nitrofurantoin-resistant subpopulations) in two E. coli blood isolates (previously whole-genome sequenced) were identified using population analysis profiling. Each isolate was known to have a loss-of-function mutation in nfsA. From each isolate, four nitrofurantoin-resistant isolates were derived at a nitrofurantoin concentration of 32 mg l−1, and a comparator isolate was obtained without any nitrofurantoin exposure. Genomes of derived isolates were sequenced on Illumina and Nanopore MinION systems. Genetic variation between isolates was determined based on genome assemblies and read mapping. Nitrofurantoin minimum inhibitory concentrations (MICs) of both blood isolates were 64 mg l−1, with MICs of major nitrofurantoin-susceptible populations varying from 4 to 8 mg l−1. Two to 99 c.f.u. per million demonstrated growth at the nitrofurantoin concentration of 32 mg l−1, which is distinct from that of a homogeneously susceptible or resistant isolate. Derived nitrofurantoin-resistant isolates had 11–66 kb deletions in chromosomal regions harbouring nfsB, and all deletions were immediately adjacent to IS1-family insertion sequences. Our findings demonstrate that the IS1-associated large-scale genetic deletion is a hitherto unrecognized mechanism of nitrofurantoin heteroresistance and could compromise UTI management. Further, frequencies of resistant subpopulations from nitrofurantoin-heteroresistant isolates may challenge conventional nitrofurantoin susceptibility testing in clinical settings.


Introduction
Nitrofurantoin is a widely used first-line antimicrobial for urinary tract infection (UTI) treatment and prophylaxis [1].Reduced nitrofurantoin susceptibility in Escherichia coli has been associated with deleterious mutations in chromosomal genes nfsA, nfsB, and ribE [2], which encode key components of the oxygen-insensitive nitroreductase system, and with acquired multidrug efflux pump OqxAB [3].Despite low prevalence of nitrofurantoin resistance in E. coli clinical isolates (<7% in Europe) [4,5], concerns over increased prevalence of nitrofurantoin resistance are growing in England, where the consumption of nitrofurantoin had increased by 41% from 2017 to 2021 [6].
We previously identified three E. coli clinical strains with apparent signs of heteroresistance to nitrofurantoin [7], a phenomenon that may easily be overlooked in clinical antimicrobial susceptibility testing.We therefore set out to test nitrofurantoin heteroresistance and explore mechanisms of this phenotype for two strains that showed similar growth patterns.

E. coli strains
Two blood strains EC0026B (SAMEA104039660) and EC0880B (SAMEA104040147) [8] were tested to confirm nitrofurantoin heteroresistance.A nitrofurantoin-resistant strain IN09 [7] and nitrofurantoin-susceptible strain ATCC 25922 [9] were included as positive and quality controls, respectively, for population analysis profiling (PAP).Illumina whole-genome sequencing (WGS) reads of parental isolates EC0026B and EC0880B have been generated previously [8] and are available under accessions ERR3142418 and ERR3142524, respectively, in the ENA.We have shown that gene nfsA in each of EC0026B and EC0880B carries a lossof-function mutation [7].

Population analysis profiling
PAP was conducted on three separate occasions for each of the four strains (Figure 1).In brief, nitrofurantoin (N7878, Sigma-Aldrich, Germany) was dissolved in dimethylformamide to produce a stock solution, which was then infused into cation-adjusted Mueller-Hinton (MH2; 90922, Sigma-Aldrich) agar to form a two-fold gradient of nitrofurantoin concentrations (4-256 mg/L).One colony of each strain was inoculated into nitrofurantoin-free MH2 broth and grown aerobically overnight at 37 °C with shaking (180 RPM).Seven serial 10-fold dilutions of each broth culture were created with phosphate buffered saline (P4417, Sigma-Aldrich).
Ten-µL aliquots of the original broth culture and serial dilutions were spread in octants of MH2 agar plates following the nitrofurantoin gradient and incubated aerobically overnight at 35 °C.
For each agar plate, colony-forming units (CFUs) were counted in the octant where colonies were well separated, and thereby nitrofurantoin MIC of each strain was determined.A maximum non-inhibitory nitrofurantoin concentration of each strain was defined as the greatest concentration in which ≥70% cells grew when compared to the growth without nitrofurantoin.

Deriving isolates from strains EC0026B and EC0880B
Following PAP, four colonies with reduced nitrofurantoin susceptibility were randomly chosen (denoted by subscript 'R', e.g., EC0026BR1) from an agar plate containing 0.5×MIC nitrofurantoin as the exposure group for each strain.Then these isolates were evenly divided into two subgroups, inoculated into MH2 broth with or without nitrofurantoin (Figure 1), and aerobically incubated overnight at 35 °C.A colony of each strain was randomly selected from a nitrofurantoin-free agar plate as a comparator (denoted by subscript 'C', e.g., EC0026BC) to the exposure group and was grown in nitrofurantoin-free MH2 broth before DNA extraction.

DNA extraction and whole-genome sequencing
For each isolate derived from a parental isolate (EC0026B or EC0880B), genomic DNA was extracted from the broth culture using Proteinase K Solution, RNase A Solution, and Cell Lysis Solution (QIAGEN, Germany), and purified using GeneJET Genomic DNA Purification Kit (ThermoFisher Scientific, UK).The mean concentration of DNA in each extract was estimated from three reads obtained with Qubit dsDNA BR Assay Kit (ThermoFisher Scientific).
Extracted DNA of each isolate was aliquoted for WGS.Short-read sequencing (101 bp, paired-end) was conducted on an Illumina HiSeq 2500 system (Illumina, USA), and long-read sequencing was conducted on a MinION R9.4.1 flow cell (Oxford Nanopore Technologies, UK) for isolates having adequate yields of DNA.For MinION sequencing, DNA libraries were prepared using Rapid Barcoding Kit SQK-RBK004 (Oxford Nanopore Technologies), and base-calling, demultiplexing, and barcode trimming were conducted with Guppy v5.0.16 (community.nanoporetech.com/downloads) and its built-in high-accuracy model.For quality control, Illumina reads (101 bp, paired-end) were trimmed for a minimal base quality of Phred Q20 (10-bp sliding window) and filtered for a minimal length of 50 bp using Trimmomatic [11]; MinION reads were filtered for a minimal per-read average quality of Q10 and minimal length of 1 kbp using NanoFilt v2.8.0 [12].

Variant identification
For each strain, the complete chromosome sequence of the comparator isolate was used as a reference to identify chromosomal genetic variation in isolates R1-R4.Since chromosomes of isolates having MinION reads were fully assembled, strain-specific alignments of these chromosomes against references were performed using Minimap2 v2.24 [20].Chromosomal mutations and structural variants were identified from alignments using paftools.js(Minimap2).

E. coli isolates and genome assemblies
Five isolates (R1-R4, and C) were derived from each of strains EC0026B and EC0880B (Figure 1), providing 10 isolates altogether for WGS.Complete genomes were assembled for nine of these isolates (Supplementary Table 2).

Genetic variation in isolates derived from strain EC0026B
Isolates EC0026BR1-R4 showed different deletions of 11-20 kbp chromosomal regions harbouring nfsB (Figure 3a), and each deletion immediately followed the right imperfect inverted repeat (IRR) of an IS1-family insertion sequence IS1A, which interrupted insertion sequence IS150 (Figure 3b).The other end of the deletion seemed variable and occurred in coding sequences.For instance, the deletion in EC0026BR4 chromosome ended at 10 bp upstream of the left imperfect inverted repeat (IRL) of another IS1A, which interrupted gene vgrG.These two IS1A elements and the relevant reference sequence in the ISFinder database mutually differed by two nucleotide substitutions in the second open reading frame (insB) of the IS1A transposase gene [27].
Twenty-four copies of IS1A were found in chromosomes of EC0026BC, EC0026BR1, and EC0026BR4, and 23 copies were found in EC0026BR2 and EC0026BR3, which is consistent with deletions illustrated in Figure 3b.These IS1A elements demonstrated 99-100% nucleotide identities and 100% coverage to the reference IS1A sequence.The explicit copy number of IS1A could not be determined in the Unicycler short-read-only assembly of parental isolates EC0026B.Nevertheless, in the assembly graph, the mean and median depths of assembled segments from which the complete IS1A was recovered suggest 20-22 copies of this element with >98% nucleotide identities to the reference sequence.
As for other chromosomal genetic variation, 14 nucleotide substitutions and two indels were found in EC0026BR1, whereas no point mutation was found in EC0026BR2, EC0026BR3, or EC0026BR4 (Table 1).No genetic variation was detected between EC0026BC and EC0026B.

Genetic variation in isolates derived from strain EC0880B
Isolates EC0880BR1-R4 showed different large deletions (24-66 kbp) of chromosomal regions harbouring nfsB (Table 2).All these deletions ended immediately upstream of two identical copies of an IS1-family insertion sequence IS1X2 (Figure 3a, c).However, the start site of each deletion was variable.This IS1X2 element showed 98% nucleotide identity and 100% coverage to its reference sequence in the ISFinder database.
Unexpectedly, an additional chromosomal 4289-bp deletion was identified in EC0880BC when compared with that of parental isolate EC0880B despite absence of any point mutations.
This deleted region was present in EC0880BR1-R4 (Table 2) and immediately followed the IRR of the same IS1X2 as those related to the deletion of nfsB (Figure 3d).This deletion caused a loss of three genes (wzi, wza, and wzb) and a 142-bp truncation of the 5' end of wzc in the conserved capsule locus wzi-wza-wzb-wzc [28].
Twelve copies of IS1X2 were found in chromosomes of EC0880BC, EC0880BR2, EC0880BR3 and EC0880BR4, with 96-98% nucleotide identities and a 100% coverage to the relevant reference sequence.The precise copy number of IS1X2 could not be determined in the Unicycler short-read-only assembly of EC0880BR1, although 100% of this insertion sequence was recovered from five connected nodes (with 96-100% nucleotide identities to the reference sequence) in the assembly graph of this genome.These five nodes had a median read depth of 11 folds, suggesting 11 copies of IS1X2, which is consistent with the identified deletion in this isolate's chromosome (Figure 3c).Similarly, 11 copies of IS1X2 were estimated from the Unicycler short-read-only assembly graph of EC0880B with ≥95% nucleotide identities.
Interestingly, the other oxygen-insensitive nitroreductase gene nfsA in chromosomes of EC0880BC, EC0880BR2, EC0880BR3, and EC0880BR4, was interrupted by the same IS1X2 element that was associated with the deletion of nfsB, and locations of these interruptions were the same as that previously observed in EC0880B [7].The same interruption of nfsA was seen in EC0880BR1, albeit the sequence of the interruptive IS1X2 could not be reliably recovered from the short-read-only assembly graph of this genome.
This frequency of less susceptible cells in each strain was greater than a proposed minimum of 10 -7 for defining heteroresistance [29].Since no genetic variation in nfsA or ribE was identified between isolates of each strain, and neither EC0026B nor EC0880B produces the multidrugefflux pump OqxAB [7], we attribute the observed nitrofurantoin heteroresistance to the IS1associated deletion of nfsB regions.This kind of deletion probably occurred in a UTI-associated E. coli ST131 strain, whose nitrofurantoin MIC increased from 8 mg/L to 128 mg/L upon intermittent in vivo exposures to a therapeutic dose of nitrofurantoin in six months [30,31].
We also show that an IS1 element, IS1X2, was associated with the truncation of capsule locus wzi-wza-wzb-wzc in EC0880BC.All deletion events observed in our study are consistent with the abortive transposition model of IS1-mediated deletion, where a DNA duplex nicks at an end of IS1 and ligates to a target site on the same molecule, causing a loss of DNA between these two breakpoints [32].Moreover, IS1X2 had interrupted nfsA in parental isolate EC0880B, and we have previously elucidated that IS1R interrupted nfsA in two nitrofurantoin-resistant isolates [7].Taken together, IS1-associated deletion and interruption of nfsA or nfsB reduces nitrofurantoin susceptibility of E. coli and generally, may occur at other loci, resulting in significant genomic and possibly phenotypic changes.As such, we emphasise the importance of monitoring the prevalence and genomic locations of IS1 using long-read sequencing.
Nitrofurantoin heteroresistance in E. coli strain K-12 MG1655 was reported earlier, with a small subpopulation growing at a maximum nitrofurantoin concentration of 8 mg/L while the majority of cells had an MIC of 4 mg/L [33].Our PAP experiment discovered similar heteroresistance in E. coli strain ATCC 25922, which showed a maximum non-inhibitory nitrofurantoin concentration of 4 mg/L and a subpopulation growing at 8 mg/L with an average frequency of 9% (Figure 2a).Nonetheless, neither the level of the reduction in nitrofurantoin susceptibility nor the nitrofurantoin MIC of each strain is as great as those of strains EC0026B and EC0880B.Because both EC0026B and EC0880B had an MIC of 64 mg/L, which is the EUCAST clinical breakpoint for nitrofurantoin resistance (>64 mg/L), these strains may grow under therapeutic or prophylactic dosing of nitrofurantoin [34,35] and gain adaptive mutations under sub-MIC exposure [36] in urinary tract or gut.Moreover, the deletion of nfsB in these strains would cause an irreversible reduction in their nitrofurantoin susceptibility, paving the way to the emergence of nitrofurantoin-resistant cells if there were to be subsequent resistance mutations, such as those in nfsA, ribE, and the marA-marB intergenic region [36,37].
Importantly however, because the loss of nfsB reduces the reproduction rate of E. coli [30], the establishment of ΔnfsB mutants may only occur under selective pressure of nitrofurantoin.
We consider IS1-associated nitrofurantoin heteroresistance as a potential threat to the management of UTI for three reasons.First, the frequency of the subpopulation in strain EC0026B and EC0880B growing at 0.5× MIC (2-99 cells per million) was 10-100 times reported frequencies of spontaneous resistance-conferring point mutations in E. coli [38,39].
An early study shows a 30-2000 folds increase in the deletion frequency when IS1 is present [40].Second, IS1 provides E. coli with adaptive fitness.Specifically, in the absence of nitrofurantoin, the nitrofurantoin-susceptible major population (the wild type) can transmit and establish colonisation or infection as homogeneously susceptible E. coli does and likely outcompete nitrofurantoin-resistant E. coli that has lost nfsA and/or nfsB; and in the presence of nitrofurantoin, the infection/colonisation may persist with the subpopulation of reduced susceptibility.Third, nitrofurantoin-heteroresistant isolates may not be detected by routine susceptibility tests in clinical settings if the less-susceptible subpopulation has a low frequency.
To determine the mechanism of nitrofurantoin heteroresistance in E. coli, our study took advantage of hybrid de novo assemblies in accurately resolving complete bacterial genomes and identifying structural variation.Further work is needed to elucidate how IS1 is associated with genetic deletions and to compare fitness between mutants and parental isolates.
Heightened surveillance for prevalence and insertion sites of IS1 in clinical E. coli isolates will help us to understand how insertion sequences contribute to bacterial evolution.

Conclusions
E. coli can be stably heteroresistant to nitrofurantoin with a potential clinical significance, and such heteroresistance may be a precursor of nitrofurantoin resistance.IS1-family elements are associated with large-scale deletion of genomic regions containing nfsB and are potentially important agents of genomic and phenotypic changes in E. coli.

Figure 1 .
Figure 1.Workflow of population analysis profiling for each E. coli strain and preparation of isolates for DNA extraction.Isolate C represents a comparator controlled for the number of passages.Isolates R1-R4, which showed reduced nitrofurantoin susceptibility when compared to the majority of cells in the same broth culture, were randomly chosen at Step 3. Nitrofurantoin (NIT) concentrations in media are noted in gold.Icons were downloaded from Bioicons (bioicons.com)under CC-BY 3.0 Licence.

Figure 2 .
Figure 2. (a) Curves of population analysis profiling for test and control E. coli strains.Test strains EC0026B and EC0880B demonstrate nitrofurantoin (NIT) heteroresistance.Control strain IN09 is known to be NIT-resistant, and ATCC 25922 is susceptible.For each strain, the proportion of colony-forming units (CFUs) growing at each NIT concentration (4-256 mg/L) was calculated by dividing the CFU with the inoculum CFU (counted from the NIT-free agar plate).Based on biological triplicates, the mean proportion of CFUs ± the standard error of the mean is shown.CFUs below the detection limit are denoted by '0' on the Y axis.(b) Growth of test and control strains on cation-adjusted Mueller-Hinton agar plates with or without NIT on the same occasion.Test strains EC0026B and EC0880B showed subpopulations growing at 32 mg/L of NIT.

Figure 3 .
Figure 3. Genetic structures of deleted regions.(a) BRIG plots comparing genome assemblies against chromosomal sequences of isolates EC0026BC and EC0880BC.Sequences were aligned using megaBLAST, and matches were filtered for a minimum nucleotide identity of 90% and a minimum query coverage of 0.6%.(b-c) Genetic structures of nfsB-carrying regions deleted from the chromosomes of EC0026BC and EC0880BC, respectively, as indicated in (a).Start and end positions of each deletion are indicated by black arrows and red labels.Asterisks indicate variants of insertion sequences against reference sequences from the ISFinder database, and grey wide arrows indicate coding sequences of hypothetical proteins.(d) Genetic structure of the capsule-encoding region in the chromosome of EC0880BR2 and the partial deletion of this region in EC0880BC.

Table 1 .
Genetic variation in isolates EC0026BR1-R4 when compared to the chromosome of 26B.C. aa, amino acids.

Table 2 .
Genetic variation in isolates EC0880BR1-R4 when compared to the chromosome of EC0880BC.Of note, the 4.3-kbp 'insertion' resulted from the deletion of this region in the chromosome of EC0880BC.