Main

Bacteriophages (or simply, phages) are an extraordinarily diverse and ubiquitous class of viruses that pose a nearly constant threat to bacteria. Phages are the most abundant biological entity on the planet, with estimates of 1031 particles that drive the daily turnover of ~20% of all bacteria in some environments1,2. Bacteria and their viral predators are locked in a perpetual coevolutionary battle, leading to the emergence of sophisticated mechanisms by which phage manipulate and exploit their hosts, and an equally diverse set of bacterial immune mechanisms collectively referred to as anti-phage defence systems3. These immunity systems include both innate mechanisms, such as restriction-modification systems, and adaptive mechanisms, such as CRISPR-Cas. Recent studies have begun to identify many new defence systems, but the full inventory probably extends well beyond what is currently defined.

Identifying additional anti-phage defence systems promises to provide new insight into the ancient coevolutionary conflict between viruses and their hosts. Recent work has found that many defence systems have homologues with similar function in eukaryotic innate immunity, indicating a potentially ancient origin and cross-kingdom conservation of many immune systems4,5,6. Additionally, previous studies of anti-phage defence have produced precision molecular tools such as CRISPR and restriction enzymes, so the discovery of new immune mechanisms may enable new tools for manipulating cells and genomes. Finally, there is growing interest in using phages to treat antibiotic-resistant bacterial infections and to manipulate microbiomes7,8,9. A more complete understanding of the diverse mechanisms by which bacteria defend themselves may be critical for these endeavours10.

Multiple groups have previously used computational methods to identify uncharacterized defence systems on the basis of the observation that anti-phage defence systems often cluster in bacterial genomes in high density, forming so-called ‘defence islands’11,12,13. However, not all defence systems may be detectable in these guilt-by-association approaches. For systems that are rare or not widely conserved, it may be difficult to detect enrichment within defence islands, and not all defence systems necessarily associate with defence islands. Additionally, candidates identified computationally must be expressed in a model laboratory organism and then tested against a panel of phages. Some systems may not work in a heterologous host or protect against the phages examined, and demonstrating that a given system provides defence in its native context is typically not tested or even possible. We reasoned that an experimental selection scheme to uncover antiviral proteins may reveal new insights into bacterial immunity, including identifying defence systems that remain uncharacterized and revealing the relative frequency of the different genomic contexts of these bacterial immune systems (that is, defence islands, mobile genetic elements (MGEs) or other contexts). To this end, we took a functional metagenomic approach to map the range of defence systems in the Escherichia coli pangenome (Fig. 1). In addition to revealing previously uncharacterized phage defence systems in E. coli, our work also provides a robust screening methodology that can now be adapted to systematically identify phage defence mechanisms in virtually any bacterial genome or metagenomic sample.

Fig. 1: Selection strategy for identifying phage defence systems.
figure 1

A fosmid library of random ~40 kb fragments of genomic DNA from 71 E. coli strains was transformed into an E. coli K12 host and then challenged with three different phages. Survivors were isolated and fragments mapped to their genome sequence. After eliminating duplicates, clones affecting adsorption, and clones harbouring restriction-modification or known defence systems, the unique fosmids corresponding to each phage selection were used to construct plasmid libraries, which were subjected to a second selection. Surviving clones were deep-sequenced and candidate defence loci pinpointed by mapping sub-library reads to genome sequences of the original fosmid inserts.

Full size image

Results

Identification of novel anti-phage defence systems

To sample the immune landscape of E. coli, we collected a diverse set of wild isolate strains from the ECOR collection as well as 19 clinical isolates14 (71 strains in total), all with available draft genome sequences (Extended Data Fig. 1). The ECOR collection is a set of strains curated to span the phylogenetic diversity of the species15. Together, the 71 strains collected encode 21,149 unique gene clusters, of which >10,000 exist in only one or two strains (Extended Data Fig. 1). From genomic DNA, we constructed a 100x-coverage library of fosmids, each harbouring a ~40 kb genome fragment, in EPI300, a derivative of E. coli K12. We used large-insert fosmids to minimize the size of the library, to include potentially large defence systems, and because the copy number is maintained around one, minimizing false positives due to overexpression.

Many anti-phage defence systems work by an abortive infection (Abi) mechanism in which an infected cell sacrifices its viability to prevent phage replication and thereby protect uninfected cells in a population16. Thus, it is impossible to directly select for clones containing an Abi-based defence system because the infected cell dies. Instead, we used a selection strategy, historically known as the ‘tab’ (T4 abortive) method17. In brief, cells harbouring the gDNA library are mixed with phage in a structured medium (soft agar) at varying concentrations of phage (Fig. 2a). At intermediate phage concentrations, individual clones from the library can grow and form small populations before encountering a phage particle. Any such micro-colonies that harbour an Abi defence system can be infected, but the initially infected cell will die without producing progeny phage, enabling the rest of the population to survive and produce a colony. Thus, our screening approach allows the identification of both Abi and conventional defence systems.

Using this general strategy, we challenged cells harbouring the fosmid library with three lytic phages: T4, λvir and T7, each representing a major class of Caudovirales, the tailed bacteriophages (Fig. 1). From each selection, we isolated approximately 90 surviving colonies and then sequenced the ends of the vector insert in each clone to identify the genomic region and strain of origin of each fragment. For each positive clone, we then measured the efficiency of plaquing (EOP)—the ratio of plaques formed by a phage plated on the positive clone to plaques formed by a phage-sensitive control strain (Fig. 2b). In some cases, these EOP assays revealed nearly complete protection even at phage titres of >109 p.f.u. ml−1. Upon initial sequencing, representative clones with this phenotype consistently encoded lipopolysaccharide (LPS) or capsule biosynthesis genes. This phenotype is consistent with a loss of adsorption that results from modification of the core surface properties of the cell (Fig. 2b and Extended Data Fig. 2a). Other positive clones produced protection accompanied by a high rate of escape plaques, a phenotype seen with many restriction-modification (RM) systems (Fig. 2b). Escape plaques can arise with many defence systems but are well known to arise at high frequencies for RM systems due to epigenetic escape18. Upon sequencing, the clones exhibiting high-frequency escape plaques indeed encoded RM systems. Thus, clones with these two phenotypes were excluded from further analyses (Methods). From 257 initial clones, 117 and 9 were eliminated as probably resulting from changes in cell surface properties affecting adsorption or RM-mediated defence, respectively. After also accounting for redundancy in the remaining 131 clones, we had 43 clones that we hypothesized to encode novel defence systems.

