A set of Bacteroidales genomes was obtained from Evans et al. (2022) and completed with genomes of publicly available bacterial strains in BEIresources.¹ Genome assemblies were downloaded from NCBI using RefSeq accession numbers as of July 2022. RefSeq entries exhibiting anomalous assembly, potential contamination, or a complete absence of the genome were removed, resulting in a total of 1,136 assemblies that were used for the in silico screening. Supplementary Table S1 provides detailed information about the analyzed strains, including taxonomy, isolation source, and accession numbers.

To screen the selected Bacteroidales genomes for the presence of putative RiPP biosynthetic gene clusters (BGCs), the BAGEL4 web-server was used (van Heel et al., 2018).² A list of the protein families or domains detected and the rules used to identify them during the screening can be found in Supplementary Table S2.

To identify promising candidate RiPP BGCs, manual curation was performed over the areas of interest (AOIs) detected by BAGEL4. The selection process involved considering the presence and arrangement of post-translational modification (PTM) enzymes, putative precursor peptides, and transporters within the AOI, as well as the occurrence of the AOI across genomes and strain isolation sources. The selected AOIs, which usually comprise over 20kbs sequence including the detected features, were manually investigated using the InterPro web-based tool (Mitchell et al., 2019)³ and protein–protein BLAST in the NCBI web service.

The BiG-MAP bioinformatic tool (Pascal Andreu et al., 2021) was employed to explore the prevalence and abundance of selected RiPP BGCs in metagenomic samples. This tool maps shotgun sequencing reads onto BGCs predicted by antiSMASH (Blin et al., 2023). For this purpose, genomes containing selected RiPP BGCs were used as an input on antiSMASH7 in loose detection strictness, which successfully detected 6 out of the subset of 9 BGCs. The resulting antiSMASH files were used as an appropriate input for BiG-MAP. Metagenomic samples downloaded from the Human Microbiome Project (HMP⁴; RRID:SCR_004919) were used to map BGCs prevalence and abundance. A total of 808 metagenomic samples from 4 different sources (205 dorsum tongue, 201 feces, 200 gingiva, and 202 posterior fornix of vagina) were randomly selected. RiPP BGC prevalence was calculated as the ratio of HMP samples with one or more units of target BGC to the total number of HMP samples for each body source. To quantify their abundance, a RPKM (Read Per Kilobase per Million mapped reads) metric was used.

The protein sequence of BGC8 YcaO (UniProt ID A0A0A2F9V0) from P. gulae COT-052 OH1451 was used to generate a Sequence Similarity Network (SSN) with the online Enzyme Function Initiative-Enzyme Similarity Tool (EFI-EST; Gerlt et al., 2015)⁵ using an E value of 5 and an alignment score (AS) of 100. The SSN for the radical S-adenosyl methionine (rSAM) enzyme (UniProt ID A0A0A2F578) was downloaded from radicalSAM.org (Oberg et al., 2022) at an AS of 45. Both SSNs were visualized using Cytoscape (Shannon et al., 2003).

The amino acid sequence of PTM enzymes used for multiple sequence analysis (MSA) was obtained using their UniProt IDs (Supplementary Table S3) and aligned in MEGA11 (Tamura et al., 2021) using MUSCLE (Edgar, 2004) on default parameters. Jalview 2.11.2.6 (Waterhouse et al., 2009) was used for MSA visualization and Logo generation.

The AlphaFold predictions (Jumper et al., 2021) for the protein 3D structures were downloaded using the enzyme’s UniProt IDs (Supplementary Table S3). When available, the X-ray structure was downloaded instead. UCSF Chimera (Pettersen et al., 2004) was used for visualization of the 3D structure and MatchMaker tool (Meng et al., 2006) for structural comparisons between proteins.

First, the 1,136 selected Bacteroidales genomes were scanned for the presence of RiPP biosynthetic gene clusters (BGCs) using BAGEL4. This analysis yielded a total of 1,340 areas of interest (AOI) related to RiPP biosynthesis (Table 1). It should be emphasized that BAGEL4 designates an AOI upon detecting homology to a specific bacteriocin structural gene and/or certain protein motifs associated with bacteriocin BGCs (van Heel et al., 2018). Therefore, AOI hits should be considered an indication of minimal potential RiPP production, rather than the actual presence of a BGC. Detailed information about the list of post-translational modification (PTM) enzyme families or domains detected and the rules used to identify them can be found in Supplementary Table S2.

