Retroelement-guided protein diversification abounds in vast lineages of Bacteria and Archaea.

Major radiations of enigmatic Bacteria and Archaea with large inventories of uncharacterized proteins are a striking feature of the Tree of Life1-5. The processes that led to functional diversity in these lineages, which may contribute to a host-dependent lifestyle, are poorly understood. Here, we show that diversity-generating retroelements (DGRs), which guide site-specific protein hypervariability6-8, are prominent features of genomically reduced organisms from the bacterial candidate phyla radiation (CPR) and as yet uncultivated phyla belonging to the DPANN (Diapherotrites, Parvarchaeota, Aenigmarchaeota, Nanoarchaeota and Nanohaloarchaea) archaeal superphylum. From reconstructed genomes we have defined monophyletic bacterial and archaeal DGR lineages that expand the known DGR range by 120% and reveal a history of horizontal retroelement transfer. Retroelement-guided diversification is further shown to be active in current CPR and DPANN populations, with an assortment of protein targets potentially involved in attachment, defence and regulation. Based on observations of DGR abundance, function and evolutionary history, we find that targeted protein diversification is a pronounced trait of CPR and DPANN phyla compared to other bacterial and archaeal phyla. This diversification mechanism may provide CPR and DPANN organisms with a versatile tool that could be used for adaptation to a dynamic, host-dependent existence.

Diverse environments host archaea and bacteria that define major lineages of predominantly uncultivated organisms: the archaeal superphylum comprising the Diapherotrites, Parvarchaeota, Aenigmarchaeota, Nanoarchaeota, and Nanohaloarchaea (DPANN – which includes the recently discovered Woesearchaeota and Pacearchaeota[1,2] phyla as well) and the bacterial candidate phyla radiation (CPR)[3,4,5]. Members of these lineages have consistently been reported to have small genomes (~0.5 – 1.5 Mbp) and some have been shown to have ultra-small cells[9,10,11]. Most DPANN and CPR genomes are missing biosynthetic pathways considered vital for autonomous growth, which points to a host-dependent lifestyle[12,13]. Despite genomic insights, little is known of the mechanisms for genetic diversification that drive either adaptation to environmental stress[14] (i.e. nutrient or energy limitation) or interactions with neighboring cells. Based on the recent identification of DGRs in two single-cell DPANN partial genomes from a subsurface environment[15], and the established role of DGRs in host-dependent bacteria and their viruses[7], we sought to address the hypothesis that organisms with minimal genomes and biosynthetic deficiencies (i.e. for nucleotides, lipids, and amino acids) belonging to CPR and DPANN phyla, are candidates for DGR utility. Diversification mechanisms in CPR and DPANN are not established, but are important, given that these radiations appear to comprise a major fraction of microbial life[3,4,5].

Among known biological mechanisms for genetic diversification, DGRs are capable of exploring the highest ceiling on coding sequence variability[16]. These retroelements are unique in that they promote rapid and targeted mutation of specific genomic loci using a reverse transcriptase (RT) that is predicted to be error-prone. They occur in genomes of bacteriophage, and in bacteria whose lifestyles are typified by parasitism, pathogenesis, and intraspecific competition[6,17,18]. Recent studies have established various forms of evidence that DGRs offer selective advantages to the host genomes that encode them[7,8,17]. First, their capacity for targeted mutation allows for diversification in a hypervariable coding scaffold, without altering conserved sequence regions that support fold stability[16]. Second, the variable proteins whose structures have been experimentally assessed to date appear to contain ligand-binding folds[16,19,20]. This unity in the structures diversified within an array of different genomes points to a common advantage in expanding ligand specificity for different viral and cellular binding proteins. Finally, DGRs appear to be conserved in multiple strains of both Legionella pneumophila and Treponema denticola and can occur within conjugative elements[16,17,21]; negative selective pressures are likely to preclude exchange and retention of DGRs across ancestral networks.

