Bacterial strains used in this study are listed in Table 1. Salmonella isolates (n=50) were previously isolated and analyzed from poultry and environmental samples by Brenner et al. (2020). Strains were paired-end sequenced on an Illumina NextSeq with a coverage of 15X. Reads were assembled with SPAdes v.3.10 and draft sequences assemblies annotated with the NCBI Prokaryotic Genome Annotation Pipeline (Bankevich et al., 2012). Salmonella sequences were deposited into Genbank under BioProject PRJNA224116.

The VTEC strains were isolated and analyzed previously in agricultural environments (irrigation water, sediment) by Nadya et al. (2016) and Falardeau et al. (2017) (Table 1). One clinical isolate from a previous foodborne outbreak, E. coli O157:H7 EDL933, was used as a reference strain because it was the first VTEC strain to be sequenced and has since been studied across many laboratories worldwide (Perna et al., 2001). VTEC strains beginning with “V-JF” were paired-end sequenced on an Illumina HiSeq with a coverage of 200X. Reads were assembled with SPAdes v.3.10 and draft sequence assemblies annotated with the NCBI Prokaryotic Genome Annotation Pipeline (PGAP). VTEC strains beginning with “SN” were subjected to paired-end sequencing on an Illumina MiSeq and assembled de novo with Spades v.3.1 (Bankevich et al., 2012). Draft sequence assemblies were annotated with Prokka v.1.10 (Seemann, 2014). VTEC sequences were deposited into Genbank under BioProject numbers PRJNA649237 and PRJNA287560.

PHASTER (PHAge Search Tool Enhanced Release) was used to identify intact prophages in the assembled draft genomes (Arndt et al., 2016). We used stringent criteria for prophage identification as described by Colavecchio et al. (2017b). Only prophages designated “intact” were used for further analysis. Hits to known phages were identified through PHASTER and subsequently confirmed with NCBI BLAST (Altschul et al., 1990).

All identified prophages were grouped if they were identified by BLAST best hit. If phages were identified by BLAST best hit and possessed 95–99.9% similarity over 50% query, they were designated variants and assigned a Roman numeral: I, II, III and so on. Nucleotide similarity was determined using the Clustal Omega alignment tool (Madeira et al., 2022), freely available at https://www.ebi.ac.uk/Tools/msa/clustalo/. Genes were annotated automatically with RAST (https://rast.nmpdr.org/) (Aziz et al., 2008).

The phylogenetic tree was constructed in R using the Ape package (Paradis and Schliep, 2019). Sequences were aligned using the ClustalW algorithm and the phylogenetic tree constructed using the Maximum-Likelihood method, employing 1,000 bootstrap replicates. The tree was imported into FigTree for additional rendering (Rambaut, 2007). Clusters were identified with ClusterPicker (Ragonnet-Cronin et al., 2013) using an inter-cluster threshold of 50% nucleotide identity (Fong et al., 2019). Comparative gene maps were constructed in R with the GeneplotR package (Guy et al., 2010) and Geneious Prime v.2022.0.2. Circular genomic visualizations and dotplot alignments were constructed with Geneious Prime v.2022.0.2 (www.geneious.com).

Bacterial CRISPR-Cas arrays were identified with CRISPRFinder (Grissa et al., 2007) using default parameters. Spacer sequences were assessed for homology to known phages using NCBI BLAST. The number of spacer sequences of VTEC and Salmonella were statistically compared with the Student's t-test in JMP version 11.1.1 (SAS Institute, Inc., Cary, NC, United States). A P-value of ≤0.05 was considered statistically significant.

Cumulatively, a total of 289 prophages were identified with bioinformatics analysis in our collection of Salmonella (n = 50) and VTEC (n = 68) strains (Supplementary Tables S1, S2). These were grouped into 130 unique (95% nucleotide identity over 50% query cover or lower) prophage sequences. Similar prophages within the same variant group were identified in diverse hosts (Supplementary Table S1) which may signify multiple phage/host interaction implications such as similar host ranges, common phage ancestry and recent or stable horizontal exchange and host acquisition (Hendrix et al., 1999). It is important to emphasize that phages were not induced experimentally, nevertheless, the intact prophages identified in this study would theoretically be capable of excision. As with all in silico prophage prediction tools, potential underrepresentation and misinterpretation of prophage diversity and abundance exists (i.e., split between two contigs in draft assemblies) and has been discussed extensively by Hurwitz et al. (2018). The prophages uncovered in this work provide a preliminary understanding of relationships in their respective hosts.