The most frequently identified AOIs in this analysis contain one or more radical S-adenosylmethionine (rSAM) genes, either alone or in combination with other genes commonly involved in RiPP BGCs (Figure 1; Table 1). AOIs containing rSAM as the sole biosynthetic enzyme are widely distributed among the different Bacteroidales families. In fact, the only families that did not exhibit rSAM containing AOIs were those underrepresented in the analysis, such as Tenuifilaceae or Paludibacteraceae. The rSAM superfamily represents one of the largest and most functionally-diverse enzyme superfamilies (Oberg et al., 2022). Through radical-mediated reactions, rSAM enzymes are able to perform complex chemical transformations and rearrangements, dehydrogenations, methylations, sulfur insertions, and modifications of DNA, RNA and peptides (Sofia et al., 2001; Shisler and Broderick, 2012). Enzymes of this superfamily have been shown to participate in the biosynthesis of various classes of RiPPs, carrying a wide array of reactions, including epimerization, thioether crosslinking, and methyltransference activity (Montalbán-López et al., 2021).

Following single rSAM AOIs, genetic regions containing both rSAM and LanC genes are particularly abundant. LanC enzymes are part of class I lanthipeptides BGCs and they install (methyl)lanthionine rings between the C_β atom of dehydrated serine or threonine residues and the sulfhydryl group of cysteine residues (Lubelski et al., 2008; Repka et al., 2017). These hybrid rSAM + LanC AOIs are notably frequent in strains of Bacteroides fragilis, Bacteroides thetaiotaomicron, Parabacteroides distasonis, and Tannerella forsythia (Supplementary Table S4). In some B. fragilis and P. distasonis strains, a combination of rSAM, LanC, and a gene containing a peptide S-glycosyltransferase (GlyS) domain is detected. rSAM genes can also be found in other combinations, including rSAM + YcaO. YcaO proteins have been shown to participate in several RiPP BGCs via azoline, thioamide, or amidine formation (Burkhart et al., 2017; Clark and Seyedsayamdost, 2022). Areas exhibiting rSAM + YcaO are particularly abundant in strains of the genus Bacteroides and Porphyromonas.

Stand-alone LanC cyclase genes were detected in several strains from different Bacteroidales families, with special prevalence in B. fragilis and T. forsythia strains. Consecutive action of LanB and LanC is required for lanthionine ring installation in class I lanthipeptides (Arnison et al., 2013) and, intriguingly, none of the screened genomes contained LanB dehydratase genes. In contrast, two AOIs exhibiting LanM-like genes were detected in Prevotella strains, either stand-alone in P. buccalis DNF00985 or in combination with a LanC-like gene in P. multisaccharivorax DSM 17128 (Figure 1; Supplementary Table S4). LanM enzymes participate in the biosynthesis of class II lanthipeptides, carrying out both dehydration and cyclization steps (Arnison et al., 2013).

Remarkably, experimentally described RiPPs were mostly absent from the in silico screening, with the exception of ComX. Five strains of B. fragilis were found to possess a protein exhibiting homology to the structural gene of the Bacillus amyloliquefaciens JRS8 quorum sensing pheromone ComX. Further exploration revealed that while the translated gene showed partial homology to the pheromone precursor, the key tryptophan residue was absent in the B. fragilis gene (data not shown). This tryptophan is a fully conserved residue whose post-translational modification results in the formation of a tricyclic structure (Okada et al., 2005). In light of this consideration and due to the unusual large size of the B. fragilis gene, ComX matches were deemed as false positives and, hence, not pursued any further.

Since the BAGEL4 analysis gave a high frequency of AOIs containing orphan enzymes (data not shown), a manual curation process of the AOIs was employed. This second round of selection focused on the presence and arrangement of PTM enzymes, putative precursor peptides, and transporters within the AOI. Additionally, the occurrence of the gene cluster across different genomes and the strain isolation source were given consideration. As a result, 9 putative interesting RiPP BGCs were selected and functionally annotated (Figure 2A). Due to their already reported abundance, it is unsurprising that 8 out 9 of the BGCs contain at least a rSAM enzyme. The PTM composition varies among the set of BGCs, as some clusters feature only rSAM enzymes (BGC1,2), whereas others exhibit a combination of rSAM and other types of PTM enzymes such as LanC (BGC3,4,5) or YcaO (BGC6,7,8). Only the last selected BGC is characterized by the absence of a rSAM and the presence of a LanM enzyme (BGC9), an enzyme involved in class II lanthipeptide biosynthesis and mostly found in Gram-positive bacteria. Strikingly, BGC3,4,5 contain LanC cyclases but lack the LanB encoded proteins, necessary for class I lanthipeptide biosynthesis. LanC cyclases are known to show highly conserved sequence features, namely a triad of residues involved in Zinc ion coordination (C-gap-CH) and additional single residues involved in active site formation (Li et al., 2006; Marsh et al., 2010). However, the putative LanC sequences identified in BGC3,4,5 do not show these motifs (Supplementary Figure S1).

