1. Genetics and Genomics
  2. Microbiology and Infectious Disease

A genome-phenome association study in native microbiomes identifies a mechanism for cytosine modification in DNA and RNA

  1. Weiwei Yang
  2. Yu-Cheng Lin
  3. William Johnson
  4. Nan Dai
  5. Romualdas Vaisvila
  6. Peter Weigele
  7. Yan-Jiun Lee
  8. Ivan R Corrêa Jr
  9. Ira Schildkraut
  10. Laurence Ettwiller  Is a corresponding author
  1. New England Biolabs, United States
  2. School of Dentistry, National Yang Ming Chiao Tung University, Taiwan
https://doi.org/10.7554/eLife.70021
Share this article
Cite this article
  1. Weiwei Yang
  2. Yu-Cheng Lin
  3. William Johnson
  4. Nan Dai
  5. Romualdas Vaisvila
  6. Peter Weigele
  7. Yan-Jiun Lee
  8. Ivan R Corrêa Jr
  9. Ira Schildkraut
  10. Laurence Ettwiller
(2021)
A genome-phenome association study in native microbiomes identifies a mechanism for cytosine modification in DNA and RNA
eLife 10:e70021.
https://doi.org/10.7554/eLife.70021
  2. Download BibTeX
  3. Download .RIS

Abstract

Shotgun metagenomic sequencing is a powerful approach to study microbiomes in an unbiased manner and of increasing relevance for identifying novel enzymatic functions. However, the potential of metagenomics to relate from microbiome composition to function has thus far been underutilized. Here, we introduce the Metagenomics Genome-Phenome Association (MetaGPA) study framework, which allows linking genetic information in metagenomes with a dedicated functional phenotype. We applied MetaGPA to identify enzymes associated with cytosine modifications in environmental samples. From the 2365 genes that met our significance criteria, we confirm known pathways for cytosine modifications and proposed novel cytosine-modifying mechanisms. Specifically, we characterized and identified a novel nucleic acid-modifying enzyme, 5-hydroxymethylcytosine carbamoyltransferase, that catalyzes the formation of a previously unknown cytosine modification, 5-carbamoyloxymethylcytosine, in DNA and RNA. Our work introduces MetaGPA as a novel and versatile tool for advancing functional metagenomics.

Editor's evaluation

This work will interest researchers who want to explore the functional potential of metagenomes. The authors present an original approach, MetaGPA, for performing enrichment analysis on cohorts of metagenomes and use it to identify a novel enzyme that can modify cytosines in DNA from natural bacteriophage populations.

https://doi.org/10.7554/eLife.70021.sa0

eLife digest

Many industrial processes, such as starch processing and oil refinement, use chemicals that cause harm to the environment. These can often be switched to more sustainable biological processes that are powered by proteins called enzymes.

Enzymes are micro-factories that speed up biochemical reactions in most living things. Communities of microorganisms (also known as microbiomes) are an amazing but often untapped resource for discovering enzymes that can be harnessed for industrial purposes. To gain a better picture of the microbes present within a population, researchers often extract and sequence the genetic material of all microorganisms in an environmental sample, also known as the metagenome. While current methods for analyzing the metagenome are good at identifying new species, they often provide limited information about the microorganism’s functional role within the community. This makes it difficult to find new enzymes that may be useful for industry.

Here, Yang, Lin et al. have developed a new technique called Metagenomics Genome-Phenome Association, or MetaGPA for short. The method works in a similar way to genome-wide association studies (GWAS) which are used to identify genes involved in human disease. However, instead of disease associated genes in humans, MetaGPA finds microbial genes that are associated with a biological process useful for biotechnology.

Like GWAS, the new approach created by Yang, Lin et al. compares two groups: the first contains microorganisms that carry out a specific process, and the second contains all organisms in the microbiome. The metagenome of each group is extracted and a computational pipeline is then applied to identify genes, including those coding for enzymes, that are found more often in the group performing the desired task.

To test the technique, Yang, Lin et al. used MetGPA to find new enzymes involved in DNA modification. Microbiome samples were collected from coastal water and sewage, and the computational pipeline was applied to discover genes that are associated with this process. Further analysis revealed that one of the identified genes codes for an enzyme that introduces a previously unknown change to DNA.

MetaGPA could be applied to other processes and microbiomes, and, if successful, may help researchers to identify more diverse enzymes than is currently available. This could scale up the discovery of new enzymes that can be used to power industrial reactions.

Introduction

Advances in next-generation sequencing technology have been reshaping metagenomics, making it possible to explore all microbes within a sample, including the overwhelming majority of those unculturable in standard growth conditions (Quince et al., 2017). From studying human gut microbiome to marine viral communities (Sunagawa et al., 2015), metagenomic analyses are being utilized across a broad spectrum of life science disciplines, contributing to novel clinical diagnoses (Chiu and Miller, 2019), antibiotics and small molecule discovery (Hu et al., 2013; Charlop-Powers et al., 2014), food safety (Cao et al., 2017), biofuel generation (Hess et al., 2011), and environment stewardship (Nesme et al., 2014; Cavicchioli et al., 2019).

Of the two primary questions metagenomic studies strive to address, profiling the taxonomic biodiversity – ‘what is out there?’ is easier to answer compared to inferring the biological functions – ‘what do they do?’. Indeed, there have been many well-established strategies to quantify taxonomic diversity such as analyzing known marker genes, binning contigs and assembling sequences into taxonomic groups or genomes (Sharpton, 2014; Quince et al., 2017). However, a large fraction of the functional potential of microbiomes remains to be discovered. Microbiomes encode a large number of evolutionarily diverse genes coding for millions of peptides and proteins for distinctive functions. The functional diversity is particularly enormous for secondary metabolites and epigenetic modifications (Sharon et al., 2014). Taking DNA modifications as an example , environmental microbiomes are a plentiful source of DNA-modifying enzymes with diverse mechanisms (Weigele and Raleigh, 2016; Hiraoka et al., 2019). In bacteriophages, DNA modifications are used to evade the restriction modification systems of the host bacterial cell. Notably, a handful of bacteriophages have been described to completely modify cytosines in their genomic DNA, for example, 5’ methylcytosine in XP12 phage (Kuo et al., 1968) and glycosylated 5’ hydroxymethylcytosine in T4 phage (Revel and Georgopoulos, 1969). Recently, additional base modifications have been discovered including 5-(2-aminoethoxy)methyluridine, 5-(2-aminoethyl)uridine, and 7-deazaguanine (Lee et al., 2018; Hutinet et al., 2019).

Current methods to harness such functional potential of microbiomes, for the most part, came from computational prediction of gene function which is based on homology search to existing databases. Nonetheless, because of the poor completeness and accuracy of microbial annotation, homology searches often fail to impute the correct functions (Huson et al., 2009; Nayfach and Pollard, 2016). Novel functions can also be found using functional screens; however, such effort entails the construction and high-throughput screening of a large number of clones which can be very time-consuming (Martinez et al., 2004). Novel approaches are therefore needed to link genetic information from microbial metagenomes to function.

For specific bacteria, determining the genetic basis of phenotypes has been addressed using genome-wide association studies (GWAS) (Falush and Bowden, 2006) combined with specific phylogenetic methods to account for the unique population structure of microbes (Collins et al., 2018). Examples of microbial GWAS have explored hundreds of isolates to identify genomic elements that are statistically associated with, for example, antibiotic resistance (Chewapreecha et al., 2014), host specificity (Sheppard et al., 2013), or virulence (Laabei et al., 2014). Nonetheless, these studies are limited to known isolates and have not yet been extended to entire complex microbial communities.

Here, we developed the Metagenomics Genome-Phenome Association (MetaGPA) framework to bridge the gap between genetic information and functional phenotype in complex microbial communities. MetaGPA is conceptually close to GWAS as it associates genotypic data with phenotypic traits. Contrasting with microbial GWAS which uses sequence variations as genetic markers, MetaGPA incorporated association analyses at the level of protein domains to reveal genes that are significantly associated with the phenotype of interest. By applying this workflow on DNA modifications as our phenotypic trait, we discovered a number of candidate enzyme families. From these candidates, we validated a novel DNA/RNA cytosine modification, 5-carbamoyloxymethylcytosine (5cmdC), and the enzyme responsible for this modification. From this example, we show that MetaGPA is a powerful and versatile method to improve metagenome functional analysis.

Results

Conceptual framework of MetaGPA studies

Like GWAS in individual species (Hirschhorn and Daly, 2005), MetaGPA requires definition of two cohorts, a ‘case’ cohort, that is, a group of organisms that share a specific phenotype under study in a given microbiome, and a ‘control’ cohort composed of all organisms within that microbiome (Figure 1 and Figure 1—figure supplement 1). While both cohorts are sequenced independently, all organisms included in a cohort are sequenced together without the need to isolate organisms. The ‘case’ group is derived from the control group after applying a selection to only retain a given phenotype.

Figure 1 with 1 supplement see all
Download asset Open asset
Schematic overview of Metagenomics Genome-Phenome Association (MetaGPA).

The core MetaGPA pipeline is applied to identify protein domains that are significantly associated with a given phenotype. The pipeline takes high-throughput sequencing datasets from two libraries representing (1) the total microbiome (control reads) and (2) a subset of the microbiome resulting from a given treatment that selects for the studied phenotype (case reads). Control and case reads are assembled into contigs and contigs are classified as enriched or depleted based on the ratio between the number of reads mapping to the contigs from the case versus control. Contigs are annotated based on homology with known protein domains and each domain is evaluated for its association with the studied phenotype based on the number of hits in the enriched compared to control contigs. In this theoretical example, the orange bar represents a protein domain that is associated with phenotype X (e.g., DNA cytosine modification, see below). We then apply additional metrics such as phylogenetic clustering, differential conservation, and co-occurrence network to refine this association and identify candidate genes for experimental validation. More details about the workflow are provided in Figure 1—figure supplement 1.

The association is computed using a computational workflow composed of the core MetaGPA pipeline that defines genetic units associated with the phenotype and further analysis tools to refine these associations. The core MetaGPA pipeline described in Figure 1 and Figure 1—figure supplement 1 takes sequencing reads from both cohorts and performs de novo assemblies into contigs. Contigs from both cohorts are combined and duplicated contigs are discarded. Reads from the case and control experiments are mapped back to the combined contig set. Each contig is either labeled as enriched or depleted from the case cohort based on the enrichment score calculated using the relative number of reads mapping to the contig. Genes identified in the contig set are annotated using homology search to known protein domains and domains that are found significantly enriched in the enriched contig set are considered to be associated with the studied phenotype. Finally, genes in the enriched contig set annotated with one or more associated domain(s) are defined as candidate genes.

Both associated domains and their corresponding candidate genes are further refined using evidence such as phylogenetic clustering and co-occurrence with other candidate domain families/genes. Phylogenetic clustering assesses for each associated domain, whether candidate genes are phylogenetically closer to each other relative to genes containing the same domain family in the depleted contigs. This analysis can be done at domain or residue resolution. Co-occurrence with other candidate genes strengthen the association by attempting to identify the entire metabolic pathway responsible for the specific phenotype under study. Taken together, this multilayer analysis effectively identifies protein domain families and their corresponding candidate genes that are truly related to the phenotypes of interest.

Application of MetaGPA for the discovery of cytosine modifying and hyper-modifying enzymes

Principle

To demonstrate the effectiveness of the MetaGPA framework, we designed a study to identify proteins associated with DNA cytosine modifications. Because the relevant phenotype (cytosine modifications) is covalently linked to the genetic material, total genomic DNA isolated from environmental sources can be used for both, the ‘control’ cohort and material for phenotypic selection.

Accordingly, a ‘case’ cohort is obtained by applying an enzymatic selection to retain only the genomic DNA containing known and unknown forms of cytosine modification. More specifically, unmodified cytosines are deaminated to uracils using the DNA cytidine deaminase apolipoprotein B mRNA editing enzyme catalytic polypeptide-like 3A (APOBEC3A) (Carpenter et al., 2012; Henry et al., 2009; Wijesinghe and Bhagwat, 2012), and subsequently excised by uracil-specific excision reagent (USER) (Bitinaite et al., 2007), resulting in fragmented DNA. Because APOBEC3A also deaminates 5-methyl-2'-deoxycytidine (5mdC) and, to a lesser degree 5-hydroxymethyl-2'-deoxycytidine (5hmdC) (Sun et al., 2021), ten-eleven translocation dioxygenase 2 (TET2) and T4 phage β-glucosyltransferase (T4-BGT) can be used to protect 5mdC and 5hmdC prior to APOBEC3A treatment (Figure 2A). Both modifications are used as internal control for MetaGPA since the enzymes that catalyzed their formation are well characterized. These enzymes, if present in the sample, are therefore expected to be associated with modifications in our MetaGPA framework.

Figure 2
Download asset Open asset
Selective sequencing of DNA containing modified cytosines.

