Microbial metabolisms in a 2.5-km-deep ecosystem created by hydraulic fracturing in shales

Hydraulic fracturing is the industry standard for extracting hydrocarbons from shale formations. Attention has been paid to the economic benefits and environmental impacts of this process, yet the biogeochemical changes induced in the deep subsurface are poorly understood. Recent single-gene investigations revealed that halotolerant microbial communities were enriched after hydraulic fracturing. Here, the reconstruction of 31 unique genomes coupled to metabolite data from the Marcellus and Utica shales revealed that many of the persisting organisms play roles in methylamine cycling, ultimately supporting methanogenesis in the deep biosphere. Fermentation of injected chemical additives also sustains long-term microbial persistence, while thiosulfate reduction could produce sulfide, contributing to reservoir souring and infrastructure corrosion. Extensive links between viruses and microbial hosts demonstrate active viral predation, which may contribute to the release of labile cellular constituents into the extracellular environment. Our analyses show that hydraulic fracturing provides the organismal and chemical inputs for colonization and persistence in the deep terrestrial subsurface.

S hale gas accounts for one-third of natural gas energy resources worldwide. It has been estimated that shale gas will provide half of the natural gas in the USA, annually, by 2040, with the Marcellus shale in the Appalachian Basin projected to produce three times more than any other formation 1 . Recovery of these hydrocarbons is dependent on hydraulic fracturing technologies, where the high-pressure injection of water and chemical additives generates extensive fractures in the shale matrix. Hydrocarbons trapped in tiny pore spaces are subsequently released and collected at the wellpad surface, together with a portion of the injected fluids that have reacted with the shale formation. The mixture of injected fluids and hydrocarbons collected is referred to as 'produced fluids'.
Microbial metabolism and growth in hydrocarbon reservoirs has both positive and negative impacts on energy recovery. Whereas stimulation of methanogens in coal beds enhances energy recovery 2 , bacterial hydrogen sulfide production ('reservoir souring') decreases profits and contributes to corrosion and the risk of environmental contamination 3 . Additionally, biomass accumulation within newly generated fractures may reduce their permeability, decreasing natural gas recovery. Despite these potential microbial impacts, little is known about the function and activity of microorganisms in hydraulically fractured shale.
Initial work by our group and others [4][5][6][7][8][9] used single marker gene analyses to identify microorganisms from several geographically distinct shale formations. These analyses showed similar halotolerant taxa in produced fluids several months after hydraulic fracturing. To assign functional roles to these organisms, we conducted metagenomic and metabolite analyses on input and produced fluids up to a year after hydraulic fracturing (HF) from two Appalachian basin shales, the Marcellus and Utica/Point Pleasant (Utica) formations. Although an earlier metagenomic study examined shaleproduced fluids 10 , the microbial communities were only sampled for nine days after HF. Here, we have reconstructed the first genomes from fractured shale, examining the microbial metabolisms sustained in these engineered, deep subsurface habitats over a period of 328 days. We provide evidence for metabolic interdependencies, and describe chemical and viral factors that control life in these economically important ecosystems. Our results show microbial degradation of chemical additives, the potential for microbially induced corrosion and the formation of biogenic methane, all of which have implications for the sustainability of energy extraction.