Additional accessory enzymes were detected in several BGCs, which have the potential to introduce PTMs. These enzymes include pseudo rSAMs (BGC3,4,5), glycosyltransferases (BGC3,9), aminotransferases (BGC7,8) and a methyltransferase (BGC9). Although no specific function has been experimentally assigned to pseudo rSAMs, InterPro reports that members of this family (IPR026418) co-occur with rSAM enzymes and a precursor peptide, as it is observed in BGC3,4,5 (Figure 2A). It has been hypothesized that these pseudo rSAM are likely working in partnership with the rSAM to introduce chemical modifications in the putative precursor peptide.

Regarding export systems, all BGCs with the exception of BGC3 and BGC5 encompass at least one gene encoding a C39 peptidase-containing ATP-binding transporter (PCAT). Transport by PCAT containing systems is a common strategy for leader peptide cleavage and export of the mature compound out of the cell (Havarstein et al., 1995; Montalbán-López et al., 2021). Most of the precursor peptides annotated exhibit residues frequently present on a GG motif (TIGR01847; Figure 2B), a signature for peptides subjected to PCAT leader cleavage and export (Arnison et al., 2013; Bobeica et al., 2019). Notably, precursor peptides associated to BGC3 and BGC5 contain several residues matching the GG motif signature despite lacking a PCAT in their BGCs.

Other Gram-negative bacteria have been shown to employ PCAT in the translocation of antimicrobial peptides. A well-known example is the class II bacteriocin colicin V (ColV) export system, which requires PCAT, a membrane fusion protein (MFP), as well as the outer membrane protein (OMP) TolC for crossing both inner and outer membranes (Gilson et al., 1990; Wu and Tai, 2004). The ColV structural gene is located within the same gene cluster as the PCAT and MFP genes, whereas TolC is encoded elsewhere (Gilson et al., 1990). In the analyzed BGCs, several genetic arrangements are found (Figure 2A). For instance, BGC4 encodes all necessary genes for peptide export in Gram-negatives, featuring PCAT, MFP, and OMP genes. BGC2 also possesses a complete set of transporter genes, although unlike BGC4, the OMP gene is not encoded in tandem with PCAT and MFP. Notably, BGC6,7,8 and 9 contain both PCAT and OMP genes, yet no MFP was identified. This raises the question whether the transport for their mature compounds is governed by a distinct mechanism or if the detection of an MFP remains to be accomplished. The remaining BGCs possess different lay-outs, including only a PCAT in BGC1, several types of transporters in BGC3, or the total absence of them in BGC5.

Besides PTM and transporter genes, RiPP BGCs usually contain other complementary genes that are involved in self-immunity and regulation of expression. Before or during manual curation, no proteins displaying homology with already described immunity proteins were detected. Nonetheless, it has been demonstrated that Bacteroidetocins, unmodified antimicrobial peptides produced by members of the Bacteroidales order, can inhibit the growth of the producer strain (Coyne et al., 2019). This observation suggests that it may be possible to produce RiPPs with antimicrobial activity in the absence of immunity proteins. On the other hand, some efflux systems that are present in selected BGCs could be involved in self-immunity as reported for the antimicrobial peptides microcin J25 (MccJ25; Bountra et al., 2017) and microcin B17 (MccB17; Garrido et al., 1988), which are produced by other Gram-negative bacteria. Immunity to MccJ25 is achieved using the same ABC exporter, opening the possibility of a dual role for PCAT exporters in selected BGCs. MccB17 and ColV exhibit dedicated transporters for immunity and, therefore, the different transporters in BGC3 could be involved in immunity.

In contrast to a notable small detection rate of self-immunity genes, several potentially regulatory proteins were detected. This includes a helix-turn-helix (HTH) domain containing gene in BGC3 and a signal transduction system in BGC5, composed of a histidine kinase (HK) and a response regulator (RR). Both types of systems have been shown to be involved in transcriptional regulation (Gallegos et al., 1993; Kuipers et al., 1998; West and Stock, 2001). Additionally, tetratricopeptide-like helical (TRP) domain-encoding genes appear in BGC2 and 9. TRP domains are involved in protein–protein interaction and, among their many biological functions, have been suggested to play a role in gene expression regulation in other bacteria after pheromone binding (Capodagli et al., 2020).