(A), Schematic of enzymatic selection. With enzymatic selection treatment, adaptor ligated products were first incubated with ten-eleven translocation dioxygenase 2 (TET2) and T4 phage β-glucosyltransferase (T4-BGT) so that 5-methyl-2'-deoxycytidine (5mdC) (m) and 5-hydroxymethyl-2'-deoxycytidine (5hmdC) (hm) in these sequences may be oxidized to 5-carboxy-2’-deoxycytidine (ca) or glucosylated to 5-β-glucosyloxymethyl-2’-deoxycytidine (gm) (Materials and methods). Unmodified cytosines were deaminated by APOBEC3A into uracils and then cut by uracil-specific excision reagent (USER), so that DNA with modified cytosines is enriched after the selection. This method is predicted to also preserve unknown forms of cytosine modification (denoted ‘x’) provided that they block C-to-U deamination. Whereas in the untreated control, only USER is added to clear up damaged DNAs generated in the previous library prep step. (B) Enrichment of modified DNA with serial dilutions. A total of 250 ng genomic DNA (from Escherichia coli) and cytosine-modified genomic DNA (from T4gt) were used for enzymatic selection or control per assay. Bar graphs show average recovery of DNA from three independent experiments ± SEM.

Effectively, unmodified DNA is depleted and the remaining material in the ‘case’ cohort should only comprise DNA containing modified cytosines resistant to APOBEC3A treatment. We hypothesize that many forms of cytosine modification, including those unknown to date, are naturally protected from deamination by APOBEC3A, and should be amenable to enrichment by our selection method. Furthermore, the experimental protocol is designed to only protect fully modified DNA, that is DNA for which cytosines are modified irrespective of sequence contexts. A number of such organisms have previously been found, notably in bacteriophages (Kuo et al., 1968; Revel and Georgopoulos, 1969). We hypothesize that this design would select for DNA-modifying enzymes with little or no sequence specificity.

Cytosine-modified DNA is retained after enzymatic selection

To demonstrate the feasibility of this approach, a mock community consisting of an equimolar amount of genomic DNA from Escherichia coli (containing unmodified cytosine, dC) and T4gT phage (containing 90–95% 5hmdC) was sheared, subjected to enzymatic selection and quantitative PCR (qPCR) was performed to quantify the remaining DNA (see Materials and methods). Compared to the original mock community, enzymatic selection resulted in a 0.5% recovery of E. coli DNA and an average 35% recovery of T4gT DNA (Figure 2B). This result shows that DNA containing modified cytosine is about 70-fold enriched compared to unmodified DNA.

To test the sensitivity and efficiency of this method, we serially diluted modified DNA and unmodified DNA to 1:3, 1:10, 1:100, and 1:1000 molar ratio and carried out the enzymatic selection. Recovery rates were calculated and compared to no-enzyme treatment control. Even at 1:1000 dilution, an average of 48.6% modified DNA was retained relative to no-enzyme control. This result showed the capability of our method to concentrate trace amounts (picogram-level) of modified DNA from a complex sample (Figure 2B).

MetaGPA association analyses on complex microbiomes identifies candidate enzymes with DNA modification potential

We collected environmental samples for the aim of identifying novel enzymes involved in cytosine modification using MetaGPA. To explore the robustness of MetGPA, we included three environmental samples from distinct sources (Figure 3A); of the three samples, two were collected at a sewage treatment plant and the other one was collected from coastal sludge.

Figure 3 with 1 supplement see all
Download asset Open asset
Screening of contigs associated with DNA modification.