Reconstruction of persisting shale genomes
Our earlier study surveyed microbial community structure in fluids from three hydraulically fractured Marcellus shale wells 5 . Five fluid samples from a single well were chosen for paired metagenomic and metabolite analyses, as these samples represented three phases of the energy extraction process. Fluid samples were collected from input materials, and at early (7 and 13 days) and late (82 and 328 day) time points (Fig. 1a) Table 3), signifying that the underlying data were well represented. We also validated that the assembled genomes reflected the microbial identities and abundances in the unassembled reads, by comparing genome bin relative abundance to reconstructed near-full-length 16S rRNA genes 11 (Supplementary Data File 1).
Consistent with our earlier taxonomic study, six halotolerant bacterial and archaeal members became enriched at later time points (82 and 328 days). We recovered six Halanaerobium, two Halomonadaceae, four Marinobacter, one Methanohalophilus, one Methanolobus and two bins from Halobacteroidaceae (Supplementary Table 2 and Supplementary Data File 2). Each of these taxa contain halotolerant and thermotolerant members. Environmental sequences closely related to 16S rRNA genes recovered here were similar to those also recovered from other hydrocarbon reservoirs or hypersaline environments (Supplementary Data File 3). For each of these persisting taxa (Fig. 1b) we obtained a representative genome that was at least 90% complete and, with the exception of Marinobacter (which contained several low-abundance Marinobacter strains), had less than 1% estimated contamination (Supplementary Table 4). The Halobacteroidaceae genus lacked closely related 16S rRNA genes (∼94% identity) and genomes (76% average nucleotide identity, ANI), suggesting this organism may be unique to shales ( Supplementary Fig. 1). Following the naming convention for near-complete (>95% sampling) genomes from metagenomic data sets 12 , we propose the genus name Candidatus Frackibacter based on the colloquial name for hydraulic fracturing, or 'fracking'. We infer that changes in membership at these later time points were due to growth of these specific taxa rather than DNA persistence in the environment, as cellular biomass increased in this well (Fig. 1b).
Emphasizing the persistence of specific taxa during energy extraction, members of the terminal community were identified from the input fluid through either identical genomes or closely related 16S rRNA genes (Fig. 1c). For instance, a Halanaerobium genome detected at both days 82 and 328 had an identical genome in the input fluid (ANI ∼99% with gene synteny) (Supplementary Table 5). Conversely, we did not detect lowerabundance members of the terminal community (for example, methanogens and Candidatus Frackibacter, <2% 16S rRNA gene  Figure 1 | Genomic sampling from hydraulically fractured Marcellus shale fluids over time. a, Non-metric multidimensional scaling ordination of 16S rRNA gene data from our previous study 5 show similar community trajectories across wells. The Marcellus fluid samples analysed here for metagenomic and metabolite analyses are indicated by green stars. b, Genome bin relative abundance at each time point, with taxa coloured according to the legend and persisting members outlined in a black box. For each time point, cell count data are overlaid with error bars representing the mean ± s.d. of technical replicates (n = 3). c, 16S rRNA gene relative abundance of persisting taxa present at low abundance in the input fluids. The star in the Halanaerobium orange bars (b,c) denotes an identical genome recovered from input, day 82 and day 328 samples. abundance) in the input fluids, probably because they were below our detection limit (Supplementary Data File 1b). This finding is the first to demonstrate that HF creates a habitat where low-abundance microorganisms are injected into the deep subsurface, bloom, and persist despite biocide addition 13 , elevated temperatures (65°C) and pressures (at least 25 MPa), and salinities that ultimately become briny.