BGC7 and 8 possess a gene containing an N-terminal aminotransferase domain (IPR004839) followed by a homeodomain (IPR009057). Aminotransferase involvement in RiPP biosynthesis has been previously described for methanobactocins (Mbns), copper-chelating peptides firstly identified in methanotrophic bacteria (Kim et al., 2004; Park et al., 2018). However, the homeodomain found in BGC7 and 8, a known DNA-binding motif, points toward their role as transcription factors. Interestingly, a similar combination of domains has been described for the Bacillus subtilis transcription factor GabR (Belitsky and Sonenshein, 2002), where the homeodomain is situated N-terminally and the aminotransferase C-terminally.

Finally, BGC5 shows an interesting cluster architecture. It presents two putative precursor peptides with a predicted long core peptide sequence, each of them followed by a putative signaling protein (Figures 2A,B). Although it appears to lack transporter sequences, BGC5 presents two rSAM, one pseudo rSAM and three LanC-like sequences, which might not perform the classical LanC activity, as previously discussed (Supplementary Figure S1). Finally, two sets of CAZymes complete the distinguishing features of BGC5. While it could be hypothesized that this cluster, and related clusters from B. fragilis species, could be a source of novel antimicrobials, no antimicrobial activity can be concluded from this analysis.

The microbiome has been found to be crucial for human health, and the production of antimicrobial peptides can have a significant impact on the diversity and function of the microbiome. To assess the potential significance of the RiPP BGCs detected in the human microbiome, the prevalence and abundance of those BGCs were studied on metagenomic samples from the Human Microbiome Project (HMP). For that purpose, a total of 808 metagenomic samples from 4 different body sources (feces, dorsum tongue, gingiva, and posterior fornix of vagina) were analyzed using BiG-MAP (Pascal Andreu et al., 2021). Out of the 9 selected BGCs, 6 were used as input. The presence of those BGCs in feces concurs with the isolation place of their producer strains (Figure 2A). As seen in Figure 3, BGC1, 2, 5 and 6 are present in almost all fecal samples. BGC7 and 8, only detected in Porphyromonas strains that were isolated from dog oral cavities, are surprisingly prevalent in human gingiva and tongue samples. As could be expected from their isolation source, the prevalence of the selected RiPP BGCs in vaginal samples is lower compared to other sources. Only BGC2, prevalent in every body source, appears in more than half of the vaginal samples analyzed.

Regarding their abundance, expressed in read per kilobase per million (RPKM), BGC8 was the most abundant cluster in both gingiva and on the tongue, while BGC2 showed the highest abundance in fecal samples. Despite being present just in around half of the vaginal samples, BGC2 is quite abundant in that body source, reaching 2.2 × 10⁴ RPKM. The opposite takes place in fecal samples with BGC1,5 and 6, whose abundance is low in spite of being quite prevalent on this body source.

Further characterization of the different components of BGC8 was undertaken due to its prevalence and abundance in oral metagenomic samples from human origin (Figure 3). The cluster is characterized by the presence of a YcaO and a rSAM (Figure 4A). Additionally, a small ORF of 61 residues is suspected to be the precursor peptide. The presence of an N-terminal bacteriocin type GG-signal sequence (TIGR01847) matches the presence of a PCAT export system (Arnison et al., 2013; Bobeica et al., 2019). Therefore, it is predicted that the leader peptide accounts for the first 27 residues.

In order to gain more insight into the YcaO enzyme function, its amino acid sequence was used to generate a sequence similarity network (SSN) with the online Enzyme Function Initiative-Enzyme Similarity tool (EFI-EST; Figure 4B). According to the SSN analysis, the YcaO enzyme from BGC8 clusters together with the YcaO from BGC6, which was experimentally characterized and named ItfB by Clark and Seyedsayamdost (2022) during the course of this study. ItfB was shown to be able to convert the backbone amide into an amidine between the last two C-terminal residues Thr and Phe, a new type of modification for the YcaO family (Clark and Seyedsayamdost, 2022). In spite of having low sequence similarity, the precursor peptides from both BGCs, A8 and A6 (ItfA), have as a final residue a Phenylalanine (Figure 4C). Further exploration of the similarity between the two YcaO enzymes showed that both enzymes share key residues for the ATP-binding motif typical of YcaO domains (IPR003776; Dunbar et al., 2014; Figure 4D). Although the sequences of these enzymes have a low percentage of identity, 36.57%, sequence-structure analysis revealed a similar protein structure for both enzymes (Figure 4E). These findings suggest that the YcaO enzyme from BGC8 may have a similar mode of action to ItfB, installing a backbone amidine on A8 between Gln33 and Phe34.