( A) Overview of the environmental samples and analysis used in this study. Three independent sequencing experiments were performed on three samples, two sewage samples (same sampling location at different times #1 and #2) and one coastal sample (Beverly, MA). Phage DNA was collected and extracted and for the sewage samples, two technical replicates (Rep 1 and 2) were performed. Analyses performed with separate or composite datasets were indicated in below blocks. (B) Normalized frequency of k-mer in sequencing reads from control (x-axis) compared to case (y-axis) groups. Each dot represents a unique 16-mer and is colored according to the resulting contig category it belonged to. The origin of the k-mer sequences were identified by exact search on the contigs or reverse complement of the contigs. K-mers found in the modified contigs were colored in orange; k-mers found in unmodified contigs were colored in blue; k-mers not found in contigs were marked as gray. See C below for definition of modified/unmodified contigs. Subsamples of randomly selected 0.1% of all possible k-mers were used for plotting. (C) Dotplot represents the enrichment scores for each non-redundant contig in the control (left) and case/enzymatic selection (right) cohorts. Dashed line separates modified from unmodified contigs based on a threshold enrichment score equaled to 3 (enrichment score = RPKM(selection) / RPKM(control)). The plot shows every contig from all three environmental samples. (D) Distributions of modified (orange) and unmodified (blue) contigs length (in kilobases, kb).

We aim at obtaining the phage-enriched fraction of each microbiome. The rationale for this selection is based on the fact that phages are known to carry a large diversity of modifications (Weigele and Raleigh, 2016) and phages have been shown to fully modify their genomes irrespective of sequencing context, which is a prerequisite for our MetaGPA selection. For this, each sample was precipitated with polyethylene glycol (PEG) followed by DNA extraction using phenol/chloroform (see Materials and methods). Libraries were made in duplicates (except for the coastal sample) using modified adapters to protect them from enzymatic degradation. Libraries were subjected to either enzymatic selection (i.e., case) or no enzymatic selection (i.e., control) prior to amplification (Figure 3A). Additionally, spiked-in genomic DNA mixture of E. coli (dC), XP12 (5mdC), and T4gt (5hmdC) were added to each sample after adapter ligation and before enzymatic treatment. Recovery of spiked-in modified DNAs was detected as expected (Figure 3—figure supplement 1A). Accordingly, we proceeded with Illumina high-throughput sequencing and obtained an average of 60 million paired-end reads per library. Analysis of the read composition reveals consistency of k-mer composition between replicates, demonstrating that our enzymatic selection for modified DNA is reproducible (Figure 3—figure supplement 1B). Reads from both replicates were combined and normalized k-mer frequency plots showed diversity of k-mer composition from different sources/samples, while highlighting a small portion of k-mers that was highly enriched in the selection libraries (Figure 3B).

Next, we apply our newly developed MetaGPA pipeline (Figure 1—figure supplement 1) to identify the domain families that are significantly associated with the case cohort. More specifically, reads from the case (selection) and control datasets from the three environmental samples were assembled into contigs. Contigs that were either too short (less than 1000 bp) or redundant were removed (Materials and methods). Reads were mapped to the contigs and the ratio between the normalized coverage in the case (selection) library (RPKM(case)) and the normalized coverage in the control library (RPKM(control)) defines the enrichment score for each contig (Materials and methods, Figure 3—figure supplement 1C). A high enrichment score suggests that the contig is derived from DNA containing modified cytosine (i.e., modified contig). We define as modified, contigs with an enrichment score equal or above 3. In total, 3901 modified contigs were identified from three DNA samples (Figure 3C and Figure 3—figure supplement 1C-D). Distribution of contigs across all length ranges was equal between modified and unmodified contigs (Figure 3D) and the proportion of modified contigs among sewage and coastal samples were comparable (about 2% in both sewage samples and 3.5% in coastal total contigs) (Figure 3—figure supplement 1D).

Annotations of protein domain families were performed using profile hidden Markov models of all protein domain families described in the Pfam database (El-Gebali et al., 2019; Materials and methods). Instances of each domain family in the modified and unmodified contigs were calculated and Fisher’s exact test was conducted to reveal domain families that were positively associated with DNA modification. Interestingly, there was a high degree of overlap between top associated domain families among different environmental samples (Figure 4—figure supplement 1A), suggesting that a group of universal protein families for DNA modification may exist regardless of the origin of the microbiome. Given this consistency, we decided to pool the data from three environmental samples together and repeat the association analysis to achieve higher statistical power. Associated domains were subsequently matched to open reading frames (ORFs) in the modified contigs to define candidate genes (Materials and methods). From this composite dataset, we identified 110 Pfam domain families that were significantly associated with modified contigs (Bonferroni-corrected p-value < 0.01) representing a total of 2365 candidate genes. To estimate the false positive rate of association, we used our two control replicate experiments for which no selection has been done. We then performed a MetGPA analysis with one replicate being the control group and the other replicate being the ‘case’ group. From the 6737 protein domains assessed in MetaGPA, none of the domains were found significantly associated with the case group (data not shown).

The resulting top associations (Figure 4A and Figure 4—figure supplement 1A, B) contain a number of domain families previously found to be involved in DNA synthesis/modification, validating that this approach can uncover relevant functional domains. For example, a subset of enzymes containing the thymidylate synthase domain (PF00303.20) have been shown to produce hydroxymethylpyrimidines (Neuhard et al., 1980). DNA ligase (PF14743.7, PF01068.22), and cytidine/deoxycytidylate deaminase (PF00383.24) are other domains that have been found in DNA-modifying enzymes (Bhattacharya et al., 1994; Subramanya et al., 1996). In addition, our analysis identified domain families that were not previously known to modify DNA or nucleotides and thus may be novel DNA-modifying enzymes or critical regulators.

Figure 4 with 2 supplements see all
Download asset Open asset
Metagenomics Genome-Phenome Association (MetaGPA) study at domain level.

(A) Identification of associated Pfam protein domains. Left, lists of domains ranked by p-values. Middle and right panels represent the normalized counts of each domain for modified or unmodified contigs, respectively. The top 20 domains significantly associated with modified contigs were colored in orange and the four domains significantly associated with unmodified contigs were colored in blue (Bonferroni-corrected p-value < 0.01). Data were normalized to total counts in each category (sums of all domain counts in modified or unmodified contigs), respectively. Domains occurring multiple times on the same contigs were counted only once. (B) Phylogenetic tree of carbamoyltransferase C-terminus domain. Orange and blue dots represent modified and unmodified contigs, respectively. Gray bars in the outer ring represent the enrichment scores of each contig with dashed lines showing the scales. Most of the carbamoyltransferase C-terminus domains from modified contigs form a distinctive phylogenetic branch. (C) Pfam neighborhood association for carbamoyltransferase C-terminus. Carbamoyltransferase C-terminus is centered in the middle and neighboring Pfams spanning 3 kb upstream and downstream are displayed as solid squares. Colored squares mark the top five Pfams co-occurred with carbamoyltransferase C-terminus in modified contigs. Predicted open reading frames (ORFs) on each contig are marked as hollow squares. Asterisk marks the contig containing the enzyme expressed and purified in Figure 6. (D) Correlation networks between the top 20 associated Pfams with modified contigs. Each node represents a Pfam. The thickness and length of edges were based on the p-values between each two nodes with thick and short edges corresponding to more significant relationships (small p-values). Only positive correlations with p-values < 0.05 are shown. The three interaction cores are circled in red.

To refine the candidate genes, we conducted phylogenetic analysis for each domain family significantly associated with modified contigs. Toward this end, all instances of a particular domain were aligned using a maximum likelihood model and the resulting phylogenetic tree was overlaid with the origin status of the contig (modified/unmodified). Particularly, several domains, including carbamoyltransferase N-terminus (PF02543.16) and C-terminus (PF16861.6) domains, exhibited a clustered pattern in which most of the sequences from modified contigs formed a distinct phylogenetic clade from the other sequences (Figure 4B and Figure 4—figure supplement 1C). These clades that are almost exclusively derived from modified contigs restating the association of the domain-of-interest with a potential differentiated phenotype of modification. Moreover, this analysis can serve as evidence for refined functional annotation and may suggest a subfamily with specific functions.

Complex DNA modifications in bacteriophages are usually carried out by multiple enzymes whose genes tend to cluster on the genome (Iyer et al., 2013). We therefore extended the analysis to study co-occurrences on the same contigs of domain families associated with modification (Figure 4C and Figure 4—figure supplement 2A-C). We found that, for example, the domains most frequently co-occurring with carbamoyltransferase C-terminus (PF16861.6) were carbamoyltransferase N-terminus (PF02543.16), thymidylate synthase (PF00303.20), phosphoribosyl-ATP pyrophosphohydrolase (PF01503.18), dCMP deaminase Zn-binding region (PF00383.24), and MafB19-like deaminase (PF14437.7) (Figure 4C). Congruously, domains belonging to the thymidylate synthase family also co-occurred with carbamoyltransferase N-terminus, phosphoribosyl-ATP pyrophosphohydrolase, dCMP deaminase Zn-binding region, and MafB19-like deaminase domains. Importantly, we found that these co-occurrences are often specific to modified contigs (Figure 4C and Figure 4—figure supplement 2A-C). For example, the co-occurence of carbamoyltransferase and thymidylate synthase domains is only found in modified contigs (Figure 4C). In unmodified contigs, carbamoyltransferase domains were flanked by genes with unrelated functions such as genes encoding for glycosyl transferases group 1 or tRNA N6-adenosine threonylcarbamoyltransferase domains.

Together, these results suggested a multi-domain network related to DNA modification. Unbiased co-occurrence analysis was therefore conducted for the top 20 associated domain families. Indeed, significant positive correlations were observed, and the results demonstrated three interaction cores centered on thymidylate synthases, DNA ligases, and carbamoyltransferases, respectively (Figure 4D).

Differential conservation of amino acids reveals key residues associated with DNA modification

Having demonstrated that MetaGPA can be successfully deployed to identify protein domains associated with a given phenotype, we next considered whether MetaGPA go as far as identifying associations at residue resolution. To identify such associations, we designed a ‘differential’ conservation score based on existing metrics (Valdar, 2002; Karlin and Brocchieri, 1996) that reflects the degree of association of individual residue with DNA modification (Materials and methods). This score ranks high the residues that show distinct conservation patterns whether they are derived from modified or unmodified DNA. Conversely, weakly conserved residues or residues conserved invariably between modified and unmodified DNA are ranked low.

We used the thymidylate synthase domain alignments to benchmark the differential conservation score. This well-characterized family of proteins involved in nucleotide biosynthesis is also ranked by MetaGPA as the most significantly associated with DNA modification (Figure 4A). The substrate specificity of thymidylate synthases is dictated by the residue at the position 177 (numbering relative to the E. coli thymidylate synthase sequence). Previous studies have demonstrated that changing this residue from asparagine (Asn) to aspartate (Asp) can switch the preference of a canonical thymidylate synthase from dUMP to dCMP resulting in the formation of 5mdCMP (Graves et al., 1992; Hardy and Nalivaika, 1992; Liu and Santi, 1992). We therefore hypothesized that position Asn177 should have a high differential conservation score with Asn found conserved in the unmodified contigs and Asp found in modified contigs.

To confirm our hypothesis, we aligned the 433 thymidylate synthase protein sequences identified in our composite dataset with the canonical E. coli thymidylate synthase which has been structurally and biochemically characterized (Stout et al., 1998). Our differential conservation score identifies position 177 and position 147 (relative to E. coli thymidylate synthase) as the top two positions (Figure 5A). Both residues are within the enzyme’s active site (Figure 5C; Stout et al., 1998). As expected, residues at positions 177 are mostly Asn for thymidylate synthase protein in the unmodified contigs suggesting that dUMP is the substrate for these enzymes. Conversely, thymidylate synthase proteins in the modified contigs have mostly Asp position 177 consistent with modified cytosines (Figure 5B).

Figure 5 with 1 supplement see all
Download asset Open asset
Key residues in the thymidylate synthase associated with DNA modification.

(A) Differential conservation score for each position in the thymidylate synthase alignment. Positions are relative to the Escherichia coli thymidylate synthase. (B) Multiple sequence alignment of 25 and 24 randomly selected thymidylate synthase sequences from the modified and unmodified contigs respectively together with the E. coli thymidylate synthase. Aligned residues at positions 147 and 177 (relative to E. coli thymidylate synthase) are colored according to the physicochemical properties of the amino acids. (C) Structure of the active site of E. coli thymidylate synthase (PDB 1BID) highlighting the H147 and N177 residues and the dUMP substrate.

The H147 forms a hydrogen bonding network through a water molecule to the keto oxygen of position 4 in dUMP (Matthews et al., 1990), suggesting a role in substrate specificity. Alternatively, H147 is part of a H-bond network with Y94, a key residue important for proton abstraction at C5. Therefore, residue occupancy at position 147 may accommodate differences between the pKa of C5 between cytosine and uracil substrates (Hong et al., 2007). Taken together, our differential conservation score identified functionally critical residues for thymidylate synthase activity and substrate specificity.

We then performed the same analysis with the top 20 associated Pfam domains from our list (Figure 4A). Among them, carbamoyltransferase sequences were aligned with the O-carbamoyltransferase family member TobZ that has been structurally and functionally characterized (Parthier et al., 2012). Our analysis identified three residues (W408, G421, and R423) with the highest differential conservation scores (Figure 5—figure supplement 1A-B). All three residues are located within 14 Å to the carbamoyl phosphate and ADP-binding pocket and may be important in defining the enzyme’s specificity toward cytosine derivatives (Figure 5—figure supplement 1C).

Altogether, these results indicate that MetaGPA is able to associate a single residue within a protein to the studied phenotype.

Characterization of a novel-modifying enzyme – 5-hydroxymethylcytosine carbamoyltransferase

Among the 110 domains that are significantly associated with DNA modifications, carbamoyltransferase domains are ranking within the top 20 most significant domains in our MetaGPA study (Figure 4A) and are exhibiting all three of our refinement metrics. Carbamoyltransferases are part of a large protein family catalyzing the addition of a carbamoyl group to various substrates but so far none of them has been shown to act on any form of cytosine, potentially revealing a new function and a new cytosine modification. We therefore sought to further explore these enzymes for their ability to modify cytosine.

The co-occurrence of carbamoyltransferase and thymidylate synthase homologs specifically in modified contigs (Figure 4C) resembles the arrangement of the T4 phages for which genes coding for the dCMP hydroxymethylase and β-glucosyltransferase co-occur on the genome (Miller et al., 2003). The T4 dCMP hydroxymethylase is homologous to thymidylate synthase and transfers a carbon from methyltetrahydrofolate (mTHF) to C5 of the pyrimidine ring producing an exocyclic methylene in the active site of the enzyme (Graves et al., 1992). However, unlike thymidylate synthase, the methylene intermediate undergoes nucleophilic attack by water producing a hydroxymethyl group. Following incorporation of 5hmC into DNA during replication, T4 β-glucosyltransferase transfers a glucose to the hydroxyl moiety of 5hmC. Thus, the pairing of a carbamoyltransferase with dCMP hydroxymethylase led us to hypothesize a novel form of DNA modification, in which the carbamoyltransferase catalyzes the transfer of a carbamoyl group to the nucleophilic hydroxyl acceptor group of 5hmdC producing 5cmdC (Figure 6A).

Figure 6 with 3 supplements see all
Download asset Open asset
Validation of the novel 5-hydroxymethylcytosine (5hmC) carbamoyltransferase.

(A) Predicted reactions. 1, carbamoylphosphate; 2, carbamoyladenylate intermediate; 3, 5-hydroxymethylcytidine; 4, 5-carbamoyloxymethylcytidine. (B) Liquid chromatography-mass spectrometry (LC-MS) trace and quantification for enzymatic reaction with a single-stranded DNA oligo containing an internal 5-hydroxymethyl-2'-deoxycytidine (5hmdC). Bar graph represents average conversion rates ± SEM from three independent experiments using three different DNA oligos. ND, not detected. (C) LC-MS trace and quantification for enzymatic reaction with 5-hydroxymethyl-2’-deoxycytidine triphosphate (5hmdCTP) . Bar graph shows average conversion rates ± SEM from three independent experiments. (D) Enzyme preference on NCN sequences. Genomic DNA from lambda (dC), XP12 (5mdC), and T4gt (5hmdC) were mixed at 1:1:1 ratio by molarity before being subjected to enzymatic selection. Modified ratio of each C in NCN motifs (with N being A, T, C, or G) was normalized to the no-enzyme control. (E) LC-MS quantification for enzymatic reaction with a single-stranded RNA oligo containing an internal 5hmC. Bar graph shows average conversion rates ± SEM from three independent experiments. (F) Quantification for enzymatic reaction with 5-hydroxymethylcytidine triphosphate (5hmCTP). Bar graph shows average conversion rates ± SEM from three independent experiments. Oligonucleotides and nucleoside triphosphate were converted to nucleosides before LC-MS analysis. A dG spike was added as internal reference for quantification of nucleoside triphosphates.

Our composite dataset contains 62 genes annotated as carbamoyltransferase for which 17 are found in the modified contigs. We selected a carbamoyltransferase candidate gene from a modified contig originally sequenced in sewage #2 sample containing both the thymidylate synthase and the carbamoyltransferase genes (Figure 6—figure supplement 1A). The ORF was cloned into pET28b vector, and the 63 kDa protein product was expressed in E. coli and purified (Figure 6—figure supplement 1B, Materials and Methods). The predicted reaction was tested by enzymatic assays and results showed that every substrate, namely carbamoyl phosphate, ATP, 5hmdC (genomic T4gt DNA was used as substrate in these experiments), and the enzyme were indispensable for the reaction (Figure 6—figure supplement 1C,D). The choice of carbamoyl phosphate and ATP substrates were guided by the enzymatic characterization of TobZ previously published (Parthier et al., 2012). The expected product was detected by liquid chromatography-mass spectrometry (LC-MS) and confirmed with corresponding fragmentation patterns (Figure 6—figure supplement 2A, B). Nearly 70% of 5hmdC were converted into 5cmdC in the denatured single-stranded T4gt genomic DNA. Interestingly, the activity of our carbamoyltransferase was several fold lower on double-stranded DNA, suggesting the preference of this enzyme for single-stranded DNA (Figure 6—figure supplement 1C). When using synthesized single-stranded DNA oligo containing a single internal 5hmdC site as substrate, the conversion rate was nearly 100% (Figure 6B). We also tested if the carbamoyltransferase could react with free deoxynucleoside triphosphate. LC-MS results demonstrated about 60% conversion of 5-hydroxymethyl-2’-deoxycytidine triphosphate (5hmdCTP) (Figure 6C). As expected, no activity was seen for 5-methyl-2’-deoxycytidine triphosphate (5mdCTP) or 5-hydroxymethyl-2’-deoxyuridine triphosphate (5hmdUTP), indicating the carbamoyltransferase is specific to 5hmdCTP. The fact that the carbamoyltransferase is active on 5hmdCTP (Figure 6C and Figure 6—figure supplement 3A) opens up the possibility that the reaction could take place before the nucleotide is incorporated into the phage DNA.

To examine if the carbamoyltransferase favors certain DNA sequences, we performed the enzymatic assay on a mixed genomic DNA library containing lambda (dC), XP12 (5mdC), and T4gt (5hmdC). Both untreated (control) and treated libraries were subjected to APOBEC3A deamination (see Materials and methods). Carbamoylation protects cytosine derivatives from deamination by APOBEC and thus, the difference in deamination rate between control and treated libraries is indicative of carbamoylation. We saw a decrease in deamination rate only in the T4gt genome indicating specific carbamoylation on 5hmdC. This result further validates the LC-MS results regarding the specificity of the enzyme (Figure 6D). Furthermore, all combinations of NCN motifs containing 5hmdC displayed comparable deamination levels compared to the control library, suggesting a general binding mechanism with no preferred context. This result was also consistent with enzymatic assays performed on free nucleotides, which further suggests that the in vivo carbamoylation reaction may take place prior to DNA replication.

Based on our finding that the carbamoyltransferase prefers single-stranded DNA, we further investigated whether RNA containing 5-hydroxymethylcytosine can be modified by the enzyme. LC-MS and fragmentation pattern confirmed that 5cmC is formed during the reaction, albeit at much lower yields (Figure 6E and Figure 6—figure supplement 3B-C). Carbamoylation of free nucleotide 5-hydroxymethylcytidine triphosphate (5hmCTP) was also detected (Figure 6F and Figure 6—figure supplement 3A). We thus concluded that this novel nucleic acid-modifying enzyme identified from MetaGPA studies acts on DNA, RNA, and nucleotide triphosphates.

Discussion

In this work, we report a novel analytic framework called MetaGPA that links functional phenotype with genetic information in metagenomes. Being cost effective and requiring reduced efforts, MetaGPA essentially rests on two sequencing libraries derived from the same starting material: one library reflects all the organisms in the community while the other library results from a selection to only contain organisms with a phenotype of interest. This experimental setup allows for candidate genes identification in less than a week.

Our MetaGPA study framework is phenotype-driven and the identification of candidate genes is agnostic to prior annotations. Theoretically, it can be performed with any case/control cohort pair as long as distinct phenotypes can be partitioned through selection of the case cohort. For example, these phenotypes may include but are not limited to DNA modifications, phage sensitivity to chemicals, and cell surface determinants such as O-antigen. Nonetheless, this partitioning is crucial for MetaGPA to succeed and therefore requires the development and optimization of the selection process for every new MetGPA application.

We used this approach to screen environmental phage metagenomes and successfully identified 110 domains representing 2365 candidate proteins, mostly enzymes that are significantly associated with DNA modifications. From these candidates, we validated a novel DNA/RNA-modifying enzyme responsible for the previously unknown 5cmdC/5cmC modification. To our knowledge, the only reported carbamoylated nucleotide is the 5-carbamoylmethyluridine, which is located in the first position of the anticodon of yeast valine tRNA34 (Yamamoto et al., 1985). Importantly, we predicted, at single residue resolution, the key positions in the enzyme associated with this DNA modification, enabling candidate prioritization and protein engineering.

Given the co-occurrence on the same contig of phage-specific genes (such as phage tail genes) with the 5-hydroxymethylcytosine carbamoyltransferase gene characterized in this study, the 5cmdC modification is most likely found in a bacteriophage. Because of the nature of the MetaGPA selection and the fact that the carbamoyltransferase show no sequence specificity, we predict that all the cytosines in this phage are modified, possibly conferring the ability to escape most of the restriction digestion system found in bacteria. This evolutionary advantage may explain the prevalence of this enzyme family in the phage fraction of all the microbiomes investigated in this study.

In these MetaGPA experiments, the sequencing libraries can be directly used as material for selection because the phenotype of interest (cytosine modification) is covalently attached to the genetic material. As such, scenarios for which the phenotype is covalently attached to the genetic material are the most straightforward applications of MetaGPA. However, for other scenarios for which phenotypes are not physically coupled with DNA or RNA, selection will have to be done while preserving the integrity of cells and viral particles to retain the link between phenotype and genotype. Thus, the limited availability of adequate selection processes that preserve such links may restrict the broader applicability of MetaGPA. While some selections can easily be adapted to MetaGPA, others may turn out to be difficult to implement.

For those applications for which a selection can be implemented, MetaGPA provides notable benefit for de novo exploratory discoveries especially with largely unknown microbial communities such as the virulome where existing knowledge is limited.

Data access

All raw and processed sequencing data generated in this study have been submitted to the NCBI Sequence Read Archive (SRA; https://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA714147.

Materials and methods

Genomic DNA

Request a detailed protocol

The E. coli K-12 MG1655, XP12, and T4gt genomic DNA used in this study were obtained from NEB.

Environmental phage collection

Request a detailed protocol

For each batch, 2–4 L of sewage or coastal seawater were used for phage collection. Large debris and bacterial cells were pelleted and removed by centrifuging at 5000× g for 30 min at 4°C. Phage particles in the supernatant were precipitated by adding PEG8000 to 10% (w/v) and NaCl to 1 M and let stand at 4°C overnight. Aggregates of phage particles were pelleted at 10,000× g for 30 min at 4°C, washed with 1 mL solution containing 10 % PEG8000 and 1 M NaCl, and resuspended in 2–4 mL of phage dilution buffer (10 mM Tris-HCl at pH 8.0, 10 mM MgCl2, 75 mM NaCl). The crude phage particle suspension was stored at 4°C for subsequent phenol-chloroform DNA extraction.

Phenol-chloroform DNA extraction

Request a detailed protocol

Two to four milliliters of crude phage suspension was divided in 400 μL aliquots. For each aliquot, phage particles were lysed at 56°C for 2 hr in 550 μL of lysis buffer (100 mM Tris-HCl at pH 8.0, 27.3 mM EDTA, 2% SDS, ~1.6 U Proteinase K [NEB #P8107]). After lysis, RNase A was added to 10 μg/mL and the reaction was incubated at 37°C for 30 min; 1× volume (~550 μL) of phenol-chloroform (Tris-HCl buffered at pH 8.0) was mixed with the lysis solution and vortexed vigorously for ~1 min, and centrifuged at 10,000× g for 5 min for phase separation. The top aqueous layer (~500 μL) was collected and mixed with 1× volume (~500 μL) of chloroform, vortex vigorously, and centrifuged for phase separation. The top aqueous layer (~450 μL) was collected; 1× volume of isopropanol was slowly added on top of the aqueous solution. Phage DNA was ‘spooled’ with a glass capillary by swirling and mixing isopropanol with the aqueous solution. The spooled DNA was washed in 70% ethanol, dried at room temperature for ~30 min, and dissolved in ~600–800 μL of TE buffer (10 mM Tris pH 7.5, 1 mM EDTA).

The phage DNA solution was further purified by ethanol precipitation. Briefly, DNA was precipitated by adding 0.1× volume of 3 M sodium acetate and 2.5× volume of 100% ethanol and incubated at −20°C overnight. Precipitated DNA was pelleted at 16,000× g for 20 min, washed twice with 1 mL of 70% ethanol, dried at room temperature, and finally dissolved in 200 μL of TE buffer for storage at −20°C. On average more than 20 μg of DNA was extracted in each batch.

Illumina library preparation

Request a detailed protocol

For each environmental sample, 1 μg of metagenomic DNA was sheared to 300 bp in 130 μL of TE buffer (10 mM Tris pH 7.5, 1 mM EDTA) using Covaris S2 Focused Ultrasonicator; 1.3 μL of 10 mg/mL RNase A (Qiagen #1007885) was added and incubated at 37°C for 30 min to remove RNA. To remove EDTA, the sheared DNA was purified with Zymo Oligo Clean & Concentrator Kit (#D4061) and eluted in 50 μL of 1 mM Tris buffer (pH 7.5).

One reaction of NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB #E7645) was used for 1 μg of input DNA, with the following modifications to the standard protocol: 5mdC Y-shaped Illumina adaptors were used to protect the adaptor from subsequent enzymatic treatment. The dCTP was replaced with 5mdCTP in the end repair reaction (5mdCTP was used instead of regular dCTP to protect end-repaired fragments from subsequent enzymatic treatment).

The DNA library was purified with 1× volume of NEBNext Sample Purification Beads (NEB #E7103) and eluted with 40 μL of 1 mM Tris buffer (pH 7.5).

For each of the two sewage DNA samples, experiments were performed in duplicate. Each one contained two pairs of replicate libraries subjected to enzymatic selection or control, respectively. The coastal sample generated only one pair: one library for enzymatic selection and one for control.

Enzymatic selection protocol

Request a detailed protocol

First to test the recovery of modified DNA with enzymatic selection, mixed genomic DNA (E. coli and T4gt) were prepared at various dilutions. A total of 250 ng mixed DNA was used per reaction; 1 µL TET2 (NEB #7120S) and 1 µL T4-BGT (NEB #M0357S) were added to the 50 µL reaction mixture containing 1× TET2 reaction buffer, 40 µM UDP-glucose, and 40 µM iron(II) sulfate hexahydrate. After 60 min incubation at 37°C, proteinase K was added at 0.4 mg/mL to inactivate the enzymes. Products were purified with Zymo Oligo Clean & Concentrator kit (#D4061) and eluted in 16 µL water. To denature double-stranded DNA, 4 µL formamide (Sigma #11814320001) was added. The 20 µL mixture was then incubated at 95°C for 10 min and immediately transferred to an ice bath. One µL APOBEC3A (NEB #E7120S) was added directly to the reaction with 10 µL of 10× APOBEC3A reaction buffer and 1 µL BSA (10 mg/mL). The reaction volume was brought up to 100 µL with water. APOBEC3A-mediated deamination was conducted at 37°C for 3 hr. Purification was performed using Zymo Oligo Clean & Concentrator kit and elution with 43 µL of water. In the final step, the reaction was incubated with 2 µL of USER (NEB #5508) in 1× CutSmart Buffer at 37°C for 15 min before final purification with Zymo Oligo Clean & Concentrator kit. Purified samples were then used for qPCR in the next step.

For each prepared phage library sample, 100 ng spiked-in genomic DNA mixture (E. coli:XP12:T4gt = 1:1:1 by molarity) were added to the library before being subjected to enzymatic selection protocol described above. The final library was eluted in 50 μL of 1 mM Tris buffer (pH 7.5).

Quantitative PCR

Request a detailed protocol

The qPCR reactions were performed with enzymatic selection or control samples using Luna Universal qPCR Master Mix (NEB #M3003S) on a Bio-Rad CFX96 real-time PCR detection system. Two µL of purified DNA (diluted 100-folds) were added per reaction. Primers used in the experiments were the following: E. coli F: 5’-TTGCTGAGTTTCACGCTTGC, E. coli R: 5’-AAAACCGCTTGTGGATTGCC, T4gt F: 5’-TCGCGAAACGGTTTTCCAAG, T4gt R: 5’-AAAGCGCTTGACCCAACAAC, XP12 F: 5’-TGCGATGTTGGATTCGTTGG, and XP12 R: 5’-ACAACCCGCCATAATGGAAC. Recovery was normalized to control using the delta-delta Ct method.

Illumina sequencing

Libraries were indexed (with NEBNext Multiplex Oligos for Illumina E7335), amplified using NEBNext Ultra II Q5 Master Mix (six cycles for control library and 12 cycles for selection library) and pooled for sequencing on an Illumina Nextseq instrument with paired end reads of 75 bp.

Reads trimming and k-mer analysis

Request a detailed protocol

Paired-end reads were downloaded as FASTQ files and trimmed with Trim Galore v0.6.4 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) using the --paired option. k-mer counting from reads was done with JELLYFISH v2.2.10 (Marçais and Kingsford, 2011) and 16-mer was chosen based on best resolution.

Computational workflow

Request a detailed protocol

The complete workflow contains the core MetaPGA pipeline and additional association analyses. The main steps included in the core pipeline are the following: first, trimmed reads are processed and mapped to contigs generated by de novo assembly. Annotation of protein domains/Pfams on each assembled contigs was also done at this step. Second, using the normalized mapped reads, the enrichment scores for every assembled contigs were calculated and modified contigs were selected based on the chosen threshold. Then, we applied statistical analysis to every protein domain/Pfam to determine if it is significantly associated with the modification phenotype. Hence, a list of candidate protein domains associated with the modification are generated. Last, MetaGPA provides three different types of association analysis (phylogenetic clustering analysis, residue differential conservation analysis, and multi-domain co-occurrence network analysis) to further guide if a custom chosen protein domain is a good candidate. Detailed methodology for each step is provided below. The scripts and detailed documentation of the pipeline are available at https://github.com/linyc74/MetaGPA, (Lin, 2021, copy archived at swh:1:rev:28c23f47fcf108b0e7b8851b92f37921358c8e8e).

Sequencing data processing

Request a detailed protocol

De novo assembly of contigs for each sample was performed with SPAdes v3.13.0 (Nurk et al., 2017) with the --meta option. We selectively reported contigs longer or equal to 1000 bp. To remove redundant contigs between case and control samples, we used CD-HIT v4.8.1 (Fu et al., 2012) nucleotide mode cd-hit-est with sequence identity threshold set to 0.95 (sequences with more than 95% similarity are considered redundant). We set this threshold due to the fact that many microorganism genomes are related. Other options used were -n 10 -d 0 M 0 T 4. The remaining non-redundant contigs were annotated with HMM-based Pfam entries (Pfam-A) using HMMER v3.3 (http://hmmer.org/). Alignment of reads onto contigs was done with BOWTIE2 v2.3.5.1 (Langmead and Salzberg, 2012) together with SAMTOOLS v1.9 (Li et al., 2009) to generate, sort, and index bam files for later analysis. ORFs on each contig were predicted with GLIMMER 3.02 (Delcher et al., 2007). We used the long-ORFs program with options --cutoff 1.15 and --linear to identify long ORFs that are very likely to contain genes. Other options used in glimmer3 program were --max_olap 100, --gene_len 110, and --threshold 30.

Contig enrichment score calculation

Request a detailed protocol

The enrichment score for each contig was calculated using the normalized mapped reads (reads per kb per million, RPKM) from selection and control as follows: enrichment score = RPKM(selection) / RPKM(control). The mapped reads counts were generated with Multicov using BEDTOOLS v2.29.2 (Quinlan and Hall, 2010). Contigs with higher enrichment score represent more mapped reads in case library relative to control library, therefore, are more likely to be associated with modification. We considered contigs with an enrichment score greater or equal to three to be modified and the rest unmodified. The calculation was done individually for three independent experiments.

Fisher’s exact test and correction to determine associated Pfams

Request a detailed protocol

The information including the number and type of Pfams on each contig was obtained with hmmsearch in the annotation step. We then separately counted the numbers of modified and unmodified contigs containing each type of Pfam. To avoid redundant counting, Pfams occurred multiple times on the same contig were counted only once. Fisher’s exact test was performed for each Pfam to identify if the count difference between the modified and unmodified contig group is significant. Because large-scale multiple testing was conducted for each Pfam, we did the Bonferroni correction to adjust the p-value. Both tests were performed in python with SciPy or Statsmodels modules.

Phylogenetic analysis

Request a detailed protocol

For each Pfam of interest, the protein sequences from contigs containing the Pfam were aligned with MUSCLE v3.8.1551 (Edgar, 2004). The resulting aligned fasta files were subjected to construct phylogenetic trees using the maximum likelihood method in the phylogenetic analysis program RAxML v8.2.12 (Stamatakis, 2014). We chose the -f a option to do rapid bootstrap analysis and the -m PROTGAMMAAUTO model to automatically determine the best protein substitution model to be used for the dataset. The parsimony trees were built with random seeds 1237. The online tool iTOL (https://itol.embl.de/) was used to visualize trees.

Co-occurrence network analysis

Request a detailed protocol

The presence-absence matrix with rows being the Pfams and columns being the contigs was generated with annotation output file from the previous step. We specifically performed co-occurrence analysis in the R package co-ocur v1.3 (Griffith et al., 2016) for the top 20 Pfams associated with modified contigs. Significant positive correlations (p-value < 0.05) were exported and the network was visualized in Cytoscape v3.8.0 (Shannon et al., 2003) with prefuse force directed layout.

Differential conservation score

Request a detailed protocol

Protein sequences were assigned to two groups according to whether they were encoded on modified or unmodified DNA. After multiple sequence alignment, positions that have less than 50% residues present were ignored. Differential conservation score was calculated at each aligned position. For each position in the alignment, intra-group similarity scores were calculated by the average of all possible ‘within-group’ pairwise similarities, while the inter-group similarity score was calculated from all possible ‘across-group’ pairwise similarities using the BLOSUM80 matrix. For a given multiple sequence alignment column, let N1 and N2 be the number of residues for the modified and unmodified groups, respectively, the two intra-group similarity scores (Imodified and Iunmodified) were defined as:

Imodified=∑i=1N1∑j§amp;gt;iN1Mai,aj×2N1(N1-1)
Iunmodified=∑i=1N2∑j§amp;gt;iN2Mai,aj×2N2(N2-1)

where Mai,aj is the value of amino acid pair ai and aj in the BLOSUM80 matrix. The inter-group similarity score (J) was defined as:

J=∑i=1N1∑j=1N2Mai,aj×1N1N2

The differential conservation score (S) was defined as the average of two intra-group similarity scores subtracted by the inter-group similarity score.

S=Imodified+Iunmodified2-J

Expression and purification of carbamoyltransferase

View detailed protocol

The following carbamoyltransferase sequence : MSDLLLTLGHNASAIAISVGDDGAAKVENAYELERLTGKKSDSAFPIDAIIALKERGMDKIDRVYVSHWSPTGRVDDLKAKYWDRSIFPPHVPVITQESMNLTHHDCHAQAAMAFAGSSFPTKDTGVLVVDGFGNLAEHLSYYRVQAGGQLHLMRRWYGYGTSLGLMYQYATSFLGLKMHEDEYKLLGYGARVATIGCDMDVLNQRIFTEAQAFLKRFRSLNSFEMSPDLAGLPAVQEKWAERFAAILDDVGFKGSSSTYEARCIVGYAVQQLLEIVIRNLFMADLPKPTNLIVTGGVAFNVELNRMLLGLIPGKLCVMPLAGDQGNALGLWAFSNRRAKLDFGDLCWGRREMTLGEPGPDTPLPDGMIVVEHDTPAVYEAIAEQLKTVGFINIVRGNMEFGPRALCNTTTLARADDRAVVEEINRINGRDTVMPFAPVVSAHEWLRYFPDASRLHRSAEFMICAVQYAPGLGEQVPGAALRTVKGLYTGRPQVYSSKYEWDSVTRILDDYGLLINTSFNVHGVPICLDLKHVVHSHQFQRERNPNVRTIVIAN* was extracted from de novo assembled contigs. The expression plasmid was synthesized from GenScript. Two 6× His-tags were co-expressed at both the N-terminus and the C-terminus of the recombinant protein using T7 Express Competent E. coli (NEB #C2566). Cells were cultured in LB media until an OD600 of 0.6 was reached and induced with 0.4 mM IPTG (Growcells #MESP-2002) for protein expression. One µM iron(II) was also added to facilitate protein folding. The induced cultures were maintained at 16°C in a shaker at 200 rpm for 23 hr. Cells were harvested by spinning down cell pellets at 3500 rpm at 4°C for 30 min. Cell pellets from 4 L culture were resuspended in 160 mL buffer A containing 20 mM Tris pH 7.5, 500 mM NaCl, 0.05% Tween-20, 20 mM imidazole, and sonicated using a Misonix S-4000 sonicator with 20 s on and 20 s off cycles until an OD260 plateau was reached. Cell lysates were spinned down at 13,000 rpm for 30 min in a pre-chilled centrifuge at 4°C. The supernatant was separated and combined with 0.2 mM PMSF (Sigma #78830); 50 mL of supernatant was loaded on AKTA (GE Healthcare Life Sciences) with 1 mL Histrap column (GE Healthcare Life Sciences) pre-equilibrated with buffer A. The column was washed with 50 mL buffer A and eluted with a gradient of buffer B containing 20 mM Tris pH 7.5, 500 mM NaCl, 0.05% Tween-20, and 500 mM imidazole. Aliquots containing concentrated proteins were pooled and diluted 1:1 with 20 mM Tris pH 7.5, 5% glycerol and 0.05% Tween-20. The diluent was reloaded on AKTA with 5 mL Hitrap Q HP column (GE Healthcare Life Sciences), followed by a wash with 35 mL buffer containing 20 mM Tris pH 7.5, 100 mM NaCl, 5% glycerol, and 0.05% Tween-20 and eluted with a gradient of buffer containing 20 mM Tris pH 7.5, 1 M NaCl, 5% glycerol, and 0.05% Tween-20. Finally, collected fractions with concentrated proteins were pooled and mixed with equal volume glycerol for storage at −20°C.

Carbamoyltransferase enzyme assay

View detailed protocol

For enzyme assay using T4gt genomic DNA as substrate, 10 min incubation at 95°C was performed to denature double-stranded DNA and the sample was immediately transferred to ice bath to prevent re-annealing. Then 0.38 nM denatured DNA was used for each 50 µL reaction with 1× NEBBuffer2.1 (NEB #B7202S), freshly prepared 10 µM iron(II) sulfate hexahydrate (Sigma #203505), freshly prepared 10 mM carbamoylphosphate and 5 mM ATP. Carbamoyltransferase was added to the reaction at 7.2 µM. The reaction mixture was incubated at 30°C for 3 hr before adding 2 µL Proteinase K to inactivate the enzyme. After 30 min incubation at 37°C with proteinase K, DNA was purified with Zymo Oligo Clean & Concentrator kit. For assays with synthesized single-stranded DNA oligos containing 5hmdC, the heat-denaturing step was omitted. Oligos were added at 1.6 µM per 50 µL reaction with the same concentration of carbamoyltransferase and other components added as listed before. Purification was performed using Norgenbiotek Oligo Clean-up and Concentrator kit (#34100). For assays with free nucleotides, 0.5 mM of the corresponding nucleotide was used per reaction. For assays with synthesized RNA oligos containing 5hmC, 1.57 µM RNA was added per reaction.

LC-MS and fragmentation analysis

View detailed protocol

Genomic DNA and synthetic oligonucleotides were digested to nucleosides by treatment with the Nucleoside Digestion Mix (NEB #M0649S) at 37°C for 3 hr. The resulting nucleoside mixtures were directly analyzed by reversed-phase LC/MS or LC-MS/MS without further purification. Nucleoside and nucleotide analyses were performed on an Agilent LC/MS System 1200 Series instrument equipped with a G1315D diode array detector and a 6120 Single Quadrupole Mass Detector operating in positive (+ESI) and negative (-ESI) electrospray ionization modes. LC was carried out on a Waters Atlantis T3 column (4.6 mm × 150 mm, 3 μm) at a flow rate of 0.5 mL/min with a gradient mobile phase consisting of 10 mM aqueous ammonium acetate (pH 4.5) and methanol. MS data acquisition was recorded in total ion chromatogram mode. LC-MS/MS was performed on an Agilent 1290 UHPLC equipped with a G4212A diode array detector and a 6490A triple quadrupole mass detector operating in the positive electrospray ionization mode (+ESI). UHPLC was performed on a Waters XSelect HSS T3 XP column (2.1 × 100 mm, 2.5 µm particle size) at a flow rate of 0.6 mL/min with a binary with a gradient mobile phase consisting of 10 mM aqueous ammonium formate (pH 4.4) and methanol. MS/MS fragmentation spectra were obtained by collision-induced dissociation in the positive product ion mode with the following parameters: gas temperature 230°C, gas flow 13 L/min, nebulizer 40 psi, sheath gas temperature 400°C, sheath gas flow 12 L/min, capillary voltage 3 kV, nozzle voltage 0 kV, and collision energy 5–65 V.

Sequence preference of carbamoyltransferase

Request a detailed protocol

Library preparation was performed except the following modifications to the standard protocol: (1) we did not perform RNase A treatment for this experiment; (2) we used pyrrolo-dC Y-shaped adaptor instead of regular adaptor so that they are protected from subsequent enzymatic treatment. For each library, 1 µg genomic DNA mixture (Lambda:XP12:T4gt = 1:1:1 by molarity) was used. After adapter ligation, DNA libraries were purified with 1× volume of NEBNext Sample Purification Beads (NEB #E7103) and eluted with 20 μL nuclease-free water. The sample was then denatured by heating to 95°C for 10 min and subjected to carbamoyltransferase reaction as described above. Carbamoyltransferase was added to the reaction at 7.2 µM for every 1 µg DNA library. Purification was performed using Zymo Oligo Clean & Concentrator kit and eluted with 16 µL water. Purified DNA samples were heated at 90°C with 4 µL formamide to generate single-stranded fragments for the deamination reaction. One µL APOBEC3A was added per reaction to both carbamoyltransferase-treated or control (untreated) samples with 10 µL of 10× APOBEC3A reaction buffer and 1 µL BSA (10 mg/mL). The reaction mixture was incubated at 37°C overnight. Final libraries were purified using Zymo Clean & Concentrator kit, indexed (with NEBNext Multiplex Oligos for Illumina E7335) and amplified with NEBNext Q5U Mater mix (NEB #M0597). Sequencing was performed on an Illumina Mi-seq instrument with pair-end reads (2 × 75 bp). Raw reads were trimmed with TrimGalore. Methylation was analyzed with Bismark v0.22.3 and plotted with RStudio v3.6.3.

Synthesis of 5hmC RNA oligonucleotide

Request a detailed protocol

Forward and reverse DNA templates were annealed at 95°C for 4 min and slowly cooled for 20 min. RNA synthesis was performed with HiScribe T7 High Yield RNA Synthesis Kit (NEB #E2040). One µg of annealed DNA template was used per reaction with 1.5 µL T7 RNA Polymerase Mix. 5hmCTP was used with the other three nucleotides ATP, UTP, and GTP at 7.5 mM each. The reaction was incubated at 37°C for 4 hr. Two µL nulease-free DNase I were added to each reaction to digest DNA templates, followed by incubation at 37°C for 15 min. Synthesized RNA was purified with Norgenbiotek Oligo and Concentrator kit and stored at −80°C.

Nucleotides and synthesized oligos

Request a detailed protocol

Single-stranded DNA oligos used in enzymatic assays were purchased from IDT. The sequences are as follows:

  • 5hmdC-1: 5'-TGTCCGATAGACT{5hmdC}TACGCA;

  • 5hmdC-2: 5'-AACTCGCCGAGGATTT{5hmdC}TAC;

  • 5hmdC-3: 5'-{Fam-AmC6}ACACCCATCACATTTACAC{5hmdC}GGGAAAGAGTTGAATGTAGAGTTGG.

The DNA templates for synthesizing RNA were purchased from IDT as follows (T7 promoter sequence was underlined):

  • Forward: 5'-GACCTAATACGACTCACTATAGGGAGTGAGAAGATGGTCTAGGTGTTTATTGGTGATGAA.

  • ComRev: 5'-TTCATCACCAATAAACACCTAGACCATCTTCTCACTCCCTATAGTGAGTCGTATTAGGTC.

5hmdCTP (D1045) and 5mdCTP (D1035) were purchased from Zymo Research. 5hmdUTP (N-2059) and 5hmCTP (N-1087) were purchased from Trilink Biotechnologies.

Code availability

Request a detailed protocol

Custom-built bioinformatics pipelines are available at https://github.com/linyc74/MetaGPA (Lin, 2021, copy archived at swh:1:rev:28c23f47fcf108b0e7b8851b92f37921358c8e8e).

Data availability

All raw and processed sequencing data generated in this study have been submitted to the NCBI Sequence Read Archive (SRA; https://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA714147.

The following data sets were generated
    1. Yang W
    (2021) NCBI Sequence Read Archive
    ID PRJNA714147. Metagenomics genotype and phenotype association analysis on DNA modification.
    https://www.ncbi.nlm.nih.gov/bioproject/PRJNA714147

References

Decision letter

  1. María Mercedes Zambrano
    Reviewing Editor; CorpoGen, Colombia
  2. Gisela Storz
    Senior Editor; National Institute of Child Health and Human Development, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "A Genome-Phenome Association study in native microbiomes identifies a mechanism for cytosine modification in DNA and RNA" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by Reviewing Editor Maria Zambrano and Michael Marletta as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this letter to help you prepare a revised submission.

Essential revisions:

The work to obtain functional information from metagenomes is both novel and promising but various issued need to be addressed to support the claim of reproducibility and broader applicability. In particular, the manuscript must provide a more detailed methodology, substantially improve the characterization of their MetaGPA approach, particularly the bioinformatics, and clearly state the limitations of this method and its capacity to be used as a general tool for other searches and metagenome functional analysis.

1) Please take into account previous work, describe how these relate to the presented work and clearly indicate the motivation for MetaGPA in the introduction. The relationship between the author's program and previous attempts at GWAS-like analyses in microorganisms are not thoroughly described in the Introduction. For example, the authors do not mention TreeWAS or other similar approaches.

2) The Methods section contains information that presumably relates to the pipeline, but there is no explicit mention of MetaGPA. There are also seemingly pipeline-related methods that are not part of the codebase. The beginning section of the Results is insufficient to explain the MetaGPA method. A reproducible example is a must.

3) An attempt was made to run the code to confirm reproducibility. The GitHub link contains no information on installation or examples of how to run the code on test or new data. The code itself contains hard-coded file paths that would make it difficult, if not infeasible, to run on another person's machine. There are also dependancies that are unstated. The code must be substantially better documented to be of utility outside this study.

4) The reasons for the use of Tet2 and BGT to modify 5mC and 5hmC are unclear. The newly discovered 5-carbamoyloxymethylcytosine is unlikely to be deaminated by A3A

and hence what is the need for Tet2 and BGT? In fact, other than 5-methylcytosine, larger modifications of cytosine would not be substrates for APOBEC3A (PMID: 28472485). The investigators could treat all the DNA with A3A, divide it into two halves, and then use the USER kit to destroy DNA that contains uracils in one-half of the samples. The DNA that survives should contain a cytosine modification that protected it against A3A. Comparison of sequences of the two populations would show that only a small fraction of DNA has survived A3A+User treatment. If the ability of A3A to deaminate 5mC is a major concern, the authors do not articulate it. However, this can be easily remedied by replacing A3A with APOBEC3B-CTD or APOBEC3G-CTD neither of which is good at deaminating 5mC.

5) The USER kit contains E. coli Exonuclease VIII. This is problematic because this enzyme will excise oxidized pyrimidines in DNA. Thus, the DNA of any organism that routinely modifies its pyrimidines in such a way that it becomes susceptible to ExoVIII, would be eliminated from the case DNA regardless of whether it contains uracils. An AP endonuclease or simple treatment with NaOH and heat would be preferable here.

6) This methodology is designed to work for cytosine modifiers that near 100% efficient and are not sequence-specific. I do not see how this methodology could isolate genes for sequence-specific cytosine modifiers. If an enzyme strongly prefers a certain sequence motif, say 5'-TpC, then all the cytosines in the VpC sequence context (V is not T) would remain unmodified. These would be deaminated by A3A and virtually all the DNA from that organism would be destroyed by the User kit. The methodology may also have difficulty isolating genes for enzymes that modify only a fraction of the cytosine bases, say 25%. Such DNA would also be destroyed during the A3A+User treatment. Such limitations of the methodology should be carefully examined.

7) How does MetaGPA handle phylogenetic resampling (i.e., dealing with the fact that genomes are related)? This is particularly important for microorganisms. It would have been preferable to see the author's method benchmarked in a similar way, if not compared against, previous approaches to similar problems.

8) The manuscript jumps quickly into the main finding of cytosine modification. What are other applications of this technique? How could one incorporate a negative control to quantify FP rates or incorporate other controls?

9) I am somewhat confused by the discussion about the putative thymidylate synthase homologs. The authors point out that several TS homologs contain a change that is equivalent to a N177D change in E. coli TS. They further note that such a change in the E. coli enzyme results in a change of substrate from dUMP to dCMP. Does this mean that these phage genes code for dCMP methyltransferases, not thymidylate synthases? If so, they may be novel enzymes, as E. coli does not contain a dCMP methyltransferase and unlike TS, the normal methyl donor for cytosine methylation is S-adenosylmethionine. The methyl donor for TS is tetrahydrofolate. It would be useful to characterize these variant TS enzymes for substrate specificity and cofactor requirements.

10) It is unclear how the investigators arrived at the substrate and co-factor requirements of the newly discovered enzyme. The classic example of a carbamoyltransferase is ornithine carbamoyltransferase which tranfers the carbamoyl moiety to a nitrogen not to an oxygen. It also appears (https://www.brenda-enzymes.org/enzyme.php?ecno=2.1.3.3) that ornithine carbamoyltransferase does not require ATP as a co-factor. With that in mind, the putative new carbamoyltransferase could have modified the N4 of cytosine without the need for ATP. Did the investigators explore this possibility? Overall, a clearer rationale needs to be provided for how the requirements for the enzyme were established.

11) Interestingly, E. coli bacteriophage Mu carries N6-carbamoylmethyladenine. Is the enzyme that performs this reaction in any way related to the newly discovered enzyme?

12) The only justification for using phage DNA as the source of new enzymes provided by the authors is "phages are known to carry a large diversity of modifications". Another obvious reason for using phage DNA is that few phage genes contain introns. However, this convenience comes with the likely drawback that many DNA modification enzymes in cellular genomes are being missed from their screen. The authors should carefully address this issue.

13) Figure 3b is unclear. What are the units of the color (scale) bars? What am I supposed to take away from this panel? What is red and what is blue (beyond the minimal description in the figure legend)?

Reviewer #2 (Recommendations for the authors):

1. Page 2- Change "who is out there?" to "what is out there?".

2. Page 6- "unmodified cytosines are deaminated to uracils using the DNA cytidine

deaminase Apolipoprotein B mRNA editing enzyme catalytic polypeptide-like 3A

(APOBEC3A) (Carpenter et al., 2012)"- This is misleading on two counts. Among the APOBEC3 subfamily of cytosine deaminases, APOBEC3A is most efficient at deaminating 5mC, in addition to C. Furthermore, the ability of APOBEC3A to cause cytosine deamination was first demonstrated in HBV genome (PMID: 19169351) and the purification of the enzyme and its biochemical characterized was first reported by a different research group than cited in the manuscript (PMID: 22798497).

3. Page 11- "DNA ligase (PF14743.7, PF01068.22), and Cytidine deaminase (PF00383.24) are other domains that have been found in DNA modifying enzymes (Subramanya et al., 1996) (Bhattacharya et al., 1994)"- Did you mean DNA-cytosine deaminase? A cytidine deaminase would convert the ribonucleoside cytidine to uridine.

4. Page 18, last paragraph- The comparison of T4 genome with the contig in which carbamoyltransferase gene is found is confusing. Why is the occurrence of dCMP hydroxymethylase and β-glucosyltransferase in T4 genome "resembles" the occurrence thymidylate synthase and carbamoyltransferase in the contig? Phages have very compact genomes and tend to aggregate genes with related functions. Furthermore, dCMP hydroxymethylase is presumably an oxidase and not a transferase. If this is what led to the prediction of the reaction of the newly discovered carbamoyltransferase, it was a really inspired guess- I say that in all sincerity.

https://doi.org/10.7554/eLife.70021.sa1

Author response

Essential revisions:

The work to obtain functional information from metagenomes is both novel and promising but various issued need to be addressed to support the claim of reproducibility and broader applicability. In particular, the manuscript must provide a more detailed methodology, substantially improve the characterization of their MetaGPA approach, particularly the bioinformatics, and clearly state the limitations of this method and its capacity to be used as a general tool for other searches and metagenome functional analysis.

The manuscript must provide a more detailed methodology, substantially improve the characterization of their MetaGPA approach, particularly the bioinformatics:

We have significantly restructured the bioinformatic pipeline (see answer to comment 3) and provided a test dataset to benchmark it. We have also re-written the first part of the manuscript describing the MetaGPA framework (see answer to comment 2).

Clearly state the limitations of this method: we have added a entire paragraph on the limitation of MetaGPA in the revised version of the discussion:

“In these MetaGPA experiments, the sequencing libraries can be directly used as material for selection because the phenotype of interest (cytosine modification) is covalently attached to the genetic material. As such, scenarios for which the phenotype is covalently attached to the genetic material are the most straightforward applications of MetaGPA. However, for other scenarios for which phenotypes are not physically coupled with DNA or RNA, selection will have to be done while preserving the integrity of cells and viral particles to retain the link between phenotype and genotype. Thus, the limited availability of adequate selection processes that preserve such links may restrict the broader applicability of MetaGPA. While some selections can easily be adapted to MetaGPA, others may turn out to be difficult to implement.”

1) Please take into account previous work, describe how these relate to the presented work and clearly indicate the motivation for MetaGPA in the introduction. The relationship between the author's program and previous attempts at GWAS-like analyses in microorganisms are not thoroughly described in the Introduction. For example, the authors do not mention TreeWAS or other similar approaches.

The introduction of the revised manuscript contains now a paragraph that describe microbial GWAS, their achievements and limitations:

“For specific bacteria, determining the genetic basis of phenotypes has be addressed using genome-wide association studies (GWAS) [PMID: 16782339] combined with specific phylogenetic methods to account for the unique population structure of microbes [PMID: 29401456]. Examples of microbial GWAS studies have explored hundreds of isolates to identify genomic elements that are statistically associated with, for example, antibiotic resistance [PMID: 25101644], host specificity [PMID: 23818615] or virulence [PMID: 24717264]. Nonetheless, these studies are limited to known isolates and have not yet been extended to entire complex microbial communities.”

2) The Methods section contains information that presumably relates to the pipeline, but there is no explicit mention of MetaGPA. There are also seemingly pipeline-related methods that are not part of the codebase. The beginning section of the Results is insufficient to explain the MetaGPA method. A reproducible example is a must.

The reviewer brought up a great point regarding what constitutes the core MetaGPA pipeline relative to the additional analysis that we performed to validate MetaGPA. The core MetaGPA pipeline lists the protein domains that are significantly associated with the defined phenotype.

We have now clearly defined what analysis is part of the core MetaGPA pipeline and added the information to the revised version of the result section:

“The association is computed using a computational workflow composed of the core MetaGPA pipeline that defines genetic units associated with the case cohort and further analysis tools to refine these associations.”

We further explain in more detail the core MetaGPA pipeline in the first part of the manuscript.

Following the description in the revised manuscript, we have bundled the core MetaGPA pipeline into a single command line that takes both, the raw sequencing reads from the control and case group, and outputs the protein domains that are significantly associated with the defined phenotype.

We also have re-organized the Method section to reflect the above changes and created a “Computational pipeline” section to provide detailed methodology in the following sections for each step. We also rewrote the legends in Figure 1 and supplementary figure 1 to improve the explanation of the conceptual and detailed workflow.

For a reproducible pipeline example, we provided a dataset for test running at https://github.com/linyc74/MetaGPA. Please see response to question (3) for codebase related improvements.

3) An attempt was made to run the code to confirm reproducibility. The GitHub link contains no information on installation or examples of how to run the code on test or new data. The code itself contains hard-coded file paths that would make it difficult, if not infeasible, to run on another person's machine. There are also dependancies that are unstated. The code must be substantially better documented to be of utility outside this study.