The DGR mechanism of mutagenic homing (Figure 1A) deploys RT to target a variable protein for diversification through an RNA intermediate[6,7,22]. Genes encoding the variable protein contain a variable repeat (VR) in close proximity to an invariant template repeat (TR); RT-induced mutation of TR-RNA adenines and cDNA replacement of VR leads to an extraordinary potential for diversification (>1020 amino acid variants)[16]. To determine whether DGRs occur in the genomes of CPR and DPANN organisms and to assess their evolutionary importance, we capitalized on the availability of a massive metagenomic dataset[2,3,5], containing numerous new draft and complete genomes. Our analyses targeted sequences of size-fractionated bacteria and archaea from an aquifer known to harbor diverse CPR and DPANN phyla. The cells were captured onto sequential 1.2 μm, 0.2 μm and 0.1 μm filters[2,3,5]. In total, we analyzed 30 metagenomic datasets, as well as six metatranscriptomes. Remarkably, we identified 1,136 non-redundant sequences that encode essential features of a DGR (i.e., an RT gene and a VR/TR pair; Figure 1A), approximately tripling the total number of DGRs that have been previously identified from more than 300 bacterial, archaeal, and viral genomes[17].

Figure 1

Prevalence of DGRs identified in groundwater metagenomes

a, Schematic of a genomic diversity-generating retroelement cassette and the mutagenic retrohoming mechanism. b, Clusters of reverse transcriptase protein sequences at >70% global pairwise alignment displayed by filter size, with the numbers of unique/shared RT clusters in parentheses. c, Distribution of DGR occurrence in reconstructed genomes, given as a fraction of 1 Mbp for each discrete interval. d, Incidence of DGRs in the archaeal DPANN superphyla, Parcubacteria (OD1) and Microgenomates (OP11), and other CPR phyla.

We determined the number of distinct DGR sequence types and examined the frequency of these sequences in the genomes of cells separated by size filtration. DGR-like RTs were grouped based on >70% amino-acid identity, generating 699 protein clusters. Only 23 clusters are unique to samples from the largest size fraction (Figure 1B); conversely, 75 are unique to the mid-size fraction, and 315 to the small size fraction. Furthermore, we identified 63 distinctive DGRs affiliated with genomes that were previously linked to cells examined by cryo-electron tomography, which passed through a 0.2 μm filter[11]. Importantly, out of the 542 genomes in which we identified DGRs, only 19 are linked to non-CPR/DPANN phyla highlighting a low frequency in organisms predicted to have larger cell sizes. Based on these results we conclude that DGRs are enriched in the genomes of ultra-small bacterial and archaeal cells, providing an entirely new niche association for this adaptive mechanism.

To estimate the genomic incidence of targeted protein diversification in these samples, we focused on DGR-containing sequences in 542 high-quality draft genomes; i.e., reconstructed bacterial and archaeal genomes, including 530 containing >75% of a set of universal single copy genes[2,23]. Notably, a prevalence of DGRs is observed in reconstructed genomes of 0.5 – 1 Mbp (Figure 1C), consistent with the hypothesis that DGRs are common to genomically-reduced phyla. Moreover, these retroelements are overrepresented (found in > 25% of genomes) in the genomes of DPANN archaea (i.e. Pacearchaeota and Woesearchaeota), Parcubacteria (i.e. Campbellbacteria, Falkowbacteria, Giovannonibateria, Jorgensenbacteria, Magasanikbacteria, Moranbacteria, and Yanofskibacteria), and Microgenomates (i.e. Beckwithbacteria and Levybacteria), while likewise abundant in other CPR phyla (Figure 1D) – providing an extraordinary association of DGRs with CPR and DPANN phyla. Among the genomes that contain a DGR, 26% of archaeal genomes and 19% of bacterial genomes encode multiple distinct cassettes. The high incidence of DGRs encoded in reduced genomes raises the question of whether mutagenic homing is a broadly available evolutionary tool among CPR and DPANN phyla (i.e. facilitating genetic adaptation to various selective pressures), and whether DGRs remain active in extant populations.

To determine whether variable proteins have been recently diversified, we assessed the pattern of heterogeneity in reads mapped to VR coding regions of DGR protein targets. The VRs were aligned with corresponding bases in cognate TRs to determine the proportion of adenine-specific mismatches leading to non-synonymous substitutions. Here, a stringent approach excluded sequences with non-synonymous changes in the VR loci that were not aligned with a TR adenine position. Overall, we detected 132 DGR sequences with pronounced levels of adenine-specific mutagenesis, unlikely to arise due to stochastic mutation (resulting in >3 non-synonymous substitutions per VR). This finding provides evidence that proteins in a subset of the recovered genomes were undergoing diversification leading up to sampling. These results are consistent with prior evidence that DGR-RT misincorporates bases in VR that correspond to TR adenines[7] (Supplementary Figure S1). Closer examination of DGRs in these genomes identified 34 sequences capable of forming stem-loops in loci proximal to their VRs (Supplementary Table S1); such cis-acting elements are important for DGR function[24]. We then searched for analogs of Avd, a low molecular weight accessory protein (14.7 kDa) required for mutagenic homing in other DGR systems[7]. This analysis uncovered a group of 72 DGRs encoded alongside conserved, small proteins (average 9.57 kDa) predicted to function analogously to Avd. Taken together, these findings suggest that the intact DGRs we identified have a capacity for driving adenine mutagenesis.