Glycine betaine and chemical additives fuel methanogenesis
Halite dissolution from the shale matrix drives large salinity increases in the produced fluids 9 , so organisms must have adaptations for tolerating a broad salinity range (Fig. 2). We recovered multiple osmoprotectant strategies from all genomes (Supplementary Table 6 and Supplementary Discussion). Our metabolite data show that, of the known osmoprotectants 14 , glycine betaine (GB) was present in the fluids, but mannitol, sorbital, ecotine and trehalose were not detected (Supplementary Table 1). Consistent as a response to salinity, GB was below detection in the input and early Marcellus shale fluids, but reached a maximum concentration at day 82 that was maintained at day 328 (Fig. 2).
Uptake and de novo synthesis of GB were features encoded in all near-complete genomes recovered over the last two time points. GB synthesis is encoded in the Methanohalophilus genome by a glycine pathway (via sarcosine and dimethylglycine intermediates) and in the Halomonadaceae and Marinobacter genomes by a pathway from choline (Supplementary Discussion). Choline, a common chemical additive in fracturing fluids, was exogenously provided in the input fluids and consumed by day 13 when Marinobacter became abundant ( Fig. 1b and Supplementary Table 1). Collectively, our paired metagenomic and metabolite findings show the production and uptake of GB is a halotolerance mechanism widely used by organisms in fractured shales.
Microbially synthesized GB, available extracellularly in the fluids, may be degraded by both of the recovered obligate fermenters (Halanaerobium, Candidatus Frackibacter) (Supplementary Discussion). Candidatus Frackibacter has two mechanisms for reducing GB. The first demethylates GB and oxidizes the methyl group via the Wood-Ljungdahl pathway, producing trimethylamine (TMA; Supplementary Fig. 2) 15 . The second pathway, also present in some shale Halananerobium genomes, uses a GB reductase (grdHI), producing TMA and acetate via a Stickland reaction 16 ( Supplementary Fig. 3). Notably, GB reduction is not widely encoded in isolated Halanaerobium genomes, being present in only 36% of the published genomes (5 of 14, http://img.jgi.doe.gov, April 2016). In the Marcellus shale, a persisting Halanaerobium (Halan-2; T82, T328) genome is the only one of our three capable of GB reduction. GB fermentation using microbially produced metabolites, rather than a dependency on input fluid chemicals, may sustain life in shales long after hydraulic fracturing.
GB fermentation yields TMA, which we infer is rapidly consumed by methanogens present at the last two time points. Each of the recovered shale methanogen genomes (Methanolobus and Methanohalophilus) has pathways for utilizing TMA, dimethylamine (DMA), monomethylamine (MMA) and methanol, but cannot use GB, hydrogen/carbon dioxide or acetate (Supplementary Discussion). In addition to microbial synthesis, HF input fluids also contain high concentrations of methylotrophic substrates (1.2 mM MMA and 331 µM methanol) that could support methanogenesis (Fig. 2). It is possible that these compounds are also assimilated as a nitrogen source (MMA) or are oxidized by Pseudomonas and Marinobacter (methanol) at earlier time points 17 . Although Methanohalophilus 16S rRNA genes have been reported from Antrim 8,9 and Burket/Geneseo 4 fractured shales, our genomic and metabolite findings identify the endogenous and exogenous sources of methyltrophic substrates, show their co-occurrence with methanogens, and confirm the metabolic pathways for methanogenesis.

Metabolisms impacting energy extraction
In addition to containing substrates that could support biogenic methane production, HF input fluids contain high concentrations  Table 1). The capacity to respire sucrose is widely encoded in our genomes (for example, Vibrio, Pseudomonas and Marinobacter), consistent with its consumption by day 7. Ethylene glycol is consumed over time, and is not detected at the last time point. This substrate is probably aerobically oxidized by Marinobacter and Vibrio at the early time points 18 (alcohol dehydrogenase), and fermented by Halanaerobium (propanediol dehydratase, acetalaldehyde dehydrogenase) later to yield ethanol, hydrogen and acetate. Candidatus Frackibacter also has the capacity to produce acetate via GB fermentation, homoacetogenesis (H 2 /CO 2 ) and sugar fermentation. Consistent with the possibility for GB, sugar and ethylene glycol fermentation at later time points, ethanol and acetate increased at day 82, when Halanaerobium, Candidatus Frackibacter and methanogens were co-enriched.
Halanaerobium are the dominant members in produced fluids from Barnett, Marcellus, Burket/Geneseo and Antrim fractured shales 4,6,7,9 . Our shale-hosted Halanaerobium genomes also have the capacity to ferment amino acids (for example, alanine; Supplementary Discussion), sucrose, fructose, glucose and maltose (Supplementary Table 7). Biofilm formation may be an important adaptation enabling the dominance of Halanaerobium across hydraulically fractured shales. Although biofilm-related genes are not detected in all surface Halanaerobium genomes 19 , these shale genomes encode genes for flagellar motility and cellular aggregation (for example, polysaccharide production and diguanylate cyclase) 20 (Supplementary Table 7).
Biofilm, organic acid and H 2 production, together with the capacity to reduce thiosulfate to sulfide (using three copies of the rhodanese-like thiosulfate:cyanide sulfur-transferase gene 3 ), implicate a role for shale Halanaerobium in steel corrosion and reservoir souring 21 . Additionally, the near-complete Halomonadaceae genome also encodes multiple thiosulfate sulfur-transferase genes, which while not previously reported in these taxa, are implicated in thiosulfate disproportionation, producing sulfide and sulfite 22 . Current microbial corrosion diagnostic practices often rely on detecting the presence of dissimilatory sulfate-reducing genes or measuring sulfate-reducing metabolic potential. However, we did not identify any sulfate-reducing genes in the Marcellus data set, suggesting the need to include alternative biological mechanisms of sulfide production for characterizing microbial corrosion potential.
Owing to the economic importance of hydrocarbons, we analysed shale-produced fluids for degradation pathways, and confirmed the presence of benzene, toluene, ethylbenzene, xylene and naphthalene (BTEX-N) in all fluid samples, while decane was detected in all but the input 5 . We failed to recover any genes for anaerobic hydrocarbon degradation, but the near-complete Marinobacter, Halomonadaceae, Pseudomonas and Vibrio genomes (input or early samples) had the capacity for aerobic hydrocarbon oxidation (Supplementary Discussion).
Of the persisting members, the Marinobacter genome encodes the capacity for alkane and BTEX-N degradation, whereas the Halomonadaceae genome lacks the first steps of these pathways but has genes for subsequent degradation ( Supplementary Fig. 4). These two taxa often co-occur in saline hydrocarbon habitats 23 ,  (Fig. 1b). Each edge represents a perfect match between a host CRISPR loci spacer sequence and a viral contig. Numbers in diamonds and circles denote the contig numbers.
including Barnett 6 and Marcellus 7 shale-produced fluids. Of our recovered genomes, these two encode the greatest metabolic versatility, enabling the use of a wide range of carbon sources (for example, acetate, lactate and hexose sugars) with O 2 or nitrate as possible electron acceptors. However, given the lack of detectable nitrate (Supplementary Table 1), we postulate that these facultative anaerobes utilize fermentative metabolisms once dissolved O 2 associated with HF has been depleted 24,25 .