We thank the reviewer’s suggestions on how to improve our code base. We have done substantial modifications such as code refactoring, packaging, and documentation, to improve the quality of our work. More specifically, we have:

1. Created a detailed, clear command-line interface with help (-h option) instructions.

2. Remove all hard-coded file paths, and make the output directory configurable in the command-line interface (-o option).

3. Make system settings (e.g. number of CPUs and max memory) configurable (-t and -m options).

4. Completely refactor the whole code base using object-oriented programming.

5. Detailed documentation in the README.md file to clearly show the following information:

a. Usage

b. Dependencies

c. Input files

d. Output files

e. Docker

6. Dockerize to completely resolve any dependency problem. Docker images could be found at https://hub.docker.com/r/linyc74/metagpa/tags.

The code is available here: https://github.com/linyc74/MetaGPA

With the aforementioned improvements, we hope future users will find it easy to apply MetaGPA to other studies.

4) The reasons for the use of Tet2 and BGT to modify 5mC and 5hmC are unclear. The newly discovered 5-carbamoyloxymethylcytosine is unlikely to be deaminated by A3A

and hence what is the need for Tet2 and BGT? In fact, other than 5-methylcytosine, larger modifications of cytosine would not be substrates for APOBEC3A (PMID: 28472485). The investigators could treat all the DNA with A3A, divide it into two halves, and then use the USER kit to destroy DNA that contains uracils in one-half of the samples. The DNA that survives should contain a cytosine modification that protected it against A3A. Comparison of sequences of the two populations would show that only a small fraction of DNA has survived A3A+User treatment. If the ability of A3A to deaminate 5mC is a major concern, the authors do not articulate it. However, this can be easily remedied by replacing A3A with APOBEC3B-CTD or APOBEC3G-CTD neither of which is good at deaminating 5mC.