To further determine if these DGRs were active contemporaneous with sampling, we analyzed six metatranscriptomes, identifying DGRs that aligned with at least one metatranscriptomic read. DGRs from Parcubacteria and Woesearchaeota appear to be disproportionately expressed, with between 13 and 94 metatranscriptomic reads that were mapped to one or more DGR feature(s) (Supplementary Table S2). In contrast to the high transcriptional levels observed for DGRs, few to no reads were mapped to protein-coding genes as putative mRNA; the majority of total metatranscriptomic reads were aligned to rRNA, tRNA, and tmRNA. In the case of each DGR cassette, reads appear to map primarily, but not exclusively, to TR regions. Strikingly, in some instances, metatranscriptomic read mapping was comparable (i.e., within 10x) between TR-RNA and other highly expressed regions (i.e. tmRNA) on the same contig (Supplementary Figure S2). These observations point to DGR expression of a stable RNA molecule in organisms affiliated with both CPR and DPANN, in a manner that is consistent with TR-specific expression in DGR systems of Trichodesmium erythraeum[18]. Taken together with evidence for adenine-mutagenesis, these findings point to a functional subset of DGRs in groundwater-associated genomes, and by their established mechanism, to protein diversification in CPR and DPANN phyla.

We next analyzed the phylogeny of RT proteins to determine whether the newly identified DGRs are closely related to known DGRs, or if the elements described in this study represent novel lineages. We find that the majority of newly discovered bacterial RTs appear only distantly related to those from known bacterial DGRs. DGR-RTs from CPR genomes almost exclusively form novel lineages (Figure 2A). Moreover, RTs identified in genomes of DPANN phyla form a monophyletic archaeal clade along with previously described DGR-RTs found in single-cell DPANN genomes[15]. Whereas representatives from archaea and CPR form separate DGR-RT lineages, we observe paraphyletic patterns at higher taxonomic resolution (Figure 2A). Notably, closely related RTs are shared between Woesearchaeota and Pacearchaetoa, and separately, groups of similar RTs appear to link Parcubacteria and Microgenomates with members of Berkelbacteria, Peregrinibacteria, and Saccharibacteria. Moreover, these DGR groups have highly similar TR nucleotide sequences, offering independent evidence of exchange by horizontal gene transfer (HGT).

Figure 2

Phylogeny of diversity-generating retroelements and radiation of novel lineages

a, Phylogeny and taxonomic association of diversifiers. Sequences obtained in this study are highlighted in black on the outer edge of the tree. Shaded slices indicate either candidate superphyla comprising groundwater organisms, or bacterial phyla and phage with previously sequenced RTs. The tree was constructed with 346 sequences and 261 alignment positions. Archaeal sequences from the DPANN superphylum are indicated in either dark red (Pacearchaeota), or light red (Woesearchaeota) shaded slices. Paraphyletic groups of species with closely related RTs are indicated with the associated range of pairwise TR sequence similarity. Symbols (hexagons and triangles) highlight DGRs identified in this study, which are found in close proximity to putative prophage or mobile elements (i.e. transposons or conjugative elements). Diversifiers that have been previously examined are: BPP, a Bordetella pertussis phage; T. denticola, Treponema denticola; L. pneumophila, Legionella pneumophila strain Corby; and ANMV-1, a marine virus of an uncultivated archaeal host. White circles indicate bootstrap values >70% for basal nodes. The scale indicates amino acid substitutions per site. b, Phylogenomic tree of CPR organisms, highlighting major lineages that contain at least one DGR. The phylogenomic tree was originally presented by Hug et al. (2016).