The rSAM enzyme of BGC8 contains a core rSAM domain (IPR007197), the unifying structural theme of the diverse rSAM superfamily (Grell et al., 2015). This domain adopts a partial triose-phosphate isomerase (TIM) barrel conformation (Nicolet and Drennan, 2004) and possesses a conserved CX₃CX₂C motif (Sofia et al., 2001). The TIM-barrel fold consists of a sixfold repeat of (βα) units, such that six parallel β-strands on the inside are covered by six α-helices on the outside. The CX₃CX₂C motif is located in a loop between strand β₁ and helix α₁ (Nicolet and Drennan, 2004). Mechanistically, the three cysteinyl sulfur groups of this motif coordinate three iron ions from a [4Fe-4S] cluster and the remaining iron anchors the SAM (Walsby et al., 2002). Additionally, an aromatic residue before the third cysteine of the CX₃CX₂C motif is suggested to be involved in SAM coordination via hydrophobic interactions (Grell et al., 2015). In Figures 5A,B, AlphaFold predicts that the rSAM enzyme of BGC8 adopts a partial TIM barrel conformation despite the lack of a complete β₆ strand. The CX₃CX₂C motif, present as “C₈₀TLNC₈₄KYC₈₇,” is located on the loop between strand β₁ and helix α₁. The Tyr in from of Cys87 is the aromatic residue suspected to be involved in SAM coordination.

Beyond their core domain, rSAM enzymes display a great diversity in structure, with functionalized N- and C-terminal extensions. N-terminal extensions involved in leader peptide recognition have been previously described for rSAM enzymes involved in peptide modification (Wieckowski et al., 2015; Grove et al., 2017). Named RiPP precursor recognition element (RRE), this N-terminal extension displays three α-helices and three β-sheets in a wHTH motif (Montalbán-López et al., 2021). The rSAM of BGC8 presents a slightly shorter extension with one α-helix and three β-sheets (Figure 5C), reminiscent of the canonical RRE described for other PTM enzymes. As a C-terminal domain, this enzyme exhibits a SPASM domain (IPR023885) that can bind two auxiliary [4Fe-4S] clusters via seven or eight cysteine residues (Haft and Basu, 2011). Enzymes containing this additional domain usually participate in post-translational modification of peptides (Grell et al., 2015). The role of the two auxiliary [4Fe-4S] clusters in RiPP modification has been suggested to be system dependent (Mendauletova et al., 2022), with an inter-Cys residue spacing varying among enzymes (Oberg et al., 2022). In Figures 5A,B, cysteine residues suspected to be involved in the binding of three [4Fe-4S] clusters- one involved in SAM coordination and two auxiliary ones- are highlighted in cyan, suggesting the role of BGC8 rSAM in peptide modification.

To elucidate whether the rSAM present in BGC8 and BGC6 (ItfD) could be catalyzing similar reactions, the new web-based resource RadicalSAM.org was used (Oberg et al., 2022). Using SSNs, RadicalSAM.org attempts to enable identification of isofunctional groups of rSAM enzymes. Exploration of the SSN generated for the rSAM superfamily showed that rSAM BGC8 and ItfD are both clustered together in Megacluster 1-1. SPASM domain containing enzymes cluster together in this Megacluster (Oberg et al., 2022), which matches the detection of such domain in rSAM BGC8 and ItfD. Nevertheless, a higher alignment score (AS) clusters the rSAM from BGC8 together with QhpD and the ranthipeptide maturases CteB and Tte1186 (Supplementary Figure S2A; Nakai et al., 2015; Bruender et al., 2016; Grove et al., 2017). Ranthipeptides subjected to CteB and Tte1186 maturases are RiPPs characterized by the presence of γ-thioether linkages between Cys and Thr residues, specifically Thr for CteA and Tte1186a. The presence of two Cys and Thr residues on the precursor peptide (Figure 4A) makes plausible that the rSAM from BGC8 performs such thioether linkages. In addition, multiple sequence alignment (Supplementary Figure S2B) and sequence-structure alignment (data not shown) shows a similar inter-Cys residue spacing between CteB and the rSAM BGC8. Notably, QhpD, CteB, and Tte1186 display 7 Cys residues on their SPASM domain, while the rSAM from BGC8 possesses a total of 8 (Supplementary Figure S2B).