The reviewer is right, larger modifications than 5mC are unlikely to be deaminated by A3A and therefore the use of Tet2 and BGT is unnecessary for the discovery of the carbamoyltransferase. This is a posteriori statement and when we designed the experiment, it was unclear whether fully C modified genomes will be present at all : we therefore took a conservative approach to include as many C modifications as possible, even the known 5mC and 5hmC. These genomes serve as a control for MetaGPA, expecting thymidylate synthases and methyltransferases to be included in the list of associated enzymes.

To clarified this point we added to the revised manuscript the underlined sentence:

“Because APOBEC3A also deaminates 5-methyl-2'-deoxycytidine (5mdC) and, to a lesser degree 5-hydroxymethyl-2'-deoxycytidine (5hmdC) (PMID: 33468551), ten-eleven translocation dioxygenase 2 (TET2) and T4 phage ß-glucosyltransferase (T4-BGT) can be used to protect 5mdC and 5hmdC prior to APOBEC3A treatment (Figure 2a). Both modifications are used as internal control for MetaGPA since the enzymes that catalyzed their formation are well characterized. These enzymes, if present in the sample, are therefore expected to be associated with modifications in our MetaGPA framework.”

Another more technical reason we used Tet2 and BGT has to do with protecting repaired DNA ends in the library construction step. Indeed, the sonication of DNA by ultrasound generates 5’ and 3’ extensions that are subsequently repaired by T4 DNA polymerase during the end-repair step (NEB Ultra II library prep kit). If regular dCTP is used, APOBEC3A will deaminate unprotected cytosines introduced in the process of end repair. To avoid this problem, we used d5mCTP instead of regular dCTP, and protected 5mC using TET2 and BGT by generating glucosylated 5hmC, which is not a substrate for APOBEC3A (se Materials and methods). As stated by the reviewer, it is possible to use APOBEC3B-CTD or APOBEC3G-CTD instead of APOBEC3A, however these two deaminases demonstrate strong substrate preference, and not commercially available (Silvas, Tania V., and Celia A. Schiffer. "APOBEC3s: DNA‐editing human cytidine deaminases." Protein Science 28, no. 9 (2019): 1552-1566). We clarified this point in the Materials and methods section:

“The dCTP was replaced with 5mdCTP in the end repair reaction (5mdCTP was used instead of regular dCTP to protect end-repaired fragments from subsequent enzymatic treatment)”.

Finally, while APOBEC3A alone will deaminate 5mC, the resulting T is not removed by USER. Instead, deaminated fragments will be sequenced but the converted sequences are uninterpretable since they are composed of only 3 nucleotides (A, T, and G).

5) The USER kit contains E. coli Exonuclease VIII. This is problematic because this enzyme will excise oxidized pyrimidines in DNA. Thus, the DNA of any organism that routinely modifies its pyrimidines in such a way that it becomes susceptible to ExoVIII, would be eliminated from the case DNA regardless of whether it contains uracils. An AP endonuclease or simple treatment with NaOH and heat would be preferable here.

In this work we used a mix which contains Endonuclease VIII. The reviewer is right to point out that Endonuclease VIII removes oxidative pyrimidines. However, many oxidative pyrimidines recognized by Endonuclease VIII are damaged bases including urea, 5, 6- dihydroxythymine, thymine glycol, 5-hydroxy-5- methylhydantoin, uracil glycol, 6-hydroxy-5, 6-dihydrothymine and methyltartronylurea. Many of these damaged bases are blocking damages for the amplifying polymerase and would not be sequenced even if they were not removed. Importantly, Endonuclease VIII does not remove 5hmC, which works as a handle for many modifications in phages. While AP endonuclease could be used instead of Endonuclease VIII in the USER mix, using a new USER mix (UDG + AP endonuclease) needs some optimization and is not commercially available. Treatment with NaOH is a good suggestion, however, we wanted to avoid harsh reaction conditions and heat treatment which could cause extra damage.

