Whole-Genome Sequencing-Based Characterization of Listeria Isolates from Produce Packinghouses and Fresh-Cut Facilities Suggests Both Persistence and Reintroduction of Fully Virulent L. monocytogenes

Listeria isolates previously obtained from environmental samples collected in zones 2 to 4 from 16 produce operations, including 13 packinghouses and 3 fresh-cut facilities, were characterized using WGS. This isolate set included 276 isolates designated “representative” (see Materials and Methods for how the representative isolates were defined). The number of isolates sequenced per operation ranged from one (operations VT-I and VT-J) to 59 (operation CU-C) (Table S1). The majority of the isolates were classified as L. monocytogenes (n = 169), followed by L. seeligeri (n = 47), L. innocua (n = 37), L. welshimeri (n = 16), L. marthii (n = 6), and L. ivanovii (n = 1) (Table S1). The 169 LM isolates represented lineage I (n = 72), lineage II (n = 66), and lineage III (n = 31), and they could further be categorized into (i) 41 unique STs, with the largest groups (i.e., >7 isolates) being ST6 (n = 12), ST1 (n = 11), ST219 (n = 11), and ST824 (n = 10) as well as (ii) 36 unique CCs, with the largest groups (i.e., >7 isolates) being CC388 (n = 15), CC4 (n = 13), CC6 (n = 12), and CC1 (n = 11) (Table S2). A total of 42 LM isolates could not be categorized into an existing ST or CC (Table S2).

In addition to the 276 representative isolates, there were 14 isolates that had been sequenced but were found to be duplicates (i.e., isolates with the same sigB allelic type [AT] that were collected from the same site on the same date). These isolates will not be discussed in the remainder of the Results section. However, these isolates, compared against their respective “representative” counterparts, were found to have as many as 4 hqSNP differences, suggesting the potential for highly related but nonidentical isolates to exist in a given sample (Table S3).

The genomes of all of the LM isolates were screened for the presence of 31 selected genes with virulence and stress response functions (Fig. 1; Table S2). Additionally, inlA was screened for the presence of premature stop codons (PMSC). While all 169 of the LM isolates carried inlAB, 9 isolates had inlA PMSCs. These 9 isolates grouped into lineage II and one of two clonal complexes: CC199 (7/9) or CC9 (2/9).

For LIPI-3, 58/169 isolates had full matches (>96% similarity over >99% of gene length) to all 8 LIPI-3 genes (Fig. 1; Table S2). In addition, isolate FSL S11-0105 had a full match to 1 gene (>99% coverage with 100% similarity), partial matches (>55% coverage with 100% similarity) to 2 genes, and no matches to the remaining 5 genes from LIPI-3 (see Table S4 for metrics). Therefore, 59/169 isolates were characterized by the full or partial presence of LIPI-3, including 53/72 of the lineage I isolates, 6/66 of the lineage II isolates, and 0/31 of the lineage III isolates (Fig. 1; Table S2). The 59 LM isolates with the full or partial presence of LIPI-3 represented (i) 13 unique STs, with the predominant groups (i.e., >7 isolates) being ST6, ST1, and ST219, and (ii) 11 unique CCs, with the predominant groups (i.e., >7 isolates) being CC4, CC6, and CC1. For LIPI-4, 42/169 isolates had full matches to all 6 LIPI-4 genes (100% similarity over >99% of gene length), with an additional 2 isolates having full or partial matches to all genes (>99% similarity over >49% of gene length). Consequently, 44/169 isolates were characterized by the full or partial presence of LIPI-4, including 32/72 of the lineage I isolates, 5/66 of the lineage II isolates, and 7/31 of the lineage III isolates (Fig. 1; Table S2). The 44 LM isolates with the full or partial presence of LIPI-4 represented (i) 7 unique STs, with the predominant groups (i.e., >7 isolates) being ST219 and ST824, and (ii) 4 unique CCs, with the predominant groups (i.e., >7 isolates) being CC4 and CC388. Overall, 39/169 LM isolates had only LIPI-3, 24/169 isolates had only LIPI-4, and 20/169 isolates had both LIPI-3 and LIPI-4 (considering both complete and partial matches) (Fig. 1; Table S2). Therefore, 86/169 LM isolates had neither LIPI-3 nor LIPI-4.