Prophages were clustered into 26 clades on the basis of 50% nucleotide similarity with 13 genomic singletons, highlighting the diversity in our subset of strains (Figure 1). Interestingly, prophages clustered together regardless of the bacterial genus, indicating that Salmonella and VTEC share genetically similar prophages with broad host ranges capable of inter-genus infectivity. Furthermore, these data suggest active dissemination of prophages amongst genetically distinct bacterial populations that may cohabitate similar niches, such as agricultural sites. We found prophages in Salmonella and VTEC that were most closely related to those of other hosts, such as Klebsiella, Haemophilus and Vibrio. The identification of genetically similar prophages amongst distinct bacteria may broaden the host range when designing future prophage therapies (e.g., engineering broad-range prophage induction agents) against foodborne pathogens (Hu et al., 2021). Furthermore, the discovery of distant hosts with closely-associated clinical outcomes may be attributed to the impact of the host environment in phage-host interactions.

We detected 13 unique prophages in 50 genomes of Salmonella (Supplementary Table S1). The abundance and type of prophage correlated with Salmonella serotype (i.e., prophage typing) (Figure 1, Thapa and Mohammed, 2019). Some prophages (NC 003315 I, NC 003315 II, NC 010393 II, NC 010393 III, NC 019488 I, NC 019488 II and NC 009542 II) were exclusively detected in strains of the same serotype (Figure 1). Polylysogeny (i.e., the carriage of more than one prophage) was also common; on average, 1.94 prophages were detected per Salmonella genome, which is lower than what has been reported (Bobay et al., 2013; Mottawea et al., 2018). Extremely divergent phages were also identified within the same genome; strain S01 contained five intact phages with an average nucleotide similarity of 51.5%. Similarly, S05 also contained five phages with an average similarity of 37.5%. These data indicate that phages with novel and diverse genetic materials are present in the genomes of our Salmonella collection. Compared to VTEC, we detected far fewer unique prophages in Salmonella. This may be due to several factors. For instance, E. coli may contain more preferable genomic integration sites and/or compatible genomic content (e.g., tRNAs, similar GC content) which has been demonstrated to favor phage infection and integration (Cardinale and Duffy, 2011). As part of the host immune response, superinfection exclusion systems also restrict superinfection by similar phages (Seed, 2015).

VTEC strain EDL933 was used in our prophage analyses as it is a well-characterized strain possessing genetically diverse prophages (Saile et al., 2016). Overall, 66 of 68 VTEC strains possessed at least one prophage. We did not identify any prophages in SN496 nor SN573, indicating these prophages may not be intact and/or defective. In a screening of 40 VTEC isolates by Zhang et al. (2020), 15.3 and 9.3% of prophages were incomplete and questionable, respectively. It currently remains unclear whether these prophages are inducible, however, Asadulghani et al. (2009) induced several defective prophages of VTEC O157 and found these phages maintained their virulence after the induction.

Clinical type strain EDL933 possessed 11 prophages, the most prophages of all strains tested. As the majority of our VTEC strains were isolated predominantly from the agricultural environment, these results suggest that these prophages may confer advantages to the host such that it allows for enhanced persistence in such environments.