To further investigate whether RTs were subject to HGT amongst CPR and DPANN organisms, we inspected genes in proximity to DGRs (i.e., +/− 10 kbp of an RT gene) for characteristics of prophage (e.g., viral proteins, terminase, integrase), or other mobile elements, such as transposons or conjugative elements. This search revealed numerous DGRs that occur in close proximity to at least one transposase gene, whereas only a single DGR could be linked to a recognizable prophage-like region (Figure 2A). The apparent capacity for horizontal transfer on mobile elements suggests that DGRs offer selective advantages in these bacteria and archaea. We additionally sought to examine monophyletic lineages from a previously constructed phylogenomic tree of CPR representatives[4] for examples of conserved DGR. Groups of related DGR RTs, and separate VP groups, appear to be conserved across the Yanofskybacteria clade of CPR phyla (Supplementary Figure S3). Retention across a broad evolutionary distance suggests that DGRs offer advantages to the genomes that encode them. It is also worth noting that, since the CPR phylogenomic tree was constructed including partial genomes, we likely have a limited view of DGR retention for certain representatives (Supplementary Figure S4).

Prior to this study, DGRs were known to occur in bacterial genomes from many phyla[7,17] belonging to various microbiomes, but were rare in the genomes of most cultured isolates[25,26,27] Our findings are in stark contrast, revealing an extraordinary radiation of novel DGR clades from newly described bacteria belonging mostly to candidate phyla (Figure 2A). Further, these elements are enriched in numerous genomes of CPR bacteria (especially Parcubacteria at >37% incidence). Importantly, DGRs were discovered across most of the major CPR lineages (Figure 2B), and not solely within a closely related subset of organisms. The findings presented herein expand upon prior evidence of two archaeal genomes that encode DGRs, revealing a multitude of diversifiers in other representatives of DPANN archaea (especially Pacearchaeota at >80% incidence, and Woesearchaeota at >28% incidence). When compared with other microbial DGRs from an array of different environments, the groundwater-associated representatives in this study account for 55% of cumulative branch length, or apparent diversity (i.e., substitutions per site), on a non-redundant tree of all representatives (Figure 2A). Remarkably, while DGRs account for an estimated 3% of RTs in previously sequenced bacterial genomes[28], they make up 57% of recognizable RTs in the analyzed metagenomes. These findings highlight DGRs as a prominent genetic feature for CPR and DPANN phyla, and reveal an exemplary biome wherein DGRs have evolved to become prominent and active agents of adaptation.

The phylogenetic distribution, horizontal exchange, and retained function of DGRs discovered here, considered alongside their established roles in cellular interaction, raise the question of whether these retroelements are recruited preferentially for adaptive evolution of proteins that enable symbiosis. To address this question we performed functional predictions for DGR variable proteins. To estimate their functional richness, we clustered variable proteins, and separately extracted VR domains (i.e., with >30% intra-cluster similarity; Supplementary Figure S5). We identified 396 variable protein clusters and 284 VR-domain clusters. The most common protein domain annotations include a conserved domain of unknown function (DUF1566), AAA+ related ATPase, and C-type lectin domain (Supplementary Figure S5). Notably, the majority of VRs are located in the C-terminus of their variable protein, which is consistent with VR placement in previously characterized DGRs[7,16]. Moreover, C-terminal VR domains appear to be predominately localized to C-type lectin (CLec-like) domains and DUF1566 domains, which were recently shown to have CLec-folds[29]. Amongst an assortment of DUF1566 domains linked to putative transmembrane proteins, we identified putative pilin structures, lipoproteins, fimbrial protein FimH, and a rearrangement hotspot (rhs) toxin, suggesting broad involvement of these proteins in cell attachment and defense (Supplementary Table S3). While most DGR variable proteins are putative single domain proteins, less than 350 amino acids in length (Supplementary Figure S6), an unexpected array of mutidomain architectures is also observed.

Through further analysis of the draft genomes from candidate phyla, we identified additional clusters of variable proteins whose domain architectures are associated with transcriptional regulation, cell-cell attachment, and signal transduction (Figure 3A). Analyses of homology identified variable proteins, which are common to multiple CPR lineages, containing AAA+ ATPase modules, pilin-like N-terminal regions fused to C-terminal CLec-like ligand-binding domains, and putative kinase-like regions. Moreover, the ATPase domains in CPR variable proteins each belong to the AAA-5 subgroup of eukaryotic-like midasins (Supplementary Table S4). Distinct classes of these variable proteins are associated with genomes representing a range of candidate phyla (Figure 3B). Whether chaperones, cell-cell attachment proteins, or signaling proteins, each of the otherwise distinct variable protein classes are likely to function as ligand-binding receptors. The general role of ligand binding has been described as a unifying attribute of both signaling and regulatory proteins serving core cellular functions in prokaryotes[30]. Our findings of HGT and adenine-specific mutations of the VR point to the selective advantage of DGRs in CPR organisms, which is perhaps related to the utility in diversifying modular ligand-binding domains – driving expansion of substrate specificity for regulation, signaling, and attachment.