The stress survival genes for which we screened included those present in the SSI-1 and SSI-2 operons. For SSI-1, 55/169 isolates had full matches (100% similarity over >99% of gene length) to all 5 SSI-1 genes, including 1/72 of the lineage I isolates, 35/66 of the lineage II isolates, and 19/31 of the lineage III isolates (Fig. 1; Table S2). The 55 isolates with the full presence of SSI-1 represented (i) 13 unique STs, with the predominant group being ST199, and (ii) 12 unique CCs, with the predominant group being CC199. However, 24 of these 55 isolates could not be assigned to an existing ST or CC. For SSI-2, 12/169 isolates had full matches to both SSI-2 genes, including 2/66 of the lineage II isolates and 10/31 of the lineage III isolates (Fig. 1; Table S2). These 12 isolates with the full presence of SSI-2 represented 3 unique STs (ST121, ST1513, ST262) and 3 unique CCs (CC121, CC262, and ST1513; ST1513 is the only ST within its respective CC). The majority of isolates with SSI-2 (8/12) could not be assigned to an existing ST or CC. No isolates had both the SSI-1 and SSI-2 operons (Fig. 1; Table S2).

The metal and detergent resistance genes screened included bcrABC, cadAC, emrE, and qacAH. Only bcrABC, a benzalkonium chloride resistance cassette, was found in the isolates studied here (Fig. 1, Table S2). Of the 169 LM isolates screened, 10 contained bcrABC. These 10 isolates all represented lineage II and were classified into one of three clonal complexes: CC155, CC199, and CC9. All 10 of the isolates with bcrABC also had the SSI-1 operon.

The WGS data were also used to perform hqSNP analyses in order to quantify the relatedness of all of the Listeria isolates characterized. While previous studies have suggested cutoffs of 20 SNPs (10) or 10 cgMLST alleles (29–31), outbreak-related isolates can differ by >50 hqSNPs (4). Therefore, 3 arbitrary hqSNP cutoffs were used to classify Listeria isolates as (i) related (<50 hqSNP differences), (ii) closely related (<20 hqSNP differences), and (iii) highly related (<10 hqSNP differences). Using the hqSNP cutoff of <50 hqSNPs allowed us to group the 276 Listeria isolates into 45 clusters (i.e., groups of isolates with <50 hqSNP differences), representing 135 isolates (Table 1) and 141 “singletons”, each of which represents a single isolate that was not related to any other isolate within <50 hqSNPs. The clusters ranged in size from 2 to 10 isolates (Table 1) and represented the species L. monocytogenes (32 clusters), L. innocua (5 clusters), L. seeligeri (5 clusters), and L. welshimeri (3 clusters). Although the majority of the clusters were comprised of isolates from a single operation, 7 clusters contained isolates from 2 operations, and 1 cluster contained isolates from 3 operations (Table 1). A total of 19 LM clusters were comprised solely of isolates obtained from the same operation on the same date but from different sites. An additional 16 clusters included isolates that were obtained from the same facility but on different dates (>60 days apart).

Table S1 provides a list of the metadata associated with each isolate, including a general site description and a cluster assignment for each isolate, if applicable. Of the 276 isolates, 27 isolates (from 10 operations) were collected from equipment frames, with 10/27 of those isolates being part of a cluster, as identified in Table 1. All 10 isolates were <10 hqSNP different from isolates from other sites in a given facility, suggesting a spread to and from equipment frames from other sites in these facilities. This includes 5/10 of these isolates that were <10 hqSNP from an isolate collected from a floor or drain site, demonstrating that Listeria found in zone 4 may be an indication that Listeria could be present in sites that are in closer proximity to food.

