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

Domain: neuroscience
Published: 2026-01-22
Editorial outcome: major_revision

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 revealed several previously undescribed PT synthesis systems, including type 4 Bacteriophage Exclusion (BREX) type 4 brx genes, which we 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 10 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.

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 N6-methyl-adenine (6mA), N4-methyl-cytosine (4mC), C5-methyl-cytosine (5mC), and 7-deazaguanine derivatives1,2, phosphorothioates (PT) are the only known modification of the sugar-phosphate backbone, with a non-bridging oxygen replaced by sulfur in Rp-specific configuration3–5 (Fig. 1A). Two PT-based restriction and modification (R-M) systems, the dnd4,6 and ssp7,8 gene clusters, have been observed in ~10% of bacteria and archaea9,10, with a third tdp system recently discovered in extremophilic bacteria11. Typically, the modification components are organized as a four-gene operon for dndBCDE and a three-gene operon for sspBCD, with an additional modification gene for dndA or sspA located adjacent to the operon, dispersed elsewhere in the genome, or replaced by a gene for a homolog such as iscS. However, the recently discovered tdpABC cluster requires only a single gene, tdpC, for PT synthesis11. Like methylation-based R-M systems8,10,12, DndACDE and SspABCD proteins catalyze PT modification on one or both strands of specific consensus sequences. For example, DndACDE confer double-stranded PTs at 5′-G AAC-3′/5′-G TTC-3′ sequences in Escherichia coli B7A and Salmonella enterica serovar Cerro 87, 5ʹ-G GCC-3′/5′-G GCC-3′ in Pseudomonas fluorescens pf0-1 and Streptomyces lividans 1326, and 5′-G ATC-3′/5′-G ATC-3′ in Hahella chejuensis KCTC 239612–14. SspABCD proteins, on the other hand, catalyze single-stranded 5′-C CA-3′ in Vibrio cyclitrophicus FF757,13. The dnd and ssp systems share some similarities, such as encoding a homolog of phosphoadenosine phosphosulfate (PAPS) reductase (DndC, SspD)7,15, a homolog of cysteine desulfurase (DndA, SspA)7,15, and a P-loop containing ATPase (DndD, SspC)7,16 (Fig. 1B). The restriction counterparts DndFGH and SspFGH sense PT modifications by poorly understood mechanisms, with only 10–15% of consensus sequences modified with PTs7,13.

Both the metabolic pathways enabling sulfur incorporation into PTs and the altered 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 gene transcription of redox homeostasis genes12. PTs also provide weak to modest protective effects in cells exposed to reactive oxygen and nitrogen species, such as peroxides12,17 and peroxynitrite18. Contrasting with this protection, PT-containing bacteria are 5-fold more sensitive to neutrophil-derived hypochlorous acid (HOCl) due to extensive DNA breaks at PTs19. 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 conditions20–22. Evidence for the presence of PT genes in bacterial strains associated with the human gut microbiome and PT dinucleotides in fecal DNA23,24 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)25, the Global Microbiome Conservancy (GMbC)26–28, and the Unified Human Gastrointestinal Genome (UHGG) collection29. Mass spectrometric analysis of PT dinucleotides in 226 of these isolates, coupled with PT-specific next-generation sequencing (PT-seq)30, led to the discovery of a previously undescribed PT system involving type 4 Bacteriophage Exclusion (BREX) genes brxPCZL. 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 previously uncharacterized PT modification systems by analyzing gene 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 analysis31 of dnd and ssp genes to find new PT modification systems. Here, we used enzyme function initiative (EFI) tools31 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 pyrophosphatase activity essential for PT biosynthesis (Fig. 1B), an activity shared by sulfur-inserting RNA modification enzymes ThiI32, MnmA33, and Ncs633, and TtcA34. 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 3120 DndC homologs and 2132 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 Data 1 and 2).

The resulting gene neighborhoods revealed several previously undescribed putative PT-modifying gene candidates. Each circle or node in Fig. 2 depicts a neighborhood containing dndC (Fig. 1A) and sspD (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 system35 (Fig. 2A, yellow nodes). Five of six BREX systems possess a BrxC/PglY ATPase, a PglX DNA methyltransferase, and BrxZ (PglZ) phosphatase. In BREX type 4 clusters containing brxP, the PglX adenine methyltransferase is replaced with the DndC S-inserting PAPS reductase domain protein35. For example, in Methanothermobacter sp. (Fig. 2A, yellow node, black outline), brxP (dndC homolog) is adjacent to a brxC ATPase gene, a brxZ 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 phosphatase, sspB nickase, and a methyltransferase. It is common to find several different defense systems in the same region due to horizontal gene transfer and genetic decay. The systems in these islands can function independently, so the methyltransferases located next to the candidate PT genes may function independently of the dnd or BREX type 4 systems.

The synthesis of PTs in the DndC-containing BREX type 4 system was subsequently validated genetically and by LC–MS, as discussed shortly.

In a second putative PT defense system, a small cluster including Bacillus cereus (Fig. 2A, light green node) 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) system36, in which the HEPN protein is an RNase toxin, but its activity is neutralized by adenylylation by the MNT antitoxin37,38. The HEPNs in the dndC clusters lack the RNase domain, but the neighboring MNT, ThiF, has the conserved motif GSX 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.

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 Cyanobacterium stanieri (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 6616 representative genomes from the Bacterial and Viral Bioinformatics Resource Center (BV-BRC; as of January 2021)39. 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 Data 3). Notably, both dndCD and sspBCD gene clusters co-occurred in the same genomes of some Cyanobacteria, such as Gloeocapsa sp., Coleofasciculus chthonoplastes, and Scytonema hofmanni, among others (Supplementary Data 3). 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 some Cyanobacteria genomes (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 epigenetics systems in the ancient photosynthetic Cyanobacteria phylum40,41.

BREX type 4 systems are homologous to Ssp proteins and catalyze PT synthesis
The gene 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 (HHpred42) of BREX family and SSP proteins shown in Fig. 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).

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 limited digest with nuclease P1 and phosphatase dephosphorylation43. 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 C, C C, and T C at 56, 228, and 98 per 106 nt, respectively (Fig. 3B, Supplementary Fig. S1A, Supplementary Data 4). As shown in Supplementary Data 4, 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 C (Prevotella sp.), A A and A C (putative Butyricimonas), T C (Parabacteroides), and A A (Bacteroides). These results clearly distinguished the BREX type 4 gene systems from the 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 Fig. 3B, loss of brxC abolished PTs, with PTs restored by complementation with in trans expression of B. salyersiae brxC (ΔbrxC brxC Bs). 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 constructed the PT pathway from B. salyersiae by cloning brxP Bs, mcrA (a PT-dependent restriction enzyme in Streptomyces coelicolor)44, and brxC Bs together in an expression plasmid or brxC Bs alone in the plasmid. These plasmids were transferred into Bacteroides thetaiotaomicron VPI, which lacks PT genes but possesses a SufS cysteine desulfurase homolog of DndA to provide a presumably complete system for PT modification (SufS, BrxP, BrxC) and restriction (McrA). The combination of brxP Bs, mcrA Bs, and brxC Bs resulted in the T C, C C, and A C dinucleotides in the same proportion as wild-type B. salyersiae. brxC Bs alone were not able to confer PTs, which is consistent with the essentiality of BrxP in PT synthesis. Furthermore, mcrA, brxZ, brxL or an unknown Bs02795 gene were not required for PT synthesis as shown in corresponding mutants (Fig. 3B). To test the generality of the PT synthesis by BrxCP, determine if brxC and sspC are interchangeable, and, since B. salyersiae and B. faecalis have different PT motifs, assess if BrxC is the component that determines the target motif, we heterologously expressed brxC Bf from B. faecalis or sspC Bo (brxC equivalent) from the human gut isolate B. ovatus CL03T12C18 in trans in the B. salyersiae brxC mutant. PTs were not detected in either strain (Fig. 3B). These results suggest species-specific complexities of the BREX system components requiring further investigation. We conclude that, in the 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 microbiome23,24, 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 the human gut microbiome, we searched for PT gene clusters in 13,663 human gut microbial genomes from the BIO-ML25 and the GMbC26–28 collection, and from isolate sequences of human gut microbes and metagenome-assembled genomes (MAGs) in the Unified Human Gastrointestinal Genome (UHGG) collection29. 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 Data 5). The essential gene sets dndCD, sspBCD, and brxPCZL were found to be present in 2.7%, 3.6%, and 1.4% of the human 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 Data 5). 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.