Figure 3

Putative functional classes of DGR variable proteins. a, Functional annotations for variable proteins identified in at least two distinct candidate phyla. Conserved domains and features based on structural homology are shown. Phyre2 confidence values are given for predictions and their closest known structures. b, Distribution of common variable proteins found in candidate phyla, labeled i – iv, based on the examples in (a).

Presumably then, variable protein genes that offer selective advantages to their genomes would be conserved in both DGR and non-DGR loci. To address this hypothesis, we identified examples of variable protein paralogues occurring both within DGR cassettes and in the absence of proximal DGR features (Supplementary Figure S7 and Supplementary Figure S8). Finding homologous variable protein genes in both DGR and non-DGR loci of a genome suggests diversification might be followed by preservation of a particular variant gene. In addition to potential advantages linked to specific cellular functions, DGRs also offer more general benefit in conferring genetic variability for minimal genomes. Targeted and localized mutagenesis could provide benefit to organisms with minimal genomes that cannot otherwise accommodate extensive variant repertoires.

Myriad selective pressures on CPR or DPANN organisms can impose a need for genetic hypervariability, and interactions with neighboring cells are likely to act as such evolutionary pressures. Based on our results, we conclude that DGRs are prevalent in ultra-small, genomically-reduced cells belonging to both CPR bacteria and DPANN archaea. This prevalence provides an indication of selective pressures that transcend the ancient divergence of CPR bacteria and DPANN archaea. As explanation, we hypothesize an enhanced utility of DGR-mediated diversification that emerges from selective pressure, balancing minimal genome size against the need for dynamic response to manage host association. Furthermore, the capacity of DGRs for accelerated protein evolution suggests a need on the part of some CPR and DPANN to rapidly evolve their symbiotic associations, which in turn suggests they may sometimes exploit intercellular associations to receive greater benefit than their host – perhaps shifting between mutualism, predation, and parasitism.

Methods

Study site and sampling

This study used data from several previously-described samples[2,3,5]. In brief, sampling was conducted within an unconfined aquifer at the Rifle Integrated Field Research Challenge (IFRC) site, which is adjacent to the Colorado River, near Rifle, Colorado, USA (39° 31′ 44.69″ N, 107° 46′ 19.71″ W). Groundwater samples were collected from three different field experiments: 6 sampling time points across the duration of acetate amendment, A-F; 4 sampling time points across the duration oxygen injection A-D; and 2 sampling time points from natural high and low oxygen conditions in the groundwater, driven by fluctuations in the water table at the site. Aquifer well CD-01 was monitored as part of a 95-day acetate amendment experiment during which acetate was added to the aquifer between August 25th and December 12th, 2011 as previously described[2]. Following this experiment, aquifer well CD-01 was monitored as part of a 132-day oxygen injection experiment where oxygen-saturated water was injected into the aquifer from August 2nd 2012 to December 12th 2012.

Aquifer well FP-101 was sampled during two specific time points characterized by high and low oxygen in the groundwater. All groundwater samples were collected from 5 m below the ground surface and cells were collected on serial 1.2, 0.2 and 0.1μm filters (Supor disc filters; Pall Corporation) towards differential sampling of small-celled organisms. Following groundwater sampling, filters were immediately frozen prior to DNA extraction, either on dry ice, or in liquid nitrogen.

Metagenomic and metatranscriptomic sequencing

As described previously[2,3,5], genomic DNA of groundwater organisms was extracted from ~1.5 g filter samples, using a PowerSoil DNA Isolation Kit (MO-BIO). Filters were cut into strips, which were then vortexed in PowerBead solution, before and after an interval of flash freezing and thawing. Following thawing, the solution was incubated for 30 minutes at 65°C while shaking. DNA was then eluted and concentrated by sodium acetate and ethanol precipitation in glycogen, followed by resuspension in 50 μl of PowerSoil elution buffer. Sequencing was performed at the Joint Genome Institute, using the Illumina HiSeq 2000 platform to generate 2×150 paired-end reads.