The analysis of WGS data by operation showed the re-isolation (detection on separate dates and at least 60 days apart) of highly related isolates (<10 hqSNP differences) in 7/16 operations (5 packinghouses and 2 fresh-cut facilities), representing 11/45 clusters. 1, 2, and 4 operations showed evidence for the re-isolation of highly related isolates representing 3, 2, and 1 clusters, respectively. 3 clusters included highly related isolates (<10 hqSNP differences) with dates spanning >1 year. These 3 clusters represented 2 fresh-cut facilities: CU-D (2 clusters) and CU-F (1 cluster). The 2 CU-D clusters (Clusters 10 and 18) (Table 1) were each comprised of isolates found during 2 sampling dates: August 2017 and September 2018. Importantly, no LM was detected during the 5 sample collection events between these 2 collection events. The isolates in Clusters 10 (3 isolates) and 18 (2 isolates) were obtained from trash cans, a trash cart, a pallet jack, and the floor. The locations of these isolation sources indicate that the isolates within Clusters 10 and 18 were dispersed throughout facility CU-D. It is possible that either or both of these strains survived in the operation over the course of the year or that these strains were reintroduced into the operation. In contrast, Cluster 7 (Facility CU-F) (Table 1), included 4 highly related LM isolates (3 to 7 hqSNPs) detected during 4 different sampling events (October 2017, December 2017, June 2018, and October 2018) (Fig. 2). The 4 isolates were from 4 sites in 3 different rooms: an overhead door, a floor/wall juncture, a forklift, and an equipment frame, suggesting a transmission of this strain between zone 2 sites and sites further removed from food contact surfaces. The isolates from the first and last sampling event (October 2017 and October 2018) were found in the same room, which was the product raw receiving area. The product was pre-washed in a neighboring packinghouse before being moved by a trailer to the raw receiving room, suggesting that isolates from this cluster could be reintroduced over time from the transfer trailer or the neighboring operation, neither of which was sampled during this study. Importantly, an investigation of the floor/wall juncture that was positive in December 2017 revealed that the wall had a covering that was used to prevent the original porous wall from getting wet. However, this covering created a sandwich in which moisture could penetrate, creating a potential niche. This was confirmed by follow-up testing done by the operation, which found positive results over repeated samplings of this site. While this wall was removed in January 2018, isolates that fell into Cluster 7 were found again in 2 subsequent sampling events, suggesting that the “Cluster 7 strain” continued to survive in this operation. Additionally, the “Cluster 7 strain” had the bcrABC cassette.

8 clusters (1, 4, 6, 13, 30, 34, 35, and 42) had highly related isolates (<10 hqSNPs) from the same operation obtained between 60 and 365 days apart. 3 of these clusters (1, 35, and 42) included highly related isolates from Packinghouse CU-C. Cluster 35 was comprised of 2 L. innocua isolates, while Cluster 42 was comprised of 2 L. seeligeri isolates. Cluster 1 isolates from Packinghouse CU-C (Fig. 3A) represented highly related LM isolates found during 2 different packing seasons and included the re-isolation of 2 highly related isolates (2 hqSNPs apart) from the same drain in April 2018 and April 2019 (subcluster 1a; isolates FSL S11-0289 and FSL S11-0472). This drain was located near the end of an approximately 29 m long, uncovered trench drain that spanned 2 rooms (i.e., the drain started in a room and then went underneath a wall to terminate in another room). That second room, which was kept at approximately 26°C (80°F) during production, was where the product was briefly dried. The third isolate within subcluster 1a (FSL S11-0553, 4 hqSNPs apart) was also from this same drain but was further “upstream” in the first room. Cluster 1 also included 3 isolates from Packinghouse CU-A that were collected in January 2018 (all classified into subcluster 1d). These findings are discussed in more detail below. The remaining 5 clusters (Clusters 4, 6, 13, 30, and 34) of highly related isolates (<10 hqSNPs) obtained from the same operation and between 60 and 365 days apart occurred in 4 packinghouses (Table 1).

Some of the isolates from Packinghouses CU-A, CU-B, and CU-C were from samples that were collected preproduction (i.e., after cleaning but before the start of production) (Table S1). Preproduction sampling was prompted by repeat positive samples being found at certain sites within the operation. Clusters 1, 3, 4, and 5 all had isolates from CU-A that were collected preproduction. Specifically, Clusters 1 and 5 included isolates from preproduction samples that were highly related (0 to 1 hqSNPs) to isolates from samples collected during the following production shift. Cluster 4 included isolates collected during production in October 2017 and January 2018. These isolates were highly related (1 hqSNP) to an isolate collected preproduction in December 2017. Most notably, Cluster 3 (Fig. 3B) had 2 highly related subclusters (3a and 3b; 5 hqSNP differences between clusters) (Fig. 3B) which each represented isolates obtained from a specific area of the facility. Each of these 2 subclusters included 3 isolates collected during production in October 2017 and 1 isolate collected preproduction in December 2017 (Fig. 3B). For each subcluster, all 4 isolates showed differences of 0 hqSNPs. These data suggest that these 2 strains (subclusters 3a and 3b) are surviving standard sanitation protocols and are persisting in the operation over time. This is supported by the fact that these isolates not only were found on separate sampling events in different months but also were found preproduction (which occurs after sanitation but before people and product begin moving through the operation). These isolates did not have the bcrABC cassette.