Validating the predicted PT modifications in gut microbiome isolates
Given the presence of genes essential for PT synthesis in ~1000 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 Data 4). In total, we identified 8 PT dinucleotides, including A A, A C, C A, C C, G A, G C, G T, and T C, and did not detect an additional two dinucleotides (C T, G G) that we found in human fecal DNA samples (Supplementary Data 4)24. Previous studies in three bacterial species with dnd and ssp genes4,43,45 revealed G A and G T at >500 per 106 nt and C C at >2000 per 106 nt, with C A, A A, A C, and T C detected at much lower levels of 1–6 per 106 nt4,43,45. The minor PT dinucleotides represent low-affinity binding sites for the DNA shape-selective Dnd proteins45. In sharp contrast, however, these minor sites in gut microbiome bacteria represent major PT modification sites, as high as 1800 per 106 nt in bacteria with brxPCZL genes (Supplementary Data 4).

The A A, A C, C A, C C, G A, G C, G T, and T 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 A and the predominance of dnd genes. The largest portion of Group 1 (32 isolates) possessed both G A and G T and harbored dndCD genes with or without dndBE, including isolates from the Bacteroidales and Clostridiales orders. Five isolates from both Bacteroidales and Clostridiales harboring the minimal dndCD gene set showed G A and G C. One isolate that possessed a solitary brxP but no dnd genes and possessed both G A and G 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 A (Supplementary Data 4).

The second group of 34 isolates was characterized by the C 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 A and G T at the detection limit of ~2 PTs per 106 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.

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 A, 12 isolates of Prevotella sp. contain C C, 3 isolates of Parabacteroides contain T C, and 9 isolates of putative Butyricimonas faecalis contain A A and A C (Fig. S1B). The B. salyersiae strain examined earlier, which contains C C, T C, and A C (Fig. 5, Supplementary Data 4).

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. It is possible but unlikely that the lack of PT dinucleotides detected in the isolates is due to a symbiont phenotype in the gut microbiome, with loss of one of the PT synthesis genes accounted for by uptake of PT intermediates from other bacteria. The PT modification pathway does not require unique small-molecule intermediates, and the modification genes are generally not interchangeable among different bacteria, so it is unlikely that gut microbes can lose one or more PT modification genes and exist as a “PT symbiont” in the gut community. Again, we cannot rule out errors in the reference genomes due to genome sequencing and sample handling.

Although based on a limited number of samples, these observations suggest the hypothesis that (1) Dnd proteins produce G A, G C, and G T, (2) Ssp proteins produce C C, and (3) Brx proteins produce everything except G A, G C, and G 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 an innovative and highly sensitive NGS sequencing technique, PT-seq30. The idea here is that the PT dinucleotides, such as G A and G T, are found in the larger consensus sequences of G ATC and 5'-G AAC-3ʹ/5ʹ-G TTC-3' in several types of bacteria5. Application of PT-seq (Supplementary Fig. S4) to the BREX type 4-containing human gut microbiome isolate B. salyersiae DSM18765, earlier found by LC-MS/MS to possess A C, C C, and T C at 6, 228, and 98 per 106 nt (Supplementary Data 4), showed 5459 PT sites: 56 at A CTC, 3881 at C CTC, and 1513 at T CTC sites, respectively (Fig. 6A, Supplementary Data 6). As observed with other dnd and ssp systems, the Brx proteins modified only a portion of the 27274 A CTC (0.2%), 22698 C CTC (17%), and 36290 T CTC (4%) total sites available, respectively. The observation of PTs at 3N CTC sites, with a strong preference for C CTC, is consistent with previous observations that Dnd proteins select their target sequences based on DNA shape rather than precise binding contacts45. We appreciate the discrepancy between LC–MS/MS and PT-seq data for the quantity of A C in the A CTC context in B. salyersiae DSM18765. We previously observed that the efficiency of iodine-induced cleavage of PT sites to strand breaks in the PT-seq method varies by 2-fold for different dinucleotides30 and that the extension efficiency of the terminal transferase used for T-tailing depends on the nucleobase sequence of the initiators, where A-ending initiators are generally less favored46–48. Together, these could account for some portion of the discrepancy between the LC–MS and PT-seq data for ApsC. However, in another BREX type 4-containing putative B. faecalis (GMbC ID 5893AJ_0218_015_F2), PT-seq revealed 9832 GA AG and 5536 GA CG (Supplementary Data 7). These sites agreed with the A A and A C dinucleotides detected by LC–MS/MS (Supplementary Data 4). It thus remains to be determined if the A CTC motif is uniquely problematic for PT-seq. Finally, PT-seq also detected sites with bistranded PTs: 1080 GA AG/CT TC and 667 GA CG/GC TC. However, we could not detect reliable signals for the corresponding T T and T G dinucleotides by LC–MS/MS, which may be due to their low abundance and to insensitive MS detection of T-containing nucleotides.

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 1474 G AGC and G CTC sites occurring among the 16,683 possible GAGC/GCTC sites, with 1154 modified on one strand and 160 modified on both strands (Fig. 6A, Supplementary Data 8). A large fraction of pileup sites occurring at NTTT and other motifs represents artifacts of poly(dT) tailing at single-strand nicks by 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 non-randomly distributed across bacterial genomes, which agrees with single-molecule real-time (SMRT) sequencing maps in Escherichia coli13. 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 Data 9). The percentage of consensus sequences modified with PTs, which is roughly 10% in genomes studied to date12,13, also varied in individual bacterial genomes, ranging from 1% in tRNA genes in Lachnospiraceae sp. to 40% in rRNA genes in B. faecalis (Supplementary Data 9). 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. 5SB, black line; Supplementary Fig. S6). Published PT-seq data for E. coli BW25113, which has a bistranded G AAC/G 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 transcription13, 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.