Active viral predation in the deep subsurface
Our results indicate that viral-mediated cell lysis is a mechanism to explain how an intracellular osmoprotectant, like GB, was detected extracellularly in fluids. We recovered 331 viral contigs including 21 closed, circular viral genomes ( Supplementary Fig. 5). A comparison of contigs across samples showed that 318 viral contigs were unique, with only 13 viral contigs shared across time points. Of the viral contigs, 86% belonged to members of the Caudovirales, tailed dsDNA viruses, with Myoviridae (44%) and Siphoviridae (26%) families predominating. Previously, only viral reads and prophage genome fragments have been reported 10,26,27 . We mined our microbial genomes for the presence of CRISPR-Cas systems (clustered regularly interspaced short palindromic repeat-CRISPR associated), which act as an acquired immunity to viruses 28 (Supplementary Table 8). CRISPR-Cas frequency estimates range from 81% of archaeal and 40% of bacterial genomes in cultivated microbes 29 , to 10% of archaeal and bacterial genomes in metagenomic data sets 30 . In contrast, 100% of the three archaeal and 84% of the 31 bacterial genome bins of the Marcellus samples had evidence of a CRISPR-Cas system, with type 1 being the most prevalent (Supplementary Table 8). In fact, all microbial genomes at the last time point had a CRISPR-Cas system, signifying that viral immunity may be an important adaption for persistence in hydraulically fractured shales.
Comparing CRISPR arrays in microbial hosts to viral contig sequences allowed us to reconstruct a history of viral encounters, and link 34 viral contigs to 11 microbial genomes ( Fig. 3 and Supplementary Fig. 5). Before our findings, the greatest number of reported CRISPR links within a data set was five, from a threeyear study in a hypersaline lake 31 . Our data showed that viral host specificity varied, with viruses linked to multiple species within a genus (for example, Halanaerobium, Marinobacter) and a single viral genome linked to two methanogen genera (Fig. 3). We observed an increase in the number of CRISPR spacers within two Halanaerobium genomes between days 82 and 328, demonstrating that adaptive viral resistance probably occurred during this time span (Supplementary Fig. 6 and Supplementary Table 8). Our metagenomic data demonstrate that viral predation and host-acquired immunity are active processes in the deep terrestrial subsurface.

Strain metabolic diversity across shales
We examined fluid metabolites collected over 302 days after HF from the Utica shale, a geographically and stratigraphically distinct Appalachian Basin formation. Despite these differences, metabolite trends in the Marcellus and Utica produced fluids were similar. For instance, methanol and ethylene glycol were detected in input fluids and salinity increased over time in both shales (Fig. 2). Unlike the Marcellus, MMA was not detected in the input but was produced over time in the Utica produced fluids, suggestive of ongoing GB production, fermentation and subsequent methanogenesis. However, due to the chemical complexity of the Utica produced fluids, we could not confirm the presence of GB.
To validate that fractured shales harbour microorganisms that produce methane from GB fermentation products, we amended Utica produced fluids with GB. The produced fluid sample was collected 96 days after HF, comparable to our Marcellus sample, where the co-occurrence of GB fermenting and methanogenic organisms was first detected (day 82). Relative to the produced fluids, the addition of GB enriched for Methanohalophilus (70%) and three Halanaerobium genomes (∼21, 3 and 0.5%) ( Supplementary  Fig. 7). The presence of a GB reductase system probably explains the changes in relative abundance within these Halanaerobium, as the dominant genome in the produced fluids lacked grdI (decreasing from 51 to 3%). This finding demonstrates the power of genomecentric metagenomics to partition local microdiversity, explaining the co-occurrence of strains with distinct functional roles.
In the enrichment, GB fermentation produced TMA, DMA and MMA in low amounts, probably due to active consumption by Methanohalophilus ( Supplementary Tables 9 and 10). Compared to the unamended control, GB addition produced 6.5 times more methane per day. Collectively, our Marcellus field and Utica laboratory data provide evidence that GB synthesis and subsequent fermentation supports biogenic CH 4 in hydraulically fractured shales.
Comparative genomics showed that the dominant Methanohalophilus and Halanaerobium near-complete genomes (Supplementary Table 2) in the Utica enrichment were closely related strains 32 to the Marcellus genomes (∼99% ANI) (Supplementary Table 5

and Supplementary Data File 2). In contrast, ANI values between the Marcellus and Utica
Methanohalophilus and other sequenced species (M. mahii and M. halophilus, both isolated from surface waters), were ∼91 and 92%, respectively. Methanohalophilus CRISPR array comparisons identified a single spacer sequence shared between the Marcellus produced fluids and Utica GB enrichment genomes; the two other non-shale-derived Methanohalophilus genomes lacked CRISPR-Cas systems. Two viral contigs also had high sequence identity (>95%), showing that these shales share genetically similar viruses. Together, our data demonstrate that environmental filtering results in populations, metabolisms and viral processes shared between these two geographically distinct fractured shale ecosystems.