On average, the incidence of prophage carriage in Salmonella was 1.94 per genome and 3.03 per VTEC genome (Figure 2). Prophage profiles clustered closely with serotype (i.e., number and type of phage correlated with bacterial serotype) (Figures 1, 2 and Supplementary Table S1), however, some prophages and their variants (e.g., NC_010392; Gifsy-1 and NC_010393; Gifsy-2) were not serotype-dependent and were identified in several strains representing different serotypes. Gifsy-1 and−2 -like phages are best known for contributing to the virulence of S. Typhimurium (Ho and Slauch, 2001), however in the present study they were identified in the genomes of S. Reading, S. Typhimurium, S. Enteritidis and VTEC. The variants differed greatly, with NC_010393 (I) possessing a genome size of ~16 kb, almost 50% smaller than variants (II) and (III). Further analysis revealed that these genes in variant (I) were mostly structural, while variants (II) and (III) also carried a slew of accessory genes for host invasion (e.g., sopE). Given that the annotation assigned mostly structural gene predictions, the carriage of variant (I) does not indicate the presence of any obvious host survival advantages but may provide an indirect mutualistic benefit if lysogeny is stably maintained. Diversity in Gifsy-1 was also greatly apparent (Figure 3) with variant (I) ~30% shorter than variant (II). An intact lysogeny cassette was not identified in this variant, although major structural and host recognition genes were identified (e.g., minor and major capsid proteins, portal proteins, tail fibers). Gifsy-1 was not only found in our subset of Salmonella strains but also in closely related to Gifsy-1 prophage in VTEC strain V-JF-010, which, to the best of our knowledge, has not been documented in E. coli before. Previous work has revealed that E. coli requires Salmonella OmpC protein for Gifsy prophage to bind (Ho and Slauch, 2001), however, upon homology analysis with NCBI BLAST no homologs were identified in strain V-JF-010. Further work is needed to elucidate the receptor of interest, and also to identify putative receptor sites in other bacterial genera, particularly in type phages such as Gifsy-like viruses that may confer virulence. The identification of such sites in taxonomically distinct species of bacteria is important for the design of prophage-engineered biocontrol strategies, an area of study still currently in its infancy (Cadieux et al., 2018).

Various genetic elements encoding an array of functions (i.e., host virulence, antimicrobial resistance) were identified in prophages of both Salmonella and VTEC (Table 2). Carriage of these virulence factors is problematic as these genes may be readily transduced throughout a bacterial population, potentially impacting host fitness such that it allows for greater survival in a variety of conditions. Phages may carry an array of genes that confer a host advantage (Colavecchio et al., 2017a; Gómez-Gómez et al., 2019; Wendling et al., 2020). For instance, prophages highly related to Gifsy-1, Gifsy-2 and Gifsy-3 are conserved within S. Typhimurium and contribute significantly to host virulence (Hiley et al., 2014). A previous study by Ross and Topp (2015) also observed the abundance of antibiotic-resistance elements in phages sourced from soil, suggesting that the agricultural environment may be a potent source of such genes.

A gene encoding putative tellurite and colicin resistance was functionally annotated in prophage SSU5 of SN141 and SN142 (Table 2). Interestingly, we also identified partitioning genes parA and parB. A high degree of homology to plasmid p4 was observed with demonstrated genomic rearrangements relative to each other (Figure 4). Plasmid p4 was originally isolated from an extended spectrum β-lactamase-producing strain of E. coli (Brolund et al., 2019; Accession: CP023851.1). Previous studies have observed the high degree of genomic similarity between SSU5 and plasmid pHCM2 (Kidgell et al., 2002; Kim et al., 2012), however, this variant of SSU5 appears to be most closely related to plasmid p4. SSU5 is 107,570 bp in length, with 130 open reading frames (Figure 5). Its putative roles as both a plasmid and phage may represent more avenues for transmission of genes that confer increased host fitness and warrants further investigation.

Tellurite resistance compounds have been found in several genera of bacteria, including E. coli, Salmonella, Shigella and Vibrio. Though uncommon, these compounds have been found in the chromosome or on plasmids (Chasteen et al., 2009). Experimentally, tellurite is used as a selective agent for the isolation of many pathogens, including E. coli O157:H7 (Turkovicova et al., 2016). Compared to the wild type strain, E. coli BL21 possessing the ter operon (terZABCDEF; responsible for conferring tellurite resistance) was able to grow in the presence of extremely high concentrations of tellurite (minimum inhibitory concentration of ~4 mM) (Turkovicova et al., 2016). In the environment, tellurium is relatively rare, however, its soluble salts such as potassium tellurite were used clinically as an antimicrobial agent in the past (Valkova et al., 2007). Additionally, it was shown previously that ter genes enhanced the ability of E. coli to survive in macrophages (Valkova et al., 2007).

