Phosphorothioate DNA modification by BREX Type 4 systems in the human gut microbiome

Phosphorothioate DNA modification by BREX Type 4 systems in the human gut microbiome | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Phosphorothioate DNA modification by BREX Type 4 systems in the human gut microbiome Yifeng Yuan , Michael S. DeMott , Shane R. Byrne , Katia Flores , Mathilde Poyet , Mathieu Groussin , Global Microbiome Conservancy , Brittany Berdy , Laurie Comstock , Eric J. Alm , View ORCID Profile Peter C. Dedon doi: https://doi.org/10.1101/2024.06.03.597175 Yifeng Yuan 1 Department of Biological Engineering, Massachusetts Institute of Technology , Cambridge, Massachusetts, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Michael S. DeMott 1 Department of Biological Engineering, Massachusetts Institute of Technology , Cambridge, Massachusetts, USA 2 Center for Environmental Health Science, Massachusetts Institute of Technology , Cambridge, Massachusetts, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Shane R. Byrne 1 Department of Biological Engineering, Massachusetts Institute of Technology , Cambridge, Massachusetts, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Katia Flores 3 Department of Microbiology, Duchossois Family Institute, University of Chicago , Chicago, IL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Mathilde Poyet 1 Department of Biological Engineering, Massachusetts Institute of Technology , Cambridge, Massachusetts, USA 4 Institute of Experimental Medicine, Kiel University , Germany 6 Global Microbiome Conservancy (), Kiel University, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site Mathieu Groussin 1 Department of Biological Engineering, Massachusetts Institute of Technology , Cambridge, Massachusetts, USA 5 Institute of Clinical and Molecular Biology, Kiel University , Germany 6 Global Microbiome Conservancy (), Kiel University, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site 6 Global Microbiome Conservancy (), Kiel University, Germany 8 Center for Microbiome Informatics and Therapeutics, Massachusetts Institute of Technology , Cambridge, MA Brittany Berdy 7 Infectious Disease and Microbiome Program, Broad Institute of MIT and Harvard , Cambridge, Massachuetts, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Laurie Comstock 3 Department of Microbiology, Duchossois Family Institute, University of Chicago , Chicago, IL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Eric J. Alm 1 Department of Biological Engineering, Massachusetts Institute of Technology , Cambridge, Massachusetts, USA 7 Infectious Disease and Microbiome Program, Broad Institute of MIT and Harvard , Cambridge, Massachuetts, USA 8 Center for Microbiome Informatics and Therapeutics, Massachusetts Institute of Technology , Cambridge, MA 9 Singapore-MIT Alliance for Research and Technology , Singapore Find this author on Google Scholar Find this author on PubMed Search for this author on this site Peter C. Dedon 1 Department of Biological Engineering, Massachusetts Institute of Technology , Cambridge, Massachusetts, USA 2 Center for Environmental Health Science, Massachusetts Institute of Technology , Cambridge, Massachusetts, USA 9 Singapore-MIT Alliance for Research and Technology , Singapore Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Peter C. Dedon For correspondence: pcdedon{at}mit.edu Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract Among dozens of microbial DNA modifications regulating gene expression and host defense, phosphorothioation (PT) is the only known backbone modification, with sulfur inserted at a non-bridging oxygen by dnd and ssp gene families. Here we explored the distribution of PT genes in 13,663 human gut microbiome genomes, finding that 6.3% possessed dnd or ssp genes predominantly in Bacillota, Bacteroidota, and Pseudomonadota. This analysis uncovered several putative new PT synthesis systems, including Type 4 Bacteriophage Exclusion (BREX) brx genes, which were genetically validated in Bacteroides salyersiae. Mass spectrometric analysis of DNA from 226 gut microbiome isolates possessing dnd , ssp , and brx genes revealed 8 PT dinucleotide settings confirmed in 6 consensus sequences by PT-specific DNA sequencing. Genomic analysis showed PT enrichment in rRNA genes and depletion at gene boundaries. These results illustrate the power of the microbiome for discovering prokaryotic epigenetics and the widespread distribution of oxidation-sensitive PTs in gut microbes. One-sentence Summary Application of informatic, mass spectrometric, and sequencing-based mapping tools to human gut bacteria revealed new phosphorothioate epigenetic systems widespread in the gut microbiome. Introduction Epigenetic modifications have been found in DNA from all domains of life. Among the variations in canonical A, T, C and G structure, such as N 6 -methyl-adenine (6mA), N 4 -methyl-cytosine (4mC), C 5 -methyl-cytosine (5mC), and 7-deazaguanine derivatives 1 , 2 , phosphorothioates (PT) are the only known modification of the sugar-phosphate backbone, with a non-bridging oxygen replaced by sulfur in R P specific configuration 3 - 5 ( Fig. 1A ). Two PT-based restriction and modification (R-M) systems, the dnd 4 , 6 and ssp 7 , 8 gene clusters, have been observed in ∼10% of bacteria and archaea 9 , 10 . Typically, the modification components dndBCDE and sspBCD are present in the form of a three- and four-gene operon, respectively. Like methylation-based R-M systems 8 , 11 , 12 , DndACDE and SspBCD proteins catalyze PT modification on one or both strands of specific consensus sequences. For example, DndACDE confer double-stranded PTs at 5′-G PS AAC-3′/5′-G PS TTC-3′ sequences in Escherichia coli B7A and Salmonella enterica serovar Cerro 87, 5ʹ-G PS GCC-3′/5′-G PS GCC-3′ in Pseudomonas fluorescens pf0-1 and Streptomyces lividans 1326, and 5′-G PS ATC-3′/5′-G PS ATC-3′ in Hahella chejuensis KCTC 2396 13 - 15 . SspABCD proteins, on the other hand, catalyze single-stranded 5′-C PS CA-3′ in Vibrio cyclitrophicus FF75 7 , 13 . The dnd and ssp systems share some similarities, such as encoding a homolog of phosphoadenosine phosphosulphate (PAPS) reductase (DndC, SspD) 7 , 16 , a homolog of cysteine desulfurase (DndA, SspA) 7 , 16 , and a P-loop containing ATPase (DndD, SspC) 7 , 17 ( Fig. 1B ). The restriction counterparts DndFGH and SspFGH sense PT modifications by poorly understood mechanisms, with only 10-15% of consensus sequences modified with PTs 7 , 13 . Download figure Open in new tab Figure 1. Microbial phosphorothioate (PT) DNA modifications. ( A ) Sulfur replaces a non-bridging phosphate oxygen in the DNA backbone in PT modifications. A limit nuclease digest of PT-containing DNA leaves PT-linked dinucleotides that can be identified and quantified by LC-MS. ( B ) Synthesis of PTs generally follows the biochemical steps performed by Dnd proteins. Both the metabolic pathways enabling sulfur incorporation into PTs and the unique chemical properties of sulfur-modified DNA have implications for interactions of PT-containing microbes with human hosts. For example, the sulfur in PTs is both readily oxidized and nucleophilic, which is proposed to provide epigenetic regulation of transcription of redox homeostasis genes 15 . PTs also provide weak to modest protective effects in cells exposed to reactive oxygen and nitrogen species, such as peroxides 15 , 18 and peroxynitrite 19 . Contrasting with this protection, PT-containing bacteria are 5-fold more sensitive to neutrophil-derived hypochlorous acid (HOCl) due to extensive DNA breaks at PTs 20 . These unusual chemical properties of PTs raise questions about how PT-containing microbes might behave in the healthy gut microbiome or be altered by inflammatory bowel disease (IBD) or other chronic inflammatory conditions 21 - 23 . Evidence for the presence of PT genes in bacterial strains associated with the human gut microbiome and PT dinucleotides in fecal DNA 24 , 25 thus motivated us to systematically analyze PT genomics in the human gut microbiome. Here we defined the landscape of PT-containing microbes in the human gut by performing a genomic analysis of 13,000 human microbiome genomes from the Broad Institute-OpenBiome Microbiome Library (BIO-ML) 26 , the Global Microbiome Conservancy (GMbC) 27 , and the Unified Human Gastrointestinal Genome (UHGG) collection 28 . Mass spectrometric analysis of PT dinucleotides in 226 of these isolates coupled with PT-specific next-generation sequencing (PT-seq) 29 led to the discovery of a new PT system involving Type 4 Bacteriophage Exclusion (BREX) genes brx PRXL . These results expand our knowledge about the diversity of PT epigenetics and lay the foundations for understanding the role of PT-containing microbes in human health and disease. Results Discovery of new PT modification systems by analyzing genome neighborhoods of dndC and sspD genes Given the power of physical clustering analyses to identify gene functions in bacteria, we first performed a comprehensive gene neighborhood analysis 30 of dnd and ssp genes to find new PT modification systems. Here we used Enzyme Function Initiative (EFI) tools 30 to first search Uniprot for DndC and SspD homologs in sequence similarity networks (SSN) followed by the EFI Genome Neighborhood Tool to identify the genomic contexts of the SSNs. DndC and SspD possess a PAPS reductase domain with ATP pyrophosphatase activity essential for PT biosynthesis ( Fig. 1B ), an activity shared by sulfur-inserting RNA modification enzymes ThiI 31 , MnmA 32 , and Ncs6 32 , and TtcA 33 . As a second criterion for PT synthesis, we required that DndC and SspD neighborhoods also possess a P-loop-containing NTPase gene. Both DndD and SspC possess P-loop-containing ATPase activity essential for PT synthesis ( Fig. 1B ). Using this approach, we retrieved 3,120 DndC homologs and 2,132 SspD homologs from Uniprot by BLAST (e-value cutoff 10 -5 ) and searched genome neighborhoods encoding both a PAPS reductase domain and a P-loop NTPase ( Supplementary Tables S1, S2 ). The resulting gene neighborhoods revealed several novel putative PT-modifying gene candidates. Each black circle or node in Figure 2 depicts a neighborhood containing dndC ( Fig. 1A ) or sspC ( Fig. 1B ), with clustering of nodes based on SSNs. Here it is clear that most of the dndC- containing gene neighborhoods fall within the largest cluster and possess additional dnd genes consistent with canonical Dnd-based PT synthesis ( Fig. 2A , nodes with red outlines). However, several of the smaller clusters contained neighborhoods with a dndC gene, a P-loop NTPase, and other known defense genes, suggesting PT-based RM systems. In one instance, dndC and the NTPase gene lie near genes for a Bacteriophage Exclusion (BREX) type 4 system 34 ( Fig. 2A , yellow nodes). Five of six BREX systems possess a BrxC/PglY ATPase, a PglX DNA methyltransferase, and PglZ phosphatase. In BREX Type 4 clusters containing dndC , the PglX adenine methyltransferase is replaced with the DndC S-inserting PAPS reductase domain protein 34 . For example, in Methanothermobacter sp. ( Fig. 2A yellow node, black outline), dndC is adjacent to an NTPase gene, a pglZ phosphatase gene, and a Lon-like protease domain-containing brxL gene typically found in BREX Type 1 and 4 systems. In Fischerella thermalis ( Fig. 2A , yellow node, red outline), the full set of dndBCDE genes cluster with brxL , brxZ (pglZ) phosphatase, sspB nickase, and a methyltransferase. The synthesis of PTs in the DndC-containing BREX type 4 system was subsequently validated genetically and by LC-MS, as discussed shortly. Download figure Open in new tab Figure 2. Gene neighborhoods analyses based on sequence similarity networks (SSNs) of DndC and SspD proteins essential for PT synthesis. SSN analysis was performed for the 3120 closest homologues of DndC (A) and 2132 homologues of SspD (B) in the UniProt database. Each node (black circle) in the network represents one DndC or SspD proteins. An edge (lines connecting nodes) is drawn between two nodes with a BLAST E -value cutoff of ≥10 –100 (alignment score of 100) in the DndC SSN or 10 – 60 (alignment score of 60) in the SspD SSN. The node outlines are colored according to the presence of other dnd or ssp genes, with red outlines denoting colocation of dndC with dndD and similarly for sspD with sspBC . Some nodes are colored according to the type of genome neighborhood structure, which are represented below for species discussed in the text. For better visualization, single nodes and clusters with a few nodes were hidden. Abbreviations: NTPase, Nucleoside triphosphatase; MTase, methyltransferase; Res, restriction enzyme; HEPN, higher eukaryotes and prokaryotes nucleotide-binding. In a second putative PT defense system, a small cluster including Bacillus cereus ( Fig. 2A , green node, black outline) pairs dndC with a putative minimal nucleotidyltransferase (MNT) and a higher eukaryotes and prokaryotes nucleotide-binding (HEPN) protein. This gene neighborhood resembles a Class II MNT-HEPN toxin-antitoxin (TA) system 35 , in which the HEPN protein is a RNase toxin, but its activity is neutralized by adenylylation by the MNT antitoxin 36 , 37 . The HEPNs in the dndC clusters lack the RNase domain, but the neighboring MNT, ThiF, has the conserved motif GSX 10 DXD of an adenylyl transferase. Coupled with the adjacent NTPase, the DndC encoded in this Class II MNT-HEPN neighborhood could confer a three-component PT-based stress response system. Finally, the dndC genes clustered with P-loop protein genes in the genome neighborhoods in two small clusters such as Sulfurimonas sp. ( Fig. 2A , blue node, black outline) and Anaerolineales bacterium ( Fig. 2A , pink node, black outline). A DUF262 gene was found in the fomer genome neighnorhood that is found in the SspE restriction component of Ssp-mediated PT systems and is involved in the PT-sensing anti-phage activity 7 . The sspD gene neighborhoods were less complicated than those involving dndC ( Fig. 2B ). Again, sspD genes were associated with BREX defense system genes, as in the small cluster containing Chloroflexota bacterium ( Fig. 2B , yellow node). This reinforces the idea that BREX type 4 genes represent a PT-based defense system. To compare the distribution of brx genes to dnd and ssp families in prokaryotes, we performed a BLASTp search with protein sequences for DndABCDE, SspABCDE and BrxPCZL as queries in 6,616 representative genomes from the Bacterial and Viral Bioinformatics Resource Center (BV-BRC; as of January 2021) 38 . Here we defined the minimal genes necessary for a functional PT synthesis system as dndCD , sspBCD , brxPCZL based on the following observations: (1) DndA/SspA are often replaced by cysteine desulfurase such as IscS; (2) DndB is a non-essential regulator; (3) DndE protein is too short to search rigorously; (4) SspE is not required for PT modification; and (5) brxPCZL are the core 4 genes that define a BREX Type 4 system, with the brxR regulator not essential. Based on gene locus information for each protein hit, we found the dndCD , sspBCD , and brxPCZL gene clusters present in 4.3%, 3.0%, and 0.6%, respectively, of the BV-BRC genomes ( Supplementary Table S3 ). Notably, both dndCD and sspBCD gene clusters co-occurred in Cyanobacteria, such as Gloeocapsa sp., Coleofasciculus chthonoplastes , and Scytonema hofmanni , among others ( Supplementary Table S3 ). This complements the gene neighborhood analysis, which revealed that dndC clusters containing brx genes, MNT-HEPN genes, and DUF3696 genes are located near dndBCD operons in Cyanobacteria such as Fischerella, Nodosilinea, and Pseudanabaena ( Fig. 2A , red-outlined yellow, green, and blue nodes in the main cluster). Similarly, 7 of 40 genomes with sspBCD operons located near brxCZL genes occurred in Cyanobacteria ( Fig. 2B , red-outlined yellow nodes in the main cluster). These observations complement the observations of Lin et al . 9 and suggest the co-evolution of the variety of PT-based epigenetic systems in the ancient photosynthetic Cyanobacteriota phylum. 39 , 40 BREX Type 4 systems are homologous to Ssp proteins and catalyze PT synthesis The genome neighborhood analyses revealed strong associations between the PAPS domain and NTPase genes essential for PT synthesis and brx genes from the BREX Type 4 family. This association was strengthened by the homology analysis (HHpred 41 ) of BREX family and SSP proteins shown in Figure 3A , which reiterates the hallmark brxPCZL genes in BREX Type 4 and the lack of the PT-defining PAPS domain and NTPase genes in BREX Type 1. For example, BrxP possesses both the PAPS reductase domain found in SspD and the DUF4007 domain found in SspB ( Fig. 3A ), which effectively renders BrxP as a SspD-SspB fusion protein. Similarly, BrxC is essentially an SspC homolog, both possessing the same DUF6079 domain ( Fig. 3A ). Download figure Open in new tab Figure 3. Homology analysis of BREX type 4 systems and evidence of PT synthesis. ( A ) The genomic organization of the BREX type 4 gene systems in Bacteroides salyersiae and putative Butyricimonas faecalis . Vibrio cyclitrophicus FF75 is included to demonstrate a typical ssp gene system and Escherichia coli HS to demonstrate a typical BREX type 1 system. Protein domains are color-coded and labeled. ( B ) The levels of PT dinucleotides in engineered B. salyersiae and Bacteroides thetaiotaomicron strains. ΔbrxC, ΔbrxC p brxC Bf , ΔbrxC p sspCBc, vector, and brxCBs are all significantly different from wilde-type (WT) by Student’s t-test, p <0.05. Based on these similarities, we tested the PT synthesis activity of proteins encoded by genes in the BREX Type 4 family. The first step was to identify PT modifications in bacterial strains known to possess BREX Type 4 genes. Here we quantified PTs as PT-linked dinucleotides by liquid chromatography-coupled triple-quadrupole mass spectrometry (LC-MS/MS) in a limit digest with nuclease P1 and phosphatase dephosphorylation 42 . In Bacteroides salyersiae DSM18765, which harbors a 7-gene brx operon but no dnd or ssp genes ( Fig. 3A ), we detected the PT dinucleotides A PS C, C PS C, and T PS C at 56, 228, and 98 per 10 6 nt, respectively ( Fig. 3B , Supplementary Fig. S1A, Supplementary Table S4 ). As shown in Supplementary Table S4 , several other bacterial isolates with a 4-5 gene complement of brx genes ( brxCLPRZ, brxCLPZ ) and insufficient sets of dnd or ssp genes also yielded C PS C ( Prevotalla sp .), A PS A and A PS C (putative Butyricimonas ), T PS C ( Parabacteroides ), and A PS A ( Bacteroides ). These results clearly distinguished the BREX Type 4 gene systems from dnd and ssp families. Building on this associative evidence, we validated PT synthesis activity by creating in-frame deletion mutants of brx genes in B. salyersiae . As shown in Figure 3B , loss of brxC abolished PTs, with PTs restored by complementation with in trans expression of brxC from an expression plasmid. Based on the similarity between BrxP and SspD, we attempted to create a brxP deletion mutant with B. salyersiae but were unable to obtain this mutant, suggesting that its deletion may be deleterious. As an alternative, we reconstructed the PT pathway from B. salyersiae by cloning brxP Bs , mcrA Bs (a PT-dependent restriction enzyme in Streptomyces coelicolor ) 43 , and brxC Bs or brxC Bs alone into the expression plasmid, which we transferred into Bacteroides thetaiotaomicron VPI, which lacks PT genes but possesses a SufS cysteine desulfurase homolog of DndA. The combination of brxP Bs , mcrA Bs and brxC Bs resulted in the T PS C, C PS C, and A PS C dinucleotides in the same proportion as wild-type B. salyersiae. While brxC Bs alone was not able to confer PTs, mcrA Bs was not required as shown in the mcrA Bs mutant ( Fig. 3B ). We conclude that, in presence of a cysteine desulfurase gene ( sufS ), brxP and brxC are the minimal set of BREX genes needed for PT synthesis. The distribution of PT modifying systems in the human gut microbiome Based on the observation of PT-containing microbes in the mouse and human gut microbiome 24 , 25 , we wondered about the presence of brx -based PT systems in gut bacteria and their relationship to dnd and ssp systems. To quantify the distribution of the PT systems in human gut microbiome, we searched for PT gene clusters in 13663 human gut microbial genomes from the BIO-ML 26 and the GMbC 27 collection, and from isolate sequences of human gut microbes and metagenome-assembled genomes (MAGs) in the Unified Human Gastrointestinal Genome (UHGG) collection 28 . As with the BV-BRC searches noted earlier, we performed a BLASTp search of proteins DndABCDE, SspABCDE and BrxPCZL as queries in these 13,663 genomes ( Supplementary Table S5 ). The essential gene sets dndCD , sspBCD , brxPCZL were found to be present in 2.7%, 3.6%, and 1.4% of the humane gut microbiome genomes, respectively ( Fig. 4 ). This represents a 1.6- and 1.2-fold reduction in dnd and ssp systems and a 2.3-fold enrichment in brx systems in the gut microbiome compared to the general BV-BRC genomes. The distributions of the three gene clusters at the genome level were nearly exclusive with few exceptions ( Supplementary Table S5 ). Among the most prevalent phyla and orders in the human gut microbiome, the PT-modifying gene clusters were mainly found in Bacteroidota (Bacteroidales), Bacillota (Clostridiales), Pseudomonadota (Enterobacterales), and, to a lesser extent, Actinomycetota (Coriobacteriales) ( Fig. 4 ). These observations raised the question of the predictive power of the genomic analyses. Download figure Open in new tab Figure 4. Phylogenetic distribution of dndCD, sspBCD, and brxPC genes in human gut microbiome genomes. The reference phylogeny was reconstructed from the concatenated alignment of 10 ribosomal proteins. The colors of the triangles in the tree show the taxa at the order level. The occurrence of gene clusters was quantitatively presented by the colors and size of circles. For better visualization, orders containing less than 5 genomes were hidden. Source data are provided as Table S5. Validating the predicted PT modifications in gut microbiome isolates Given the presence of genes essential for PT synthesis in ∼1,000 out of 13,663 gut microbiome isolates, we next used LC-MS to validate the presence of PT dinucleotides in DNA in a collection of 226 bacterial isolates possessing dndCD , sspBCD , brxPCZL gene neighborhoods or individual PT synthesis genes not predicted to lead to PTs ( Supplementary Table S4 ). In total, we identified 8 PT dinucleotides, including A PS A, A PS C, C PS A, C PS C, G PS A, G PS C, G PS T, and T PS C, and did not detect an additional two dinucleotides (C PS T, G PS G) that we found in human fecal DNA samples ( Supplementary Table S4 ) 25 . Previous studies in three bacterial species with dnd and ssp genes 4 , 42 , 44 revealed G PS A and G PS T at >500 per 10 6 nt and C PS C at >2000 per 10 6 nt, with C PS A, A PS A, A PS C, and T PS C detected much lower levels of 1-6 per 10 6 nt 4 , 42 , 44 . The minor PT dinucleotides represent low-affinity binding sites for the DNA shape-selective Dnd proteins 44 . In sharp contrast, however, these minor sites in gut microbiome bacteria represent major PT modification sites, as high as 1,800 per 10 6 nt in bacteria with brxPCZL genes ( Supplementary Table S4 ). The A PS A, A PS C, C PS A, C PS C, G PS A, G PS C, G PS T, and T PS C dinucleotides were distributed unevenly among the isolates based on the type of PT modification system, with four groups emerging from the analyses ( Fig. 4 ). The first group of 39 isolates was characterized by G PS A and a predominance of dnd genes. The largest portion of Group 1 (32 isolates) possessed both G PS A and G PS T and harbored dndCD genes with or without dndBE , including isolates from Bacteroidales and Clostridiales orders. Five isolates from both Bacteroidales and Clostridiales harboring the minimal dndCD gene set showed G PS A and G PS C. One isolate that possessed a solitary brxP but no dnd genes and possessed both G PS A and G PS T. Given the presence of different dinucleotides for the ssp and brx gene families, this latter observation could be explained by genome sequencing errors or a sample mix up during regrowth, either resolved by re-sequencing. A single Blautia species harboring dnd genes showed only G PS A ( Supplementary Table S4 ). The second group of 34 isolates was characterized by the C PS C dinucleotide and ssp genes, with all but 2 harboring sspBCD ± sspE ( Fig. 5 ). One isolate harbored dndCD genes that are not proximal to the sspBCD operon and showed LC-MS evidence of low levels of G PS A and G PS T at the detection limit of ∼2 PTs per 10 6 nt ( Supplementary Fig. S2 ), though these signals await high-resolution mass spectrometric validation. One isolate in this group lacks sspBC genes and one lacks all essential sspBCD genes. Again, the consistency of dinucleotide distributions suggests that sspBCD or brxPCZL genes are present in the genomes and perhaps missed due to insufficient genome sequencing coverage. Download figure Open in new tab Figure 5. PT dinucleotides and PT consensus sequences in human gut microbiome isolates. The numbers of isolates were analyzed are listed on left. The circles represent PT dinucleotides quantified by LC-MS (left) and the presence of corresponding genes (right). The PT-modified consensus motifs were characterized in one representative isolate from each group using PT-seq. t.b.d., to be determined. *, the C PS AG motif was characterized using metagenomics PT-seq in another on-going study. **, G PS A and G PS T were detected using LC-MS QQQ but the exact mass was not verified by Orbitrap. The third group of 30 isolates carried the BREX Type 4 genes brxPCZL and the most diverse sets of PT dinucleotides. Four isolates of Parabacteroides contain C PS A, 12 isolates of Prevotella sp. contain C PS C, 3 isolates of Parabacteroides contain T PS C, and 9 isolates of putative Butyricimonas faecalis contain A PS A and A PS C (Fig. S1B). The B. salyersiae strain examined earlier, which contains C PS C, T PS C, and A PS C ( Fig. 5 , Supplementary Table S4 ). The fourth group of 122 bacteria lacked reliable LC-MS signals for PT dinucleotides. The majority (83) lacked the minimal sets of dndCD , sspBCD or brxPC gene clusters. For example, Bacteroides dorei CL03T12C01 harbors a PAPS reductase-encoding gene next to an ATPase without an adjacent dndD ( Supplementary Fig. S3 ). Several isolates harbor dndD and dptH next to a methylase and a restriction enzyme but lack dndC ( Supplementary Fig. S3 ). The absence of PT dinucleotides in these isolates agrees with the requirement for a PAPS reductase domain-containing protein and an NTPase. However, several (10) carry the minimal set of dndCD or sspBCD genes ( Fig. 4 ), so PT dinucleotides were expected. Again, we cannot rule out errors in the reference genomes due to genome sequencing and sample handling. These observations lead to a general conclusion that (1) Dnd proteins produce G PS A, G PS C, and G PS T, (2) Ssp proteins produce C PS C, and (3) Brx proteins produce everything except G PS A, G PS C, and G PS T. Understanding the basis for these differences requires knowledge about the longer consensus sequences containing the PT dinucleotides in the gut microbes. Novel PT consensus sequences in gut microbiome isolates Here we defined novel larger PT consensus sequences, especially for the BREX Type 4 system, in gut microbiome isolates using a new and highly sensitive NGS sequencing technique, PT-seq 29 . The idea here is that the PT dinucleotides such as G PS A and G PS T, are found in the larger consensus sequences of G PS ATC and G PS AAC-3ʹ/5ʹ-G PS TTC in several types of bacteria 5 . Application of PT-seq ( Supplementary Fig. S4 ) to the BREX Type 4-containing human gut microbiome isolate B. salyersiae DSM18765, earlier found to possess A PS C, C PS C, and T PS C at 62, 228, and 98 per 10 6 nt ( Supplementary Table S4 ), showed 5459 PT sites: 56 at A PS CTC, 3881 at C PS CTC, and 1513 at T PS CTC sites, respectively ( Fig. 6A , Supplementary Table S6 ). As observed with other dnd and ssp systems, the Brx proteins modified only a portion of the 27274 ACTC (0.2%), 22698 CCTC (17%), and 36290 TCTC (4%) total sites available, respectively. The observation of PTs at 3 NCTC sites, with a strong preference for CCTC, is consistent with previous observations that Dnd proteins select their target sequences based on DNA shape rather than precise binding contacts 44 . In another BREX Type 4-containing putative B. faecalis (GMbC ID 5893AJ_0218_015_F2), PT-seq revealed 9832 GA PS AG and 5536 GA PS CG ( Supplementary Table S7 ). These sites agreed with the A PS A and A PS C dinucleotides detected by LC-MS/MS ( Supplementary Table S4 ). PT-seq also detected sites with bistranded PTs: 1080 GA PS AG/GT PS TG and 667 GA PS CG/CT PS GC. However, we could not detect reliable signals for the corresponding T PS T and T PS G dinucleotides by LC-MS/MS, which may be due to their low abundance and to insensitive MS detection of T-containing nucleotides. Download figure Open in new tab Figure 6. Biogeographical maps of PTs in bacteria with BREX type 4 systems. ( A ) Pie charts depict the number of single- or bi-stranded PT-modified consensus sequences. The structural diagrams depict these motifs. ( B ) Analysis of the distribution of PT consensus sequences and modified sites in 1kb upstream and downstream regions in B. salyersiae ( left ) and putative B. faecalis ( right ). The number of total motif sites in the sense strand ( upper ) or both strands ( lower ) are represented in pink. The PT modified sites are represented in purple. The fraction of PT-modified motif sites is represented by the black line. Application of the original PT-seq protocol without biotin labeling and capturing to a dnd system-containing human gut microbiome isolate, Lachnospiraceae sp. (GMbC ID 2807EA_1118_063_H5), revealed 1,474 G PS AGC and G PS CTC sites occurring among the 16,683 possible GAGC/GCTC sites, with 1,154 modified on one strand and 160 modified on both strands ( Fig. 6A , Supplementary Table S8 ). That a large fraction of pileup sites occurred at NTTT and other sites represent artifacts of poly(dT) tailing of single strand nicks using terminal transferase. The uneven genomic distribution of PTs In addition to defining the consensus sequences for BREX Type 4 and other PT synthesis systems, PT-seq also revealed that PTs are nonrandomly distributed across bacterial genomes, which agrees with single-molecule real-time (SMRT) sequencing maps in Escherichia coli 13 . For example, the percentage of gene classes and intergenic regions containing PTs in Lachnospiraceae sp., B. salyersiae, and B. faecalis varied significantly: 6% of intergenic regions, 13% of tRNA genes, and 58% of coding sequences (CDS) ( Supplementary Table S9 ). The percentage of consensus sequences modified with PTs, which is roughly 10% in genomes studied to date 13 , 15 , also varied in individual bacterial genomes, ranging from 1% in tRNA genes in Lachnospiraceae sp. to 40% in rRNA genes in B. faecalis ( Supplementary Table S9 ). Interestingly, the consensus sequences recognized by PT synthesis proteins were underrepresented at the boundaries of CDS ( Supplementary Fig. S5 ). We further assessed the distribution of PTs within and between CDSs by quantifying PT-modified and unmodified consensus sequences in the sense strand of gene bodies and in regions 1 Kbp upstream and downstream of all CDSs. For both B. faecalis and B. salyersiae, this analysis revealed that the percentage of PT-modified consensus sequences was lower at the boundaries of coding regions ( Fig. 5B , black line; Supplementary Fig. S6 ). Published PT-seq data for E. coli BW25113, which has a bistranded G PS AAC/G PS TTC PT consensus, also showed this biased distribution of PTs ( Supplementary Fig. S7A ). Thus, both the consensus sequences for PTs and the proportion modified with PT were underrepresented on either side of CDSs. While PTs have been proposed to interfere with transcription 13 , the percentage of PT-modified consensus sequences did not correlate with the level of transcription activity ( Supplementary Fig. S7B ), nor did the number of PTs in possible promoter regions ( Supplementary Fig. S7C ). The basis for these biased distributions remains to be defined. Discussion Here we applied systems-level informatic, mass spectrometric, and sequencing tools to discover and characterize new PT modification systems in bacteria, with a focus on gut microbes. A rigorous neighborhood analysis of the key PT synthesis genes in dnd and ssp systems led to the discovery of several potential new PT modification systems, including MNT-HEPN, DUF262, and BREX Type 4 gene gene families. Genetic and bioanalytical validation of the BREX Type 4 system revealed the essentiality of brxPC genes for PT synthesis in bacteria with brxPCZL clusters. With the discovery of the dnd system in 2005 3 and the ssp system in 2020 7 , the brxPCZL cluster represents the third established PT modification system. The extent of distribution of the dnd , ssp , and brx gene clusters was assessed first among 6,616 representative prokaryotic genomes and then among 13,663 gut microbiome isolates. The difference in the frequency of dndCD , sspBCD , brxPCZL gene clusters in the BV-BRC general set of bacteria (4.3%, 3.0%, and 0.6%, respectively) and the gut microbiome isolates (2.7%, 3.6%, and 1.4%, respectively) suggests a bias for the BREX Type 4 systems in the gut environment. The total of 7.7% of gut microbiome isolates possessing PT-synthesis genes is consistent with the estimate of 5-10% of stool microbes possessing PT modifications, as we observed using mass spectrometric measurement of PT dinucleotides in human fecal DNA 25 . With the discovery of BREX Type 4 as the third PT synthesis gene family, one feature of PT modifications appears to be universal for PT epigenetics: partial modification of available consensus sequences in individual genomes. In agreement with previous sequencing studies with bacteria containing Dnd proteins 13 , 45 , PT-seq analysis revealed that the BREX Type 4 proteins also only partially modify their four-nucleotide consensus sequences, ranging from 0.2% at ACTC to 17% at CCTC in B. salyersiae ( Fig. 5A , Supplementary Table S6 ). The BREX Type system also points to another emerging feature of PT modification systems: sequence-selectivity of the dnd, ssp , and brx gene families. Previous 7 , 8 , 46 and present studies ( Fig. 4 ) consistently show that ssp systems catalyze C PS C dinucleotides mainly in CCA motifs, while the dnd systems insert PTs mainly in GN motifs: G PS A, G PS C, G PS G 42 , and G PS T dinucleotides. The latter occur mainly as bistranded modifications of GNNC motifs. We previously showed that Dnd proteins select modification sites based on DNA shape, with GAAC, GTTC, and GATC all sharing similar shapes and all being modified with PTs at G PS A and G PS T in S. enterica 44 . The fact that ssp is selective for CCA suggests a similar shape selectivity. BREX Type 4-dependent PT systems are more complicated, which may reflect the hybrid nature of the PT-synthesizing gene families in different bacteria. The five types of PT nucleotides produced by Brx proteins (C PS A, C PS C, A PS C, A PS A, T PS C) share C PS C with ssp systems but lack the G PS N unique to dnd systems ( Fig. 5A ). This may be partly explained by the sequence similarities between SspBCD and BrxPC proteins. The shape selectivity argument is again supported with brxPCZL in B. salyersiae DSM18765, with the ACTC, CCTC, and TCTC sites in an NCTC motif all modified to differing extents ( Fig. 5A ) and in B. faecalis with GAAG and GACG motifs. Our studies revealed the complexity of interactions among the three PT systems, with widespread combinations of the various components of all three systems and coexistence of two or more gene clusters in the same organism ( Fig. 1 , Supplementary Table S3 ). These evolutionary co-occurrences likely reflect a fitness advantage, such as that observed for the presence of Dnd and Ssp together, which provides complementary and synergistic protection against temperate and lytic phages as well as phage induction 46 . The co-occurrence of combinations of all three established PT systems and two putative systems in Cyanobacteria supports the hypothesis of an evolutionary divergence that may originate from a sulfur-based metabolism in ancient Cyanobacteria ancestor 9 . All three PT systems contain a gene encoding a PAPS reductase domain protein that has been proposed to involve in the initial sulfur mobilization step 5 . The observation of a biased toward brx -based PT systems and the unique biochemical properties of PTs raise questions about the role of PT epigenetics in the human gut microbiome in health and disease. Methods Bacterial strains and growth conditions Bacteria revived from BIO-ML and GMbC were cultured as previously decribed 26 . Parabacteroides spp. and Bacteroides spp., E. faecalis , and BIO-ML isolates were grown on ASF plates (Becton Dickinson) and BHIS plates 47 , BHI plates (Becton Dickinson), and Brucella agar, (Catalog # AS-141, Anaerobe Systems, Morgan Hill, CA), respectively. Culture manipulations were performed at 37 °C in an anaerobic chamber with an atmosphere of 80% nitrogen, 5% carbon dioxide, and 5% hydrogen. Growth on agar plates was harvested with phosphate buffered saline (Catalog # 10010-023, Gibco, Paisley, PA). The bacterial suspensions were pelleted by centrifugation at 10,000 rpm for 15 min at ambient temperature. The supernatant was removed and the pellets were stored at -20 °C. Bacterial stains and plasmids used and generated for BREX 4 system deletion and reconstruction are listed in Supplementary Table S10 . Bacteroidales strains were grown in basal liquid medium 48 or on BHIS plates 47 . Antibiotics used for selection include erythromycin (10 µg/mL), gentamicin (200 µg/mL), anhydrotetracycline (50 ng/mL). E. coli S17 λ pir was grown in LB broth or plates with carbenicillin (100 µg/mL) added for selection. Creation of deletion mutants and complementing clones Internal non-polar deletion mutants of brx genes HMPREF1532_02792, HMPREF1532_02793, HMPREF1532_02794, HMPREF1532_02795, and HMPREF1532_02796 were constructed by amplifying DNA upstream and downstream of each gene using the primers listed in Supplementary Table S11 . These flanking pieces were cloned into BamHI-digested pLBG13 47 using NEBuilder (New England BioLabs) and transformed into E. coli S17 λ pir. PCR confirmed plasmids were sequenced (Plasmidsaurus) to confirm correct error-free PCR amplification. The correct construct was conjugally transferred from E. coli into B. salyersiae and cointegrates were selected on gentamycin/erythromycin. Double recombination cross-outs were selected on BHIS plates with anhydrotetracycline (aTC, 50 ng/ml) and screened via PCR for mutant genotype. Genes expressed in trans in B. thetaiotaomicron and B. salyersiae Δ2793 were PCR amplified and cloned into BamHI-digested pFD340 49 using NEBuilder. Transformants were PCR screened and the plasmid was sequenced. Plasmids were conjugally transferred from E. coli to B. salyersaie and transconjugants were selected on gentamycin/erythromycin plates. Sequence Similarity Networks Sequence similarity networks (SSNs) were generated by submitting the sequences of E.coli B7A DndC (GeneBank AIF62362.1) and V. cyclitrophicus FF75 SspD (NCBI Refseq WP_016789110.1) to the EFI-EST webtool 30 using the BLAST option (e-value cutoff 10 -5 , maximum number of sequences retrieved 5000). The initial SSN was generated with an alignment score cutoff set such that each connection (edge) represents a sequence identity of approximately 40%. Sequences that share 100% sequence identity were grouped into a single node. More stringent SSNs were created by increasing the alignment score cutoff in small increments (usually by 5–10). This process was continued until each cluster is estimated to be isofunctional, in which the nodes represent enzymes that catalyze the same reaction. Isofunctionality was determined by mapping the known dnd clusters and following their movements as the alignment cutoff increased until the dnd clusters fell into subclusters. The network were visualized in alignment score weighted Prefuse Force-Directed Layout using Cytoscape 3.9 50 . The resulted SSNs were submitted to EFI-EST webtool 30 for genome neighborhood analysis with default setting. Nodes were highlighted in SSNs if the neighboring Pfam families, with maximal median distance of no more than 4 and minimal co-occurrence rate larger than 0.2, have NTPase activity, NTP binding, or DNA binding activities. Sequence analyses For sequence analyses, the BLAST tools 51 (with e-value cutoff of 10 -10 and query coverage of 25%) and HHpred 41 , and the resources of BV-BRC 38 were routinely used. In total, the data from 20279 bacterial and archaeal genomes were retrieved from the BIO-ML 26 , GMbC 27 , Unified Human Gastrointestinal Genome (UHGG) 28 , and BV-BRC representative bacteria as collected in Jan 2021. The queries used in this study were listed in Supplementary Table S12 , including IscS (GeneBank AIF64277.1), DndB (GeneBank AIF62361.1), DndC (GeneBank AIF62362.1), DndD (GeneBank AIF62363.1), DndE (GeneBank AIF62364.1) in E.coli B7A and SspA (Refseq WP_016789103.1), SspB (Refseq WP_022570853.1), SspC (Refseq WP_016789109.1), SspD (Refseq WP_016789110.1), SspE (Refseq WP_016789111.1), and BrxP (GMbC tag OIFBNFKG_02855), BrxC (GMbC tag OIFBNFKG_02856), BrxZ (GMbC tag OIFBNFKG_02857), BrxL (GMbC tag OIFBNFKG_02858). The genes dndCD , sspBCD , and brxPCZL were considered to be present when all genes were adjacent. Phylogeny tree We first built a concatenated alignment of 10 nearly universal and single-copy ribosomal protein families. We used Diamond v0.8.22 52 (with parameters blastx -more-sensitive -e 0.000001 -id 35 -query-cover 80) to BLAST all proteomes in out collection against the RiboDB database v1.4.1 53 of bacterial ribosomal protein genes. We excluded proteins bL17, bS16, bS21, uL22, uS3 and uS4, as they were not sufficiently distributed across all genomes. In each RiboDB gene family, we excluded genomes that contained gene duplicates. Then, we aligned all protein families individually with MUSCLE v5.1 54 (with default setting). We filtered out misaligned sites using BMGE v1.12 55 (with parameters -t AA -g 0.95 -m BLOSUM30) and concatenated all individual alignments using Seaview v4.7 56 . The phylogenomic tree was reconstructed using FastTree v2.1.10 57 (with parameters -lg -gamma) and visualized and modified in iTOL 58 . DNA isolation The cells were resuspended in 500 μL of phosphate buffered saline upon receiving and re-pelleted by centrifugation at 6,000 g for 5 min at 4 °C. Genomic DNA was isolated using E.Z.N.A Bacterial DNA kit (Catalog # D3350-02), with the bead beating option (speed of 4 m/s, 45 s “on”, 1 min rest, 3 cycles). DNA was eluted with 150-200 μL of RNase free water and stored at -80 °C. Digestion of DNA for LC-MS/MS analysis of PT dinucleotides DNA (20 μg, 78 μL) was incubated with Nuclease P1 (1.5 U, 3 μL, US Biological) in 30 mM ammonium acetate pH 5.3 and 0.5 mM ZnCl 2 (90 μL total reaction volume) for 2 h at 55 °C. The reaction mixture was diluted with Tris-HCl (100 mM final concentration, pH 8.0, 9 μL) and incubated with calf intestinal alkaline phosphatase (51 U, 3 μL, Sigma) for 2 h at 37 °C. Enzymes were removed by passing the mixture through a VWR 10 kDa spin filter with centrifugation at 12,000 g for 12 min. The solution was lyophilized to dryness and resuspended in H 2 O (50 μL). LC-MS/MS analysis of DNA PT dinucleotides Synthetic PT DNA dinucleotides or nuclease P1 hydrolyzed DNA were analyzed by LC-MS/MS on an Agilent 1290 series HPLC system equipped with a Synergi Fusion RP column (2.5 μm particle size, 100 Å pore size, 100 mm length, 2 mm inner diameter) and a DAD. The HPLC was coupled to an Agilent 6490 triple quadrupole mass spectrometer. The column was eluted at 0.35 mL/min at 35 °C with a linear gradient of 3-9% acetonitrile in 97% solvent A (5 mM ammonium acetate pH 5.3) over 15 min. The column was rinsed with 95% acetonitrile in solvent A for 1 min, and then the initial conditions were regenerated by rinsing the column with 97% solvent A for 3 min. Canonical deoxyribonucleosides that eluted from the column were quantified by their 260 nm absorbance with the DAD. PT-containing dinucleotides were identified and quantified by tandem quadrupole mass spectrometry with electrospray ionization operated with the following parameters: N2 temperature, 200 °C; N2 flow rate, 14 L/min; nebulizer pressure, 20 psi; capillary voltage, 1800 V; and fragmentor voltage, 380 V. For product identification, the mass spectrometer was operated in positive ion multiple reaction monitoring mode using the conditions tabulated in Supplementary Table S13 . High-resolution mass spectrometry Synthetic PT dinucleotides (2 pmol per 10 μL injection) or Nuclease P1 hydrolyzed RNA (4 μg per 10 μL injection) were analyzed on a Dionex Ultimate 3000 UHPLC system equipped with a Synergi Fusion RP column (2.5 μm particle size, 100 Å pore size, 100 mm length, 2 mm inner diameter). The HPLC was coupled to a Thermo Fisher Q Exactive Hybrid Quadrupole-Orbitrap mass spectrometer. The column was eluted at 0.35 mL/min at 35 °C with a linear gradient of 3-9% acetonitrile in 97% solvent A (5 mM ammonium acetate pH 5.3) over 15 min. The column was rinsed with 95% acetonitrile in solvent A for 1 min, and then the initial conditions were regenerated by rinsing the column with 97% solvent A for 3 min. High resolution mass spectra for the PT-containing dinucleotides were obtained by hybrid quadrupole-Orbitrap mass spectrometry with the following parameters: sheath gas flow rate, 50 L/min; aux gas flow rate, 15 L/min; sweep gas flow rate, 3 L/min; spray voltage, 4.20 kV; and capillary temperature, 275 °C. For product identification, the mass spectrometer was operated in positive ion targeted single ion monitoring mode using the conditions tabulated in Supplementary Table S14 . PT-seq library preparation PT-seq applied on Lachnospiraceae sp. was adapted from previously described version 45 ( Supplementary Fig. S4B ). 10 ug of genomic DNA was diluted in 500 μl of ddH2O in 2ml Eppie tube on ice and fragmented by probe sonication for 5 cycles (2 min ON with 50% duty cycle and ∼20 output control and 1 min OFF). The samples were subjected to SpeedVac concentration for ∼2.5 h at ∼450 mTorr to the final volume of ∼40 μL. Blocking of pre-existing strand-break sites was achieved by 4 cycles of denaturation-dephosphorylation-blocking. Each cycle started by denaturating at 94 °C for 2 min and immediately cooling down on ice for 2 min. The initial dephosphorylation reaction was in a mixture (50 μL) containing 5 μl of terminal transferase buffer (NEB Tdt Reaction Buffer, Catalog # M0315S), 1 μL of shrimp alkaline phosphatase (rSAP, NEB Catalog # M0371S), and 5 μg of fragmented DNA, with incubation at 37°C for 30min to remove phosphate at 3′ end of the strand-breaks. The phosphatase was then inactivated by heating at 65 °C for 10 min. After cooling, 1 μL of Tdt Reaction Buffer, 6 μL of CoCl 2 (0.25 mM), 2 μL of ddNTPs (2 mM each, TriLink) and 1 μL of terminal transferase (20 units, NEB Catalog # M0315S) was added to the reaction with incubation at 37 °C for 1 h to block any pre-existing strand-break sites. Blocking cycles were repeated 3 times, in which fresh reagents were added in each cycle. The blocked DNA was purified using a DNA cleanup kit (Zymo Catalog # 11-304C). For iodine cleavage, 32 μL of the blocked DNA was incubated with 4 μL of 500 mM Tris-HCl pH 9.0 and 4 μL of iodine solution (50 mM, Fluka, Catalog # 318981-100) at room temperature for 5 min. Then, the reaction product was purified using two DyeEX columns (QIAGEN Catalog # 63206) to remove salts and iodine. The purified DNA was denatured by incubating at 94 °C for 2 min and cooling on ice for 2 min and then incubated with 5 μL of NEB rCutsmart buffer and 1 μl of rSAP (50 μL reaction) to remove 3′-phosphates arising from iodine cleavage. After incubation at 37°C for 30 min and 65 °C for another 10 min, the product was incubated with 1 μL of dTTP (1 mM, NEB Catalog # N)443S), 1 μL of Tdt buffer, 6 μL of CoCl 2 , 1 μL of Tdt and 1 μL of H 2 O. After incubation at 37 °C for 45 min, the dTTP was removed by DyeEX columns. PT-seq applied on B. faecalis and B. salyersiae was adapted by adding biotin labeling and streptavidin bead-capturing ( Supplementary Fig. S4A ), to resolve the inadequate poly(dT) tailing of iodine cleaved nicks. DNA (10 μg) was subjected to blocking of pre-existing strand-break sites, iodine cleavage, and dTTP labeling as described above. DNA (60 μL) was then incubated with 7.8 uL of Tdt buffer, 7.8 uL of CoCl 2 , 1 μL of ddUTP-biotin (1 mM, Jena Bioscience NU-1619-BIOX-S), 1 μL of Tdt and 1 μL of H 2 O at 37 °C for 1 h to terminate the T-tails by ddUTP-biotin. After cleaned with DyeEX columns, the DNA was diluted in 500 μL of H 2 O and fragmented by probe sonication as described above. The DNA fragments was mixed with 10 μL of streptavbidin coated meganetic beads (NEB Catalog # S1421S) and 500 μL of binding buffer (5 mM Tris-HCl, pH7.5, 1 M NaCl, 0.5 mM EDTA) and incubated on a shaker at ambient temperature for 1 h. The beads were pull down by a magnetic and washed 3 times with 100 μL of binding buffer. After discarding the supernatant, the beads were resuspended in 20 μL of H 2 O. The purified product and captured beads in upgraded ICDS method were subjected to Illumina library preparation with SMART ChIp-seq kit (Takara, Catalog # 634865) by following the manufacturer’s protocol. The final step of PCR was performed using the Illumina primers provided in ChIp-seq kit and 12 cycles were used for amplification. The PCR product with unique sequencing barcode was submitted to Illumina MiSeq and NovaSeq instrument for 150 bp paired-end sequencing. Data analysis As illustrated in the workflow ( Fig. S4 ), the data analysis started with triming of Adapters. Adapters were removed using bbduk from BBtools (sourceforge.net/projects/bbmap/), (with parameters ktrim=r k=18 hdist=2 hdist2=1 rcomp=f mink=8 qtrim=r trimq=30 for R1, ktrim=r k=18 mink=8 hdist=1 rcomp=f qtrim=r trimq=30 for R2). T-tails were removed using bbduk with parameters ktrim=r k=15 hdist=1 rcomp=f mink=8 for R1 and ktrim=l k=15 hdist=1 rcomp=f mink=8 for R2. PhiX were removed using bbduk with parameters k=31 hdist=1. Trimmed reads were aligned to the corresponding genome using Bowtie 2 with setting “sensitive”. The bam files were cleaned using SMARTcleaner 59 and split into two strands using samtools 60 . The coverage of each strand was calculated separately using bedtools 61 with the genomecov -d -5 option. Then, custom shell scripts were used for pileup calling. Briefly, all read start positions were recorded for both strands separately. The read pileup depth at each position was represented by the number of read starting (5ʹ-end) at each position. The 13 nt sequences centered at positions with depth ≥1 were retrieved using samtools 60 . The consensus motifs were analyzed with incrementing depth, typically from 1 to 100 with step of 10, using MEME 62 with parameters -dna -objfun classic -nmotifs 5 - mod zoops -evt 0.05 -minw 3 -maxw 6 -markov_order 0 -nostatus -oc. The read starting sites mapping within less than 3 bps were collapsed into the centremost consensus motif sites. Regions where starting positions mapped within 4 bp on opposite strands and within reverse complementary sequences were considered as double stranded modifications. Data availability Custom scripts for processing the sequencing data are described in Methods and are available upon request. Sequencing data have been deposited in NCBI SRA database under BioProject ID PRJNA1006039. Acknowledgements We thank Susan Weir and Katya Moniz at the Openbiome for providing us with bacterial isolates. This work was supported by National Institutes of Health (R01 ES031576), by a NIEHS Training Grant in Environmental Toxicology T32-ES007020 (S.R.B), and by funding for the GMbC from the MIT Center for Microbiome Therapeutics and the Neil and Anna Rasmussen Family Foundation. Footnotes https://www.ncbi.nlm.nih.gov/sra References 1. ↵ Sanchez-Romero , M.A. & Casadesus , J. The bacterial epigenome . Nat Rev Microbiol 18 , 7 – 20 ( 2020 ). OpenUrl CrossRef PubMed 2. ↵ Thiaville , J.J. et al. Novel genomic island modifies DNA with 7-deazaguanine derivatives . Proc Natl Acad Sci U S A 113 , E1452 – 1459 ( 2016 ). OpenUrl Abstract / FREE Full Text 3. ↵ Zhou , X. et al. A novel DNA modification by sulphur . Mol Microbiol 57 , 1428 – 1438 ( 2005 ). OpenUrl CrossRef PubMed Web of Science 4. ↵ Wang , L. et al. Phosphorothioation of DNA in bacteria by dnd genes . Nat Chem Biol 3 , 709 – 710 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 5. ↵ Wang , L. , Jiang , S. , Deng , Z. , Dedon , P.C. & Chen , S . DNA phosphorothioate modification—a new multi-functional epigenetic system in bacteria . FEMS Micro Rev 43 , 109 – 122 ( 2019 ). OpenUrl CrossRef PubMed 6. ↵ Xu , T. , Yao , F. , Zhou , X. , Deng , Z. & You , D . A novel host-specific restriction system associated with DNA backbone S-modification in Salmonella . Nucleic Acids Res 38 , 7133 – 7141 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 7. ↵ Xiong , X. et al. SspABCD-SspE is a phosphorothioation-sensing bacterial defence system with broad anti-phage activities . Nat Microbiol 5 , 917 – 928 ( 2020 ). OpenUrl 8. ↵ Wang , S. et al. SspABCD-SspFGH Constitutes a New Type of DNA Phosphorothioate-Based Bacterial Defense System . mBio 12 ( 2021 ). 9. ↵ Jian , H. et al. The origin and impeded dissemination of the DNA phosphorothioation system in prokaryotes . Nat Commun 12 , 6382 ( 2021 ). OpenUrl 10. ↵ Xiong , L. et al. A new type of DNA phosphorothioation-based antiviral system in archaea . Nat Commun 10 , 1 – 11 ( 2019 ). OpenUrl CrossRef PubMed 11. ↵ Xiong , L. et al. A new type of DNA phosphorothioation-based antiviral system in archaea . Nat Commun 10 , 1688 ( 2019 ). OpenUrl CrossRef 12. ↵ Tong , T. et al. Occurrence, evolution, and functions of DNA phosphorothioate epigenetics in bacteria . Proc Natl Acad Sci U S A 115 , E2988 – E2996 ( 2018 ). OpenUrl Abstract / FREE Full Text 13. ↵ Cao , B. et al. Genomic mapping of phosphorothioates reveals partial modification of short consensus sequences . Nat Commun 5 , 3951 ( 2014 ). OpenUrl CrossRef PubMed 14. Chen , C. et al. Convergence of DNA methylation and phosphorothioation epigenetics in bacterial genomes . Proc Natl Acad Sci U S A 114 , 4501 – 4506 ( 2017 ). OpenUrl Abstract / FREE Full Text 15. ↵ Tong , T. et al. Occurrence, evolution, and functions of DNA phosphorothioate epigenetics in bacteria . Proc Natl Acad Sci U S A 115 , E2988 – E2996 ( 2018 ). OpenUrl Abstract / FREE Full Text 16. ↵ You , D.L. , Wang , L.R. , Yao , F. , Zhou , X.F. & Deng , Z.X . A novel DNA modification by sulfur: DndA is a NifS-like cysteine desulfurase capable of assembling DndC as an iron-sulfur cluster protein in Streptomyces lividans . Biochemistry 46 , 6126 – 6133 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 17. ↵ Yao , F. , Xu , T. , Zhou , X. , Deng , Z. & You , D . Functional analysis of spfD gene involved in DNA phosphorothioation in Pseudomonas fluorescens Pf0-1 . FEBS Lett 583 , 729 – 733 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 18. ↵ Dai , D. et al. DNA phosphorothioate modification plays a role in peroxides resistance in Streptomyces lividans . Front Micro 7 , 1 – 13 ( 2016 ). OpenUrl 19. ↵ Huang , Q. et al. Defense Mechanism of Phosphorothioated DNA under Peroxynitrite-Mediated Oxidative Stress . ACS Chem Biol 15 , 2558 – 2567 ( 2020 ). OpenUrl 20. ↵ Kellner , S. et al. Oxidation of phosphorothioate DNA modifications leads to lethal genomic instability . Nat Chem Biol 13 , 888 – 894 ( 2017 ). OpenUrl CrossRef 21. ↵ Bhattacharyya , A. , Chattopadhyay , R. , Mitra , S. & Crowe , S.E . Oxidative stress: an essential factor in the pathogenesis of gastrointestinal mucosal diseases . Physiol Rev 94 , 329 – 354 ( 2014 ). OpenUrl CrossRef PubMed 22. Zhu , S. et al. Development of Methods Derived from Iodine-Induced Specific Cleavage for Identification and Quantitation of DNA Phosphorothioate Modifications . Biomolecules 10 ( 2020 ). 23. ↵ Mangerich , A. et al. Infection-induced colitis in mice causes dynamic and tissue-specific changes in stress response and DNA damage leading to colon cancer . Proc Natl Acad Sci U S A 109 , E1820 – 1829 ( 2012 ). OpenUrl Abstract / FREE Full Text 24. ↵ Sun , Y. et al. DNA Phosphorothioate Modifications Are Widely Distributed in the Human Microbiome . Biomolecules 10 ( 2020 ). 25. ↵ Byrne , S ., et al. Temporal dynamics and metagenomics of phosphorothioate epigenomes in the human gut microbiome bioRxiv ( 2024 ). 26. ↵ Poyet , M. et al. A library of human gut bacterial isolates paired with longitudinal multiomics data enables mechanistic microbiome research . Nat Med 25 , 1442 – 1452 ( 2019 ). OpenUrl PubMed 27. ↵ Groussin , M. et al. Elevated rates of horizontal gene transfer in the industrialized human microbiome . Cell 184 , 2053 – 2067 e2018 ( 2021 ). OpenUrl 28. ↵ Almeida , A. et al. A unified catalog of 204,938 reference genomes from the human gut microbiome . Nat Biotechnol 39 , 105 – 114 ( 2021 ). OpenUrl PubMed 29. ↵ Yuan , Y. , DeMott , M. , Byrne , S. & Dedon , P. PT-seq for highly sensitive metagenomic mapping of phosphorothioate DNA modifications. bioRxiv ( 2024 ). 30. ↵ Oberg , N. , Zallot , R. & Gerlt , J.A. EFI-EST , EFI-GNT, and EFI-CGFP: Enzyme Function Initiative (EFI) Web Resource for Genomic Enzymology Tools . J Mol Biol 435 , 168018 ( 2023 ). 31. ↵ Kambampati , R. & Lauhon , C.T . Evidence for the transfer of sulfane sulfur from IscS to ThiI during the in vitro biosynthesis of 4-thiouridine in Escherichia coli tRNA . J Biol Chem 275 , 10727 – 10730 ( 2000 ). OpenUrl Abstract / FREE Full Text 32. ↵ Shigi , N . Biosynthesis and functions of sulfur modifications in tRNA . Front Genet 5 , 67 ( 2014 ). 33. ↵ Bouvier , D. et al. TtcA a new tRNA-thioltransferase with an Fe-S cluster . Nucleic Acids Res 42 , 7960 – 7970 ( 2014 ). OpenUrl CrossRef PubMed 34. ↵ Goldfarb , T. et al. BREX is a novel phage resistance system widespread in microbial genomes . EMBO J 34 , 169 – 183 ( 2015 ). OpenUrl Abstract / FREE Full Text 35. ↵ Yao , J. et al. Identification and characterization of a HEPN-MNT family type II toxin-antitoxin in Shewanella oneidensis . Microb Biotechnol 8 , 961 – 973 ( 2015 ). OpenUrl 36. ↵ Songailiene , I. et al. HEPN-MNT Toxin-Antitoxin System: The HEPN Ribonuclease Is Neutralized by OligoAMPylation . Mol Cell 80 , 955 – 970 e957 ( 2020 ). OpenUrl 37. ↵ Yao , J. et al. Novel polyadenylylation-dependent neutralization mechanism of the HEPN/MNT toxin/antitoxin system . Nucleic Acids Res 48 , 11054 – 11067 ( 2020 ). OpenUrl 38. ↵ Olson , R.D. et al. Introducing the Bacterial and Viral Bioinformatics Resource Center (BV-BRC): a resource combining PATRIC, IRD and ViPR . Nucleic Acids Res 51 , D678 – D689 ( 2023 ). OpenUrl 39. ↵ Schirrmeister , B.E. , Gugger , M. & Donoghue , P.C . Cyanobacteria and the Great Oxidation Event: evidence from genes and fossils . Palaeontology 58 , 769 – 785 ( 2015 ). OpenUrl CrossRef GeoRef PubMed 40. ↵ Luo , G. et al. Rapid oxygenation of Earth’s atmosphere 2.33 billion years ago . Sci Adv 2 , e1600134 ( 2016 ). OpenUrl FREE Full Text 41. ↵ Zimmermann , L. et al. A Completely Reimplemented MPI Bioinformatics Toolkit with a New HHpred Server at its Core . J Mol Biol 430 , 2237 – 2243 ( 2018 ). OpenUrl CrossRef PubMed 42. ↵ Wang , L. et al. DNA phosphorothioation is widespread and quantized in bacterial genomes . Proc Natl Acad Sci U S A 108 , 2963 – 2968 ( 2011 ). OpenUrl Abstract / FREE Full Text 43. ↵ Liu , G. et al. Cleavage of phosphorothioated DNA and methylated DNA by the type IV restriction endonuclease ScoMcrA . PLoS Genet 6 , e1001253 ( 2010 ). OpenUrl CrossRef PubMed 44. ↵ Wu , X. et al. Epigenetic competition reveals density-dependent regulation and target site plasticity of phosphorothioate epigenetics in bacteria . Proc Natl Acad Sci U S A 117 , 14322 – 14330 ( 2020 ). OpenUrl Abstract / FREE Full Text 45. ↵ Cao , B. et al. Nick-seq for single-nucleotide resolution genomic maps of DNA modifications and damage . Nucleic Acids Res ( 2020 ). 46. ↵ Jiang , S. et al. A DNA phosphorothioation-based Dnd defense system provides resistance against various phages and is compatible with the Ssp defense system . mBio 14 , e0093323 ( 2023 ). OpenUrl 47. ↵ Garcia-Bayona , L. et al. Nanaerobic growth enables direct visualization of dynamic cellular processes in human gut symbionts . Proc Natl Acad Sci U S A 117 , 24484 – 24493 ( 2020 ). OpenUrl Abstract / FREE Full Text 48. ↵ Pantosti , A. , Tzianabos , A.O. , Onderdonk , A.B. & Kasper , D.L . Immunochemical characterization of two surface polysaccharides of Bacteroides fragilis . Infect Immun 59 , 2075 – 2082 ( 1991 ). OpenUrl Abstract / FREE Full Text 49. ↵ Smith , C.J. & Callihan , D.R . Analysis of rRNA restriction fragment length polymorphisms from Bacteroides spp. and Bacteroides fragilis isolates associated with diarrhea in humans and animals . J Clin Microbiol 30 , 806 – 812 ( 1992 ). OpenUrl Abstract / FREE Full Text 50. ↵ Shannon , P. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks . Genome Res 13 , 2498 – 2504 ( 2003 ). OpenUrl Abstract / FREE Full Text 51. ↵ Altschul , S.F. , Gish , W. , Miller , W. , Myers , E.W. & Lipman , D.J . Basic local alignment search tool . J Mol Biol 215 , 403 – 410 ( 1990 ). OpenUrl CrossRef PubMed Web of Science 52. ↵ Buchfink , B. , Reuter , K. & Drost , H.G . Sensitive protein alignments at tree-of-life scale using DIAMOND . Nat Methods 18 , 366 – 368 ( 2021 ). OpenUrl CrossRef PubMed 53. ↵ Jauffrit , F. et al. RiboDB Database: A Comprehensive Resource for Prokaryotic Systematics . Mol Biol Evol 33 , 2170 – 2172 ( 2016 ). OpenUrl CrossRef PubMed 54. ↵ Edgar , R.C . MUSCLE: multiple sequence alignment with high accuracy and high throughput . Nucleic Acids Res 32 , 1792 – 1797 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 55. ↵ Criscuolo , A. & Gribaldo , S . BMGE (Block Mapping and Gathering with Entropy): a new software for selection of phylogenetic informative regions from multiple sequence alignments . BMC Evol Biol 10 , 210 ( 2010 ). 56. ↵ Gouy , M. , Guindon , S. & Gascuel , O . SeaView version 4: A multiplatform graphical user interface for sequence alignment and phylogenetic tree building . Mol Biol Evol 27 , 221 – 224 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 57. ↵ Price , M.N. , Dehal , P.S. & Arkin , A.P . FastTree: computing large minimum evolution trees with profiles instead of a distance matrix . Mol Biol Evol 26 , 1641 – 1650 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 58. ↵ Letunic , I. & Bork , P . Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation . Nucleic Acids Res 49 , W293 – W296 ( 2021 ). OpenUrl CrossRef PubMed 59. ↵ Zhao , D. & Zheng , D . SMARTcleaner: identify and clean off-target signals in SMART ChIP-seq analysis . BMC Bioinformatics 19 , 544 ( 2018 ). 60. ↵ Li , H. et al. The Sequence Alignment/Map format and SAMtools . Bioinformatics 25 , 2078 – 2079 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 61. ↵ Quinlan , A.R. & Hall , I.M . BEDTools: a flexible suite of utilities for comparing genomic features . Bioinformatics 26 , 841 – 842 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 62. ↵ Bailey , T.L. , Johnson , J. , Grant , C.E. & Noble , W.S . The MEME Suite . Nucleic Acids Res 43 , W39 – W49 ( 2015 ). OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted June 03, 2024. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Phosphorothioate DNA modification by BREX Type 4 systems in the human gut microbiome Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Phosphorothioate DNA modification by BREX Type 4 systems in the human gut microbiome Yifeng Yuan , Michael S. DeMott , Shane R. Byrne , Katia Flores , Mathilde Poyet , Mathieu Groussin , Global Microbiome Conservancy , Brittany Berdy , Laurie Comstock , Eric J. Alm , Peter C. Dedon bioRxiv 2024.06.03.597175; doi: https://doi.org/10.1101/2024.06.03.597175 Share This Article: Copy Citation Tools Phosphorothioate DNA modification by BREX Type 4 systems in the human gut microbiome Yifeng Yuan , Michael S. DeMott , Shane R. Byrne , Katia Flores , Mathilde Poyet , Mathieu Groussin , Global Microbiome Conservancy , Brittany Berdy , Laurie Comstock , Eric J. Alm , Peter C. Dedon bioRxiv 2024.06.03.597175; doi: https://doi.org/10.1101/2024.06.03.597175 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Microbiology Subject Areas All Articles Animal Behavior and Cognition (7652) Biochemistry (17752) Bioengineering (13936) Bioinformatics (42084) Biophysics (21501) Cancer Biology (18655) Cell Biology (25586) Clinical Trials (138) Developmental Biology (13410) Ecology (19949) Epidemiology (2067) Evolutionary Biology (24378) Genetics (15639) Genomics (22562) Immunology (17779) Microbiology (40505) Molecular Biology (17219) Neuroscience (88825) Paleontology (667) Pathology (2845) Pharmacology and Toxicology (4840) Physiology (7666) Plant Biology (15182) Scientific Communication and Education (2048) Synthetic Biology (4305) Systems Biology (9840) Zoology (2274)