Among all 45 hqSNP clusters, 8 represented isolates from multiple operations (Clusters 1, 2, 4, 5, 23, 24, 38, and 40), with 2 clusters (Clusters 2 and 40) having isolates from separate operations with <20 hqSNP differences and 1 additional cluster (Cluster 38) having isolates from separate operations with <10 hqSNP differences. Among the 8 clusters with isolates from multiple operations, 6 represented LM clusters. Only 1 LM cluster (Cluster 2) with isolates from multiple operations had isolates from separate operations with <20 hqSNP differences. More specifically, subcluster 2b included 6 isolates from Packinghouse CU-B (all obtained in June 2018) and 1 isolate from Packinghouse CU-C (obtained in December 2017), and these differed from the CU-B isolates by as few as 17 hqSNPs (Fig. 3C). This represented the least hqSNP difference between LM isolates from 2 operations for this study. These 2 packinghouses were located in the same region of the United States and packed the same commodity. These findings suggest a potential common upstream source, as opposed to these strains persisting in each operation over time. A third operation (VT-A) also had 3 isolates that were part of this cluster (subcluster 2c), These isolates were collected in July and August 2017 and differed from all of the other isolates in cluster 2 by >25 hqSNPs. The remaining 5 LM clusters with isolates from multiple operations ranged from 30 to 49 hqSNP differences between isolates from separate operations, potentially representing Listeria that are more broadly distributed in the environment.

The 2 non-LM clusters with isolates from multiple operations both represented L. seeligeri, and both contained isolates from separate operations with <20 hqSNP differences (Table 1). Cluster 40 had 2 isolates (1 from Packinghouse CU-A and 1 from Fresh-cut Facility CU-F) that were 19 hqSNP apart. These two operations were located in the same region of the United States. Cluster 38 (Table 1) included 4 isolates with 3 identical L. seeligeri isolates (0 hqSNP difference) isolated from samples collected in 2 packinghouses (VT-A and VT-B) on the same day. Sample collection was performed by different personnel associated with each operation (Fig. 4), and the 4 isolates within this cluster were collected from 4 different sites. Although the authors acknowledge the possibility of cross-contamination within the lab, the positive-control used throughout the study was a strain of LM, whereas Cluster 38 was comprised of L. seeligeri isolates. In addition, all of the negative-controls were negative throughout the study. Moreover, the samples from the two packinghouses in this instance were processed by two different people (i.e., one person processed all of the samples collected from the first packinghouse, and a different person processed the samples from the second packinghouse). Further investigation revealed that in addition to being geographically close (<20 miles), it was not uncommon for these two packinghouses to employ the same workers, who resided at the same housing site for temporary, nonimmigrant workers, for different shifts (e.g., one employee may help with packing in the morning at one operation and with sanitation in the evening at the second operation), suggesting that shared employees could have played a role in cross-contamination between the operations. This possibility is supported by an L. seeligeri isolate being detected in a sample collected from the employee breakroom of one of the packinghouses (Packinghouse VT-B, isolate FSL S10-3508), although the isolate did not cluster with any of the other detected isolates.

Overall, our data suggest that a large proportion of the LM isolates obtained from the United States produce operations included in this study have the ability to cause human disease, as supported by (i) the classification of a considerable number of isolates into hypervirulent and outbreak-associated CCs, (ii) the frequent presence of LIPI-3 and LIPI-4, and (iii) the infrequent presence of strains with virulence attenuated inlA PMSCs. More specifically, the most frequent CC (CC388), which represented 15 of the 169 LM isolates in this study, was associated with a pork-related outbreak that sickened more than 200 people in Spain in 2019 (32, 33). The 3 next most frequent CCs from this study represented CC4 (n = 13), CC6 (n = 12), and CC1 (n = 11), which have also been associated with human illness (15, 34–38) and have all been suggested to represent hypervirulent LM (15). In addition to other outbreaks (38, 39), CC4 has been associated with a 2013 outbreak in Switzerland that was linked to salad (39). We also identified two CC7 isolates, and this CC was associated with the 2011 United States cantaloupe outbreak (39), indicating that CCs that have previously been associated with produce related outbreaks continue to be found in produce processing facilities.

Only 9/169 (5%) of the LM isolates in this study were found to have inlA PMSCs, which have been shown to result in virulence attenuation of LM isolates due to their low invasion efficiency into intestinal epithelial cells (40, 41). These isolates were categorized into two CCs that have historically been isolated from food sources rather than from clinical sources: CC199 and CC9 (15, 42). More specifically, a previous study of 6,633 isolates from France, including 2,584 clinical isolates, classified CC9 as a food-associated clone that rarely causes human disease (15). A different study of 300 LM isolates (117 clinical) from 5 continents included 6 isolates from CC199, which were all from food or environmental sources (42). Our finding that only 5% of the LM isolates in our study had inlA PMSCs is in contrast to previous studies that suggest that inlA PMSCs are common among isolates from RTE foods, processing plants, and retail environments (40, 41, 43–45). For example, Van Stelten et al. found that 45% of 502 isolates from RTE foods (i.e., bagged salads, fresh soft cheeses, soft ripened cheeses, smoked seafood, seafood, and deli salads and meats) carried an inlA PMSC, compared to 5% of human clinical isolates (n = 507) (43).

