A query of the PulseNet database revealed 356 isolates related to the outbreak strain that had <15 core-genome multilocus sequence typing (MLST; cgMLST) allele differences (Table 1; Appendix 1 Table 1) (6). Of those, 302 isolates corresponded to human cases associated with 6 outbreaks spanning 3 years; dates of isolation ranged from September 27, 2016, to January 3, 2020. An additional 54 isolates were either clinical isolates not associated with a recognized outbreak (n = 14) or from environmental (n = 20), food (n = 8), or animal (n = 12) samples. Seven-gene MLST and Manning clade typing revealed all isolates were sequence type (ST) 11 and belonged to Manning clade 2 (Appendix 2). In silico PCR of the Shiga toxin (stx) genes revealed that all but 2 isolates contained stx2a, whereas 2 remaining isolates had no detectable stx genes. We generated a closed-reference genome, 2019C-3201 (Strain: PNUSAE020169; BioSample: SAMN10432148), using PacBio Sequel technology (https://www.pacb.com) and assembled with Flye version 2.6 (7). The sequence data assembled into a single complete chromosomal contig and 3 plasmids (Table 2).

We selected a subset of 245 isolates for further genomic analysis to more evenly sample across outbreaks and to reduce computational demands. Isolates were characterized by core genome MLST implemented in BioNumerics 7.6 (6) and high-quality single-nucleotide polymorphism (SNP; hqSNP) methods using Lyve-SET version 1.1.4f (9), using the chromosomal sequence of 2019C-3201 as a reference and the Lyve-SET presets for E. coli. Overall, hqSNP was more discriminatory, differentiating isolates by a median of 10 pairwise hqSNPs (0–39 SNPs), whereas cgMLST differentiated isolates by a median of 2 allele differences (0–8 alleles) (Table 1). This finding was foreseeable because hqSNP does not depend on a predefined scheme; therefore, intergenic SNPs between loci, multiple SNP differences within a given locus, or SNPs in loci not included in the cgMLST schema can result (9).

Figure 1

Figure 1. Tip-dated maximum clade credibility tree of 245 isolates of reoccurring Escherichia coli O157:H7 strain REPEXH02 linked to leafy greens–associated outbreaks, 2016–2019, generated in BEAST2 (https://www.beast2.org). Tips...

Time-tree analysis using BEAST version 2.6.3 (10) revealed the divergence of this strain into 2 clades that last shared a common ancestor around late 2015 (median December 19, 2015; 95% highest posterior density interval December 7, 2014–July 10, 2016) (Figure 1). After outbreak D in 2016, sequences corresponding to a given outbreak belonged to 1 of 2 clades; outbreaks B2 and C were associated with clade 1, and outbreaks A, B1, and B3 were associated with clade 2. Of note, outbreak A was traced to romaine lettuce from Salinas Valley, whereas traceback and sampling in outbreak B2 linked some illnesses to romaine lettuce from Santa Maria, California (https://www.fda.gov/food/outbreaks-foodborne-illness/investigation-summary-factors-potentially-contributing-contamination-romaine-lettuce-implicated-fall; https://www.fda.gov/food/outbreaks-foodborne-illness/outbreak-investigation-e-coli-romaine-salinas-california-november-2019). Lettuce from Salinas was not considered a source of any illnesses in outbreak B2. Environmental sampling in Santa Maria in 2019 yielded isolates clustering closely with outbreak B2 in the time tree.

We analyzed the closed reference sequence of 2019C-3201 using Prokka version 1.8 to enable SNP annotation (11). We examined output from Lyve-SET to determine the SNPs differentiating the 2 clades in our phylogenetic analysis (Appendix 1 Table 2). This work confirms a previous study reporting a nonsense mutation in the arsR gene, an arsenical resistance operon repressor (12). All clade 1 isolates in this study possess a G→A mutation resulting in a premature stop codon. This mutation could decrease the activity of this repressor and lead to constitutive expression of this operon. Agricultural soils and water sources can contain increased arsenic levels because of natural processes, industrial sources, or agricultural uses of arsenic, such as application of arsenic-containing herbicides, pesticides, or animal drugs (13). This mutation could provide an ecologic advantage in environments containing high levels of arsenic. This finding underscores the potential need to routinely screen enteric bacterial strains for heavy metal resistance determinants, as well as to consider heavy metal levels in soil as part of traceback investigations.

Figure 2

Figure 2. Annotated plasmids of reference genome 2019C-3201 of Escherichia coli O157:H7 containing clade-specific genomic features. A) p2019C-3201_1 annotated with prokka version 1.14.5 (yellow annotations) (11). Mapped list...

We further characterized isolates through assembly and annotation using Shovill-SPAdes version 1.0.9 and Prokka version 1.14.5 (11) and subsequent analysis in Roary version 3.11.2 (14) and scoary version 1.6.16 (15) to identify differences in the pangenome among isolates. We compared differentially distributed genes with the reference genome using BLASTn (https://blast.ncbi.nlm.nih.gov/Blast.cgi) to identify feature location (chromosome/plasmid). Roary/scoary analysis revealed a subset of clade 1 isolates with additional genomic content. A total of 156 genomic features had >90 sensitivity and >90 specificity to this subset of clade 1. Of those, 87 (56%) are on plasmid p2019C-3201_1, and 69 (44%) are on p2019C-3201_2 (Figure 2; Appendix 1 Tables 3, 4). Prokka-annotated features associated with p2019C-3201_1 (Figure 2; Appendix 1 Table 3) were predominantly genes encoding hypothetical proteins with unknown functions and common plasmid-associated genes. Annotated features associated with p2019C-3201_2 (Figure 2; Appendix 1 Table 4) were predominantly associated with conjugation and span a large portion of that plasmid. Additional work is necessary to characterize the role of these plasmids in clade 1. When visualizing the distribution of these clade 1–specific features alongside the maximum-clade credibility tree (Figure 1; Appendix 1 Table 5), it appears those features were acquired after clade 1 and clade 2 diverged. Given the geographic distribution of isolates, these features might be a result of adaptation to a particular niche or environment.

In summary, a specific strain of E. coli O157:H7 associated with leafy greens has been the source of ongoing enteric illness since late 2016. This strain is estimated to have emerged in late 2015 and consists of 2 clades with different geographic distributions, 1 of which has notable genomic features. After this analysis, an additional outbreak associated with this strain was detected in late 2020 in which a reported 40 infections occurred in 19 states; 20 persons were hospitalized, and 4 developed hemolytic uremic syndrome (https://www.cdc.gov/ecoli/2020/o157h7-10-20b/index.html). After that outbreak, no further outbreaks have been detected, and only a single clinical isolate associated with this strain has been identified by PulseNet. The Centers for Disease Control and Prevention has classified this strain as a reoccurring, emerging, or persistent (REP) strain (https://www.cdc.gov/ncezid/dfwed/outbreak-response/rep-strains.html) with the designation REPEXH02. REP strains represent a new paradigm in enteric molecular surveillance, distinct from discrete outbreaks where numerous cases occur in a relatively short time frame. Detailed genomic characterization of additional REP strains, using the types of approaches outlined in this study, is necessary to elucidate factors contributing to their emergence and persistence in specific environments.