Conclusion
Resolving genomes from Marcellus and Utica produced fluids unveiled microbial metabolisms, adaptations and viral predation resistance mechanisms in fractured shales. From 16S rRNA gene analyses we could not have predicted the role Halanaerobium strains play in fermenting GB and HF chemical additives such as ethylene glycol, nor would we have associated the Halmonadaceae with detrimental sulfide production. Our genomic analyses show that closely related strains are niche differentiated. For instance, GB addition selected for the only Halanaerobium genome with GB reduction capacity. We identified the metabolic capabilities of Candidatus Frackibacter, unique to fractured shales, which can also ferment GB. Our metagenomic data revealed a possible role for viruses in the top-down (predation and lysis) and bottom-up (release of cellular contents; for example, GB) control of microbial communities in fractured shale. Notably, unlike earlier studies, all host genomes recovered at the last time point contained a CRISPR-Cas system. We also identified active host responses to viral predation, including new spacer incorporation over time. Together, our viral findings demonstrate the probable importance of CRISPR-Cas-mediated immunity for microbial persistence in fractured shales.
Here, we show that hydraulic fracturing provides the organisms, chemistry and physical space to support microbial ecosystems in ∼2,500-m-deep shales (Fig. 4). Ultimately, our metagenomic and metabolite results indicate that adaptation to high salinity, metabolism in the absence of oxidized electron acceptors, and viral predation are potential controlling factors mediating longterm microbial metabolism during energy extraction from fractured shales. This study highlights the resilience of microbial life to adapt to, and colonize, a habitat structured by physical and chemical features very different from their origin.