It would be nice, in a subsequent study, to compare contig enrichment with and without the use of Endonuclease VIII using a more promiscuous amplification polymerase. In this version of MetaGPA, we took a conservative approach with the aim to identify the subset of modifications for which the reaction steps can be easily deduced (for experimental validation). Having too many types of modified bases included in the sample cohort can confound the identification of the relevant ones.

6) This methodology is designed to work for cytosine modifiers that near 100% efficient and are not sequence-specific. I do not see how this methodology could isolate genes for sequence-specific cytosine modifiers. If an enzyme strongly prefers a certain sequence motif, say 5'-TpC, then all the cytosines in the VpC sequence context (V is not T) would remain unmodified. These would be deaminated by A3A and virtually all the DNA from that organism would be destroyed by the User kit. The methodology may also have difficulty isolating genes for enzymes that modify only a fraction of the cytosine bases, say 25%. Such DNA would also be destroyed during the A3A+User treatment. Such limitations of the methodology should be carefully examined.

The reviewer is correct: the current procedure only selects DNA for which near 100% of the cytosines are modified. The resulting enzymes are therefore more likely to show no sequence specificity which is a desired property for further applications of the enzymes. We agree with the reviewer that the scope and the rationale for targeting such non-specific enzymes has not been addressed in the manuscript. We have therefore added the following sentences to the revised manuscript:

“the experimental procedure is designed to only protect fully modified DNA, that is DNA for which cytosines are modified, irrespective of sequence contexts. A number of such organisms have previously been found, notably in bacteriophages (Kuo et al., 1968) (Revel and Georgopoulos, 1969). We hypothesize that this design would select for DNA modifying enzymes with little or no sequence specificity.”

7) How does MetaGPA handle phylogenetic resampling (i.e., dealing with the fact that genomes are related)? This is particularly important for microorganisms. It would have been preferable to see the author's method benchmarked in a similar way, if not compared against, previous approaches to similar problems.

The reviewer brought up a great point that genome sequence similarities need to be carefully considered in metagenome analysis.

First, we used metaSPAdes for metagenome assembly for which population contigs (e.i, contigs derived from multiple related species) are generated. In brief, this tool uses more aggressive settings and introduces an additional metagenomics-specific decision rule to optimize assembly within a metagenome. The detailed algorithm and benchmarking against other metagenome assemblers is published [PMID: 28298430].

We also used CD-HIT to further reduce sequence redundancy of assembled contigs. Considering the highly related sequences in microorganisms, we set the sequence identity threshold at 95%.

Finally we used RAxML which is a standard phylogenetic analytical tool [PMID: 24451623]. We enabled rapid bootstrap analysis mode and allowed for automatic determination of sufficient bootstrapping replicates (with option -f a -# autoMRE) for resampling.

8) The manuscript jumps quickly into the main finding of cytosine modification. What are other applications of this technique? How could one incorporate a negative control to quantify FP rates or incorporate other controls?

This is a multi-layered question. Here are the answer to each individual points:

The manuscript jumps quickly into the main finding of cytosine modification.

This is a good point. Before describing the application of MetaGPA on cytosine modification, we have now provided an expanded explanation of the MetaGPA pipeline and associated analysis for the reader to better understand the general concept. These expanded explanations are reflected in the manuscript main text, figure 1 and the figure 1 legend.

What are other applications of this technique?

We have a section in the discussion, describing other possible applications of MetaGPA:

“Our MetaGPA study framework is phenotype-driven and the identification of candidate genes is agnostic to prior annotations. Theoretically, it can be performed with any case/control cohort pair as long as distinct phenotypes can be partitioned through selection of the case cohort. For example, these phenotypes may include but are not limited to DNA modifications, phage sensitivity to chemicals and cell surface determinants such as O-antigen. Nonetheless, this partitioning is crucial for MetaGPA to succeed and therefore requires the development and optimization of the selection process for every new MetGPA application.”

How could one incorporate a negative control to quantify FP rates or incorporate other controls?

This is a great question raised by the reviewer. Negative and positive controls to estimate the selection procedure has been done. More specifically, we had incorporated in the samples a mixture of spiked-in genomic DNA consisting of E. coli (dC), XP12 (5mdC) and T4gt (5hmdC). Recovery of spiked-in modified DNAs were as expected (see manuscript).

A negative control to evaluate MetaGPA would be a microbiome for which no selection has been applied. MetaGPA requires that the control and case groups are derived from the same microbiome, thus the negative control should be a subset of the control group without selection. We therefore use the two replicate experiments done on the control group and perform a MetGPA analysis with one replicate being the control group while the other replicate being the “case” group. None of the 6737 protein domains are found significantly associated with one of the control replicates when compared to the other control replicate, indicating that the FP rate is low. We have incorporated this analysis to the revised manuscript:

“To estimate the false positive rate of association, we used our two control replicate experiments for which no selection has been done. We then perform a MetGPA analysis with one replicate being the control group and the other replicate being the “case” group. from the 6737 protein domains assessed in MetaGPA, none of the domains are found significantly associated with the case group (data not shown).”

9) I am somewhat confused by the discussion about the putative thymidylate synthase homologs. The authors point out that several TS homologs contain a change that is equivalent to a N177D change in E. coli TS. They further note that such a change in the E. coli enzyme results in a change of substrate from dUMP to dCMP. Does this mean that these phage genes code for dCMP methyltransferases, not thymidylate synthases? If so, they may be novel enzymes, as E. coli does not contain a dCMP methyltransferase and unlike TS, the normal methyl donor for cytosine methylation is S-adenosylmethionine. The methyl donor for TS is tetrahydrofolate. It would be useful to characterize these variant TS enzymes for substrate specificity and cofactor requirements.

We agree that the explanation in the manuscript is confusing. We have therefore added the following paragraph to the revised result section of the manuscript:

“The substrate specificity of thymidylate synthases is dictated by the residue at the position 177 (numbering relative to the E. coli thymidylate synthase sequence). Previous studies have demonstrated that changing this residue from asparagine to aspartate can switch the preference of a canonical thymidylate synthase from dUMP to dCMP resulting in the formation of 5mdCMP (Hardy and Nalivaika, 1992) (Graves et al., 1992) (Liu and Santi, 1992). We therefore hypothesised that position 177 should have a high differential conservation score with Asn found conserved in the unmodified contigs and Asp found conserved in modified contigs.”

More detailed information: Thymidylate synthase belongs to a family of proteins that alkylate pyrimidines at C5 by transfer of a single carbon unit from methyltetrahydrofolate (mTHF). The thymidylate synthase family encompasses four enzymatic activities: canonical thymidylate synthase (all domains of life), dUMP hydroxymethylase (Bacillus phage SPO1 and others), dCMP hydroxymethylase (bacteriophage T4 and others), and dCMP methylase (Xanthomonas phage Xp12). Sequences encoding these enzymes share a high degree of sequence and structural similarity. For example, the crystal structures of E. coli thymidylate synthase and it's related T4 dCMP hydroxymethylase have ~1 Å RMSD. The transfer of a carbon from mTHF to C5 results in a product intermediate that contains an exocyclic methylene (i.e. a double bond between the transferred carbon and C5 of the pyrimidine ring). For canonical exocyclic methylene can be reduced to the methyl form by hydride transfer from THF, forming dihydrofolate (DHF). But for some enzymes, the exocyclic methylene undergoes nucleophilic attack from water instead, resulting in the formation of a hydroxymethyl moiety at C5. The substrate specificity of thymidylate synthase homologs is dictated by the residue at the homologous position Asn177 (numbering relative to the E. coli sequence, or 221 for thymidylate synthase of Lactobacillus casei). For the phage T4 dCMP hydroxymethylase, this position is occupied by an aspartate residue. Researchers have demonstrated that changing this residue from asparagine to aspartate can switch the preference of a canonical thymidylate synthase from dUMP to dCMP resulting in the formation of 5mdCMP. Similarly, mutating this homologous position within the T4 dCMP hydroxymethylase from aspartate to asparagine results in an enzyme that prefers dUMP and produces 5hmdUMP. The residues governing formation of a hydroxymethyl versus a methyl group at C5 are less clear. In the case of the viral metagenome sequences encoding a DNA O-carbamoyltransferase identified here, the function of the thymidylate synthase homolog is to provide the hydroxyl "handle" for the attachment of the carbamoyl group.

References:

1. Graves, K. L., Butler, M. M., and Hardy, L. W. (1992) Roles of Cys148 and Asp179 in catalysis by deoxycytidylate hydroxymethylase from bacteriophage T4 examined by site-directed mutagenesis. Biochemistry-us. 31, 10315–10321

2. Liu, L., and Santi, D. V. (1992) Mutation of asparagine 229 to aspartate in thymidylate synthase converts the enzyme to a deoxycytidylate methylase. Biochemistry-us. 31, 5100–5104

3. Hardy, L. W., and Nalivaika, E. (1992) Asn177 in Escherichia coli thymidylate synthase is a major determinant of pyrimidine specificity. Proc National Acad Sci. 89, 9725–9729

4. Agarwalla, S., LaPorte, S., Liu, L., Finer-Moore, J., Stroud, R. M., and Santi, D. V. (1997) A Novel dCMP Methylase by Engineering Thymidylate Synthase †. Biochemistry-us. 36, 15909–15917

5. Liu, L., and Santi, D. V. (1993) Exclusion of 2’-deoxycytidine 5’-monophosphate by asparagine 229 of thymidylate synthase. Biochemistry-us. 32, 9263–9267

10) It is unclear how the investigators arrived at the substrate and co-factor requirements of the newly discovered enzyme. The classic example of a carbamoyltransferase is ornithine carbamoyltransferase which tranfers the carbamoyl moiety to a nitrogen not to an oxygen. It also appears (https://www.brenda-enzymes.org/enzyme.php?ecno=2.1.3.3) that ornithine carbamoyltransferase does not require ATP as a co-factor. With that in mind, the putative new carbamoyltransferase could have modified the N4 of cytosine without the need for ATP. Did the investigators explore this possibility? Overall, a clearer rationale needs to be provided for how the requirements for the enzyme were established.

Those are excellent points. The enzyme requires 3 substrates: ATP, carbamoyl phosphate and 5hmdC. We have explained throughout the manuscript the reason for 5hmdC (see also answer to reviewer 2 point 4). The choice of ATP and carbamoyl phosphate substrates were guided by the enzymatic characterization of TobZ previously published by Parthier et al., 2012. In this publication, the authors used ATP and carbamoyl phosphate as substrates for TobZ, a member of the O-carbamoyltransferase that has been enzymatically and structurally characterized and has strong homology with our enzyme. We have clarified this point in the revised manuscript (underlined text):

“The predicted reaction was tested by enzymatic assays and results showed that every substrate, namely carbamoyl phosphate, ATP, 5hmdC (genomic T4gt DNA was used as substrate in these experiments) and the enzyme were indispensable for the reaction (Supplementary Figure 5c-d). The choice of carbamoyl phosphate and ATP substrates were guided by the enzymatic characterization of TobZ previously published (Parthier et al., 2012).”

Regarding the requirement of ATP, we have performed the reaction with and without ATP. The result showed that ATP is a required cofactor (see Supplementary Figure 5c).

Finally regarding the possibility of a reaction on the N4 of cytosine without the need for ATP, we found that neither the predicted mass of an N-carbamoylated 5hmC, nor an N-carbamoylated cytidine matched the modified base produced in our reconstituted in vitro reaction.

11) Interestingly, E. coli bacteriophage Mu carries N6-carbamoylmethyladenine. Is the enzyme that performs this reaction in any way related to the newly discovered enzyme?

Following the reviewer’s comment, we explored whether N6-carbamoylmethyladenine in Mu phage (Methylcarbamoylase mom) is related to carbamyltranferase.

Methylcarbamoylase mom protein belongs to the GCN5-related N-acetyltranferase (GNAT) family. It is suggested that the Mom enzyme uses a coenzyme-A carrier to donate a formamide moiety to an m6A substrate. This is different from a carbamoyltranferase. We did protein sequence comparison with protein BLAST and HMMER. The results showed no similarity between the enzyme we found and Methylcarbamoylase mom (uniprot p06018). HMMER search against the pfam database (Pfam-A.hmm) for the two protein sequences has no overlapping pfam hit. Together, these results indicate that these two enzymes are mechanistically and evolutionarily unrelated.

12) The only justification for using phage DNA as the source of new enzymes provided by the authors is "phages are known to carry a large diversity of modifications". Another obvious reason for using phage DNA is that few phage genes contain introns. However, this convenience comes with the likely drawback that many DNA modification enzymes in cellular genomes are being missed from their screen. The authors should carefully address this issue.

The reviewer is right that we did not carefully explain the underlying reasons for having the current design of MetaGPA enriching for phages. Indeed, phages are known to carry a large diversity of modifications but importantly (and we have omitted to state this point in the manuscript) phages tend to fully modify their genome to evade any restriction systems. Other organisms tend to have modifications only in certain sequence contexts. Because our experimental procedure is designed to only select for fully modified cytosines (also addressed in essential revision 6), we hypothesized that the phage fraction would be the most likely fraction of the microbiome to contain such DNA.