It is interesting to note that the ter operon is functionally diverse and may encode for a variety of other bacterial phenotypes. The ter genes play a coordinated role in stress resistance and may offer resistance against a broad spectrum of agents, including colicins (small antimicrobial peptides (i.e., bacteriocins) produced by E. coli to kill non-host E. coli cells) (Jin et al., 2018) and phages by mounting a restriction or suicidal action upon phage infection (Anantharaman et al., 2012). Given the increasing demand for antibiotic alternatives in sectors such as food processing and clinical medicine, the dissemination of phage resistance genes in bacterial hosts can be problematic as it may limit the efficacy of phage-based treatments for host elimination (e.g., direct phage application). Because SSU5 may be able to act intracellularly as a plasmid and a phage, multiple avenues for dissemination throughout a bacterial population exist and may therefore amplify the spread of critically important genes such as those encoding antimicrobial resistance. Such phages may help to drive evolution of bacterial populations in certain environments.

The CRISPR loci is exploited for various analyses, such as microbial typing and tracking (Dion et al., 2021). Since spacer regions are hypervariable and provides historical information on phage resistance, we screened the spacer regions of identified CRISPR arrays with the aim to elucidate and characterize possible associations between spacer elements, prophage lysogeny and strain-specific differences in our dataset. All strains contained at least six unique spacers (Figure 6 and Supplementary Table S3). Salmonella contained significantly more spacers (an average of 27 regions) than VTEC (an average of 19 regions) (p <0.05). Previous work has revealed that in over >600 strains of Salmonella representing four serotypes, an average of 16 Class 1 and 20 Class 2 systems were identified (Shariat et al., 2015). More recent work also elucidated the diversity of spacers among various serotypes of Salmonella and identified a high number of spacers (440 and 330 unique spacers within 2221 and 2211 of CRISPR 1 and CRISPR2 arrays, respectively) in all Salmonella strains analyzed (Kushwaha et al., 2020). Many of the core CRISPR-associated (Cas) proteins that comprise the effector complex have evolved in line with the arms race against evolving phages (Koonin and Makarova, 2019). Thus, on the basis of Cas proteins and array organization, CRISPR defense systems have been categorized into two groups, Class 1 and Class 2 that are respectively subdivided into types I, III, IV and types II, V, VI (Makarova and Koonin, 2015). In a separate report surveying 100 strains of E. coli, 745 unique spacers were identified (Díez-Villaseñor et al., 2010). This is in stark contrast to a previous study by Touchon and Rocha (2010), where they observed no more than three CRISPRs in 51 genomes of Salmonella and Escherichia. These conflicting data may be because of several factors, not the least being differences in serotype, pathogenicity, and source of isolation of the strains tested. Zeng et al. (2017) observed that compared to clinical isolates (n = 17 spacers), food isolates of Cronobacter sakazakii possessed a significantly greater number of spacers in their CRISPR loci (n = 30). Correspondingly, foodborne isolates possessed less prophages (n = 2.81) compared to clinical isolates (n = 4.15). Enhanced evolution in CRISPR loci is important to consider when developing strategies such as phage cocktails for pathogen control because a key criterion for successful biocontrol is that resistance to phages is absent or delayed (Fong et al., 2019). It should be noted, however, that our bacterial strains of agricultural origin are draft genomes and may potentially underrepresent the true number of spacer regions. Nevertheless, our findings provide a preliminary basis for further probing of such genomes.

Spacers with significant homology to the identified prophages in this study were not observed, and we did not identify any self-targeting spacers. Surprisingly, only a small subset of strains (14%) encoded spacers with homology to known prophages, and these were only identified in VTEC (Supplementary Table S3). It should be noted that the scarcity of sequenced prophages in public databases may have supported this observation (Zeng et al., 2017). We also did not observe a significant correlation between spacer and prophage abundance in our strains (Figure 7). In silico work by Wang et al. (2020) also did not discern an association between these two elements in Bifidobacterium pseudocatenulatum. Despite their high degree of conservation, the role of CRISPR-Cas systems in prophage immunity has been questioned previously (Shariat et al., 2015), which has led to the hypothesis that these systems in Salmonella and VTEC may have alternative functions yet to be clarified (Touchon and Rocha, 2010; Shariat et al., 2015).