Fig. 2: Identification of phage defence systems.
figure 2

a, Left: schematic of the ‘tab’ selection method. At intermediate concentrations of phage, tab selection facilitates the survival of cells with either abortive infection or direct defences. Right: examples of T4 selection plates for cells containing the fosmid library or an empty vector control. b, Left: tenfold dilutions of λvir phages on lawns of a sample of 15 positive clones from the λvir screen. Multiple phenotypes were observed, including reduction of plaquing with individual escape plaques indicative of a restriction-modification system; no lysis at any concentration of phage typically reflecting a loss of adsorption; or reduction of plaquing, generally indicative of a phage defence system. Right: examples of fosmid inserts corresponding to exemplar phenotypes in b (left), with relevant genes coloured. c, Examples of read coverage (100 bp moving average) from deep sequencing of sub-libraries generated from positive fosmid clones, with maxima delineating defence system candidates. Genes were coloured or shaded as indicated at the top. d, Summary counts of defence systems identified.

Full size image

To pinpoint possible defence systems, we pooled the three sets of remaining fosmids corresponding to each phage and generated random 6–12 kb fragments. These fragments were then sub-cloned into a low-copy plasmid vector to create three high-coverage sub-libraries. Cells harbouring these high-coverage sub-libraries were then selected for resistance to their respective phage, and plasmids from positive clones were sequenced by Nanopore long-read sequencing. The sub-library clones that survived selection each contained the originally selected defence system, flanked by random lengths of adjacent DNA from the original fosmid insert. Thus, when reads were mapped back to the genome fragment in a given fosmid, the coverage maxima typically delineated the boundaries of each candidate defence system (Fig. 2c). In some cases, this identified previously characterized systems, including type III and IV RM systems and an Old-family endonuclease, or non-defence genes that account for phage resistance, such as malI, a regulator of the lambda receptor LamB (Supplementary Table 1). After discarding such cases, there were 21 unique candidate defence systems, with 10, 6 and 5 systems from the selections for T4, λvir and T7, respectively (Fig. 2d and Supplementary Table 2). Each candidate system was provisionally named with a PD prefix, for phage defence, followed by the phage used for selection and a unique number.

Validation of candidate defence systems

To validate these novel defence systems, we cloned each candidate open reading frame (ORF) or operon into a low-copy vector under the control of its native promoter in wild-type MG1655 (Supplementary Table 2). We confirmed that each system did not affect phage adsorption (Extended Data Fig. 2b) and then challenged each with a panel of 10 diverse phages (Fig. 3a). Each candidate system was confirmed to substantially reduce the EOP for the phage originally used to select the system, and often others. In some cases, a given defence system did not change the EOP of a phage but instead produced smaller plaques. Although most defence systems were relatively specific, protecting against only a few phages, some systems provided relatively broad protection, such as PD-λ-5, which affected EOP or plaque size for all but one of the phages tested. The 10 systems selected to defend against T4 also generally protected against the other, related T-even phages, T2 and T6. Most systems protected most strongly against the phage originally used to select it, but with some exceptions. For instance, PD-λ-5, PD-T7-1 and PD-T7-3 protected more strongly against the T-even phages than against the λvir and T7 phages used to identify them. The fact that these systems provided robust defence against T4 but were not identified in the T4 selection indicates that our screen was not saturating and that the systems identified represent only a subset of the defence systems in the original 71 E. coli isolates used.

Fig. 3: Summary and annotation of 21 previously uncharacterized defence system loci.
figure 3

a, Left: each defence system was cloned into a low-copy plasmid with its native promoter and the EOP tested for a panel of ten phages. Darker colours indicate a higher level of protection. Systems leading to smaller plaque sizes are noted with an ‘S’, and systems that protect via an Abi mechanism are indicated with an asterisk. Right: for each defence system identified, the operon structure and predicted domain composition of each component are shown. Shaded regions correspond to domain predictions using HHpred, summarized by association to the PFAM clan, with short descriptions at the top. TM, transmembrane domain; HTH, helix-turn-helix. b, Bacterial growth in the presence of phage at MOIs of 0, 0.05 or 5. Robust growth at MOI 0.05 but not at MOI 5 indicates an Abi mechanism. Lines represent the mean of three technical replicates, with shaded regions indicating s.d. (see Extended Data Fig. 3 for extended MOI data). c, Plaquing of T4 on E. coli isolates ECOR22 and ECOR65 or the isogenic defence system deletions. Dilutions were done on two different plates and images combined for presentation. d, Plaquing of T4 on strains harbouring the indicated defence system or isogenic site-directed mutants of predicted domains. Asterisks indicate approximate location of mutations made. e, Instances of homologues of defence systems by bacterial class, sorted by number of instances, descending from left to right. Known systems are listed in bold for comparison to newly identified systems.

Full size image

We then sought to classify whether each system functioned via Abi or provided direct immunity such that infected cells could survive an infection. Because Abi systems require killing of infected cells, these systems typically only provide defence at a low multiplicity of infection (MOI) in which bacteria outnumber phages, whereas direct immunity provides defence to the infected cell and thus allows comparable, although not identical, rates of growth at MOIs above and below 1.0. We thus tested the growth of strains harbouring each defence system infected at an MOI of 0.05 and 5. Of the 21 systems tested, 9 provided direct inhibition of phage infection, producing comparable protection at both MOIs; the other 12 probably use an Abi mechanism of protection with stronger protection at the lower MOI (Fig. 3b and Extended Data 3a). These results validate the ability of our screening strategy to detect Abi defences and underscore the notion that Abi systems contribute considerably to E. coli immunity.