Methods
Bacterial strains and growth conditions
Bacteria revived from BIO-ML and GMbC were cultured as previously decribed25. Parabacteroides spp. and Bacteroides spp., E. faecalis, and BIO-ML isolates were grown on ASF plates (Becton Dickinson) and BHIS plates52, 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 2.5-3.5% hydrogen, 10% carbon dioxide, and the remainder nitrogen. Growth on agar plates was harvested with phosphate-buffered saline (Catalog #10010-023, Gibco, Paisley, PA). The bacterial suspensions were pelleted by centrifugation at 4000 × g 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 BREX4 system deletion and reconstruction are listed in Supplementary Table S1. Bacteroidales strains were grown in basal liquid medium53 or on BHIS plates52. Antibiotics used for selection include erythromycin (10 µg/mL), gentamicin (200 µg/mL), and 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, including mcrA (HMPREF1532_02792), brxC (HMPREF1532_02793), brxZ (HMPREF1532_02794), brxL (HMPREF1532_02795), and Bs02796 (HMPREF1532_02796), were constructed by amplifying DNA upstream and downstream of each gene using the synthetic primers listed in Supplementary Table S2. These flanking pieces were cloned into BamHI-digested pLBG1352 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 ΔbrxC were PCR amplified and cloned into BamHI-digested pFD34054 using NEBuilder v2.5. Transformants were PCR screened, and the plasmid was sequenced. Plasmids were conjugally transferred from E. coli to B. salyersiae 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 web tool v2021_0331 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 ~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 was visualized in alignment score weighted Prefuse Force-Directed Layout using Cytoscape v3.955. The resulting SSNs were submitted to EFI-EST web tool31 for genome neighborhood analysis with default settings. Nodes were highlighted in SSNs if the neighboring Pfam families, with a maximal median distance of no more than 4 and a 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 v2.956 (with E-value cutoff of 10-10 and query coverage of 25%) and HHpred v2020071742, and the resources of BV-BRC vbeta39 were routinely used. In total, the data from 20,279 bacterial and archaeal genomes were retrieved from the BIO-ML25, GMbC26–28, Unified Human Gastrointestinal Genome (UHGG)29, and BV-BRC representative bacteria as collected in January 2021. The queries used in this study were listed in Supplementary Data 10, 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.2257 (with parameters blastx -more-sensitive -e 0.000001 -id 35 -query-cover 80) to BLAST all proteomes in our collection against the RiboDB database v1.4.158 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.159 (with default settings). We filtered out misaligned sites using BMGE v1.1260 (with parameters -t AA -g 0.95 -m BLOSUM30) and concatenated all individual alignments using Seaview v4.761. The phylogenomic tree was reconstructed using FastTree v2.1.1062 (with parameters -lg -gamma) and visualized and modified in iTOL v563.

DNA isolation
The cells were resuspended in 500 µL of phosphate-buffered saline upon receiving and re-pelleted by centrifugation at 6000 × 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 ZnCl2 (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 H2O (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 S3. No statistical method was used to predetermine sample size. Comparisons between Bacteroides WT and engineered strains were based on three biological replicates analyzed using a t-test. For the remaining species, only one sample was analyzed owing to limited material availability, and the most recent technical replicate was included in this study.

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 S4.

PT-seq library preparation
PT modifications were mapped in the genomes of Lachnospiraceae sp., B. faecalis, and B. salyersiae by PT-seq as described elsewhere29 (Supplementary Fig. S4). DNA (10 µg) was diluted to 30.5 µL in filtered H2O and subjected to four blocking cycles, each consisting of the following: (a) DNA was heated at 94°C for 2 min to denature and immediately cooled on ice for 2 min. (b) rSAP (1 µL) and rCutSmart buffer (3.5 µL) were added, the reaction was incubated for 30 min at 37°C, and the phosphatase was heat deactivated for 10 min at 65°C. (c) Four ddNTPs (1 µL, 2 mM each), TdT reaction buffer (1.5 µL), CoCl2 (5 µL 0.25 mM), terminal transferase (1 µL, 20 units), and filtered H2O (3.5 µL) were added and the reaction incubated for 60 min at 37°C. Fresh reagents were added to each subsequent cycle at specific steps: rSAP (1 µL) and rCutSmart buffer (0.3 µL) were added at Step b, while ddNTPs (0.3 µL each), Tdt reaction buffer (0.3 µL), CoCl2 (0.3 µL), and terminal transferase (1 µL) were added at Step c. After all cycles were complete, excess unincorporated ddNTPs were removed by treatment with a DNA Clean & Concentrator kit (Zymo #11-304C). For iodine cleavage of PTs, the blocked DNA (32 µL) was incubated with 500 mM Tris–HCl pH 9.0 (4 µL) and iodine solution (4 µL, 5 mM) (Fluka #318981-100) for 10 min at 65°C and cooled to 4°C at a speed of 0.1°C/s. The reaction was then split into two 20 µL aliquots, and each filtered using a single DyeEX column (QIAGEN #63206) to remove salts and iodine, followed by recombining eluents. Filtered DNA was again denatured by heating at 94°C for 2 min and then cooling on ice for 2 min. Thereafter, the reaction was adjusted with rCutSmart buffer (5 µL), rSAP (1 µL), and filtered H2O to a final volume of 50 µL, and incubated for 30 min at 37°C to remove 3ʹ-terminal phosphates arising from iodine cleavage. The phosphatase enzyme was then inactivated for 10 min at 65°C. The reaction was further adjusted with dTTP (1 µL, 1 mM) (NEB #N0443S), TdT reaction buffer (1 µL), CoCl2 (6 µL, 0.25 mM), filtered H2O (1 µL), and terminal transferase (1 µL, 20 units), and incubated for 45 min at 37°C to add dT-tails to the 3’-terminal ends created specifically by iodine at PT sites. The reaction was split into three 20 µL aliquots, excess dTTP was removed using three DyeEX columns, and then eluents were recombined. The reaction (60 µL) was then adjusted with TdT reaction buffer (7.8 µL), CoCl2 (7.8 µL), ddUTP-biotin (1 µL, 1 mM) (Jena Bioscience #NU-1619-BIOX-S), filtered H2O (0.4 µL), and terminal transferase (1 µL) and incubated for 60 min at 37°C to terminate the dT-tails with a conjugated biotin moiety. After cleaning with four DyeEX columns and recombining eluents, the processed DNA was diluted in filtered H2O (500 µL) and fragmented by probe sonication as described above. Thereafter, DNA fragments were mixed with Hydrophilic Streptavidin Magnetic Beads (5 µL) (NEB #S1421S) and binding buffer (500 µL) (5 mM Tris–HCl, pH 7.5, 1 M NaCl, 0.5 mM EDTA) and incubated for 60 min on a Nutator at room temperature. The beads were subjected to a pull-down process using a magnetic separation rack (Dynal #MPC-S) and washed three times with binding buffer (100 µL). After discarding the supernatant, beads were resuspended in filtered H2O (20 µL). Finally, the enriched DNA that was captured on the beads was directly subjected to Illumina library preparation using a SMART ChIP-seq kit (Takara #634865) following the manufacturer’s protocol. The final PCR step was performed using Illumina primers provided in the kit, and 12 cycles were run for amplification. The PCR product, with unique sequencing barcodes, was submitted to the MIT BioMicro Center for next-generation sequencing. Paired-end sequencing (150-bp) was performed on an Illumina MiSeq instrument, using a Custom Read 2 sequencing primer that was supplied in the kit (poly(dA) followed by universal read 2 primer).

Data analysis
As illustrated in the workflow (Supplementary Fig. S4B), the data analysis started with trimming adapters. Adapters were removed using bbduk from BBtools v35.85 (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 v2.4.5 with the setting “sensitive”. The bamfiles were cleaned using SMARTcleaner v1.064 and split into two strands using samtools v1.19.265. The coverage of each strand was calculated separately using bedtools v2.30.066 with the genomecov -d -5 option. Then, custom 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 reads starting (5ʹ-end) at each position. The 13 nt sequences centered at positions with depth ≥ 1 were retrieved using samtools65. The consensus motifs were analyzed with incrementing depth, typically from 1 to 100 with a step of 10, using MEME v5.3.367 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 < 3 bps collapsed into the centermost consensus motif sites. Regions where starting position were mapped within 4 bp on opposite strands and within reverse complementary sequences were considered as double-stranded modifications.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this work, Yuan et al. discovered a novel type of PT modification system involving BrxP and BrxC, validated its function, identified the modification sites, and quantified the frequency. As the only known backbone modification in DNA, PT modifications are widespread in bacteria and not only have unique chemical properties but also play roles in host defense. In addition, PT modifications have been developed into new biotechniques such as nucleic acid detection and site-directed RNA editing. The authors and collaborators elucidated the chemical structure of PT modifications in previous studies, opening the window to a new epigenetic era. They subsequently uncovered the ssp system, a PT modification-based bacterial defense system, and the antioxidation mechanism. PT modifications can be specifically recognized by sulfur-binding domains, with some able to recognize all sequence patterns, even those not discovered in nature at that time. The discovery of the third type of PT modification system in this study significantly expands the range of PT consensus sequences, offering a perfect explanation for the recognition issue of PT modifications.\nThis work expands our understanding of epigenetic diversity and provides a foundation for epigenetics in gut microbiota and its association with health and disease.\nOverall it is very well written manuscript and the experiments appear to have been conducted properly. I would recommend acceptance. Here are some minor issues that I recommend the authors to address:\n1. Line 353, For the 122 bacterial strains where PT dinucleotides were not detected, the authors suggest that this could be due to the majority of strains lacking the essential dndCD, sspBCD, or brxPC gene clusters, and errors in the reference genomes cannot be ruled out. Is it possible that gut microbiota can achieve PT modification in a symbiotic environment but lose this ability when cultured independently?\n2. Figure 2B, in the figure legend, the authors claimed that an edge (lines connecting nodes) is drawn between two nodes that have a BLAST E-value cutoff of ≥10–100. To my knowledge, a lower E value indicates higher homology, so is there any mistake in the “cutoff of ≥10–100”?\n3. Line 276, “13,663”. Line 283, “humane” should be “human”", "reviewer_2": "The authors explore the distribution of phosphorothioation genes in a large number of bacteria present in the human gut microbiome, finding that 6.3% of these organisms possess the genes necessary to produce PT modification in their genomes. They demonstrate PT modification is indeed present in these strains and in total DNA isolated from human gut microbiome samples. They identify several putative new PT synthesis systems through gene neighborhood analyses and characterize PT modification in one of these, a Type 4 BREX system. They show PT modification requires only BrxP and BrxC, and occurs partially at relatively frequent 3+ basepair motifs which differ for the several BREX Type 4 systems investigated. The finding that brxP + brxC are the minimal set of genes required for PT modification in the BREX 4 system is important. This differs from other BREX Types, such as Type I, where brxA, brxB, brxC, brxX and brxZ are all required for host modification (methylation). This work is significant as a step toward identifying what role PT epigenetics may play in the human gut microbiome in health and disease.\nThe data presented generally supports the conclusions made. The analysis on the distribution of PT modifications relative to coding sequences is less convincing as it relies on examining 1 kb regions upstream and downstream of coding regions (CDS) rather than actual coding and non-coding sequences. Given the dense coding in bacteria, most 1kb regions adjacent to a CDS will themselves contain coding sequences for other gene(s), so a more precise analysis here seems warranted.\nThe methodology of identifying PT systems by identifying homologs of the core PT genes using BLAST then clustering PT systems based on neighborhood gene analysis is sound and is applied to an impressive number of human gut bacteria for which genome sequences are available. The methodology can be readily reproduced (as the set of genome sequences grows).\nThe potential for microbes containing PT to contribute to either human health or disease through their role in the human gut microbiome is worth investigating!\nSpecific Comments:\nLine 91-92: States “DndACDE and SspBCD are three and four gene operons, respectively” but then line 94 refers to a four gene DndACDE system while line 98 also refers to a four gene SspABCD system: just may be a bit confusing to the reader.\nLine 127: the genes of BREX type IV systems are brx PCZL, not brxRXZL (brxR is a regulator gene in BREX Type I systems while brxX is the DNA methyltransferase that is replaced by the brxP PT gene in Type IV systems.\nFigure 2: suggest: be consistent with BREX names: either use PglZ or BrxZ, not sometimes one or the other (they are the same).\nFor Methanothermobacter sp: Is the DndC gene a close homolog of the BrxP gene in the original BREX paper? If so, include the BrxP designation. Likewise, is the “NTPase” here a close homolog of the Type IV BrxC gene in the original BREX paper? If so, also call it BrxC as well as NTPase.\nFigure 2 and accompanying text: Suggest stating that phage defense systems often occur together in ‘islands’, so it is common to see multiple different defense systems in the same neighborhood even though they function independently of one another. Systems can be partial, due to decay or incomplete horizontal transfer. For example, the MTase gene in the Fischerella strain here is most likely independent of the BREX-Dnd system.\nFigure 2: it is a bit difficult to tell in which clusters the strains shown in the genome neighborhood box reside, at least from the colored circles in my printed copy. Is there a way to help the reader more easily see which cluster contains each strain example ?\nFigure 2: The system in Fischerella thermalis seems to represent a variant on the Brex Type IV architecture, in that here there are a full set of Dnd genes (BCDE), whereas the originally described BREX Type IV system only has a single Dnd/PT gene called BrxP (in Goldfarb 2014 BREX paper). This is an interesting new wrinkle in BREX systems (maybe something like Type 4A and Type 4B?)\nFigure 2B: The Chloroflexota example contains an intact Brex Type I system adjacent to the ssp genes. This is likely a DNA-methylation based BREX system that happens to reside next to ssp genes in this defense island, rather than a PT-based BREX system (based on the presence of the BrxX methyltransferase gene in the proper context.\nThe Cyanobacterium ATCC 29140 system here might be PT based (with SspB,C,D and second B gene) functioning as the (PT) modifying and DNA targeting role, although it could be a partial Brex Type I system that is unrelated but residing next to an Ssp system in a defense island context. This would need genetic validation (beyond what is needed for this manuscript).\nLine 332: Missing genes could potentially be located on a low copy plasmid that did not make it into the genome sequence assembly.\nLine 355ff: Is there any reason BREX systems would exclude GpsA, GpsC and GpsT from their recognition repertoire? Or is there just a small sample size of 3 characterized Type 4 BREX systems to date, so maybe these were simply not seen (yet)?\nLine 367-8: Why is the ApsC dinucleotide seen at so much higher frequency in the LC-MS data compared to the PT-seq data? What does this reveal about PT-seq motif discovery, if anything?\nLine 380: seems like it should be 1080 5’-GApsAG-3’/5’-CTpsTC-3’ and 667 5’-GApsCG-3’/5’-CGpsTC-3’ ??? Convention is 5’-3’ for both DNA strands.\nLine 393ff: uneven genomic distribution of PTs\nGiven that bacterial genomes are mostly coding sequences, and can even have slightly overlapping genes, how relevant is the “1kb upstream and downstream” from any given gene as an identifier of non-coding DNA, as these regions likely are also gene coding regions, just for other genes?\nThis analysis is not my area of expertise, but it seemed unconvincing to me.\nFor each given motif in the genome, say CpsCTC, are any locations never PT modified, while other locations are nearly always modified, or is there partial PT modification at every CCTC motif? Perhaps some comment on the level of modification at the various CCTC sites and whether any differences might be functional or are simply due to stochastic variation or perhaps variation in accessibility to the BrxP/BrxC modifying complex would be helpful.\nAgree that basis for biased distributions remains to be defined.\nLine 410: correct Fig. 5B to Supplementary Fig S5B\nLine 443 paragragh:\nIt seems only a very few BREX Type 4 systems have been characterized for recognition motifs (3 strains in Figure 5), so it seems premature to suggest BREX Type 4 do not recognize GpsN: just seems too little data to date. It will be interesting to identify the DNA binding pocket for BrxP(C) and SspBCD complexes to see how these specify their target sequence motifs.\n(Hopefully someone is working on CryoEM for the complex with DNA)\nPerhaps comment on the overall frequency of PT modification in host genome and how partial modification at frequent sites might be sufficient to identify self from non-self invasive DNAs that are targeted for restriction. It is interesting that 17% modification at CCTC means there is a host-protective PT modification on average every 1500 bp (1 four base motif every 256 bp, divided by .17 modification frequency = 1 PT mod every 1506 bp on average). This level of protective host modification is similar to Type I BREX methylation, since Type I BREX generally have a 6 bp recognition motif that is 100% methylated, or approximately 1 modification every 4000bp on average. Perhaps the PT-modification BREX systems just need to sense the presence of some PT modifications at frequent binding sites within a larger region of DNA rather than modification at nearly all of less frequent but highly-specific binding sites as in the BREX Type I systems?\nBREX defense systems may be well adapted to use distributed partial modification as a means of protecting the host from the effects of the defense system. Their requirement for forming a complex of multiple protein partners to accomplish their defensive role may make them able to tolerate partial PT modification without triggering restrictive complex formation when they encounter unmodified individual target site motifs in the host DNA, while still enabling restrictive complex formation upon encountering completely unmodified invasive DNAs.", "reviewer_3": "In the work \"Phosphorothioate DNA modification by BREX Type 4 systems in the human gut microbiome.\" Yuan et al. searched publicly available human microbiome metagenomes to find the genomic neighborhoods surrounding homologs of DndC and SspD, enzymes catalyzing the transfer of sulfur to the DNA backbone. In so doing they identified DndC and SspD homologs not only in the expected context of Dnd and Ssp gene clusters, but also within Type IV BREX systems, suggesting previously a previously unrecognized pathway of PT modification. Genetic knockouts within culturable representatives of the hits obtained in their search, as well as heterologous expression of genes in non-PT modifying hosts confirmed that the BrxC and BrxP of Type IV BREX systems, were sufficient to produce PT modifications of a variety of sequence contexts. Additionally, they performed PT-seq of select strains expressing Type IV BREX systems, as well as previously uncharacterized strains encoding Dnd and Ssp gene clusters, expanding the known repertoire of sequence motifs containing a PT backbone modification.\nThe work is significant in that it not only expands the known pathways of PT modification and the sequences targeted for sulfuration, but reveals an entirely new role of PT modification in the host-dependent modification of self DNA in a BREX system, functioning analogously methylation of dA at N6. I believe the evidence presented support the claims stated in the work, with proper interpretation of the results and sufficient methods for others to reproduce the work.\nI think there is sufficiently new ground broken here to demand additional experiments, and it is sufficient for publication pending minor revisions, but there are a few areas which, if addressed, would elevate the work. First, would be demonstration of anti-phage activity by the Type 4 BREX systems characterized in this work, i.e. biological function. Perhaps this function was already demonstrated by in earlier work or by other groups, and if it was cited, I might have missed it in the text. Second would be some preliminary characterization of the BrxP, even a large language model predicted structure, to account for its unexpected functions, i.e. sulfur transferase, target recognition, and protein-protein interactions. Thirdly, how does this work square with the recent indentification of a minimal PT-synthesizing enzyme encoded in bacterial and archaeal extremophiles: TdpC (see An T, Tan Q, Jiang L, Liu L, Jiang X, Liu L, Chang X, Tian X, Deng Z, Gao S, Wang L, Chen S. A DNA phosphorothioation pathway via adenylated intermediate modulates Tdp machinery. Nat Chem Biol. 2025 Aug;21(8):1160-1170. doi: 10.1038/s41589-024-01832-w. Epub 2025 Jan 16. PMID: 39820821.) Are there any similarities with this system and the systems uncovered in this work? How might what is known about TdpC affect hypotheses and results for the minimal sets of genes needed for PT-modification revealed in this work?\nSome other suggested edits:\n• lines 91 and 93- is the four gene cluster dndBCDE or is it DndACDE?\n• line 106- would \"altered\" be more appropriate than \"unique\"\n• Figure 2- scale bars might be helpful so as to compare homologs between different systems, some of which, though they have the same gene name, might have additional domains revealed by the different lengths.\n• is Type 4 capitalized or is it simply type 4? e.g. line 170 versus line 204\n• line 235 limit should be limited\n• line 283 human instead of humane\n• Figure 5 labels of Group 1, Group 2, Group 3, Group 4 is a bit confusing. What was the criterion by which they were grouped? Are they clusters?\n• lines 389-391 is not a complete sentence. Minimally, represent could be changed to represents"}, "author_response": "Response to reviewer comments for NCOMMS-25-42857-T “Phosphorothioate DNA modification by BREX type 4 systems in the human gut microbiome” by Yuan et al.\nWe thank the reviewers for their constructive and thoughtful comments. We have revised the manuscript to address all the concerns, and the resulting work is much stronger. The reviewer comments are noted in black font, and our responses are noted in blue font. Revised text in the manuscript and supplementary information is noted in red font.\nReviewer #1\nIn this work, Yuan et al. discovered a novel type of PT modification system involving BrxP and BrxC, validated its function, identified the modification sites, and quantified the frequency. As the only known backbone modification in DNA, PT modifications are widespread in bacteria and not only have unique chemical properties but also play roles in host defense. In addition, PT modifications have been developed into new biotechniques such as nucleic acid detection and site-directed RNA editing. The authors and collaborators elucidated the chemical structure of PT modifications in previous studies, opening the window to a new epigenetic era. They subsequently uncovered the ssp system, a PT modification-based bacterial defense system, and the antioxidation mechanism. PT modifications can be specifically recognized by sulfur-binding domains, with some able to recognize all sequence patterns, even those not discovered in nature at that time. The discovery of the third type of PT modification system in this study significantly expands the range of PT consensus sequences, offering a perfect explanation for the recognition issue of PT modifications.\nThis work expands our understanding of epigenetic diversity and provides a foundation for epigenetics in gut microbiota and its association with health and disease.\nOverall it is very well written manuscript and the experiments appear to have been conducted properly. I would recommend acceptance. Here are some minor issues that I recommend the authors to address:\n1. Line 353, For the 122 bacterial strains where PT dinucleotides were not detected, the authors suggest that this could be due to the majority of strains lacking the essential dndCD, sspBCD, or brxPC gene clusters, and errors in the reference genomes cannot be ruled out. Is it possible that gut microbiota can achieve PT modification in a symbiotic environment but lose this ability when cultured independently?\nThe reviewer makes a good point here. While some bacterial symbionts have lost genes for enzymes catalyzing specific steps in multistep RNA modifications (e.g., losing the ability to make the queuosine precursor queuine), the same is unlikely to be true for phosphorothioate DNA modifications. The PT modification pathway does not require unique small-molecule intermediates and the modification genes are generally not interchangeable among different bacteria, so it is unlikely that gut microbes can lose one or more PT modification genes and exist as a “PT symbiont” in the gut community. We thus believe that the 122 strains lacking detectable PT dinucleotides most likely lack PT synthesis genes or have errors in their reference genomes. We have addressed this possibility in the revised manuscript on lines 370-371.\n2. Figure 2B, in the figure legend, the authors claimed that an edge (lines connecting nodes) is drawn between two nodes that have a BLAST E-value cutoff of ≥10–100. To my knowledge, a lower E value indicates higher homology, so is there any mistake in the “cutoff of ≥10–100”?\nWe thank the reviewer for pointing out this typographical error. We have corrected it in the revised manuscript on lines 761-762 in the legend for Figure 2B.\n3. Line 276, “13,663”. Line 283, “humane” should be “human”\nWe thank the reviewer for pointing out these errors, which have been corrected in the revised manuscript.\nReviewer #2\nThe authors explore the distribution of phosphorothioation genes in a large number of bacteria present in the human gut microbiome, finding that 6.3% of these organisms possess the genes necessary to produce PT modification in their genomes. They demonstrate PT modification is indeed present in these strains and in total DNA isolated from human gut microbiome samples. They identify several putative new PT synthesis systems through gene neighborhood analyses and characterize PT modification in one of these, a Type 4 BREX system. They show PT modification requires only BrxP and BrxC, and occurs partially at relatively frequent 3+ basepair motifs which differ for the several BREX Type 4 systems investigated. The finding that brxP + brxC are the minimal set of genes required for PT modification in the BREX 4 system is important. This differs from other BREX Types, such as Type I, where brxA, brxB, brxC, brxX and brxZ are all required for host modification (methylation). This work is significant as a step toward identifying what role PT epigenetics may play in the human gut microbiome in health and disease.\nThe data presented generally supports the conclusions made. The analysis on the distribution of PT modifications relative to coding sequences is less convincing as it relies on examining 1 kb regions upstream and downstream of coding regions (CDS) rather than actual coding and non-coding sequences. Given the dense coding in bacteria, most 1kb regions adjacent to a CDS will themselves contain coding sequences for other gene(s), so a more precise analysis here seems warranted.\nThe methodology of identifying PT systems by identifying homologs of the core PT genes using BLAST then clustering PT systems based on neighborhood gene analysis is sound and is applied to an impressive number of human gut bacteria for which genome sequences are available. The methodology can be readily reproduced (as the set of genome sequences grows).\nThe potential for microbes containing PT to contribute to either human health or disease through their role in the human gut microbiome is worth investigating!\nSpecific Comments:\nLine 91-92: States “DndACDE and SspBCD are three and four gene operons, respectively” but then line 94 refers to a four gene DndACDE system while line 98 also refers to a four gene SspABCD system: just may be a bit confusing to the reader.\nWe apologize for the confusion here. We have revised the manuscript to clarify structure of the dnd and ssp modification systems (lines 92-93): “Typically, the modification components are organized as a four-gene operon for dndBCDE and a three-gene operon for sspBCD, with an additional modification gene for dndA or sspA located adjacent to the operon, dispersed elsewhere in the genome, or replaced by a gene for a homolog such as iscS.”\nLine 127: the genes of BREX type IV systems are brx PCZL, not brxRXZL (brxR is a regulator gene in BREX Type I systems while brxX is the DNA methyltransferase that is replaced by the brxP PT gene in Type IV systems.\nWe apologize for this error, which has been corrected in the revised manuscript on line 132.\nFigure 2: suggest: be consistent with BREX names: either use PglZ or BrxZ, not sometimes one or the other (they are the same).\nWe thank the reviewer for pointing out the inconsistency. We have revised the manuscript as suggested in Figure 2 and lines 166 to 167.\nFor Methanothermobacter sp: Is the DndC gene a close homolog of the BrxP gene in the original BREX paper? If so, include the BrxP designation. Likewise, is the “NTPase” here a close homolog of the Type IV BrxC gene in the original BREX paper? If so, also call it BrxC as well as NTPase.\nWe appreciate the reviewer’s suggestions. The homology between BrxP and DndC proteins had been mentioned and BrxC had been referred to as a DUF499 domain-containing ATP-binding protein in the original BREX paper. The manuscript has been revised as suggested in line 169.\nFigure 2 and accompanying text: Suggest stating that phage defense systems often occur together in ‘islands’, so it is common to see multiple different defense systems in the same neighborhood even though they function independently of one another. Systems can be partial, due to decay or incomplete horizontal transfer. For example, the MTase gene in the Fischerella strain here is most likely independent of the BREX-Dnd system.\nWe appreciate the suggestion. The manuscript has been revised as suggested in lines 176 to 179.\nFigure 2: it is a bit difficult to tell in which clusters the strains shown in the genome neighborhood box reside, at least from the colored circles in my printed copy. Is there a way to help the reader more easily see which cluster contains each strain example?\nWe appreciate the reviewer’s suggestion. We revised Figure 2 by labeling the nodes corresponding to the various illustrative bacterial strains with a two-letter label. We have revised the legend to Figure 2 to more clearly describe the node color schemes and neighborhood examples.\nFigure 2: The system in Fischerella thermalis seems to represent a variant on the Brex Type IV architecture, in that here there are a full set of Dnd genes (BCDE), whereas the originally described BREX Type IV system only has a single Dnd/PT gene called BrxP (in Goldfarb 2014 BREX paper). This is an interesting new wrinkle in BREX systems (maybe something like Type 4A and Type 4B?)\nWe agree that the PT-catalyzing BREX systems that we observed represent a variant on the BREX Type IV system, which raises interesting phylogenetic questions. For example, we found that 7 strains harboring the DndBCD-BrxCZL type gene islands and 7 out of 40 strains harboring the SspBCDB-BrxCZL type (or similar) gene islands are Cyanobacteria and 3 of them are Deinococcota. Some phylogenetic analyses suggest that Deinococcota and Cyanobacteria are both relatively close to the root of the bacterial tree of life as early-diverging lineages (PMID: 31792218), which agrees with the argument that dnd systems originated in ancient Cyanobacteria after the Great Oxygenation Event (PMID: 34737280). Similarly, BREX systems have also been found in branches near the root of the bacterial tree and deep branches (PMID: 25452498). It will require a comprehensive investigation, both bioinformatic and experimental, to see if we can define BREX Type 4 variants and identify the relationships among dnd, ssp, BREX typ4 and other potential PT systems at the time of their origin. We have added this comment to the Discussion section in lines 511-517.\nFigure 2B: The Chloroflexota example contains an intact Brex Type I system adjacent to the ssp genes. This is likely a DNA-methylation based BREX system that happens to reside next to ssp genes in this defense island, rather than a PT-based BREX system (based on the presence of the BrxX methyltransferase gene in the proper context.\nThe Cyanobacterium ATCC 29140 system here might be PT based (with SspB,C,D and second B gene) functioning as the (PT) modifying and DNA targeting role, although it could be a partial Brex Type I system that is unrelated but residing next to an Ssp system in a defense island context. This would need genetic validation (beyond what is needed for this manuscript).\nWe agree with the reviewer. The Chloroflexota example looks to result from gene loss in an adjacent BREX type 1-Ssp system. The genomic organization of the noted Cyanobacterium appears to be an Ssp system next to an incomplete BREX type 1 system. The sequence of the “BrxC-like” protein clustered with sspB and brxLZ in the Cyanobacterium example shares more similarity with SspC than BrxC from a typical BREX type 1 system. A sequence comparison and phylogenetic analysis of BrxC, the “BrxC-like” protein, and SspC sequences may determine whether the “BrxC-like” protein is a variation of BrxC (type 1) or a variation of SspC. We agree that this analysis is beyond the scope of the present manuscript, but we will pursue this in future studies. Based on these arguments, we renamed the “BrxC-like” protein as “BrxC-like ATPase” in Figure 2.\nLine 332: Missing genes could potentially be located on a low copy plasmid that did not make it into the genome sequence assembly.\nWe agree with the reviewer about this alternative explanation.\nLine 355: Is there any reason BREX systems would exclude GpsA, GpsC and GpsT from their recognition repertoire? Or is there just a small sample size of 3 characterized Type 4 BREX systems to date, so maybe these were simply not seen (yet)?\nWe agree with the reviewer that the limited sample size (a couple of hundred isolates) could have missed dnd, ssp or BREX systems that confer other dinucleotide patterns. We have added a note on lines 370-371 about this point. Our observations were based on our previously published correlation between PT consensus sequences identified by metagenomic PT sequencing and LC-MS analysis of PT dinucleotides in human fecal DNA (PMID: 40128848). It might not be surprising that the BREX-based PT modifications have different motifs than dnd and ssp systems given the observed differences between dnd and ssp themselves. However, further broad analyses may reveal variation of the PT motifs for all three systems. It will be interesting to identify the basis for the PT sequence selectivity for all PT modification systems.\nLine 367-8: Why is the ApsC dinucleotide seen at so much higher frequency in the LC-MS data compared to the PT-seq data? What does this reveal about PT-seq motif discovery, if anything?\nThe reviewer raises a good point that we did not address in the manuscript. We had previously observed that the efficiency of iodine-induced cleavage of PT sites to strand breaks in the PT-seq method varies by 2-fold for different dinucleotides (PMID: 38895297) and that the extension efficiency of the terminal transferase used for T-tailing can strongly depends on the nucleobase sequence of the initiators, where A-ending initiators are generally less favored (PMID: 34156829, 39988321, and 30367148). Together, these could account for the discrepancy between the LC-MS and PT-seq data for ApsCTC. However, another ApsC context, GAPSCG, showed strong agreement between PT-seq and LC-MS, so there is likely something unique about ApsCTC as a PT-seq target. Further, we have not observed such large discrepancies for other sequence contexts but future experience will inform about the need to revise the PT-seq method to obviate biases. We have added a comment about this in the Results section in lines 391-398.\nLine 380: seems like it should be 1080 5’-GApsAG-3’/5’-CTpsTC-3’ and 667 5’-GApsCG-3’/5’-CGpsTC-3’ ??? Convention is 5’-3’ for both DNA strands.\nThe typographical error has been corrected in line 402.\nLine 393: uneven genomic distribution of PTs\nGiven that bacterial genomes are mostly coding sequences, and can even have slightly overlapping genes, how relevant is the “1kb upstream and downstream” from any given gene as an identifier of non-coding DNA, as these regions likely are also gene coding regions, just for other genes?\nThis analysis is not my area of expertise, but it seemed unconvincing to me.\nWe appreciate the reviewer concerns here. We point out that the 1 kb upstream frame ensured the coverage of regulatory elements, but that only 20 nucleotides in the boundary regions proved to be strongly biased. Further, excluding polycistronic operons did not affect the observed biases. We are therefore confident on the conclusions based on the 1 kb frame.\nFor each given motif in the genome, say CpsCTC, are any locations never PT modified, while other locations are nearly always modified, or is there partial PT modification at every CCTC motif? Perhaps some comment on the level of modification at the various CCTC sites and whether any differences might be functional or are simply due to stochastic variation or perhaps variation in accessibility to the BrxP/BrxC modifying complex would be helpful.\nAgree that basis for biased distributions remains to be defined.\nWe appreciate the reviewer’s astute observation here, which was addressed in the third paragraph of the Discussion in the original manuscript. Partial modification of all available sequence motifs, hemi-modification within a sequence motif, cell replication-dependent shifts in the locations of modified motifs, and active maintenance of a constant density of PT modifications have all been observed in PT systems identified to date (PMID: 28604692, 30289455, 24899568, 32518115), despite the presence of active restriction components (PMID: 24899568). While some portion of the hemi-modification and partial modification could result from naturally occurring oxidation and repair or replacement, the unusual behavior of all PT systems to date appears to be inherent to the mechanisms of modification target selection and surveillance for restriction. This manuscript has been revised to emphasize these issues in lines 474-481.\nLine 410: correct Fig. 5B to Supplementary Fig S5B\nThis error has been revised in line 433.\nLine 443 paragraph:\nIt seems only a very few BREX Type 4 systems have been characterized for recognition motifs (3 strains in Figure 5), so it seems premature to suggest BREX Type 4 do not recognize GpsN: just seems too little data to date. It will be interesting to identify the DNA binding pocket for BrxP(C) and SspBCD complexes to see how these specify their target sequence motifs.\n(Hopefully someone is working on CryoEM for the complex with DNA)\nWe agree with the reviewer about drawing firm conclusions about BREX-based PT modification sequence motifs based on the limited data in the present studies. We revised the manuscript to address the reviewer’s concern in lines 370-371 and 499-500.\nPerhaps comment on the overall frequency of PT modification in host genome and how partial modification at distributed frequent sites might be sufficient to identify self from non-self invasive DNAs that are targeted for restriction. It is interesting that 17% modification at CCTC means there is a host-protective PT modification on average every 1500 bp (1 four base motif every 256 bp, divided by .17 modification frequency = 1 PT mod every 1506 bp on average). This level of protective host modification is similar to Type I BREX methylation, since Type I BREX generally have a 6 bp recognition motif that is 100% methylated, or approximately 1 modification every 4000bp on average. Perhaps the PT-modification BREX systems just need to sense the presence of some PT modifications at frequent binding sites within a larger region of DNA rather than modification at nearly all of less frequent but highly-specific binding sites as in the BREX Type I systems?\nBREX defense systems may be well adapted to use distributed partial modification as a means of protecting the host from the effects of the defense system. Their requirement for forming a complex of multiple protein partners to accomplish their defensive role may make them able to tolerate partial PT modification without triggering restrictive complex formation when they encounter unmodified individual target site motifs in the host DNA, while still enabling restrictive complex formation upon encountering completely unmodified invasive DNAs.\nWe appreciate the reviewer’s comment here, which parallels the previous comment about partial modification. We have revised the manuscript in lines 474-481 to discuss this issue. The basis for this unusual modification and restriction system will be the subject of major future studies.\nReviewer #3\nIn the work \"Phosphorothioate DNA modification by BREX Type 4 systems in the human gut microbiome.\" Yuan et al. searched publicly available human microbiome metagenomes to find the genomic neighborhoods surrounding homologs of DndC and SspD, enzymes catalyzing the transfer of sulfur to the DNA backbone. In so doing they identified DndC and SspD homologs not only in the expected context of Dnd and Ssp gene clusters, but also within Type IV BREX systems, suggesting previously a previously unrecognized pathway of PT modification. Genetic knockouts within culturable representatives of the hits obtained in their search, as well as heterologous expression of genes in non-PT modifying hosts confirmed that the BrxC and BrxP of Type IV BREX systems, were sufficient to produce PT modifications of a variety of sequence contexts. Additionally, they performed PT-seq of select strains expressing Type IV BREX systems, as well as previously uncharacterized strains encoding Dnd and Ssp gene clusters, expanding the known repertoire of sequence motifs containing a PT backbone modification.\nThe work is significant in that it not only expands the known pathways of PT modification and the sequences targeted for sulfuration, but reveals an entirely new role of PT modification in the host-dependent modification of self DNA in a BREX system, functioning analogously methylation of dA at N6. I believe the evidence presented support the claims stated in the work, with proper interpretation of the results and sufficient methods for others to reproduce the work.\nI think there is sufficiently new ground broken here to demand additional experiments, and it is sufficient for publication pending minor revisions, but there are a few areas which, if addressed, would elevate the work.\n• First, would be demonstration of anti-phage activity by the Type 4 BREX systems characterized in this work, i.e. biological function. Perhaps this function was already demonstrated by in earlier work or by other groups, and if it was cited, I might have missed it in the text.\nWe are also interested in the question of anti-phage activity of the BREX type 4 systems. We cited the paper describing this activity but did not specifically mention the anti-phage activity. However, we do not know if the variant type 4 systems that we observed also have anti-phage activity, which will be the subject of future studies. We have revised the Discussion in lines 454-456 to note the issue of anti-phage activity.\n• Second would be some preliminary characterization of the BrxP, even a large language model predicted structure, to account for its unexpected functions, i.e. sulfur transferase, target recognition, and protein-protein interactions.\nWe share the reviewer’s interest in defining the biochemical function of BrxP. A full biochemical characterization is beyond the scope of the present manuscript and will be the subject of future publications.\n• Thirdly, how does this work square with the recent indentification of a minimal PT-synthesizing enzyme encoded in bacterial and archaeal extremophiles: TdpC (see An T, Tan Q, Jiang L, Liu L, Jiang X, Liu L, Chang X, Tian X, Deng Z, Gao S, Wang L, Chen S. A DNA phosphorothioation pathway via adenylated intermediate modulates Tdp machinery. Nat Chem Biol. 2025 Aug;21(8):1160-1170. doi: 10.1038/s41589-024-01832-w. Epub 2025 Jan 16. PMID: 39820821.) Are there any similarities with this system and the systems uncovered in this work? How might what is known about TdpC affect hypotheses and results for the minimal sets of genes needed for PT-modification revealed in this work?\nWe agree with the reviewer that this recent publication is an exciting and relevant discovery, which we had not noted in the original manuscript. While we revisiting all the analyses performed here with the new TdpABC gene information would be beyond the scope of the present studies, the new information expands the realm of PT epigenetics. We have revised the Introduction in lines 91 and 95-96 and the Discussion in lines 453 and 487 to note the tdpABC gene system.\n• lines 91 and 93- is the four gene cluster dndBCDE or is it DndACDE?\nWe have corrected this to dndBCDE in the revised manuscript in lines 92.\n• line 106- would \"altered\" be more appropriate than \"unique\".\nWe appreciate the suggestion. The manuscript has been revised in line 111.\n• Figure 2- scale bars might be helpful so as to compare homologs between different systems, some of which, though they have the same gene name, might have additional domains revealed by the different lengths.\nWe agree with the reviewer and have revised Figure 2 with the change.\n• is Type 4 capitalized or is it simply type 4? e.g. line 170 versus line 204\nThe manuscript has been revised with ‘type 4’ in lowercase consistently throughout.\n• line 235 limit should be limited.\nThe manuscript has been revised as suggested.\n• line 283 human instead of humane\nThe manuscript has been revised as suggested.\n• Figure 5 labels of Group 1, Group 2, Group 3, Group 4 is a bit confusing. What was the criterion by which they were grouped? Are they clusters?\nThe isolates are grouped by both PT dinucleotides detected and their PT genes. Group 1, GpsA – dnd system; Group 2, CpsC – ssp system; Group 3, PT dinucleotide other than GpsA; Group 4, no PT dinucleotide detected. The genes for each PT system are gene clusters.\n• lines 389-391 is not a complete sentence. Minimally, represent could be changed to represents\nThe manuscript has been revised as suggested in line 412-414.\nReviewer #3 (Remarks on code availability):\nI'm not very proficient in shell scripts (.sh) or Python scripts (.py) but the repository is a collection of scripts dealing with processing the data obtained following Illumina sequencing of libraries prepared from input DNA treated by iodine resulting in cleavage at phosphorothioate linkages followed by 3' extension and biotin tailing for enrichments. The scripts were deposited for an earlier publication on the PT-seq method: Yuan Y, DeMott MS, Byrne SR; Global Microbiome Conservancy; Dedon PC. PT-seq: A method for metagenomic analysis of phosphorothioate epigenetics in complex microbial communities. bioRxiv [Preprint]. 2024 Jun 3:2024.06.03.597111. doi: 10.1101/2024.06.03.597111. PMID: 38895297; PMCID: PMC11185561.\nThe collection of scripts assumes a certain level of proficiency at the command line and does not provide information regarding package dependencies, installation, and execution. The scripts themselves, however, do contain comment lines that likely would be useful to an experienced bioinformatician.\nA closer reading of the Methods section also gives sufficient guidance, in my view.\nWe appreciate the reviewer’s comments. We updated the repository by (1) adding more comments and explanations of usage; (2) adding software versions; and (3) more details in README file."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have very successfully addressed the questions raised in the initial reviews and have made appropriate modifications (pun intended) to the manuscript. I can enthusiastically recommend publication of this study. Well done.", "reviewer_3": "Thank you for addressing my concerns. The work is, in my view, suitable for publication without any further changes. Congrats and looking forward to sharing the paper with colleagues when it is published."}, "author_response": ""}]

Decision letter