DNA Phosphorothioate Modifications Are Widely Distributed in the Human Microbiome

The DNA phosphorothioate (PT) modification existing in many prokaryotes, including bacterial pathogens and commensals, confers multiple characteristics, including restricting gene transfer, influencing the global transcriptional response, and reducing fitness during exposure to chemical mediators of inflammation. While PT-containing bacteria have been investigated in a variety of environments, they have not been studied in the human microbiome. Here, we investigated the distribution of PT-harboring strains and verified their existence in the human microbiome. We found over 2000 PT gene-containing strains distributed in different body sites, especially in the gastrointestinal tract. PT-modifying genes are preferentially distributed within several genera, including Pseudomonas, Clostridioides, and Escherichia, with phylogenic diversities. We also assessed the PT modification patterns and found six new PT-linked dinucleotides (CpsG, CpsT, ApsG, TpsG, GpsC, ApsT) in human fecal DNA. To further investigate the PT in the human gut microbiome, we analyzed the abundance of PT-modifying genes and quantified the PT-linked dinucleotides in the fecal DNA. These results confirmed that human microbiome is a rich reservoir for PT-containing microbes and contains a wide variety of PT modification patterns.


Introduction
Phosphorothioate (PT) DNA modifications, in which the nonbridging oxygen in the phosphate backbone is replaced by sulfur, are widespread among bacteria and archaea [1]. The PT-modifying gene cluster, dndA-E, has been found in more than 1300 sequenced genomes [2] and confers cells with 5 -G ps AAC-3 /5 -G ps TTC-3 or 5 -G ps GCC-3 consensus sequences [3]. Recently, another PT-modifying gene cluster has been reported, sspA-D, which confers cells with 5 -C ps CA-3 on the single strand [4]. The DndA protein, which can be functionally substituted by an IscS (a cysteine desulfurase located elsewhere in the genome) [5,6], transfers sulfur into the Fe-S cluster of DndC [5,7]. SspA and SspD may share the same initial sulfur mobilization pathway with DndA and DndC, respectively [4]. Meanwhile, SspC shows ATPase activity and may provide energy in a manner similar to DndD [4]. Though PT modifications are usually introduced in a sequence-specific manner, only 12%-14% of target motifs are modified [3,8]. In some bacteria, PT modifications are involved in a restriction-modification (R-M) system with gene cassette dndF-H/sspE or an antiviral system with gene cassette pbeA-D [9]. However, nearly half of the PT-containing strains lack dndF-H/sspE and pbeA-D, indicating that the PT modifications possess other functions [4,10]. A recent transcriptomic and metabolomic analysis showed that PT modifications contribute to the cellular redox state [11], while, at the same time, PT modifications confer sensitivity to the HOCl produced by neutrophils [12]. To date, little is known about the presence of PT-containing bacteria in the human microbiome. Here, we performed an informatic analysis of sequenced microbiome-related genomes and a mass spectrometric analysis of fecal DNA to show the widespread presence of PT-modifying genes and PT-linked dinucleotides in the human microbiome.