Several of the isolates included in this study had full or partial matches to the pathogenicity islands LIPI-3 (59/169) or LIPI-4 (44/169), including 20/169 isolates with matches to both. Both LIPI-3 and LIPI-4 have been shown to be associated with hypervirulence (15, 46). Overall, LIPI-3 was found in 35% of the LM isolates characterized in our study. By comparison, Chen et al. (12), found that 25% of the 102 LM isolates (recovered from 27,389 United States refrigerated RTE food samples) carried LIPI-3, while Kim et al. (14) found that 37% of 121 LM isolates (recovered from milk, milk filters, and milking equipment on bovine dairy farms) had LIPI-3. However, Hurley et al. (13) reported that only 10% of 100 LM lineage I and II isolates from food processing environments carried LIPI-3. More specifically, LIPI-3 was found among 74%, 9%, and 0% of the lineage I, II, and III LM isolates, respectively, that were characterized here. While LIPI-3 is generally thought to be restricted to lineage I (29), previous studies (13, 14) have found at least one lineage II isolate that had partial matches to LIPI-3. The frequencies of LIPI-3 per lineage found in our study (74%, 9%, and 0% of lineages I, II, and III, respectively) were higher than or consistent with those reported in previous studies. Kim et al. (14) found that 73%, 2%, and 0% of lineage I, II, and III isolates had LIPI-3, respectively. Hurley et al. (13) found that 39% and 1% of lineage I and II isolates had LIPI-3, respectively. Chen et al. (12) found that 51% of the lineage I isolates, but none of the lineage II or III isolates, had LIPI-3.

Overall, LIPI-4 was found in 26% of the LM isolates characterized in our study. This is higher than the percentages reported by Kim et al. (14), Chen et al. (12), and Hurley et al. (13), who found that 17%, 15%, and 1% of isolates had LIPI-4, respectively. More specifically, our study identified LIPI-4 in 44%, 8%, and 23% of lineage I, II, and III isolates, respectively. Kim et al. (14) found LIPI-4 among 32%, 0%, and 33% (1 isolate) of lineage I, II and III isolates, respectively. Chen et al. (12) found LIPI-4 among 30.6% of lineage I isolates and none of the lineage II or III isolates in their study. Hurley et al. (13) found LIPI-4 in 4% of lineage I isolates and in no lineage II isolates. The fact that we identified LIPI-4 among lineage II isolates is surprising and may require further follow-up studies (e.g., long range sequencing) to confirm the presence of LIPI-4, determine the genomic location of LIPI-4, and probe for possible lateral transfer events that may have introduced LIPI-4 into lineage II strains.

There was a surprisingly high proportion of isolates in this study (31/169) that were categorized as lineage III, which has been historically found to be underrepresented among isolates from human clinical cases (47). However, lineage III has also been considered to be underrepresented among isolates from foods and instead appears to be more commonly associated with food-production animal sources, particularly ruminants (47, 48), with a recent report also indicating a high prevalence of lineage III isolates among isolates collected from natural environments across the United States (49). This could suggest that a considerable proportion of produce-associated isolates could come from an animal source or a yet to be identified source of lineage III isolates that is common to farm animals and produce operations.

In total, there were slightly fewer isolates that had either SSI-1 (55/169) or SSI-2 (12/169), compared to the proportions of isolates with LIPI-3 or LIPI-4. Various studies have shown a diverse range of occurrence of SSI-1 among LM isolates characterized, ranging from 33% to 70% (12–14, 50–53). When comparing within lineages, our data showed that 1%, 53%, and 61% of lineage I, II, and III isolates had SSI-1, respectively. In comparison, Chen et al. (12) found SSI-1 in 35%, 80%, and 100% of lineage I, II, and III isolates, respectively, while Kim et al. (14) found SSI-1 in 47%, 43%, and 33% of lineage I, II, and III isolates, respectively. The 33% SSI-1 prevalence among the isolates characterized here is lower than what was reported by Chen et al. (12), Hurley et al. (13), and Kim et al. (14), who found that 57%, 51%, and 45% of isolates had SSI-1, respectively. For SSI-2, our data showed that 0%, 3%, and 32% of lineage I, II, and III isolates had SSI-2, respectively. Fewer studies appear to screen isolates for the SSI-2 operon, with most studies showing few isolates (0 to 5%) having the operon (12, 13, 52, 53) and Hurley et al. (13) showing 12% of isolates (16% of lineage II isolates) having SSI-2. Overall, compared to previous studies, our study found a lower frequency of LM isolates with SSI-1 (12–14, 50–53) but a higher frequency of LM isolates with SSI-2 (12, 13, 52, 53), possibly due to the fact that the isolates characterized here included a larger number of lineage III isolates.