One advantage to our screen is that we have the strains from which these defence systems originated, in contrast to computational screens that have identified many defence systems from species that have sequenced genomes but are not readily available. We were thus able to delete candidate systems in their originating strains and test whether, in their native context, they protect against phage. Specifically, we tested the native role of PD-T4-2 and PD-T4-9, which originate from ECOR65 and ECOR22, respectively, strains to which T4 can adsorb but not infect. Deleting each system dramatically increased the plaquing efficiency of T4, demonstrating that these systems provide defence in both their original, native context and when introduced into E. coli K12 (Fig. 3c). We also asked whether our defence systems work in E. coli strains other than MG1655 by testing four candidate defence systems in ECOR13 and E. coli C, which are natively susceptible to the three phages. All four systems provided protection, indicating that the function of the systems identified is not strictly dependent on strain background (Extended Data Fig. 3b).

In total, 26 of 32 proteins in the 21 systems identified here were annotated in GenBank as either ‘hypothetical protein’ or as containing domains of unknown function and had no primary sequence homology to any characterized anti-phage defence system. To more sensitively characterize each protein, we used HHpred to detect even remote similarity to PFAM domains19. This did not reveal any homology to a known system, and in the majority of cases, most of each protein remained uncharacterized. We were, however, able to detect potential similarity to some motifs or domains characteristic of defence systems, for example, nucleic-acid binding or cleavage domains (Fig. 3a and Supplementary Table 3). In most cases, this was limited to small regions hinting at enzymatic function, but not a mechanism of activation or specific targets.

Remote homology detection revealed several intriguing features uncharacteristic of known defence systems. These included (1) similarity to a ribosome-dependent ribonuclease (RelE) in conjunction with a phage-sheath-like domain, (2) a three-gene operon encoding an exotoxin A-like domain known to participate in bacterial virulence although not phage defence, an unknown protein and a SecB-like chaperone, (3) a putative membrane-anchored protein with a central coiled-coil domain (DUF4041) and a C-terminal DNA binding/cleavage domain, (4) a beta propeller fused to a DNA binding/cleavage domain, (5) a P4 phage β-like protein, (6) a lipoprotein, (7) a zinc-finger-like domain fused with a C-terminal domain of unknown function that belongs to an extended family of ribonucleases, (8) a CoiA domain and (9) DUF6575.

Several defence systems showed similarity to domains that are less frequently associated with defence systems compared with domains such as nucleases and helicases. These less frequent domains included a peptidase, a eukaryotic-like Ser/Thr kinase20, a NAD+-binding Sir2 homologue21 and a GIY-YIG nuclease22. Four of the systems identified contained a component with similarity to a toxin of toxin-antitoxin (TA) systems, but either none of these systems were found in existing TA databases or they encoded additional uncharacterized components aside from simply toxin and antitoxin. Collectively, our results reveal a striking diversity of proteins involved in bacterial defence and highlight the vast, unexplored landscape of antiviral immunity in bacteria.

In total, 11 defence systems had a component suggestive of DNA binding or cleavage activity (often including distant similarity to the restriction endonuclease-like fold PD(D/E)XK), although this similarity was, as noted above, typically restricted to a small region or motif. These putative DNases are unlikely to be part of restriction-modification systems as only one included a predicted methylase. Notably, 7 of 11 of these provided direct immunity (not Abi), suggesting non-self, nucleic-acid-targeting activity in potentially novel ways. Remote similarity to HEPN motifs or domains were found in 6 of the 21 systems. These domains are also present in the ribonucleases associated with CRISPR-Cas and toxin-antitoxin systems, supporting the notion that they are common, versatile components of defence systems23. In four cases, we mutagenized key conserved residues in the predicted domains and found that they were essential to phage protection, suggesting that the remote domain predictions feature in the function of these systems (Fig. 3d).

To assess the conservation of the 21 systems identified, we investigated the phylogenetic distribution of each (Fig. 3e). Homologues of each system (that is, those that encode all components) were found in other gamma-proteobacteria, with 16 and 18 having homologues in alpha- and beta-proteobacteria, respectively. More than half were also abundant in Firmicutes, Actinobacteria, Bacteroidetes and Spirochaetia, suggesting that many of the systems represent new, widespread classes of anti-phage defence systems. Seven of the systems we identified were restricted to proteobacteria, and three were exclusively found in γ-proteobacteria. Thus, the phage immune landscape of E. coli is composed of both widespread and clade-specific systems.

Mobile elements dominate the defence system landscape in E. coli

Previous searches for proteins enriched in defence islands identified and validated 38 novel defence systems12,13. None of the 21 systems identified here are homologues of those systems, and only one component (of PD-T7-2) resembles (32% identity) a protein of a previously validated multi-component system. PD-T4-8 has a DUF4263 domain in common with the Shedu defence system but is not homologous to the validated Bacillus cereus system. The computational approach in ref. 13 also identified 7,472 protein families enriched in defence islands that have yet to be validated. Only 14 of 32 proteins identified here have homology to those, and often with <35% identity over limited regions of the proteins (Supplementary Table 4). These observations suggest that our experimental-based selection may uncover different types of defence systems than can be found computationally by searching for enrichment in defence islands.

To further probe this idea, we analysed the native genomic context of the 21 systems identified here. We found that 12 of the 21 systems were located in intact prophages (Extended Data Fig. 4). Seven of these systems were found within P2-like prophages, with four located in the same position of the P2 genome, directly between the genes encoding the P2 replication endonuclease and portal protein (Fig. 4a). This location has been previously found to harbour the defence systems tin and old, the former being also identified in our screen against λvir24 (Supplementary Table 1). More recently, this location was found to encode a wide array of previously uncharacterized defence systems25. We also observed a second defence-enriched locus in P2-like phages, which contained three of the systems discovered here, and the previously identified defence gene fun in the P2 reference genome24 (Fig. 4a). This hotspot was not noted in the recent analysis of P2 defences. Notably, 47 of our 71 E. coli strains together encode 63 P2 portal proteins, with 110 unique proteins found in the adjacent hotspot. Only 13 of these 63 hotspots contain previously known defence systems, and most genes are annotated as hypothetical. These findings suggest that not only do P2 prophages encode a rich diversity of anti-phage proteins25 but that P2 defence hotspots make up a substantial fraction of the immune landscape in E. coli (Fig. 4b).