Methods
Sample collection and fluid chemistry. Our earlier study characterized some of the geochemistry and conducted 16S rRNA gene surveys from fluids collected (June 2012 to May 2013) from three Marcellus gas wells located in Pennsylvania, USA 5 . Hydraulic fracturing (input, noted as T0) and shale-produced fluids were collected from well heads (days 3-14) and gas-fluid separators (49, 82 and 328 days), with fluids from well 1 used for more detailed metagenomics and NMR metabolite analyses here. For our Utica samples, injected fluids and produced fluids from gasfluid separators 4,8,9 were collected between July 2014 and May 2015 from an oil-gas well in Ohio, USA. The gas-fluid separators at the Marcellus and Utica sites had a capacity of ∼5,560 l, approximately half gas and half produced fluids (2,780 l). Flow rates ranged from ∼380,000 l per day at early time points to ∼190,000 l per day at later time points (Marcellus day 328), with an estimated maximum residence time of 8 h on the days of sampling. Additionally, 16S rRNA sequences from key taxa we genomically sample here were either sampled here at earlier time points directly from the well head (for example, Halomonas and Marinobacter), or also recovered from other shale produced fluids where samples were collected exclusively from the well head (for example, Halanaerobium, Methanolobus and Methanohalophilus 9 ).
Hydraulic fracturing included the injection of freshwater amended with chemicals, proppant and the addition of 20% recycled produced fluids for the Marcellus (not Utica). Notably, unlike the Marcellus wells, the Utica well was shut in for 86 days after fracturing, before initiating fluid and hydrocarbon collection at the surface (denoted on Fig. 2). Input and produced fluid samples (1 l) were collected in sterile bottles filled to capacity. As described elsewhere 5 , ethoxylated surfactants and hydrocarbons were assessed in Marcellus fluids using liquid chromatography quadrapole time of flight with electro-spray ionization (Agilent Technology) and gas chromotography (Hewlett Packard), respectively. Conductivity and pH were measured on unfiltered fluids in the field using Orion star probes, while fluid dissolved anions (F − , acetate, formate, Cl − , Br − , NO 3 − , SO 4 2− , oxalate, S 2 O 3 2− and PO 4 3− ) were analysed using a Thermoscientific Dionex ICS-2100 ion chromatograph with the exception of the Utica input fluid. Samples for ion chromatography were diluted by a factor of 10-100, due to the high salinity. The Utica input fluid had high viscosity, which precluded analysis by ion chromatography. Conductivity and Cl − were not measured for this input sample. NO 3 − and SO 4 2− in the Utica input were measured on unfiltered samples with a HACH DR/890 portable colorimeter using cadmium reduction method 8171 and turbidimetric method 8051, respectively. Fluid samples for major and trace cations (Na, Mg, K, Ca, Si, Sr, Ba, Li, Mn and Fe) were acidified immediately after filtration to ∼0.5% with nitric acid and then analysed using a Perkin Elmer Optima 4300DV inductively coupled plasma optical emission spectrometer. The charge balance discrepancy for the input samples is probably due to unmeasured cations (for example, ammonium and organic additives in the fracking fluids), and has been documented in other studies 33 Figure 4 | Interconnected metabolisms catalysed by persisting microorganisms in hydraulically fractured shales. a, HF input fluids from both Marcellus and Utica shales contain substrates that sustain microbial metabolism. Parentheses indicate metabolites detected in one shale. b, Microorganisms in shales adapt to high salinities by producing and using osmoprotectants such as GB (red circles), which can be released into fluids by viral lysis. c, Marinobacter and Halomonadaceae have the potential to aerobically oxidize hydrocarbons and respire sugars using nitrate and oxygen as electron acceptors. d, Candidatus Frackibacter and Halanaerobium ferment GB, yielding trimethylamine, which supports methanogenesis by Methanohalophilus and Methanolobus (blue box). Methylamines and methanol in the input fluids can also support methanogenesis (yellow box).
internal standard. All NMR spectra were collected using a Varian Direct Drive 600 MHz NMR spectrometer equipped with a 5 mm triple resonance salt-tolerant cold probe. The 1D 1 H NMR spectra of all samples were processed, assigned and analysed using Chenomx NMR Suite 8.1 with quantification based on spectral intensities relative to the internal standard. Candidate metabolites present in each of the complex mixtures were determined by matching the chemical shift, J-coupling and intensity information of experimental NMR signals against the NMR signals of standard metabolites in the Chenomx library. The 1D 1 H spectra were collected following standard Chenomx data collection guidelines 34 , using a 1D nuclear Overhauser effect spectroscopy (NOESY) presaturation experiment with 65,536 complex points and at least 512 scans at 298 K. Additionally, 2D spectra (including 1 H-13 C heteronuclear single-quantum correlation spectroscopy (HSQC), Metagenomic sequencing and assembly. For genomic sample collection, 300-1,000 ml samples were concentrated onto 0.22-µm pore size polyethersulfone (PES) filters (Millipore, Fisher Scientific). Viruses were probably obtained on filters by flocculation with iron, which precipitated during the filtering process when the samples were first exposed to oxygen 34 . Total nucleic acids were extracted from the filter using the PowerSoil DNA Isolation kit (MoBio) for Marcellus fluids and a modified phenol chloroform nucleic extraction 35 for Utica fluids and enrichment cultures. Total cells with intact membranes were enumerated from unfiltered fluid samples for calibrated gate ranges on a Guava EasyCyte flow cytometer (EMD Millipore). Briefly, samples were fixed with 1% glutaraldehyde and stained with 0.1% SYBR Gold (Life Technologies), and quantified via flow-cytometry. For each time point, technical triplicates were measured and the data reported in Fig. 1 represent the mean ± s.d. (n = 3).
For the Marcellus input and produced fluids, Illumina HiSeq 2000 libraries were prepared using the Nugen Ovation Ultralow Library System following the manufacturer's instructions. Genomic DNA was sheared by sonication, and fragments were end-repaired. Sequencing adapters were ligated and library fragments were amplified with ∼8-10 cycles of PCR before Pippin Prep size selection, library quantification and validation. Libraries were sequenced on the Illumina HiSeq platform and paired-end reads of 113 cycles were collected. Fastq files were generated using CASSAVA 1.8.2. Similar protocols were used for the Utica fluid GB enrichment culture, where sequencing was conducted on an Illumina HiSeq 2500 platform using a Kapa Hyper Prep library system with five cycles of PCR before solid-phase reversible immobilization (SPRI) size selection.
All metagenomics methods and scripts contributing to analyses in this manuscript are included in Supplementary Data File 4. Briefly, Illumina sequences from each of the five samples (input, T7, T13, T82 and T328) were first trimmed from both the 5′ and 3′ ends using Sickle (https://github.com/najoshi/sickle), then each sample was assembled individually using IDBA-UD (refs 36,37) with default parameters. Scaffold coverage was calculated by mapping reads back to the assemblies using Bowtie2 (ref. 37). Given the dominance and high strain variation in some samples, highly abundant genomes (>400×) often failed to assemble. Using an approach outlined in ref. 38, subassemblies were performed to reconstruct the dominant genomes in the day 13 and day 82 samples, using 10 and 8% of the reads, respectively. Results from the subassemblies are included (Supplementary Table 2).
Metagenomic annotation and genomic binning. All scaffolds ≥5 kb (≥1 kb for subassemblies, Methanohalophilus-1, and the GB enrichment culture) were included when binning genomes from the metagenomic assembly. Scaffolds were annotated as described previously 36,37 by predicting open reading frames using MetaProdigal 39 . Sequences were compared using USEARCH 40 to KEGG, UniRef90 and InterproScan 41 with single and reverse best hit (RBH) matches greater than 60 bits reported. The collection of annotations for a protein were ranked: reciprocal best BLAST hits (RBH) with a bit score >350 were given the highest (A) rank, followed by reciprocal best blast hit to Uniref with a bit score >350 (B rank), blast hits to KEGG with a bit score >60 (C rank), and UniRef90 with a bit score greater than 60 (C rank). The next rank represents proteins that only had InterproScan matches (D rank). The lowest (E) rank comprises the hypothetical proteins, with only a prediction from Prodigal but a bit score of <60. Complete annotation files for all contigs >1,000 are available for download from https://chimera.asc.ohio-state.edu/daly_et_al_nature.html.
Within each sample we obtained the genome resolved 'bins' using a combination of phylogenetic signal, coverage and GC content 36,37 . For each bin, genome completion was estimated based on the presence of core gene sets (highly conserved genes that occur in single copy) for Bacteria (31 genes) and Archaea (104 genes) using Amphora2 (ref. 42; Supplementary Table 4). Overages (gene copies >1 per bin) indicating potential misbins, along with GC and phylogeny, were used to manually remove potential contamination from the bins.
To illustrate the microbial similarities shared between three Marcellus wells located in close proximity, we used 16S rRNA gene membership and abundance data from our earlier 454 amplicon study 5 to generate a non-metric multidimensional scaling ordination using Bray-Curtis distances. The ordination had a stress of 0.12, indicating that the matrix data were well represented by the ordination (Fig. 1a). Samples selected for metagenomics reflect the changing community over time and are denoted by a star in Fig. 1a. The relative abundance of each assembled genome in a sample was calculated as a proportion of the summed average coverage of the binned contigs in each sample. The relative abundance of the taxa in the input fluid that become dominant at later time points (Fig. 1c) was based on the normalized relative abundance of reconstructed 16S rRNA genes using EMIRGE (ref. 11), as genomic bins were not recovered for all these taxa. To verify the accuracy of our binned genomes, near-full-length ribosomal 16S rRNA gene sequences were reconstructed from unassembled Illumina reads from Marcellus fluids and an Utica fluid GB enrichment culture using EMIRGE (Supplementary Data File 1) 11 . To reconstruct 16S rRNA gene sequences we followed the protocol with trimmed paired-end reads where both reads were at least 20 nucleotides used as inputs and 50 iterations. EMIRGE sequences were chimaera checked before phylogenetic gene analyses.
For some genomes that lacked high strain resolution, such as Methanohalophilus-1, we confirmed the manual binning using an emergent selforganizing map (ESOM) using both the metagenomic data and isolate genomes from Marinobacter, Methanohalophilus, Methanolobus, Halanaerobium and Halomonas isolated species, as described previously 36,43 . Tetranucleotide frequencies were calculated for ≥5 kb fragments, with the number of tetranucleotides in each fragment normalized on the basis of the total number of observations in all fragments, with these values robust Z-transformed. The resulting matrix was used to train an ESOM for 30 epochs using scripts previously reported (https://github.com/ tetramerFreqs/Binning). The ESOM was visualized using the Databionic ESOM Tools software.
Taxonomic placement of the genome bins relied on the phylogenetic analyses of 16S rRNA and/or ribosomal proteins. To determine if the same genome was present in different time points, we calculated ANI values (http://enve-omics.ce.gatech.edu/ani/) using a two-way ANI, with ≥99% ANI considered an initial cutoff for identical genomes through time. For genera with ≥99% ANI values (Arcobacter, Halanaerobium and Idiomarina), we then aligned contigs to examine synteny using the progressive Mauve aligner 44 in Geneious R8 (ref. 45). If a clustered regularly interspaced short palindromic repeat (CRISPR) array was present in each genome being compared over time, the contigs with CRISPR arrays were preferentially chosen for alignment, as CRISPR arrays are hyper-variable regions and are dynamic at short timescales due to new spacer incorporation. Several high-quality bins (with a representative of each taxa persisting in later time points) were selected for manual curation and genome finishing.
Viral genome binning, CRISPR identification and links to microbial hosts. Viral contigs were identified through annotations by including contigs ≥5 kb with viral structural genes (for example, capsid proteins, tail proteins and terminases), contigs containing a high number of proteins with no known homology, and using Metavir 2 (ref. 46) comparison to the viral RefSeq database. Circular contigs, indicating complete viral genomes, were determined using two methods: (1) analysis in the Metavir 2 (ref. 46) software by identifying identical k-mers at the two ends of the sequence and (2) manually examining SAM files generated by Bowtie2 (ref. 47) for paired reads present at the two ends of the sequence at the appropriate coverage level. Similar contigs shared across time points, or between Utica and Marcellus shales, were identified by comparing contigs using the criteria of >95% ANI over 80% of the contig length, analogous to the clustering in ref. 48. Contigs within clusters were then individually aligned and manually inspected to confirm identical contig sequences.
Crass was used to identify CRISPR repeat and spacer sequences 49 , and was run both on individual sample reads and combined reads from all samples. To identify the microbial hosts of viruses and determine if viral predation was ongoing in the deep shale, we used BLASTn with an E-value cutoff of 1e-8 to identify contigs with repeat and spacers. First, we matched repeats to genomic bins to identify host CRISPR loci. Next, within each identified host CRISPR loci, we matched the spacers identified by Crass to viral contigs, identifying the viruses (or highly similar viruses) the host had encountered previously. All matches were manually confirmed as perfect matches by aligning sequences in Geneious R8, and we used the CRISPR Recognition Tool plugin (CRT, version 1.2) to confirm CRISPR loci in genomic bins. Spacer links between host genomic bins and viral contigs were used to construct a network (Fig. 3) in Cytoscape (version 3.1.0). Each spacer match between a host CRISPR loci and a virus represents an edge in the network; nodes represent hosts (ovals) or viruses (diamonds and circles). Multiple edges indicate cases where a host CRISPR loci had multiple spacer links to the same viral contig, and these are represented by multiple (curved) edges in the network. For identical genomes present across samples (for example, Halanaerobium-1) where the viral-host links were only detected in genomes from some samples, the genome in other samples (oval nodes without edges) was also included. These links were probably not detected due to low estimated genome completion in some samples. Based on a recent paper 50 , we classified CRISPR-Cas modules by manually examining the operon architectures of annotated contigs. Although we cannot determine whether a particular viral contig was in a virulent/lytic state at the time of sampling, in the input and day 7 samples, the contigs with the highest read coverage, identified as viral by the above criteria, had coverage several-fold higher than the most abundant bacterial or archaeal contig (input, 5.1-fold higher; T7, 4.0-fold higher), suggesting these viruses were virulent at the time of sampling.
Phylogenetic and metabolic analyses. The 16S rRNA genes recovered from the five nearest neighbours to each EMIRGE sequence from our Marcellus metagenomic samples were obtained from SILVA (release 123) 51 , anchored with cultivated representatives (Supplementary Date File 1). 16S rRNA gene sequences were aligned in Geneious R8 using MUSCLE, a phylogenetic tree was constructed with RAxML 7.2.8 (GTR Gamma nucleotide model, 999 bootstrap replicates), and relative abundance data over time were graphed using iTOL. For the S3 protein tree, aminoacid sequences were pulled from the Marcellus and Utica GB enrichment genomic bins and were augmented with sequences mined from NCBI and JGI IMG databases (Supplementary Data File 2). Sequences were aligned using MUSCLE version 3.8.31 and run through ProtPipeliner, a python script developed in-house for generation of phylogenetic trees (https://github.com/lmsolden/protpipeliner). A maximum likelihood phylogeny for the alignment of S3 ribosomal proteins was conducted using RAxML version 8.3.1 under the LG+α+γ model of evolution with 100 bootstrap replicates and visualized in iTOL. The 16S rRNA genes recovered for the key terminal taxa in Fig. 1 were also blasted to NCBI to show the relationship between key terminal taxa in this study and environmental sequences. Sequences were trimmed to the V3-V4 region, and the top ten non-redundant hits from NCBI were included for analysis. 16S rRNA gene sequences were aligned in Geneious R8 using MUSCLE, and a phylogenetic tree was constructed with RAxML 7.2.8 (GTR Gamma nucleotide model, 999 bootstrap replicates) (Supplementary Data File 3).
Metabolic profiling was largely conducted by manual analyses. For the osmoprotectants inventory, the annotated gene lists were searched by name, KEGG number and E.C. number with positive records saved to files that were manually inspected to remove misidentified genes. The results were compared to the same functions in available genomes from the same genus in the IMG database on 15 August 2015 (http://img.jgi.doe.gov/). For key functional genes we used both a list and homology-based approach to help annotate genes. The latter is important, as many methylamine cycling genes were incorrectly annotated or not included on scaffolds >5,000 bp. Putative GrdEGI/PrdA were identified from the Utica GB enrichment by blasting known homologues capable of glycine/sarcosine/GB/proline reduction. The MttB and GrdEGI/PrdA trees were constructed similarly to above by aligning manually curated amino sequences in MUSCLE with RAxML version 8.3.1 using the LG model of evolution with 100 bootstrap replicates. For key genes in the metabolic processes outlined in the main text, we confirmed protein structure and functional prediction via modelling, catalytic or structural residues, or phylogenetic analyses.
GB enrichment culture from Utica produced fluids. NMR metabolite analyses were performed on Utica produced fluids that had been filtered and stored at −80°C since the time of collection, after 483 days of anoxic (100% N 2 headspace) incubation without amendment, and before and after amendment to stimulate GB utilizing and methane-producing organisms. The GB enrichment consisted of 30% anoxic, incubated, produced Utica fluid (day 96) and 70% sterile modified DSMZ 479 media dispensed in Balch tubes sealed with butyl rubber stoppers and aluminium crimps under an atmosphere of N 2 /CO 2 (80:20, vol/vol). Before mixing with produced fluids, the modified DSMZ medium (per litre) included 87 g sodium chloride, 1.5 g potassium chloride, 6.0 g magnesium chloride, 0.4 g calcium chloride, 1.0 g ammonium chloride, 2.0 g yeast extract, 2.0 g trypticase peptone, 0.2 g coenzyme M, 0.2 g sodium sulfide, 4.0 g sodium bicarbonate and brought to a pH of 7.1 using 1 mM NaOH. After autoclaving, this medium was supplemented with 9.7 mM GB and a trace element solution (DSMZ 141). As a no-donor GB control, 30:70 produced fluid:medium was established in parallel that lacked GB substrate amendment. GB enrichment and control cultures were monitored for methane production over 32 days using a Shimadzu gas chromatograph equipped with a thermal conductivity detector (TCD) using helium as a carrier gas at 100°C. At 32 days after initial enrichment, DNA was extracted for 16S rRNA gene bacterial and archaeal clone libraries (data not included) and metagenomics. Analyses were conducted as described above, with the exception that 1% of the reads were used to recover the dominant Halanaerobium genomic bin. All scaffolds in the assembly, regardless of length, were searched for homologues for MttB, GrdI, and methanogenesis pathways, CRISPR arrays were only analysed for the binned scaffolds (nucleotide ≥1 kb).