Multiple Sequence Alignments
All publicly available genomes (including the complete and partial genomes) were downloaded from the online databases of NCBI (National Center for Biotechnology Information) (https://www. ncbi.nlm.nih.gov/genome/), EMBL (European Molecular Biology Laboratory) (https://www.ebi.ac.uk/), and HMP (Human Microbiome Project) (https://www.hmpdacc.org/). We used the MakeDB of multigeneblast (software version 1.1.13) [13] to convert the original databases into the GBK format. The dnd genes of Salmonella enterica serovar Cerro 87 (CP008925: 3477655...3481641) and the ssp genes of Vibrio cyclitrophicus FF75 (NZ_ATLT01000001: 2194844...2200061) were used as PT-modifying component queries [4]. The DndF-H cluster of S. enterica serovar Cerro 87 (CP008925: 3467758 . . . 3475796) was used as the restriction component query. The PbeA-D cluster of Haloterrigena jeotgali A29 (CP031303: 460843 . . . 466000) was used as the antiviral component query. The amino acid translation of each gene sequence within the query gene cluster is searched against the selected GBK (Genebank) database, yielding a dataset of BLAST (Basic Local Alignment Search Tool) hits [13]. The BLAST hits are then mapped to their parent nucleotide scaffolds, based on the information from the database [13]. The nucleotide scaffolds are then sorted according to their empirical similarity scores with the query gene cluster [13]. The number of blast hits per gene to be mapped is 250. The weight of synteny conversion in hit sorting is 0.5. The minimal sequence coverage of blast hits is 25. The minimal identity of blast hits is 30%. The maximum distance between genes in locus is 20 kb. These parameters above were set up for multiple sequence alignments in multigeneblast (software version 1.1.13) [13].

Gene Abundance Calculation
A thenon-redundant human gut microbial gene catalog was constructed in our previous study [15]. In brief, the reads from the metagenomes were de novo assembled into contigs. Genes were predicted from the contigs and merged into non-redundant genes based on sequence similarity. The abundance of the genes was obtained by mapping the reads on the catalog, obtaining the gene-length normalized base counts, and adjusting the sequencing depth with a resampling procedure. All the genes were clustered into CAGs based on their abundance data using the canopy-based algorithm with default parameters. CAGs with more than 700 genes were regarded as bacterial CAGs for further analysis. The CAG abundance profiles were calculated as the sample-wise median gene abundance, essentially as described elsewhere [16].

Fecal DNA Preparation
All subjects gave their informed consent for inclusion before they participated in the study [15]. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of the School of Life Sciences and Biotechnology, Shanghai Jiao Tong University (No. 2012-016). Feces samples were obtained and immediately frozen on dry ice upon collection and stored at −80 • C until further analysis [15]. DNA extraction from the human fecal samples was conducted as previously described [15] and purified by the QIAamp DNA mini kit (51304, QIAGEN, Germany).

Detection of PT-Linked Dinucleotides
All the fecal DNA samples were hydrolyzed with 4 units of nuclease P1 (Sigma, St. Louis, MO, USA) and subsequently dephosphorylated by 10 units of alkaline phosphatase (Fermentas), essentially as described elsewhere [17]. For qualitation, the digested DNA samples were pre-purified by reversed-phase high-performance liquid chromatography on a ThermoHypersil Gold Aq column (250 × 4.6; SN:1292941W) at a flow rate of 0.8 mL/min with the following parameters: solvent A: water with 8 mM NH 4 OAC; B: Acetonitrile; gradient: 3% B for 17 min; 40% B for 23 min; 100% B for 10 min; 3% B for 10 min. The pre-purified samples were dried and re-suspended in 50 µL of deionized water for analysis by the Agilent 6410 Triple Quad liquid chromatography mass spectrometer, as previously described [18].
For quantification, the digested DNA samples were purified by ultrafiltration, dried, and resuspended in 40 µL of deionized water. Then, the mixture containing PT-dinucleotides was resolved on an Agilent ZORBAX SB-C18 column (2.1 × 150 mm, 3.5 µm bead size) with a flow rate of 0.3 mL/min and the following parameters: column temperature: 25 • C; solvent A: 0.1% formic acid in H 2 O; solvent B: 0.1% formic acid in acetonitrile; gradient: 4% B for 5 min, 4% to 15% B over 15 min, 15 to 20% B for 5 min, and 20 to 100% B for 5 min. The high-performance liquid chromatography column was coupled to an Agilent G6470A Triple Quadrupole mass spectrometer with an electrospray ionization source in positive mode with the following parameters: gas flow, 10 L/min; nebulizer pressure, 30 psi; drying gas temperature, 325 • C; and capillary voltage, 3000 V. Multiple reaction monitoring modes were used for the detection of ions derived from the precursor ions, with all the instrument parameters optimized for maximal sensitivity (retention time in min, precursor ion m/z, fagmentor voltage, product ion m/z and collision energy for qualitation, product ion m/z and collision energy for quantification): d(A ps A), 10

PT Genes Are Widely Distributed in the Human Microbiome
We searched for PT-modifying genes in over 40,000 bacterial genomes under the human microbiome category from websites (https://www.ncbi.nlm.nih.gov/genome/, https://www.ebi.ac.uk/, and https: //www.hmpdacc.org/). Multiple sequence alignments were performed using the genomes for strains harboring PT-modifying genes. A total of 2621 strains harboring PT-modifying genes (dataset) were obtained from the results queried by the PT-modifying genes. The amino acid sequences of DndC/SspD were used to construct a phylogenetic tree ( Figure 1). The strains in the same genus are clustered into different branches according to sequence similarity ( Figure 1). The bubbles represent collapsed clades and each bubble clustered strains with sequence identities over 60%. These bubbles were used to show the distribution of homogeneous sequences due to the strong correlation between modifying types and amino acid sequences [17]. The PT-linked dinucleotides of G ps A/G ps T, G ps G, and C ps C were separated on relatively independent clades with sequence identities over 30%. The boundaries of these clades were defined by the brackets (Figure 1). It was reported that strains with the same PT-linked dinucleotides were presumed to be closer on the phylogenetic tree [11]. Intriguingly, over 60% of the strains were clustered on distinct branches far from the clades marked with PT-linked dinucleotides ( Figure 1). Thus, we speculated that unknown modification patterns could exist in the human microbiome. 4 same genus are clustered into different branches according to sequence similarity ( Figure 1). The bubbles represent collapsed clades and each bubble clustered strains with sequence identities over 60%. These bubbles were used to show the distribution of homogeneous sequences due to the strong correlation between modifying types and amino acid sequences [17]. The PT-linked dinucleotides of GpsA/GpsT, GpsG, and CpsC were separated on relatively independent clades with sequence identities over 30%. The boundaries of these clades were defined by the brackets (Figure 1). It was reported that strains with the same PT-linked dinucleotides were presumed to be closer on the phylogenetic tree [11]. Intriguingly, over 60% of the strains were clustered on distinct branches far from the clades marked with PT-linked dinucleotides ( Figure 1). Thus, we speculated that unknown modification patterns could exist in the human microbiome.
We speculated that PT-harboring species spread across the entire human body due to the distribution of PT genes in diverse species ( Figure 1). Thus, we collected bacterial genomes directly from distinct anatomical sites on the human body (https://www.hmpdacc.org/) to perform multiple sequence alignments. Although the total number of genomes is limited, we found diverse species harboring PT in four body parts (Table 1). Notably, there are species differences in the distribution of dnd genes and ssp genes (Table 1). Intriguingly, the PT-modifying gene-harboring strains in the genus Prevotella were found only in the oral and urogenital systems (Table 1), despite the relatively abundant distribution of Prevotella species in the human intestinal tract [19]. All the data above ( Figure 1, Table 1) confirm that strains harboring PT-modifying genes are widely distributed in the human body.  We speculated that PT-harboring species spread across the entire human body due to the distribution of PT genes in diverse species (Figure 1). Thus, we collected bacterial genomes directly from distinct anatomical sites on the human body (https://www.hmpdacc.org/) to perform multiple sequence alignments. Although the total number of genomes is limited, we found diverse species harboring PT in four body parts (Table 1). Notably, there are species differences in the distribution of dnd genes and ssp genes (Table 1). Intriguingly, the PT-modifying gene-harboring strains in the genus Prevotella were found only in the oral and urogenital systems (Table 1), despite the relatively abundant distribution of Prevotella species in the human intestinal tract [19]. All the data above ( Figure 1, Table 1) confirm that strains harboring PT-modifying genes are widely distributed in the human body.

Species Difference of PT-Related Genes
It has been reported that PT gene clusters can be spread by horizontal transfer [17,20], and we find that the distribution of PT-modifying genes varied among genera. For instance, over 90% of the Mycobacterium strains harbor dnd genes rather than ssp genes ( Figure 2). On the other hand, over 60% of the Pseudomonas strains harbor ssp genes rather than dnd genes ( Figure 2). Notably, dnd genes are widely distributed in Mycobacterium, Salmonella, and Escherichia rather than the other genera ( Figure 2). To further investigate the role of PT, we performed a multiple sequence alignment for dndF-H and PbeA-D. Notably, the dndF-H gene cassette was "missing" in 90% of strains of Mycobacterium, which means that most Mycobacterium strains lack the R component ( Figure 3). On the contrary, the dndF-H gene cassette was widely distributed in Escherichia strains. Meanwhile, the pbeA-D gene cassette existed in 10% of the Mycobacterium and Vibrio, which suggests that most human microorganisms lack the antiviral component ( Figure 3). We speculated that PT may play other roles in these species without both dndF-H and PbeA-D. It seems that there exist undiscovered PT-dependent genes in their genomes. On the contrary, over 60% of the Bacillus strains harbor the R component without any modifying genes (Figure 3), which means these strains harbor undiscovered modifying systems to protect their DNA against nickase [3].

Atypical PT Gene Clusters of Different Genera
The organization of the PT-related genes varied among the genera (Figure 4). In 146 strains of E. coli and 20 strains of Bacillus cereus, the genes surrounding the PT-modifying genes encoded homologs of dndA (cysteine desulfurase) [5], dndB (negative transcriptional regulator) [18], and dndD (possibly related to DNA structure alteration) [21,22]. For instance, the identities are 23.98% between DndB and its homolog from E. coli Ecol_316 (NCBI: CP018957). Genomic islands may facilitate the horizontal transfer of PT systems and result in multiple copies of genes for E. coli and B. cereus [23]. Several strains from Bacteroidales (Prevotella and Bacteroides) showed a rare dndEi in PT-modifying

Atypical PT Gene Clusters of Different Genera
The organization of the PT-related genes varied among the genera (Figure 4). In 146 strains of E. coli and 20 strains of Bacillus cereus, the genes surrounding the PT-modifying genes encoded homologs of dndA (cysteine desulfurase) [5], dndB (negative transcriptional regulator) [18], and dndD (possibly related to DNA structure alteration) [21,22]. For instance, the identities are 23.98% between DndB and its homolog from E. coli Ecol_316 (NCBI: CP018957). Genomic islands may facilitate the horizontal transfer of PT systems and result in multiple copies of genes for E. coli and B. cereus [23]. Several strains from Bacteroidales (Prevotella and Bacteroides) showed a rare dndEi in PT-modifying

Atypical PT Gene Clusters of Different Genera
The organization of the PT-related genes varied among the genera (Figure 4). In 146 strains of E. coli and 20 strains of Bacillus cereus, the genes surrounding the PT-modifying genes encoded homologs of dndA (cysteine desulfurase) [5], dndB (negative transcriptional regulator) [18], and dndD (possibly related to DNA structure alteration) [21,22]. For instance, the identities are 23.98% between DndB and its homolog from E. coli Ecol_316 (NCBI: CP018957). Genomic islands may facilitate the horizontal transfer of PT systems and result in multiple copies of genes for E. coli and B. cereus [23]. Several strains from Bacteroidales (Prevotella and Bacteroides) showed a rare dndEi in PT-modifying genes that have an additional DNA helicase domain compared to canonical DndE [24]. In 211 Acinetobacter strains, the histidine kinase genes are adjacent to the SspE, which is the R component of the ssp system [4]. We speculate that the histidine kinase genes (Figure 4) involve in a two-component system to regulate the expression of SspE. Moreover, it was previously reported that there are complex interactions between DNA methylation and PT [25]. Intriguingly, 72% of PT-modifying genes in Mycobacterium strains lack both dndF-H and pbeA-D ( Figure 3) and are adjacent to methyltransferase genes (Figure 4). This suggests that there is a shared mechanism linking the R system between DNA methylation and the PT [25] in Mycobacterium strains. These results indicate that PT plays various roles in different species.
Biomolecules 2020, 10, x FOR PEER REVIEW 7 of 12 7 genes that have an additional DNA helicase domain compared to canonical DndE [24]. In 211 Acinetobacter strains, the histidine kinase genes are adjacent to the SspE, which is the R component of the ssp system [4]. We speculate that the histidine kinase genes (Figure 4) involve in a two-component system to regulate the expression of SspE. Moreover, it was previously reported that there are complex interactions between DNA methylation and PT [25]. Intriguingly, 72% of PT-modifying genes in Mycobacterium strains lack both dndF-H and pbeA-D ( Figure 3) and are adjacent to methyltransferase genes (Figure 4). This suggests that there is a shared mechanism linking the R system between DNA methylation and the PT [25] in Mycobacterium strains. These results indicate that PT plays various roles in different species.

Detection of PT-Linked Dinucleotides in Human Fecal DNA
For DNA methylation, restriction-modification systems are classified into four types depending on multiple factors, including subunits and target motifs [26]. In bacteria, the sequence specificity of DNA methylation depends on the target recognition domains (TRDs) of the modifying complex, and the variation in the TRDs impacts the sequence specificity upon modification [27]. For PT, although the TRDs have not been demonstrated, there is a strong correlation between the PT sequence contexts and modifying complex sequences [17]. Based on the variation in the dndC/sspD genes ( Figure S1) and the phylogenic distribution of PT-linked dinucleotides (Figure 1), it remains a possibility that there are undiscovered modifying patterns, including different motifs and PT-linked dinucleotides. Exploiting the fact that most strains possessing PT-modifying genes were assigned to the gut microbiome (Table 1), we initiated a search for new PT consensus sequences by analyzing the diversity of PT-linked dinucleotides in human fecal DNA. Notably, six new PT-linked dinucleotides (CpsG, CpsT, ApsG, TpsG, GpsC, ApsT) were detected from human fecal DNA (Table 2, Figure S2). Additionally, two PT nucleotides (CpsA, TpsA) previously discovered in mutants were detected, which were previously discovered in E. coli DH10B expressing PT-modifying genes from S. enterica 87 [17].

Detection of PT-Linked Dinucleotides in Human Fecal DNA
For DNA methylation, restriction-modification systems are classified into four types depending on multiple factors, including subunits and target motifs [26]. In bacteria, the sequence specificity of DNA methylation depends on the target recognition domains (TRDs) of the modifying complex, and the variation in the TRDs impacts the sequence specificity upon modification [27]. For PT, although the TRDs have not been demonstrated, there is a strong correlation between the PT sequence contexts and modifying complex sequences [17]. Based on the variation in the dndC/sspD genes ( Figure S1) and the phylogenic distribution of PT-linked dinucleotides (Figure 1), it remains a possibility that there are undiscovered modifying patterns, including different motifs and PT-linked dinucleotides. Exploiting the fact that most strains possessing PT-modifying genes were assigned to the gut microbiome (Table 1), we initiated a search for new PT consensus sequences by analyzing the diversity of PT-linked dinucleotides in human fecal DNA. Notably, six new PT-linked dinucleotides (C ps G, C ps T, A ps G, T ps G, G ps C, A ps T) were detected from human fecal DNA (Table 2, Figure S2). Additionally, two PT nucleotides (C ps A, T ps A) previously discovered in mutants were detected, which were previously discovered in E. coli DH10B expressing PT-modifying genes from S. enterica 87 [17]. Table 2. PT-linked dinucleotides in the human fecal DNA. Precursor  Ion   Product  Ion  z1  z2  z3  1488 1489 1494 1493

PT-Linked Dinucleotides
Blank means not detected; "Y" means sufficient levels to be quantified; "N" means close to detection limit; "*" means retention time slightly shifted from expected values; "**" means retention time shift.

The Abundance of PT Modifications Varied Among Individuals
Although we could detect PT-linked dinucleotides in the fecal DNA, the abundance of PT genes and the modification types remain unknown in the gut microbiome without quantitative analysis. Meanwhile, we speculated that the distribution of the dnd/ssp genes and modification patterns varied among individuals. Thus, we performed a quantitative analysis of the PT genes and modifications in human fecal DNA. Based on the metagenomic sequencing performed on 109 human fecal DNA samples, a gene catalog with~2 million non-redundant human gut microbial genes was constructed [15]. By BLAST alignment, 209 potential PT genes with at least 30% identity and 80% coverage compared with the PT references were identified. We were able to assemble 118 high-quality draft genomes of prevalent gut bacteria from these metagenomic data. Among these bacterial genomes, nine harbored PT genes: Bacteroides plebeius CAG00296, Bacteroidesplebeius CAG0079, Klebsiella sp. CAG00146, and Faecalibacterium prausnitzii CAG00158 were queried by dnd genes, while Bacteroides sp. CAG00020, Eubacterium ventriosu CAG00166, Clostridium bolteae CAG00012, Clostridiales bacterium CAG00239, and Clostridiales bacterium CAG00048 were queried by ssp genes. To further investigate the PT landscapes in the human gut microbiome, we analyzed the strains harboring PT-modifying genes (dataset) in the 14 individuals ( Figure 5A). Then, we quantified the PT-linked dinucleotides in these fecal samples and found seven different PT-dinucleotides ( Figure 5B). 8 CAG00146, and Faecalibacterium prausnitzii CAG00158 were queried by dnd genes, while Bacteroides sp. CAG00020, Eubacterium ventriosu CAG00166, Clostridium bolteae CAG00012, Clostridiales bacterium CAG00239, and Clostridiales bacterium CAG00048 were queried by ssp genes. To further investigate the PT landscapes in the human gut microbiome, we analyzed the strains harboring PTmodifying genes (dataset) in the 14 individuals ( Figure 5A). Then, we quantified the PT-linked dinucleotides in these fecal samples and found seven different PT-dinucleotides ( Figure 5B).

Discussion
Over 80% of the strains harboring PT-modifying genes could be assigned to Actinobacteria, Proteobacteria, and Bacteroidetes [11]. Notably, our survey in the human microbiome shows that PTmodifying genes are distributed in these phyla (Figure 1). The phylogenic tree (Figure 1) was used to demonstrate the distribution of homogeneous sequences rather than the genera distribution, though it exhibits the horizontal transfer of PT-modifying genes, which is consistent with the previous conclusions [28]. Although we could not demonstrate DndC/SspD-containing TRDs, the DndC/SspD may play the role of site recognition according to their function of inserting sulfur into the sugarphosphate backbone of DNA [4]. Thus, we built the tree by DndC/SspD and assumed that the human microbiome contained new PT-linked dinucleotides due to the phylogenic distribution of GpsA/GpsT, GpsG, and GpsC ( Figure 1). Fortunately, we discovered new PTs in human fecal DNA. We speculated that 5 μg of fecal DNA was not enough to detect multiple PT-linked dinucleotides. Thus, we tested

Discussion
Over 80% of the strains harboring PT-modifying genes could be assigned to Actinobacteria, Proteobacteria, and Bacteroidetes [11]. Notably, our survey in the human microbiome shows that PT-modifying genes are distributed in these phyla (Figure 1). The phylogenic tree (Figure 1) was used to demonstrate the distribution of homogeneous sequences rather than the genera distribution, though it exhibits the horizontal transfer of PT-modifying genes, which is consistent with the previous conclusions [28]. Although we could not demonstrate DndC/SspD-containing TRDs, the DndC/SspD may play the role of site recognition according to their function of inserting sulfur into the sugar-phosphate backbone of DNA [4]. Thus, we built the tree by DndC/SspD and assumed that the human microbiome contained new PT-linked dinucleotides due to the phylogenic distribution of G ps A/G ps T, G ps G, and G ps C (Figure 1). Fortunately, we discovered new PTs in human fecal DNA. We speculated that 5 µg of fecal DNA was not enough to detect multiple PT-linked dinucleotides. Thus, we tested 30 µg of fecal DNA with a control group of 1493. However, we found no more new PT-linked dinucleotides compared to 5 µg of fecal DNA. Additionally, the variety of PT-linked dinucleotides varied among the fecal DNA samples, which suggested that the PT-containing strains carried by each person could be classified through metadata. As shown in Table 2, z1 and z3 both carried G ps T, assigned to species queried by dnd genes (Figure 1); z2 and z3 both carried C ps C and A ps C, assigned to species queried by ssp genes (Figure 1); z2 and z3 both carried A ps T, assigned to unknown species. The qualitation and quantification of PT-linked dinucleotides confirmed the existence of PT in the human gut microbiome with a strong correlation between the PT-modifying gene abundance and the modified dinucleotide yield. Thus, the abundance of PT modifications varied among individuals due to the different quantities of PT-harboring bacteria.
Intriguingly, our data showed that the PT-modifying genes seemed widely distributed in different pathogens or opportunistic pathogens, such as the strains that belonged to A. baumannii, P. aeruginosa, C. difficile, and M. abscessus (Dataset). The different bioinformatic approaches including this work confirmed that PT-modifying genes widely spread in human microbiome-associated strains, especially pathogenic candidates [11,20,23]. The previous study reported that nearly half of the US clinical isolates of M. abscessus exhibited the PT phenotype during pulsed field gel electrophoresis [29]. It was also reported that PT-modifying genes were frequently found in pathogenic E. coli [28]. One of the possible causes is that PT modifications were presumed to enhance the antioxidative ability of these bacteria against environmental stress [14,30]. In contrast, halogens could induce PT-linked DNA breaks, which implicated that PT reduced the fitness for bacterial pathogens during human infections [12]. Thus, it remains a possibility that the microenvironment or drugs may intervene in the abundance of bacteria with PT in the human microbiome, while there are few related reports. Additionally, there are few reports about how PT-modifying genes spread in human pathogens or opportunistic pathogens.

Conclusions
PT confers bacteria with multiple physiological functions, making them a specific microbial population sharing unusual characteristics [10][11][12]. However, our understanding of PT-containing bacteria tends to be restricted to isolated strains rather than a community in the ecosystem. This study, to our knowledge, is the first time to systematically investigate the wide distribution of PT-containing microbes in the human microbiome. Our results firstly reveal the landscape of PT modifications in the gut microbiome. This work will guide future research on the transmission of PT among microbial communities in the micro-ecosystem of the human body. With the rapid growth of the human microbiome database, more PT-containing bacteria could be found to analyze the routine of PT transmission and discover the unknown motive force behind this phenomenon with effects on human health.