Fig. 4: Prophages and MGEs are major sources of defence systems.
figure 4

a, Hotspots of previously uncharacterized defence systems. Left: the native genome context of seven defence systems identified here, showing the boundaries of P2-like prophages in the genome from which they originated. Genes are colour-coded as indicated at the bottom. Right: two defence systems were identified in the accessory region of an ICE-like element within the indicated genome. Homologous elements from other bacterial genomes contain known and putative defence systems in the same location. b, All identified instances of P2 defence hotspot #1 in our 73-strain E. coli collection. White genes are flanking conserved P2 genes. Each colour of gene within the hotspot represents a protein cluster (30% identity). All grey genes belong to a lone cluster. Double slashes denote the end of a contig. c, Number of defence and prophage-associated genes ±10 kb from system homologues. The scatterplots indicate, for each homologous system, the numbers of prophage and defence-associated genes within ±10 kb. Examples in c represent systems that were found outside of prophages in our genome collection. For all 21 systems, see Extended Data Fig. 5. d, Same as c but for systems we found in prophages. e, Examples of the ±10 kb context for PD-T7-1 homologues in a defence island or prophage. f, Distribution of the number of nearby defence-domain-containing genes in homologues of systems commonly found in prophages (>10% homologue-containing regions with 8+ prophage genes in proximity) or not; n = 12 and 9 systems, respectively. Boxes indicate bounds of the distribution as median ± quartiles, and limits exclude outliers. g, Linear least squares regression for total nearby prophage gene and defence-domain-containing genes for each system. Pearson r = −0.442, P = 0.045; error indicates 95% CI.

Full size image

Five defence systems were found associated with other types of prophages or their remnants, including P4 satellite prophages or related integrases, Mu-like phages and lambdoid phages (Extended Data Fig. 4). We also observed a defence-enriched locus within an integrative and conjugative-like element (ICE), containing two systems we identified—PD-T4-4 and PD-T7-2 (Fig. 4a). This element is not widely distributed in our E. coli strain collection, but similar elements in other genomes encode known defence systems at this locus, as well as hypothetical proteins that may also function in phage defence (Fig. 4a). PD-T7-4 and its homologues often overlap an integrase gene, while PD-T4-5 was identified on a plasmid. The remaining 4 systems did not appear to reside within active MGEs, but each had a nearby integrase gene suggesting that they may be part of decaying MGEs (Extended Data Fig. 4).

Our study supports previous findings that in addition to defence islands, prophages and other MGEs are a rich reservoir of defence systems26. However, these categories are not mutually exclusive, as some defence islands may be carried on or derived from MGEs27. To more systematically document the different genomic contexts for the systems identified here, we collected, for each of the 21 systems, all homologues in a set of 844,603 publicly available bacterial genomes. We classified the genes within 10 kb upstream and downstream of each homologue as either defence-related, prophage-related or neither (Methods). We then tabulated the number of defence- and prophage-related genes flanking each homologue (Fig. 4c,d and Extended Data Fig. 5).

We detected two distinct patterns. For systems that we found outside of prophages in our strains, the homologues were also typically not prophage-associated (Fig. 4c and Extended Data Fig. 5) and were often near several other genes encoding defence-associated domains. Thus, these systems do appear in defence islands, even though they were not previously detected as enriched in them. For the systems we identified within E. coli prophages, some of their homologues were also found in prophages, as evidenced by dozens of flanking prophage-related genes (Fig. 4d and Extended Data Fig. 5). These prophage-associated homologues were typically near 1–2 defence-associated genes but rarely more than 2, suggesting that some systems reside in small defence islands or clusters within a prophage, as with the P2 hotspots (Fig. 4d)25.

Notably, there are homologues of each system we identified that can be found in defence islands (Fig. 4e), some more rarely than others, indicating that they do not require a prophage context to function. In aggregate, we observed an inverse correlation between the number of neighboring defence-associated genes for homologues of systems commonly found in prophages versus those that are not (Fig. 4f,g). This highlights the constraint on how many defence systems can be carried by prophages, or within a given hotspot, due to size limitations in DNA that can be packaged. We suspect that this is one reason why many of these genes are enriched in prophages as compared with defence islands.

Finally, we compared defence island and prophage enrichment (Methods) between systems discovered here and those previously predicted computationally and validated experimentally12,13. We found that our experimentally selected systems were on average less frequently associated with known defence genes but more frequently associated with prophage genes (Extended Data Fig. 6). These analyses suggest that defence-island enrichment methods may be less sensitive in identifying defence systems frequently found in prophage.

Previously uncharacterized toxin-antitoxin-like systems

Toxin-antitoxin (TA) systems are typically composed of a protein toxin that can arrest cell growth but is normally neutralized by a cognate, co-expressed antitoxin. TA systems are extremely prevalent in bacterial genomes and MGEs, but their functions remain poorly understood28. A handful of TA systems have been found to provide anti-phage defence through an Abi mechanism29,30,31. Our selection yielded four different systems that were recognizable as TA-like in nature. These encoded gene products with sequence similarity to toxic proteins, mostly featured multiple components, and provided Abi defence (PD-T4-5, PD-T4-7, PD-T4-9 and PD-λ-2). A fifth (PD-T4-10) facilitated Abi defence and had two overlapping ORFs, reminiscent of many TA systems. As noted above, none of these systems were previously annotated as TA systems, so we sought to validate the three featuring multiple components: PD-T4-10, PD-λ-2 and PD-T4-9 (Fig. 5a).

Fig. 5: Previously uncharacterized toxin-antitoxin-derived defence systems.
figure 5