Our study found that only a few isolates contained any of the selected genes we screened for from the “metal and detergent resistance” category (see Materials and Methods). The only genes detected in our isolates were those representing the bcrABC resistance cassette, which has been shown to confer reduced sensitivity to a quaternary ammonium (“quat”) compound called benzalkonium chloride, which is commonly used for sanitation in food production environments. bcrABC was found in 10/169 (6%) of the LM isolates studied here, including all 4 isolates from Cluster 7, which appear to have persisted in facility CU-F. The 10 isolates with bcrABC found here were all from lineage II (10/66, 15%), and the proportion of isolates with bcrABC is lower than what was found by Chen et al. (12), who found this cassette in 10/49 (20%) and 35/51 (69%) of lineage I and lineage II isolates, respectively. For the produce-specific isolates within that study (12), 2/14 (14%) and 9/14 (64%) of the lineage I and II isolates had bcrABC, respectively. A study investigating 100 LM isolates from three meat and vegetable processing facilities found that 19% of the isolates had the bcrABC cassette (13). A different study that characterized 15 produce-associated LM isolates in the United Kingdom found that 2/15 (13%) isolates had bcrABC. Our data suggest that bcrABC presence may be less common in isolates found in the produce-associated operations studied here, compared to previous studies. This could at least be partially due to the fact that a considerable proportion of the LM isolates characterized here (i.e., 91%) were obtained from packinghouses, which may be less likely to use quaternary ammonium compounds. Additionally, it is important to note that cadAC, emrE, and qacAH were not detected in any of the isolates studied here. Overall, our findings are consistent with those reported in previous studies (54, 55) that have not identified a strong association between the presence of specific “persistence” genes. While the findings to date could suggest that the establishment of persistence may include a strong element of chance (meaning persistence is likely to occur when an appropriate strain is introduced into a facility location that represents a potential niche where Listeria would not be removed by sanitation), further studies that use even larger isolate sets than those described here and tools such as genome-wide association to identify new genetic markers that are putatively associated with persistence would be valuable. In addition to the large sample sizes needed for these types of studies, a continued challenge with these types of studies will be classifying isolates as truly “sporadic”. Isolates may be misclassified as sporadic if they persist in locations that are difficult to sample or are not sampled for other reasons.

In addition to 141 isolates that did not fall into any hqSNP cluster, we also found that 19/45 clusters in this study (representing 42 isolates) were comprised of isolates from a single operation obtained on a single date but from different sites. Hence, the majority of Listeria or LM positive sites appear to be due to sporadic contamination, with some representing short-term Listeria spread within an operation, with contamination apparently controlled via standard cleaning and sanitation practices that were in place. However, our data showed the re-isolation (detection on separate dates that are at least 60 days apart) of highly related isolates (<10 hqSNP differences) in 7/16 operations (5 packinghouses and 2 fresh-cut facilities), suggesting that persistent contamination (or reintroduction, as discussed further below) is still frequent among United States produce operations (Table S5). Additionally, our data indicate that persistent Listeria contamination can occur in both packinghouses and fresh-cut facilities. This is consistent with previous studies that found that a significant proportion of food-associated operations show evidence for LM or Listeria persistence. For example, in a study of 9 small cheese processing facilities, 7 facilities showed evidence for Listeria spp. persistence (20). Similarly, in a study of 30 deli operations, 12 showed evidence for persistence (18). More specifically, we identified 6 LM clusters that had isolation dates spanning >1 year, providing evidence for long-term persistence. This includes 3 clusters (representing 2 fresh-cut facilities) that had isolates from a single operation that were highly related (<10 hqSNPs) and detected >1 year apart. Importantly, these 2 fresh-cut facilities were operated continuously, whereas the packinghouses in this study, from which data were collected over >1 year, were operated seasonally. Therefore, the packinghouses all had an “off-season” in which equipment was down and able to be disassembled, cleaned, and dried for an extended period. This could potentially explain why we did not detect highly related isolates among those collected >1 year apart in any of the participating packinghouse operations, which is consistent with the findings of a previous study that found low LM prevalence and no evidence for persistence in 2 seasonally operated crawfish processing facilities (56). Overall, our findings are also consistent with prior observations that LM strains may survive in operations for extended periods of time, as supported by a study that showed a processing plant that had a single strain persisting over at least 12 years (57). Interestingly, the fact that a single produce processing facility included 3 distinct LM strains that showed evidence for persistence (Packinghouse CU-C had 3 clusters of highly related isolates found between 60 and 365 [exclusive] days apart) suggests that some facilities may be more prone to allowing for the establishment of persistence. Based on observations of this operation, persistence may occur due to equipment with poor sanitary design or infrequent sanitation procedures, and these observations are in concordance with a previous review, which found that one of the two most common risk factors for persistence mentioned in the literature was equipment cleanability (58). However, we did not formally assess here which specific factors may have the greatest impact on the likelihood of persistence occurring in a given facility. While future studies on risk factors for persistence will be valuable, they will be challenging, as a large number of facilities would need to be enrolled.