We have included this point in the revised manuscript (see underlined text):

“We aim at obtaining the phage-enriched fraction of each microbiome. The rationale for this selection is based on the fact that phages are known to carry a large diversity of modifications (Weigele and Raleigh, 2016) and phages have been shown to fully modify their genomes irrespectively of sequencing context which is a prerequisite for our MetaGPA selection.”

13) Figure 3b is unclear. What are the units of the color (scale) bars? What am I supposed to take away from this panel? What is red and what is blue (beyond the minimal description in the figure legend)?

We realized Figure 3b is less informative for the point we want to convey. We have now redone the kmer analysis, this time annotating each kmer in one of the 3 following categories: [1] kmer found in the modified contigs (in orange), [2] kmer found in the unmodified contigs (in blue) and [3] kmer that did not result in assembled contigs (grey). In the revised manuscript, we replace figure 3b accordingly and added the new legend:

“b, Normalized frequency of k-mer in sequencing reads from control (x-axis) compared to case (y-axis) groups. Each dot represents a unique 16-mer and is colored according to the resulting contig category it belonged to. Orange: kmer at modified contigs; blue: kmer at unmodified contigs; grey: unassembled kmers. See c below for definition of modified/unmodified contigs. Subsamples of randomly selected 0.1% of all possible kmers were used for plotting.”

Reviewer #2 (Recommendations for the authors):

1. Page 2- Change "who is out there?" to "what is out there?".

We have now changed this sentence to “what is out there?”

2. Page 6- "unmodified cytosines are deaminated to uracils using the DNA cytidine

deaminase Apolipoprotein B mRNA editing enzyme catalytic polypeptide-like 3A

(APOBEC3A) (Carpenter et al., 2012)"- This is misleading on two counts. Among the APOBEC3 subfamily of cytosine deaminases, APOBEC3A is most efficient at deaminating 5mC, in addition to C. Furthermore, the ability of APOBEC3A to cause cytosine deamination was first demonstrated in HBV genome (PMID: 19169351) and the purification of the enzyme and its biochemical characterized was first reported by a different research group than cited in the manuscript (PMID: 22798497).

We have clarified the fact that APOBEC3A also deaminates 5mC in the revised manuscript (underlining text) and added the two references:

“unmodified cytosines are deaminated to uracils using the DNA cytidine deaminase Apolipoprotein B mRNA editing enzyme catalytic polypeptide-like 3A (APOBEC3A) (PMID: 19169351, PMID: 22798497, Carpenter et al., 2012), and subsequently excised by Uracil-Specific Excision Reagent (USER) (Bitinaite et al., 2007), resulting in fragmented DNA. Because APOBEC3A also deaminates 5-methyl-2'-deoxycytidine (5mdC) and, to a lesser degree 5-hydroxymethyl-2'-deoxycytidine (5hmdC) (PMID: 33468551), ten-eleven translocation dioxygenase 2 (TET2) and T4 phage ß-glucosyltransferase (T4-BGT) can be used to protect 5mdC and 5hmdC prior to APOBEC3A treatment (Figure 2a)”

3. Page 11- "DNA ligase (PF14743.7, PF01068.22), and Cytidine deaminase (PF00383.24) are other domains that have been found in DNA modifying enzymes (Subramanya et al., 1996) (Bhattacharya et al., 1994)"- Did you mean DNA-cytosine deaminase? A cytidine deaminase would convert the ribonucleoside cytidine to uridine.

We would like to thank the reviewer for spotting this typo. We meant the protein family PF00383.24 (cytidine and deoxycytidylate deaminase zinc-binding region). We corrected it in the manuscript: For example, a subset of enzymes containing the thymidylate synthase domain (PF00303.20) have been shown to produce hydroxymethylpyrimidines (Neuhard et al., 1980). DNA ligase (PF14743.7, PF01068.22), and Cytidine/deoxycytidylate deaminase (PF00383.24) are other domains that have been found in DNA modifying enzymes (Subramanya et al., 1996) (Bhattacharya et al., 1994).

4. Page 18, last paragraph- The comparison of T4 genome with the contig in which carbamoyltransferase gene is found is confusing. Why is the occurrence of dCMP hydroxymethylase and β-glucosyltransferase in T4 genome "resembles" the occurrence thymidylate synthase and carbamoyltransferase in the contig? Phages have very compact genomes and tend to aggregate genes with related functions. Furthermore, dCMP hydroxymethylase is presumably an oxidase and not a transferase. If this is what led to the prediction of the reaction of the newly discovered carbamoyltransferase, it was a really inspired guess- I say that in all sincerity.

We have not properly explained how the T4 phage gene configuration led us to hypothesize the mechanism of action of the carbamoyltransferase. We have therefore added to the revised version of the manuscript the following paragraph:

“The co-occurrence of carbamoyltransferase and thymidylate synthase homologs specifically in modified contigs (Figure 4c) resembles the arrangement of the T4 phages for which genes coding for the dCMP hydroxymethylase and β-glucosyltransferase co-occur on the genome (Miller et al., 2003). The T4 dCMP hydroxymethylase is homologous to thymidylate synthase and transfers a carbon from methyltetrahydrofolate (mTHF) to C5 of the pyrimidine ring producing an exocyclic methylene in the active site of the enzyme (Graves et al. 1992). However, unlike thymidylate synthase, the methylene intermediate undergoes nucelophilic attack by water producing a hydroxymethyl group. Following incorporation of 5hmC into DNA during replication, T4 β-glucosyltransferase transfers a glucose to the hydroxyl moiety of 5hmC. Thus, the pairing of a carbamoyltransferase with dCMP hydroxymethylase led us to hypothesize a novel form of DNA modification, in which the carbamoyltransferase catalyzes the transfer of a carbamoyl group to the nucleophilic hydroxyl acceptor group of 5hmdC producing 5-carbamoyloxymethylcytosine (5cmdC) (Figure 6a).”

https://doi.org/10.7554/eLife.70021.sa2

Article and author information

Author details

  1. Weiwei Yang

    New England Biolabs, Ipswich, United States
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – original draft, Writing – review and editing
    Contributed equally with
    Yu-Cheng Lin
    Competing interests
    is an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8836-9018

  2. Yu-Cheng Lin

    1. New England Biolabs, Ipswich, United States
    2. School of Dentistry, National Yang Ming Chiao Tung University, Taipei, Taiwan
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review and editing
    Contributed equally with
    Weiwei Yang
    Competing interests
    was an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4787-2565

  3. William Johnson

    New England Biolabs, Ipswich, United States
    Present address
    Department of Molecular Biology and Microbiology, Tufts University, School of Medicine, Boston, United States
    Contribution
    Investigation, Methodology
    Competing interests
    was an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents

  4. Nan Dai

    New England Biolabs, Ipswich, United States
    Contribution
    Investigation, Methodology, Validation, Writing – review and editing
    Competing interests
    is an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents

  5. Romualdas Vaisvila

    New England Biolabs, Ipswich, United States
    Contribution
    Conceptualization, Investigation, Methodology, Writing – review and editing
    Competing interests
    is an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents

  6. Peter Weigele

    New England Biolabs, Ipswich, United States
    Contribution
    Conceptualization, Supervision, Writing – review and editing
    Competing interests
    is an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents

  7. Yan-Jiun Lee

    New England Biolabs, Ipswich, United States
    Contribution
    Formal analysis, Investigation, Methodology, Writing – review and editing
    Competing interests
    is an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents

  8. Ivan R Corrêa Jr

    New England Biolabs, Ipswich, United States
    Contribution
    Conceptualization, Data curation, Supervision, Writing – review and editing
    Competing interests
    is an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3169-6878

  9. Ira Schildkraut

    New England Biolabs, Ipswich, United States
    Contribution
    Investigation, Methodology, Validation
    Competing interests
    is an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents

  10. Laurence Ettwiller

    New England Biolabs, Ipswich, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Supervision, Validation, Writing – original draft, Writing – review and editing
    For correspondence
    laurence.ettwiller@gmail.com
    Competing interests
    is an employee of New England Biolabs Inc, a manufacturer of restriction enzymes and molecular reagents
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3957-6539

Funding

New England Biolabs (no data)

  • Weiwei Yang
  • Yu-Cheng Lin
  • William Johnson
  • Nan Dai
  • Romualdas Vaisvila
  • Peter Weigele
  • Yan-Jiun Lee
  • Ivan R Corrêa Jr
  • Ira Schildkraut
  • Laurence Ettwiller

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors would like to thank Barry Cohen and the junior team (Grace, Alina, Camille, and Jeanne) for environment sample collections, Lana Saleh for useful discussions and providing the RNA oligos, Vladimir Potapov, Zhiyi Sun, and Tamas Vincze for analysis and IT support and the sequencing core. John Buswell for providing the Pyrrolo dC adaptor. This work was supported by New England Biolabs.

Senior Editor

  1. Gisela Storz, National Institute of Child Health and Human Development, United States

Reviewing Editor

  1. María Mercedes Zambrano, CorpoGen, Colombia

Publication history

  1. Preprint posted: March 23, 2021 (view preprint)
  2. Received: May 4, 2021
  3. Accepted: November 5, 2021
  4. Accepted Manuscript published: November 8, 2021 (version 1)
  5. Version of Record published: December 14, 2021 (version 2)

Copyright

© 2021, Yang et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,395
    Page views
  • 182
    Downloads
  • 2
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Weiwei Yang
  2. Yu-Cheng Lin
  3. William Johnson
  4. Nan Dai
  5. Romualdas Vaisvila
  6. Peter Weigele
  7. Yan-Jiun Lee
  8. Ivan R Corrêa Jr
  9. Ira Schildkraut
  10. Laurence Ettwiller
(2021)
A genome-phenome association study in native microbiomes identifies a mechanism for cytosine modification in DNA and RNA
eLife 10:e70021.
https://doi.org/10.7554/eLife.70021

Categories and tags

Further reading

    1. Genetics and Genomics
    2. Neuroscience

    Tau polarizes an aging transcriptional signature to excitatory neurons and glia

    Timothy Wu, Jennifer M Deger ... Joshua M Shulman
    Research Article Updated

    Aging is a major risk factor for Alzheimer’s disease (AD), and cell-type vulnerability underlies its characteristic clinical manifestations. We have performed longitudinal, single-cell RNA-sequencing in Drosophila with pan-neuronal expression of human tau, which forms AD neurofibrillary tangle pathology. Whereas tau- and aging-induced gene expression strongly overlap (93%), they differ in the affected cell types. In contrast to the broad impact of aging, tau-triggered changes are strongly polarized to excitatory neurons and glia. Further, tau can either activate or suppress innate immune gene expression signatures in a cell-type-specific manner. Integration of cellular abundance and gene expression pinpoints nuclear factor kappa B signaling in neurons as a marker for cellular vulnerability. We also highlight the conservation of cell-type-specific transcriptional patterns between Drosophila and human postmortem brain tissue. Overall, our results create a resource for dissection of dynamic, age-dependent gene expression changes at cellular resolution in a genetically tractable model of tauopathy.

    1. Evolutionary Biology
    2. Genetics and Genomics

    Native American genetic ancestry and pigmentation allele contributions to skin color in a caribbean population

    Khai C Ang, Victor A Canfield ... Keith C Cheng
    Research Article

    Our interest in the genetic basis of skin color variation between populations led us to seek a Native American population with genetically African admixture but low frequency of European light skin alleles. Analysis of 458 genomes from individuals residing in the Kalinago territory of the Commonwealth of Dominica showed approximately 55% Native American, 32% African, and 12% European genetic ancestry, the highest Native American genetic ancestry among Caribbean populations to date. Skin pigmentation ranged from 20 to 80 melanin units, averaging 46. Three albino individuals were determined to be homozygous for a causative multi-nucleotide polymorphism OCA2NW273KV contained within a haplotype of African origin; its allele frequency was 0.03 and single allele effect size was -8 melanin units. Derived allele frequencies of SLC24A5A111T and SLC45A2L374F were 0.14 and 0.06, with single allele effect sizes of -6 and -4, respectively. Native American genetic ancestry by itself reduced pigmentation by more than 20 melanin units (range 24 - 29). The responsible hypopigmenting genetic variants remain to be identified, since none of the published polymorphisms predicted in prior literature to affect skin color in Native Americans caused detectable hypopigmentation in the Kalinago.

    1. Genetics and Genomics

    Unraveling the influences of sequence and position on yeast uORF activity using massively parallel reporter systems and machine learning

    Gemma E May, Christina Akirtava ... Joel McManus
    Research Article Updated

    Upstream open-reading frames (uORFs) are potent cis-acting regulators of mRNA translation and nonsense-mediated decay (NMD). While both AUG- and non-AUG initiated uORFs are ubiquitous in ribosome profiling studies, few uORFs have been experimentally tested. Consequently, the relative influences of sequence, structural, and positional features on uORF activity have not been determined. We quantified thousands of yeast uORFs using massively parallel reporter assays in wildtype and ∆upf1 yeast. While nearly all AUG uORFs were robust repressors, most non-AUG uORFs had relatively weak impacts on expression. Machine learning regression modeling revealed that both uORF sequences and locations within transcript leaders predict their effect on gene expression. Indeed, alternative transcription start sites highly influenced uORF activity. These results define the scope of natural uORF activity, identify features associated with translational repression and NMD, and suggest that the locations of uORFs in transcript leaders are nearly as predictive as uORF sequences.