a, Schematics of PD-T4-9, PD-T4-10 and PD-λ-2 defence system operons and their domain predictions. Representative homologues of the systems are shown in their genomic contexts and indicate conservation and order of the system components. Blue, putative antitoxin; red, putative toxin; green, accessory factor. b,d,f,i, Each component or pair of components indicated was expressed (+) or not (−) from an inducible promoter and assayed for viable colony-forming units in tenfold serial dilutions. c,e, Plaquing assays for the phage indicated on cells harbouring an empty vector or a vector containing a given defence system with all components (WT) or lacking the component indicated. g,j, Plaquing of phages on TAC-containing strains expressing a second copy of the chaperone component to varying levels during infection. h, Schematic and T2 plaquing assay of the mqsRAC TAC system.

Full size image

For PD-T4-10, neither of the two proteins had a predicted function. We expressed each component from an inducible promoter and found that the second component, PD-T4-10B, was toxic. This toxicity was fully neutralized when PD-T4-10A was co-expressed (Fig. 5b). PD-T4-10A could not be deleted, consistent with it being an antitoxin, whereas deletion of the toxin PD-T4-10B abolished resistance to T4 infection (Fig. 5c). Thus, this system comprises a novel, bona fide TA system that provides strong protection against T-even phages.

PD-λ-2 features three components. The first is similar to HigB toxins that inhibit translation, the second has a HigA antitoxin-like domain (Xre-family helix-turn-helix) fused to a C-terminal peptidase domain and the third component is related to a P4 phage antitoxin. An Xre-peptidase fusion, co-expressed with an upstream toxin-like gene, has been observed in sequence-based TA searches but the function of these loci is unknown32. Unexpectedly, we found that overexpressing the second component, PD-λ-2B, alone was toxic, with toxicity rescued by co-expression with PD-λ-2A (Fig. 5d). We also confirmed that PD-λ-2C, although not required to neutralize the toxin, is required for defence against λvir (Fig. 5e).

PD-T4-9 also contains a third component, a SecB-like chaperone, suggesting that it is related to an enigmatic class of TA systems called toxin-antitoxin-chaperone (TAC) systems, so we renamed this system CmdTAC for chaperone-mediated defence. The antitoxins of TAC systems are typically homologous to canonical antitoxins but feature an unstructured extension of their C termini called the chaperone addiction domain (ChAD)33. In the absence of the cognate SecB-like chaperone, the ChAD renders the antitoxin prone to aggregation and proteolytic degradation, thereby freeing the toxin33. Inducing the expression of CmdT, a previously uncharacterized toxin that has an ADP-ribosyltransferase-like domain, was toxic. Co-expression with the presumed antitoxin, CmdA, only marginally improved viability. However, co-expression with both the putative antitoxin and chaperone components completely restored viability (Fig. 5f). Together, these results suggest that CmdTAC is a novel TAC system that protects against phage.

To further characterize the role of the chaperone, CmdC, we overproduced it during T4 infection of CmdTAC+ cells. Interestingly, an oversupply of CmdC abolished phage protection by CmdTAC, suggesting that destruction of CmdC, or sequestration of CmdC from the complex may allow this TA system to activate in response to phage infection (Fig. 5g). The chaperone normally promotes neutralization of CmdT by CmdA but, following infection, could be depleted or sequestered by a phage product, leading to liberation of the toxin and abortive infection. Providing additional CmdC prevents the loss or complete sequestration of chaperone, thereby preventing release of CmdT from CmdA. Additional work is needed to further dissect the mechanistic nature of chaperone-phage interactions.

We hypothesized that TAC systems may be a broad class of phage defence systems. To test whether other TAC systems can protect E. coli MG1655 against phage infection, we cloned and tested an MqsRAC system from E. coli C496_10. Although completely unrelated in toxin and antitoxin sequence to CmdTAC, MqsRAC is a canonical TAC system, has been characterized in Mycobacteria and includes a SecB-like chaperone homologous to CmdC33. This system conferred robust protection against T2 (Fig. 5h) but not T4 as with cmdTAC. Similar to CmdTAC, toxicity of MqsR could only be rescued by expressing MqsA and MqsC (Fig. 5i), and inducing additional MqsC in MqsRAC+ cells inhibited phage protection (Fig. 5j). Thus, our work indicates that TAC systems may be a widespread and diverse new class of phage defence system.

Discussion

Our work indicates that a large reservoir of diverse, previously unknown phage defence genes is distributed across the E. coli pangenome. Like many bacteria, there is tremendous variability in the ‘accessory’ genomes of different strains of E. coli. Many of these accessory genes are probably associated with phage defence, as recently suggested for marine Vibrio species34. Although efforts to find new defence systems on the basis of proximity to known systems have proven fruitful, our work reveals that many phage protective systems remain unidentified.

Our functional screening, which is agnostic to the genomic context of defence systems, indicates that many systems are not commonly or detectably enriched among known defence islands. Indeed, none of the 21 systems we identified were previously reported to provide phage defence in previous discovery efforts based on defence island enrichment. Moreover, only three of the systems we identified were within 10 kb of other known defence systems in our E. coli genomes. However, notably, 15 of the 21 systems we identified were present in apparently active or recently active MGEs (prophages, ICEs or plasmids) (Extended Data Fig. 4), with the other 6 located in regions suggestive of MGEs in the late stages of decay. Homologues of the systems discovered here are sometimes present in defence islands (Figs. 4d,e and 5a), but these associations are often relatively rare.

It is well established that MGEs contribute to antiviral defence in bacteria26. By providing a glimpse into the relative contributions of MGEs and defence islands to immune system context, our results support the notion that active prophages and other MGEs are probably the primary reservoirs of anti-phage defence systems in E. coli. This idea is consistent with studies that have identified diverse anti-phage systems in P2 and P4 prophages24,25 and another recent study revealing how the transfer of defence-containing ICEs drove the emergence of phage resistance in clinical isolates of V. cholerae34. The defence systems found in functional MGEs probably help these elements to protect themselves and host resources by preventing infection by other phages27. However, given genome size and packaging constraints, there is a limit to the number of defence systems carried by a given prophage. Such constraints may help explain why increased prophage association was correlated with lower defence gene association (Fig. 4g). Future investigation should work towards further uncovering the vast reservoir of anti-phage elements carried on MGEs.