As reported previously[3], prior to RNA extraction, genomic DNA removal and cleaning was done using RNase-Free DNase Set kit (Qiagen) and Mini RNeasy kit (Qiagen). Next, RNA extractions were performed using Invitrogen TRIzol Reagent (Invitrogen). Sample aliquots were analyzed before sequencing using the Agilent 2100 Bioanalyzer to assess the quality of purified RNA. The cDNA sequence preparation was performed using a SOLiD Total RNA-Seq kit (Applied Biosystems). Samples were sequenced on the SOLiD 5500XL platform, at the DOE Environmental Molecular Sciences Laboratory, a facility of the Pacific Northwest National Laboratory. Initial genome sequence mapping was conducted using LifeScope software (version 2.5; SOLiD) with default parameters; additional metatranscriptomic read mapping details are given below.

Assembly, annotation, and binning of metagenomic sequences

This study involved the analysis of sequences that were previously preprocessed, assembled, and binned[2,3,5]. Briefly, sequencing reads were filtered for quality using Sickle software version 1.33 (https://github.com/najoshi/sickle), with default parameters. Next, sequences were assembled with IDBA_UD, using default parameters for paired-end reads[31]. Only assembled scaffolds exceeding 5kb were used for downstream annotation and binning steps. Open reading frame (ORF) annotation was performed using Prodigal[32] with the metagenome mode setting. In the present study, functional annotations for genes were determined using hmmsearch v3.1[33] against an in-house HMM database constructed based on KEGG orthology, while sequences for DGR variable proteins and putative reverse transcriptases were also compared with the Uniprot database using pHMMER[33]. Here, variable proteins were also analyzed for homology to known protein structures using Phyre2[34]. Read mapping was conducted using Bowtie2[35]. Previously, scaffolds were binned to specific organisms by using coverage across the samples, phylogenetic identity, and GC content, both automatically with the ABAWACA algorithm[3] and manually using ggKbase (http://ggkbase.berkeley.edu/). ABAWACA is an algorithm that generates genome bins based on scaffold taxonomic affiliations, time-series abundance patterns, and nucleotide frequencies (https://github.com/CK7/abawaca). Genome bins generated by ABAWACA were manually inspected within ggKbase. Binning purity was confirmed using an Emergent Self-Organizing Map (ESOM)[36,37]. Each bin was previously assigned a genome phylogeny if it met the following criteria: 1) determined to be high quality[2,3,5] (i.e. scaffolds could not be further separated into distinct bins); 2) contains at least 75% of the conserved, single-copy genes found broadly across bacterial genomes, or separately across archaeal genomes[38]. For draft genomes that contain less than 75% of conserved single copy genes (12 genomes examined in this study), taxonomy was determined from ribosomal protein and 16S rRNA phylogeny as previously described[2].

Identification, annotation, and clustering analysis of DGR Features

The following methods were carried out in this study. Several genomic features that are encoded by all DGRs – namely a reverse transcriptase gene and VR/TR pairs – can be used as diagnostic indicators for in silico identification of these retroelements[7,28]. To this end, we used a consensus sequence of aligned DGR-RT protein sequences from Treponema denticola, Bordetella pertussis phage, Legionella pneumophila, archaeal virus ANMV-1, and uncultivated nanoarchaeota, to conduct a tblastn search for RT-like hits in assembled scaffolds from groundwater metagenomes (Supplementary Figure S9). Next, a custom python script was used to identify near-repeats within 10 kb of the putative RT gene, by applying a sliding window (200 bp windows; 50 bp step) blastall search with the following parameters: -word_size 8 -reward 1 -penalty -1 -evalue 1e-5 -gapopen 6 -gapextend 6 -perc_identity 50. This output was filtered for near repeats: >60bp pairs, which contain >5 adenine-specific mismatches and no more than one non-adenine mismatch (i.e. putative VR/TR pairs). Given that adenines of AAY codons are selectively targeted by DGRs[16], ORFs were only identified as variable protein genes where the majority of TR adenines correspond to the first two positions of a given codon in VR.

Putative stem-loop encoding regions (100bp downstream of VR) were extracted from DGR cassette nucleotide sequences, and regions capable of forming a stem-loop were identified using MFold[39]. Translated reverse transcriptase sequences were clustered using CD-HIT[40] with a global alignment (identity) threshold of 70%. Variable protein sequences were clustered using H-CD-HIT with three iterative rounds of clustering and identity thresholds at 90%, 60%, and 30%. VR domains, including 50 flanking amino acids, were extracted from the variable protein sequences and separately clustered using H-CD-HIT with the same three-iteration identity settings as above.

Clustering was also conducted to assess DGR occurrence in genomes by putative cell size. DGR-like RTs were clustered using H-CD-HIT[40] as above, resulting in 699 protein clusters, which were then inspected for representatives from individual metagenomic libraries for each filter pore size (i.e., 1.2 μm, 0.2 μm, & 0.1 μm). It should be noted that the cells in the smallest size fraction passed through a filter commonly used for filter-sterilization. Larger organisms (e.g., Spirochaetes), and viruses, are able to pass through 0.2 μm filters, thus highlighting the need to carefully assess the phylogenetic affiliation of DGR-containing contigs linked to smaller filter fractions. To address the concern that filtration is an imperfect method for size exclusion, we examined genomes from one sample that was previously used for cryo-electron tomography to quantify the size of cells belonging to CPR bacteria Parcubacteria (OD1 superphylum) and Microgenomates (OP11 superphylum)[11].

DGR readmapping and metatranscriptomic analysis

We analyzed the variability of TR adenine sites leading to non-synonymous versus synonymous substitutions, by assessing reads mapped to the corresponding VR sequence on assembled scaffolds. First, assembly errors were inspected for each scaffold using stringent criteria for each basecall: any regions that were not supported with paired reads with at most one mismatch were replaced with Ns; errors are re-assembled using stringent mapping for one read in a pair; scaffolds were split if insert (Ns) coverage is zero. For the adenine-variability analysis, only uninterrupted VR and TR sequences were used (i.e. without Ns). The scaffold sequence of aligned TR and VR regions was inspected in-frame with respect to the variable protein-encoding gene’s stop codon. Reads mapped to each VR codon were analyzed for non-synonymous or synonymous substitutions at TR adenines. To compare with variability that could have resulted from stochastic mutation, non-synonymous and synonymous substitutions were also tabulated for non-adenine positions in TR.

Stringent metatranscriptomic read mapping to DGR regions was performed using Bowtie2[35], whereby only matching metatranscriptomic reads, with at most a single mismatch, were mapped to the scaffold. Regions without an apparent ORF, but with >10 reads mapped, were searched against the rFAM database[41], to identify potential RNA-encoding sections of DGR-containing scaffolds. We determined metatranscriptomic read coverage for DGR features separately (variable protein; RT; VR; TR), in addition to calculating coverage for the whole DGR cassette.

Phylogenetic analyses

To compare DGR representatives that were derived from groundwater metagenomes with DGR-like RTs from other bacterial and archaeal genomes, we sought a phylogenetic reconstruction for RT protein sequences. A consensus sequence from an alignment of previously studied DGRs[15,17] was used to search the NCBI-nr protein database for additional DGR-like hits (blastp; e-value <10−20). Next, the returned hits were individually used towards blast searches against NCBI-nr, to obtain up to 20 top hits for each DGR-like representative (e-value <10−20). Before alignment and tree construction, we performed clustering on the redundant list of RT sequences using CD-Hit46 with an intracluster global alignment threshold of 90% and default parameters. Representatives that were assessed as DGR-like (i.e. included in the DGR-specific RT tree) exhibited a monophyletic association with known DGR representatives, branching separately from Group-II intron-associated RTs, as previously shown[15,42]. Groundwater-associated and other DGR-like representatives were aligned to a hidden markov model for the reverse transcriptase protein family (PF00078) using HMMalign[43]. A phylogenetic tree of the RT alignment was constructed using FastTree2 with the WAG model and CAT approximation.

Variable protein domains, including midasin-like ATPase, and separately, DUF1566-like, were aligned using ClustalW and manually inspected to remove ambiguously aligned sites. Trees for variable protein clusters were constructed in Geneious v 8.1.4 (Biomatters), using PhyML[44] with the model LG + G, while branch support was determined with 100 bootstrap replicates. All trees for this study were visualized using FigTree (v1.4.2).

Data Availability

The data supporting the results of this study are available within the paper and its Supplementary Information and Supplementary Data Files. Assembled metagenomic sequences are available in the NCBI BioProject database under the accession numbers PRJNA268032, PRJNA273161, PRJNA288027, and KY476664-KY476802.