While, as discussed above, the repeat isolation of closely related Listeria often is interpreted as providing evidence for persistence, these types of findings could also be due to reintroduction from outside or upstream sources that are persistently contaminated with a given Listeria strain. Interestingly, we found that 7 of the 45 hqSNP clusters in this study were comprised of isolates from 2 operations. An additional cluster (Cluster 2) included isolates from 3 operations, with all isolates within this cluster differing by <30 hqSNPs. Our findings that 8 clusters contained isolates from at least 2 operations demonstrate the potential for closely related isolates to be collected at different operations and provides evidence of the introduction of specific Listeria strains into multiple facilities from a common source. More specifically, 1 cluster showed LM isolates that were as few as 17 hqSNPs apart and were collected from 2 packinghouses. This could be due to a common source of raw materials obtained from facilities or fields that harbor this strain or could represent a widespread presence of isolates representing this given hqSNP cluster in the environment. Interestingly, we also identified L. seeligeri isolates that had 0 hqSNP differences in 2 different packinghouses that shared employees, despite being separately owned and operated. This provided a potential mechanism for cross-contamination between facilities. In concordance with our findings, other studies have also identified closely related Listeria spp. isolates from different facilities, although it is important to note that different SNP-based data analysis approaches may not be directly comparable (59). For example, 1 study showed that LM isolates with <10 SNP differences were isolated from multiple delis in separate states in the United States (22), suggesting the introduction of closely related LM into multiple facilities, likely from an upstream source, such as a common supplier. Another study, which investigated isolates from a cold-smoked salmon facility, showed that an isolate from a different cold-smoked salmon facility was within 11 to 23 SNPs of the other isolates within the cluster (8). Other studies have shown instances in which closely related isolates have been associated with separate operations, such as an LM strain that was tied to two ice cream production facilities, one of which purchased ingredients from the other (25). Our findings further support the importance of using WGS data in combination with metadata to help differentiate between re-contamination and persistence. For example, the repeat isolation of closely related Listeria isolates after sanitation and preoperation in locations with limited traffic during off hours (e.g., production rooms) supports persistence, while the repeat isolation of closely related Listeria isolates only during operations and from sites close to potential introduction routes (e.g., a receiving dock), more likely indicates re-contamination. In addition, advanced WGS data analysis, including the construction of tip-dated phylogenies (see Harrand et al. for an example [60]) can provide information on the most recent common ancestor (MRCA) of closely related isolates. In this case, an MRCA that predates the construction date of a facility could also suggest a reintroduction rather than persistence (particularly if supported by metadata). Importantly, if closely related isolates are found in different facilities, this may also suggest contamination from higher-up in the food supply chain (e.g., agricultural water, fields, field equipment) or from common employees, for example. Re-introduction may be more likely in supply chains where no kill steps (e.g., heat treatment) are applied (e.g., fresh produce), which would facilitate survival throughout the supply chain. Our study specifically shows that the availability of larger WGS data sets that are comprised of isolates from multiple facilities, along with detailed root cause and epidemiological approaches, will help to differentiate persistence from reintroduction or cross-contamination. This will help facilities more rapidly address the true root causes of contamination events.