Of the 32 proteins in 21 systems found here, 13 feature domains never before reported to function in phage defence. These include a protein with distant similarity to RelE and a region of a phage tail sheath, an Exotoxin A-like domain previously only shown to function in bacterial virulence, a SecB-like chaperone, DUF6575, DUF4041, a CoiA competence-related domain, a β-propeller fold, an ImmA/IrrE peptidase, a HigB toxin and three proteins with no ascertainable similarity to deposited domains. In addition, ten proteins contained large regions with no predicted domains (Fig. 3a).

Although some regions of the proteins identified have distant similarity to known nuclease motifs, these domains are found here in new or unusual contexts and associations, which raises fascinating questions for future investigation. For example, one features a 7-bladed beta propeller with a separate N-terminal Mrr-family nuclease domain. The β-propeller fold consists of separate modules that adopt a disc-like, circularly arranged structure with a central channel that can accommodate many substrates including protein and DNA35. Determining how this domain aids in activation or target specificity of this system opens many avenues of future discovery. We also uncovered five single-protein systems that exclusively contain a putative DNase domain. Although these domains are found in defence systems such as RM and CRISPR, how orphan DNase-like proteins sense and respond to phage infection, especially to provide direct defence, is unclear. As all but one is encoded without a DNA methylase, they probably do not distinguish phage and host DNA in the same way as RM systems. Additionally, these DNase-containing systems protected against T4 and T7, which are intrinsically resistant to most RM systems36,37. Type-IV RM systems are single-component defences known to target the modified cytosine of T-even phages38, but none of the five proteins discussed here resemble these systems. How those that are specific to T7 might target T7 DNA, which is not modified, is unclear. The existence of unconventional, phage-targeting, nucleic-acid degradation systems underscores a knowledge gap in the molecular mechanisms of viral resistance and self/non-self recognition.

We also identified and validated three TA systems as phage defence elements including cmdTAC, prompting the discovery that unrelated TAC systems such as mqsRAC also function in antiviral immunity. MqsRA is a well-studied TA system but has no documented role in phage defence; our results suggest that bacteria have co-opted this TA system by addition of a chaperone-dependent antitoxin to activate the toxin in response to infection. Our findings support the notion that TA systems play a central role in phage defence29.

Finally, we note that four systems identified here probably belong to classes of previously described systems but have diverged quite significantly such that they share little significant sequence homology, reinforcing the extreme divergence and adaptations that typify many immune systems. These divergent systems include PD-T7-2 in which the second protein is similar to HerA of Sir2/HerA13, PD-T4-8 whose central domain bears distant resemblance to Shedu12, PD-λ-5 which appears to be a highly compacted prophage version of an RM+Abi system and PD-T4-5—a plasmid-borne gene that is a distant relative of AbiF/D.

Our screening methodology enables the experimental discovery of anti-phage defence genes and has several powerful features. First, we can return to the strain of origin for a given system and demonstrate that it provides defence in its native context. Second, there is a built-in pairing of defence systems to the phage they defend against, whereas with computational studies, the phage(s) a given system defends against must be subsequently identified. Third, in the pipeline described here, the source DNA comes from other E. coli strains, which probably minimizes false negatives that can arise from producing candidates in a heterologous host. Finally, our experimental approach is not limited to genes that are detectably enriched in defence islands. As noted, only three of the systems we identified were natively associated with obvious defence islands, and some also do not appear to have many close homologues in defence islands. Some defence systems may not associate with defence islands, while some may have arisen too recently or not be widespread enough to detect an association. Indeed, some of the systems identified here show a relatively limited phylogenetic distribution. However, the phage defence capabilities of bacteria probably include both broadly conserved and clade-specific systems adapted to the unique biology of a given organism and its phages.

A limitation of our method is that the initial selection does not distinguish between clones that lead to cell surface modifications that block adsorption and those providing intracellular defence systems; these must be distinguished in downstream steps. Further, at present, our method probably does not reach saturation as some systems are strongly overrepresented among the positive clones, which probably hinders the identification of others. The reason for this uneven distribution is unclear, but efforts are underway to increase throughput and sensitivity of the method.

The methodology developed here can be powerfully extended in several ways. First, genomic DNA from other sources, including metagenomic DNA, could be used as input material. From just 71 strains, we identified 21 new defence systems, suggesting that fertile ground remains to be discovered both in and beyond E. coli. Second, the panel of phages tested here was limited to three and could easily be expanded, particularly given the enormous diversity of phages. Finally, with only small modifications, any transformable bacterium could be used as the host strain. Further identification and characterization of bacterial immune systems promise to shed new light on the ancient arms race between bacteria and their viral predators, and may also have practical applications, providing the foundation for precise molecular tools and helping to inform future efforts to develop phage as therapeutics.

Methods

Bacteria and phage growth and culture conditions

Cultures were routinely grown at 37 °C in LB medium unless otherwise stated. Phage stocks were propagated on MG1655, filtered through a 0.2 µm filter and stored at 4 °C. Select ECOR strains were obtained from the Thomas S. Whittam STEC Center at Michigan State University, and UMB isolates were obtained from Alan J. Wolfe at Loyola University Chicago (Supplementary Table 5). Other strains, plasmids and primers + synthetic gene fragments are listed in Supplementary Tables 6–8, respectively.

Library construction

Genomic DNA was collected from pooled overnight cultures of each E. coli isolate using the PureLink genomic DNA mini kit (Invitrogen). From this sample, a fosmid library was constructed by Rx Biosciences using the CopyControl fosmid library production kit (Lucigen) according to the manufacturer’s protocol. Plasmid sub-libraries were constructed first by extracting fosmid DNA from select positive clones using the ZR plasmid miniprep kit (Zymo Research). Equimolar pooled fosmids were sheared to an average of 8 kb using g-TUBEs (Covaris). Sheared DNA was end-repaired and 5’-phosphorylated using the End-It DNA end-repair kit (Lucigen) and then purified using the DNA Clean and Concentrator kit (Zymo Research). The plasmid vector was prepared by PCR, and blunt-ended fragments were ligated to the plasmid using T4 DNA Ligase (NEB) for 2 h at room temperature. The ligation reaction was electroporated into MegaX DH10B T1R Electrocomp cells (Thermo Fisher) and selected on LB with 50 µg ml−1 kanamycin.

Defence system selection

We used a variant of the previously described ‘tab’ selection procedure to select for fosmids that provide resistance to phage infection17. A heavy inoculum (>30 µl) of a high-titre library freezer stock or empty vector freezer stock was inoculated into 5 ml LB containing 20 µg ml−1 chloramphenicol and grown to stationary phase at 37 °C (approximately 4 h, optical density (OD)600 = 2–3). Cultures were adjusted to OD600 = 1.0, and 0.1 ml (~8 × 107 cells) was pipetted onto one side of 3–6 empty 15 cm Petri dishes. A tenfold dilution series of phage stock was prepared, and 0.1 ml of each dilution was pipetted onto the empty plates containing the bacterial cultures (onto a separate area of the plate, preventing mixing). Molten LB (20 ml) 0.5% agar was added to the plate and briefly mixed to disperse bacteria and phage. Plates were incubated at 37 °C overnight, except for one T4 screen which was conducted at room temperature. Bacterial colonies were picked from the plate containing the phage dilution that produced the largest difference in number of colonies between the control and library samples and then streaked onto fresh plates to isolate single colonies. To test phage resistance phenotypes, single colonies were cultured overnight, and 30 µl of culture were mixed with 5 ml of molten LB 0.5% agar and 30 µg ml−1 chloramphenicol in 8-well rectangular dishes. Serial dilutions of phage were spotted onto the solidified culture media and incubated overnight at 37 °C. Fosmid DNA was isolated from positive clones, and the end of each fosmid insert was Sanger sequenced.

With some phage-resistant clones, we observed no lysis even at very high phage concentrations (Fig. 2b). All sequenced clones with this phenotype contained either LPS or capsule biosynthesis genes. In our experience, no intracellular defence system completely prevents visible lysis at extremely high phage concentrations, whereas changes in the phage receptor or cell surface can, so we suspected that these clones disrupted phage adsorption. Similarly, with regard to T7, genes for capsule biosynthesis survived selection but did not display any difference from the control in plaquing efficiency. All clones showing a complete lack of lysis, or no change in EOP and encoding capsule genes, were discarded. Twelve positive clones (all from the λvir selection) produced a high number of discrete escape plaques, consistent with RM systems. Strains that displayed this escape phenotype and whose fosmids contained an identifiable RM system were discarded from further analysis.

Sub-libraries were screened as described above using a variation that allowed bulk collecting of all positive colonies directly from the screening plate (similar to the previously published gro screen)39. In this variation, instead of molten LB 0.5% agar, bacteria and phage samples were spread on the surface of LB 1.2% agar using glass beads.

Long-read sequencing and defence system identification

After cells were collected in bulk from the screening plates, total plasmid DNA was extracted and linearized by digestion with the restriction enzyme NdeI, EagI or FsoI. For Oxford Nanopore sequencing, linearized samples were characterized on a FemtoPulse (Agilent Technology) to confirm integrity, and high-quality samples were indexed by native barcoding (ONT kits EXP-NBD104/114) with supporting reagents from New England Biolabs. Libraries were prepared using the LSK-109 chemistry, and samples were run on either a MinION (R9) or PromethION (R9.4) flowcell. Basecalling was done using built-in ONT tools. Processed reads were aligned to public reference genomes of the source organisms, or the relevant portions of the genomes contained in the fosmid inserts, using the Minimap2 plugin within the Geneious Prime 2020.2.4 software suite. Fosmids that could not be mapped to their genomes due to contig gaps were also sequenced identically, and de novo assembly was conducted using Flye40, also in Geneious. Candidate defence systems were predicted to be any gene or operon residing under the coverage maxima. In the few cases in which the result was ambiguous, candidates were cloned after prioritization by features including domain prediction, location in defence hotspots, hypothetical proteins and by general comparative genomic investigation. Multi-component systems (operons) were predicted by ORF proximity, promoter prediction and gene co-occurrence in homologues.

Strain construction

Defence system cloning was performed using Gibson Assembly of PCR products containing predicted defence systems and their predicted promoters into a destination vector lacking an upstream promoter. Assembled plasmids were transformed into MG1655 using the transformation and storage solution (TSS) method41. MG1655 with a deletion of the region containing mrr was used as the host strain for PD-λ-5. Site-directed mutagenesis was conducted by PCR using outward facing primers containing the desired mutation and with compatible overlapping regions. Amplification of the wild-type template plasmid was cycled 18 times, and the reaction was chemically transformed into DH5α cells. In-frame deletions of defence systems were constructed by transforming a temperature-sensitive plasmid expressing λ-red recombinase into the target strain. Oligos with overlapping regions to the genome targeted for deletion were used to amplify a kanR-resistance marker. The amplicon was then electroporated into the target strain induced to express λ-red, and recombinants were selected on kanamycin. The recombinase plasmid was then cured from the target strain by growth at 37 °C.

Efficiency of plaquing assays

50 µl of overnight culture were mixed with 3 ml LB 0.5% agar and overlayed onto plates containing LB 1.2% agar and appropriate antibiotics. Two microlitres of phage from a tenfold serial dilution of stocks were pipetted onto the surface of the overlay plate. Spots were allowed to dry and incubated at 37 °C until plaques were visible. Plaques were then enumerated and EOP was measured as total plaque-forming units on the experimental strain divided by plaque-forming units on the control WT strain. Often, individual plaques were not distinguishable, that is, no viable phage were produced in an infection, resulting in a lysis zone but no discrete plaques. In such an event, samples were counted as having one plaque on the last dilution that showed lysis. For EOP assays with TAC systems, chaperone expression was titrated by overlaying cultures on media with increasing concentrations of arabinose before spotting phage dilutions.

Bacteriophage adsorption assay