Listeria isolates were obtained from two previous studies (61, 62) that investigated Listeria prevalence in 13 produce packinghouses (VT-A, VT-B, VT-C, VT-D, VT-E, VT-F, VT-G, VT-H, VT-I, VT-J, CU-A, CU-B, CU-C) and 3 fresh-cut facilities (CU-D, CU-E, CU-F) (61, 62). These two studies collected a total of 4,152 environmental samples, of which 217 were positive for Listeria. From these 217 Listeria-positive samples, a total of 679 Listeria isolates were collected. Additionally, a third study collected 156 environmental samples from 3 of the operations that participated in the Sullivan et al. (62) study (i.e., operations CU-A, CU-C, CU-E) approximately 1 year after the original study, resulting in 21 Listeria-positive samples and 155 Listeria isolates. Only LM isolates from the third study were included in the study reported here. To identify unique representative isolates from each sample, sigB allelic typing was performed on all of the isolates as previously described (63). One isolate of each sigB allelic type that was present in a given sample was considered “representative”. Therefore, a single positive sample could have more than one representative isolate if these isolates showed distinct sigB allelic types. Additionally, different representative isolates could have the same sigB allelic type, as long as they originated from different samples. These “representative” isolates (n = 276) were characterized via whole-genome sequencing, as detailed below. In certain instances, such as when preliminary sequencing data were unavailable at the time of strain selection, the whole-genome sequencing of an isolate that was ultimately designated representative was not performed. Operation code designations for each isolate are consistent with those of the original studies. However, to differentiate between studies, the operations studied by Sullivan et al. were given a “CU” prefix (i.e., Cornell University), and the operations studied by Estrada et al. were given a “VT” prefix (i.e., Virginia Tech) (61, 62).

The total DNA was extracted from representative isolates for whole-genome sequencing as previously described (22). Sequencing was performed using a HiSeq 2500 (Illumina, Inc., San Diego, CA, United States) with a maximum read length of 2 × 150 bp. Due to the sensitivity of the research reported here and to protect the data privacy of the participating operations, the sequences were not uploaded to NCBI. The assemblies, however, are publicly available on the Cornell University eCommons repository: https://doi.org/10.7298/74sp-fg52 (64). Genomes were assembled as previously described (22). Briefly, Trimmomatic (v 0.39) was used to trim and filter the raw reads (65), which were then evaluated based on quality using FastQC (v 0.11) (66). The trimmed reads were assembled using SPAdes (v 3.13) (67) with k-mer sizes of 21, 33, 55, and 77 bp. All of the isolates, including the non-LM Listeria species, were sorted into groups based on sigB allelic type before performing a reference-free SNP analysis using kSNP3 (v 3.1) (68) to determine the clusters for the hqSNP analysis. Using the results of the kSNP3 analysis, isolates that were <100 SNPs apart were then analyzed using the CFSAN SNP pipeline (v1.0.1) (69) to determine high-quality single nucleotide polymorphisms (hqSNPs). The reference assembly for the United States Food and Drug Administration CFSAN SNP pipeline analysis was selected to represent the isolate with the maximum Q value within each cluster. The Q value was obtained via the following equation.

If, after the hqSNP analysis, clusters had isolates that were >50 hqSNPs away from all of the other isolates, the isolate was removed, and the cluster was reanalyzed. Therefore, if a cluster had 2 or more subclusters that were >50 hqSNPs apart, the cluster was split into multiple subclusters for analysis. Maximum likelihood phylogenetic trees based on the hqSNP data for select clusters were created using RAxML (v 8.2.12) (70) with 1,000 bootstraps and the GTRCAT nucleotide substitution model.

To assess whether virulence and stress-response genes were present in the genomes sequenced here, allele sequences were downloaded from the Pasteur Institute′s BIGs-LM database from the following schemes: (i) Virulence (all 8 genes in LIPI-3 and all 6 genes in LIPI-4), (ii) Metal & Detergent Resistance (bcrABC, cadAC, emrE, and qacAH), and (iii) Stress Islands (all 5 and 2 genes in SSI-1 and SSI-2, respectively) (29, 71–73). Additionally, the classical 7-gene MLST scheme (abcZ, bglA, cat, dapE, dat, idh, ihkA) was downloaded from the Pasteur Institute’s BIGs-LM database for the in silico MLST of the isolates using the WGS data. A database was created of the seven MLST genes, which were used to assign each LM isolate to a sequence type (ST) and a clonal complex (CC). The inlA and inlB alleles were also downloaded from the Pasteur Institute’s BIGs-LM database, and a BLAST database was created to determine whether each gene was present or absent. MEGA X (v 10.1.7) (74) was used to identify isolates with PMSCs in inlA. The presence or absence of genes and the presence of inlA PMSCs were annotated, using iTOL (v1.0) (75), in a phylogenetic tree that was generated from the kSNP3 analysis of all of the LM isolates.