The method is adapted from ref. 42. Overnight bacterial cultures were diluted 1:100 and grown to OD600 = 0.5. Cultures were infected at an MOI of 0.1. Samples were then incubated at 37 °C (T4 and λvir) or 25 °C (T7) for 15, 25 or 15 min for T4, λvir and T7, respectively. Samples (500 µl) were then added to a tube of ice-cold chloroform and vortexed, and unadsorbed phage were enumerated by the top agar overlay method using a susceptible indicator strain. Percent adsorption was determined relative to a simultaneous mock control experiment that contained growth medium but no host cells.

Abortive infection assays

Overnight cultures were normalized to OD600 = 1.0 and then diluted 100-fold. Diluted cultures (150 µl) were dispensed in a flat-bottomed 96-well plate. Phage dilutions (10 µl) were added to each well such that the MOI varied from 50 to 0.005. Wells were then overlayed with 20 µl mineral oil, and plates were covered with a breathable membrane. Plates were incubated at 37 °C in a Biotek Synergy H1 microplate reader. OD600 was measured every 15 min. Three technical replicates were conducted for each strain.

Toxicity assays

Strains containing plasmids with inducible promoters were grown overnight at 37 °C in LB under repressing conditions (LB or LB + 0.2% glucose). Cultures were washed in LB, and tenfold serial dilutions were spotted on LB agar with and without inducer (0.2% arabinose, 200 µg ml−1 vanillate, or 100 µg ml−1 anhydrotetracycline). Plates were incubated overnight at 37 °C.

Bioinformatic analyses

Sanger sequences of fosmid ends were mapped to their strains of origin using BLASTn43, followed by manual inspection. General remote domain prediction was done using the HHpred online web server (https://toolkit.tuebingen.mpg.de/) or locally (HHblits and HHsearch) against Pfam A domains (v. 35.0)19,44. To label domains in Fig. 3, we used the top HHsearch hit for each independent region of the protein. If there were many good matches, the bounds of the predicted domain were taken from the top hit, while the label was chosen on the basis of the Pfam clan to which the top hits belonged. The only exception was PD-λ-5, for which the top hit, ‘methyltransferase’, was chosen as we deemed it more descriptive than the Pfam clan designation. Investigation of the P2 defence hotspot was conducted by identifying homologues of the P2 portal protein using BLASTp, extracting the surrounding genes, clustering to 30% protein identity and visualizing using Clinker and clustermap.js45. Pangenome analyses were performed by annotating Genbank assemblies with Prokka46, followed by analysis with Roary47. The phylogenetic tree was generated using FastTree48 on the core genome alignment produced by Roary, using a generalized time-reversible model. We used DefenseFinder49 on our protein sequences to determine whether homology to known systems could be detected. Other software that were instrumental for routine genome analyses were PATRIC web server, DNAFeaturesViewer and Mauve50,51,52.

To assess whether defence systems were potential toxin-antitoxin systems, we used WU-BLAST 2.0 to search against TADB v2.0 (https://bioinfo-mml.sjtu.edu.cn/cgi-bin/TADB2/nph-blast-TADB.pl)53.

Identification of defence system homologues and genomic context analysis

For each defence system, we searched for homologues of each individual component using blastp against all bacterial proteins in the NCBI non-redundant (nr) protein database using the following parameters: -evalue 0.00001 -qcov_hsp_perc 80 for single-gene systems and -evalue 0.00001 for multi-gene systems. The NCBI nr database was downloaded for local use in March 2021. All instances of the homologues identified in the nr search were then located within all full bacterial genomes (n = 844,603) downloaded from Genbank and RefSeq (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/) in April 2021. For multi-gene systems, the system was only considered present in a given genome if all components of the system were present in the same genomic region. The only exception was CmdTAC in which otherwise clearly homologous systems were widely variable in the antitoxin sequence. For this system, homologues were required to have a CmdT homologue, a second downstream ORF and a SecB-like chaperone as the third component.

The local genomic context of a defence system homologue was defined as all coding sequences ±10 kb of the system (or to the end of a contig if less than 10 kb). All coding sequences within this local context were searched for defence-related and prophage-related domains using HMMER3 hmmscan54 with E-value cut-offs of 10−5 and 10−15 for the defence-related and prophage-related domain searches, respectively. For the defence-related domain search, sequences were searched against defence-related pfam and COG domains identified and used in refs. 11,12,13. We considered each gene flanking a given homologue separately, even if multiple adjacent genes were part of a single multi-component defence system. For the prophage-related domain search, sequences were searched against all pVOG55 domains available as of May 2021 when the pVOG database was downloaded. For scatterplots and marginal histograms in Fig. 4 and Extended Data Fig. 5, any regions with <10 coding sequences (that is, located on short contigs) were excluded. In native context schematics, prophage genes were predicted by annotation, pVOG analysis and BLASTp against the ACLAME database56.

To compare defence and prophage context between our systems and those that were identified computationally and subsequently validated in refs. 12,13 (Extended Data Fig. 6), we identified homologues of each system and their flanking genes as described above. The collected flanking proteins were then clustered using the function cluster within MMseqs257 with the following parameters: –cluster-mode 1–min-seq-id 0.9. Each resulting cluster was called as defence- or prophage-related if at least 90% of the proteins within the cluster contained the same defence- or prophage-related domain(s), respectively. This clustering helps to control for overrepresentation of closely related sequences. Defence and prophage enrichment for a given system was then calculated as the number of defence- or prophage-domain-containing clusters divided by the total number of clusters.

Taxonomy analysis

The taxonomic distribution of each system was defined by the system’s presence across the downloaded bacterial genomes with the same parameters as described above. For a given genome with a defence system present, the NCBI taxid was extracted and translated to major bacterial classes using taxon kit58. For comparison, we also examined the taxonomic distribution of the following known systems: type I-IV RM systems, EcoKI, EcoRI, EcoPI and McrBC, respectively; P2 old, AAD03309.1; Cas9, WP_032462936.1; Zorya I, system containing BV17222.1; Zorya II, system containing ACA79490.1; ToxN, WP_000675353.1; Kiwa, system containing AEZ43441.1; Druantia, system containing ERA40829.1.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.