Produced water samples were collected from 14 hydraulically fractured wells in North Dakota actively producing oil from the Bakken Formation or the underlying Three Forks Formations at the time of sampling in May 2018 (Figure 1). As production data from these formations is typically reported together, we refer to the sampled region as the Bakken Shale or Bakken region throughout this manuscript. Each produced water sample was obtained from the associated three-phase separator, which separates hydrocarbons from the remaining fluids. Separators are often the closest available sampling point to the wellhead and have been used in previous studies (Cluff et al., 2014; Mohan et al., 2014; Lipus et al., 2017). Samples were collected in unused plastic carboys that were pre-rinsed with sample waters and allowed between 1 and 3 h to settle into distinguishable oil and water phases in the sealed container. A portion of the water phase was also collected in a sterile 1 L nalgene bottle and immediately placed on dry ice until arrival in the lab for microbiology analysis. Dry ice was replenished as needed, such that samples remained frozen until they arrived in the lab. Upon arrival in the lab, samples were stored according to previously described guidelines at −20°C until microbial processing (Lipus et al., 2017, 2018, 2019). The remaining water phases were then passed through glass wool cartridges under gravity and then a 0.45 μm inline filter (EnviroTech, Salt Lake City, UT, United States), minimizing oxygen exposure in a closed loop. This filtered water was collected and measured for pH using a Horiba multimeter (Horiba, Edison, NJ, United States) and titrated for alkalinity using a Hach digital titrator (Hach, Loveland, CO, United States). Samples were collected with zero-headspace in Shimadzu TOC vials and transported on ice for total organic carbon (TOC) analysis. A 30 mL volume of the filtered water was acidified to 2% by volume using nitric acid and transported on ice for inductively coupled plasma optical emission spectroscopy (ICP-OES). Finally, a 15 mL volume of the filtered water was passed through a second 0.22 μm sterile polyethersulfone filter (Millipore, Inc.) and shipped on ice for ion chromatography (IC) measurement.

Major cations and anions were detected using ion chromatography (IC) on a ThermoFisher (Thermofisher, Waltham, MA, United States) ICS-5000+ with AS11-HC column for anion quantification and CS16 column for cation quantification. All samples were run in triplicate and the standard error of IC measurements reported here was less than 3%. Additionally, every 10–20 samples, a cation/anion control sample (Sigma Aldrich, St Louis, MO, United States) with known certified concentrations was added during measurements and all control samples demonstrated an accuracy within 95–105%.

Trace metals were analyzed using inductively coupled plasma optical emission spectroscopy (ICP-OES) on an Optima 7300 DV (Perkin Elmer, Waltham, MA, United States), which is a dual-view spectrometer with solid state SCD detectors. U.S. EPA Method 6010D was employed for analysis with one duplicate, one standard recovery, and one spike recovery. Sample solutions were nebulized using a glass Seaspray concentric nebulizer and a glass Baffled Cyclonic spray chamber (Glass Expansion, Pocasset, MA, United States) using 5 ppm yttrium as an online internal standard. Calibration standards were purchased from Inorganic Ventures (Christiansburg, VA, United States) and are traceable to NIST standard reference materials.

Total organic carbon (TOC) was analyzed with a Shimadzu TOC Analyzer (Model: TOC-LCSN) for non-purgeable organic carbon (NPOC) and total inorganic carbon (TIC). Average values from 3 to 5 replicates were reported with <2% precision (%RSD), and quality control samples for every 10–12 samples demonstrated consistent accuracy within 95–105%.

We extracted DNA from the samples using a modified version of the DNeasy Powersoil kit (Qiagen, Hilden, Germany). Produced water samples were filtered through a 0.2 micron polyethersulfone membrane filter (Pall Corporation, Port Washington, NY, United States). The filter with the collected biomass was then transferred to a 1.5 mL microcentrifuge tube incubated with 1 mL TE buffer and 10 μl of 20 mg/mL lysozyme (Sigma Aldrich) at 37°C for 30 min. Samples were exposed to 10 min of bead beating on a vortexer followed by extraction using the recommended guidelines in the manufacturer’s protocol. DNA was eluted from spin filters in 100 μL nuclease free water. Four kit blanks were also concurrently extracted to confirm that no contamination occurred.

The bacterial abundance in produced water samples was measured using qPCR with 16S rRNA gene primers (F: GTGSTGCAYGGYTGTCGTCA; R: ACGTCRTCCMCACCTTCCTC) designed by Maeda et al. (2003) with an expected amplicon size of 146 base pairs. qPCR reactions were run in triplicate on a Magnetic Induction Cycler (MIC) (Bio Molecular Systems, Upper Coomera, Australia). Each reaction contained 2X SensiFAST SYBR No-Rox master mix (Bioline, London, United Kingdom), 400 nM forward primer, 400 nM reverse primer, and 1 μL of template DNA for a total reaction of 20 μL. qPCR conditions consisted of a polymerase activation step at 95°C for 2 min followed by 40 amplification cycles each consisting of: denaturation at 95°C for 5 s, annealing at 62°C for 5 s, and an extension step at 72°C for 1 s. Standard curves were generated using gBlocks Gene Fragments (Integrated DNA Technologies, Coralville, IA, United States) (Supplementary Data 1 and Supplementary Figure 1) and negative control samples were included in each amplification assay.

Extracted DNA and kit blanks were amplified using universal primers targeting the V4 region of the 16S rRNA gene, as previously described (Caporaso et al., 2011, 2012). In order to obtain enough PCR product, we pooled 8 technical replicates per sample. The pooled PCR products were cleaned using AMPure or SPRIselect beads (Beckman Coulter, Pasadena, CA), quantified using the Qubit dsDNA High Sensitivity Assay Kit (Life Technologies, Carlsbad, CA, United States), and visualized on a Bioanalyzer using the High Sensitivity DNA Kit (Agilent, Santa Clara, CA, United States). Negative PCR controls were also amplified, quantified, and visualized in order to confirm that no contamination occurred. In total, we attempted to amplify 4 kit extraction blanks and 6 negative PCR controls. After quantification and visualization, only one amplified kit extraction blank and one negative PCR control contained amplicons. These samples, identified as Blank 1 and 2 in the Supplementary Data 1, were sequenced with our experimental 16S rRNA libraries. Purified 16S rRNA libraries were pooled, diluted to a concentration of 2 nM, and denatured using fresh 0.2 M NaOH. Libraries were further diluted according to the manufacturer’s instructions and sequenced on an Illumina Miseq (Illumina, San Diego, CA, United States) using a 300 cycle V2 Nano kit.

16S rRNA gene sequences were analyzed using QIIME2 version 2019.10 (Bolyen et al., 2019). Sequences were imported as EMPSingleEndSequences and demultiplexed using the demux emp-single command. DADA2 (Callahan et al., 2016) was used to filter, denoise, and remove chimeras from the demultiplexed sequencing data. For this command, we utilized default settings and set the truncation length to 250 base pairs. The classify-sklearn (Pedregosa et al., 2011) command was used to classify representative sequences identified through DADA2 using a pre-trained Naive Bayes classifier trained on Silva 132 99% OTUs (Quast et al., 2012; Yilmaz et al., 2014) from the 515F/806R region and any sequences identified as chloroplast or mitochondria were filtered from the data.

Data generated by Qiime2 was imported into R for further analysis (R Core Team, 2013). We utilized the vegan package (Oksanen et al., 2012) to calculate the Chao1 index, Shannon index, and Bray–Curtis dissimilarity values. We also used the vegan package (Oksanen et al., 2012) to complete non-metric multidimensional scaling (NMDS), calculate Analysis of Similarities (ANOSIM), complete the Mantel test, and to fit environmental parameters onto the generated NMDS ordination plot. All NMDS plots were constructed with k = 3 dimensions. NMDS, ANOSIM, and the Mantel test were completed using Bray–Curtis distance measurements calculated after sequence libraries were resampled to the depth of the sample with the fewest sequences from this experiment (1942 sequences) (Supplementary Table 1). The Mantel test also relied on the Euclidean distance of the geochemical data, which was calculated in Base R. The significance of the ANOSIM analysis, the Mantel test, and the environmental vectors was based on 999 permutations of the grouping, geochemical, and environmental data, respectively. Finally, we also used the RColorBrewer package to generate colorblind accessible palettes when necessary (Neuwirth, 2014).

Produced water samples were obtained from 14 wells in the Bakken region (Figure 1). The pH was slightly below circumneutral with a range of 6.1–6.7 (Table 1), which is similar to previously reported results (Wang et al., 2019). The exception was well 6, which had a pH of 5.8 (Table 1). Alkalinity varied from low to high, with the majority of wells ranging from 186 to 740 mg/L as CaCO₃. Exceptions were wells 5, 6, 7, and 10, which ranged from 1440 to 2000 mg/L as CaCO₃ (Table 1). Although alkalinity of 2000 mg/L and higher have been reported in this region (Blondes et al., 2019; Chang et al., 2019), the variation is notable, as previous work involving these wells observed only low alkalinity measurements (Lipus et al., 2018). TOC was measured in the form of non-purgeable organic carbon (NPOC), which ranged from 63 to 543 mg as C/L (Table 1), and is similar to previously reported DOC values for the Bakken region (Wang et al., 2019).

The TDS concentrations ranged between 135,000 mg/L and 357,000 mg/L (Table 1) and were characterized by 86.1% wt to 88.6% wt Na⁺ and Cl^– (Table 2). These TDS concentrations are similar to those previously reported within the Bakken Shale (Thyne and Brady, 2016; Lipus et al., 2018; Wang et al., 2019) and reinforce that produced water from the Bakken region has higher TDS concentrations than samples collected from other shale plays, including the Marcellus and Barnett (Akob et al., 2015; Lipus et al., 2017; Wang et al., 2019). These high salinities are likely derived from post-depositional halite and anhydrite dissolution processes (Bachu and Hitchon, 1996; Iampen and Rostron, 2000). TDS measurements in conjunction with major cation and anion concentrations can either reflect the formation brine (Iampen and Rostron, 2000; Engle et al., 2016) or fluid mixing processes related to hydraulic fracturing at unconventional wells (Haluszczak et al., 2013; Vengosh et al., 2017). Samples in this study display a positive linear association (R² = 0.9892) between TDS and chloride (Figure 2A), a trend that is seen in other Bakken produced waters (Thyne and Brady, 2016) and is expected since chloride is the primary anion in this study’s samples. This trendline includes well 11, suggesting that this water’s TDS and chloride concentrations are still compositionally similar to other samples in this study. This is notable, as produced water from well 11 exhibited a lower TDS concentration (135,000 mg/L) compared to the other samples, which had an average of 287,000 mg/L (Table 1). TDS levels from well 11 also fall within the range of recorded Bakken produced water measurements when data reported in the USGS PWDB were plotted along with samples from this study (Blondes et al., 2019) (Supplementary Figure 2A). The low TDS of well 11 may be due to its spatial heterogeneities within the reservoir; Bakken produced waters have high spatial TDS variability, as there is minimal fluid flow in certain regions due to the formation’s aquitard characteristics (Bachu and Hitchon, 1996; Hitchon, 1996).

Sulfate is a typical secondary anion in produced waters and increases in sulfate concentration can reflect drilling activity, either due to injected chemicals associated with hydraulic fracturing operations (Paukert Vankeuren et al., 2017) or due to injected water reacting with sulfur and sulfate bearing minerals (Osselin et al., 2019). Although anhydrite content is negligible, trace amounts of pyrite have been observed in the Bakken Formation (Pitman et al., 2001; Ashu, 2014). When TDS is plotted against sulfate and the youngest and most spatially unique sample (well 18) is excluded there is a negative linear relationship (R² = 0.853) (Figure 2B). Although most wells were hydraulically fractured more than 4 years ago (Table 1), well 18 is 1 month removed from a completed hydraulic fracturing procedure. Therefore, the elevated sulfate concentration in well 18 is either the result of spatial heterogeneity or recent hydraulic fracturing. When TDS and sulfate for other Bakken produced waters are plotted (Supplementary Figure 2B), there is no clear linear association (R² = 0.0811) between the parameters. This suggests samples from this study are exhibiting a localized trend. Despite this, the young age of well 18 makes it difficult to conclude whether it falls off the linear trendline (Figure 2B) due to its location or the recent hydraulic fracturing operation.

Previous research has demonstrated that evaluating produced water chemistry using compositional data analysis can remove inaccurate correlations (Engle and Rowan, 2013; Engle and Blondes, 2014; Blondes et al., 2016). When they are applied to produced waters, compositional data analysis techniques include using isometric log-ratio transformations to understand whether concentration variations indicate water-rock or fluid mixing processes. Molar ratios of four analytes (calcium, sulfate, sodium, and chloride) were transformed using isometric log-ratios and were plotted (Figure 2C) for Bakken produced waters. The region near the origin represents equilibrium conditions with the evaporite minerals halite and sulfate. Samples from this study plotted within the range of USGS PWDB samples, but they plotted on the high end of the calcium-sulfate ratios (Figure 2C). Previous work studying produced waters has suggested similarly elevated calcium-sulfate ratios may be the result of sulfate reduction (Engle et al., 2016). This suggests that sulfate reduction was more prevalent in samples from this study than other samples from the Bakken region.

In order to determine the total microbial abundance in the evaluated Bakken region produced waters, qPCR analyses were conducted on the 14 collected samples. The median microbial load was found to be 1.15E+2 16S rRNA gene copies/mL of produced water. The microbial abundances across all produced water samples ranged from 2.15E+1 to 3.06E+5 16S rRNA gene copies/mL. The microbial load in well 11 was significantly higher than in the other samples and appeared to be an outlier (Table 1 and Supplementary Figure 3). However, it falls within the microbial loads recorded in the Marcellus Shale and Denver-Julesburg Basin, which range from 1.00E+5 to 2.10E+8 and 1.0E+3.1 to 1.0E+8.5, respectively (Lipus et al., 2017; Hull et al., 2018).

To assess the microbial community composition, we conducted 16S rRNA gene amplicon sequencing for 11 of the 14 produced water samples. No DNA could be amplified for the remaining three samples, potentially due to low overall biomass. Additionally, 10 blanks were prepared, two of which had measurable DNA and were therefore sequenced. 16S rRNA amplicon sequencing yielded a total of 675,581 raw sequences, of which 289,457 remained after quality control filtering (Supplementary Table 1). We also concurrently analyzed previously published sequences (Lipus et al., 2018) generated from samples collected at the same sites and within the same region using the same bioinformatic pipeline (Supplementary Table 1).

Alpha diversity was determined by counting the total number of ASVs in a sample and calculating the Shannon richness, Chao1 index, and Pielou’s evenness metric for each site (Table 3, Supplementary Table 1, and Supplementary Figure 4). Each sample had an average of 58 unique ASVs per 1942 sequences. Chao1 indices ranged from 37 to 103, Shannon indices ranged from 2.42 to 3.68, and Pielou’s evenness measurements ranged from 0.64 to 0.89 per 1942 sequences across all produced water samples (Table 3, Supplementary Table 1, and Supplementary Figure 4). As a whole, our samples were significantly more rich, diverse, and even when compared to samples previously collected from the same wells by Lipus et al. (2018) (Supplementary Table 1 and Supplementary Figure 4). Our samples also had a greater spread of Chao1 indices and a smaller spread of both Shannon indices and Pielou’s evenness measurements (Supplementary Figure 4). The notable exception is that samples collected in January 2015 by Lipus et al. (2018) had comparable levels of evenness and a much smaller spread for all alpha diversity metrics (Supplementary Figure 4). One possible reason for this trend is that after several years of production, well conditions select for a microbial community with a high diversity of adapted microorganisms.

The analyzed produced water samples were dominated by microbiota from 19 bacterial families, each of which represents ≥ 5% of sequences from any one sample. Twelve of the dominant families have been previously isolated and cultured from hydrocarbon reservoirs and 5 have been previously identified through sequencing hydrocarbon-associated environmental samples (Korenblum et al., 2012; Tatar, 2018). These dominant families comprise between 76.24% and 97.08% of the total microbial community in each well (Figure 3A, Supplementary Table 3, and Supplementary Figure 5). The families Halanaerobiaceae, Methylomonaceae, Xanthobacteraceae, Burkholderiaceae, Desulfohalobiaceae, and Pseudomonadaceae were especially abundant, with an average relative abundance of 57.02% and a range between 4.43 and 83.24% (Figure 3A, Supplementary Table 2, and Supplementary Figure 5). The one exception is produced water from well 10, which was dominated by the Clostridiales Family XI and Shewanellaceae, representing 35.80 and 31.32%, respectively, of the total microbial community (Figure 3A, Supplementary Table 3, and Supplementary Figure 5). Shewanellaceae has previously been detected in saline environments, oil field wastes, and other hydrocarbon environments (Nealson and Scott, 2006), suggesting an adaptability to conditions similar to unconventional reservoirs. It was previously found in produced water from the Bakken region (Lipus et al., 2018), although only at low abundances, suggesting adaptation to well conditions over time. No other wells contained such high abundances of Family XI and Shewanellaceae. However, Family XI and Shewanellaceae were present in the remaining produced water samples with an average relative abundance of 2.17 and 1.58%, respectively (Figure 3A, Supplementary Table 3, and Supplementary Figure 5).

Two other identified bacterial families include Desulfobacteraceae and Desulfomicrobiaceae, which have an average relative abundance of 1.13 and 1.09% across all produced water samples, respectively (Figure 3A, Supplementary Table 3, and Supplementary Figure 5). Both families were especially dominant in well 4, constituting 5.72 and 6.54%, respectively, of the total relative abundance (Figure 3A, Supplementary Table 3, and Supplementary Figure 5). Desulfobacteraceae also represented 5.6% of the relative abundance in well 18 (Figure 3A, Supplementary Table 3, and Supplementary Figure 5). Desulfobacteraceae and Desulfomicrobiaceae are putative sulfate reducers and may directly contribute to well souring, corrosion, and/or other operational failures. Bacillaceae were found in every produced water sample and accounted for 0.01–8.67% of the total community within each well (Figure 3A, Supplementary Table 3, and Supplementary Figure 5). This is especially notable as samples previously collected at the same sites were dominated by Bacillales (Supplementary Figure 5) (Lipus et al., 2018). Finally, Veillonellaceae (Order: Selenomonadales) were detected in 7 of the 11 produced water samples and constituted 16.17% of the microbiota present in well 17 (Figure 3A, Supplementary Table 3, and Supplementary Figure 5). The remaining dominant bacterial families include Dysgonomonadaceae, Marinobacteraceae, Anaerolineaceae, Rhodocyclaceae, Arcobacteraceae, Clostridiaceae 1, and Enterobacteriaceae (Figure 3A, Supplementary Table 2, and Supplementary Figure 5). None of these taxa were found in every produced water sample and they were identified at average relative abundances of 0.88–3.50% across all wells. Archaea were also present in six of the 11 produced water samples and contributed a maximum of 0.89% of the relative abundance in any one well. The archaeal orders present in our samples included Methanobacteriales and Methanomicrobiales, which are common methanogens previously implicated in hydrocarbon degradation (Jiménez et al., 2016; Toth and Gieg, 2018).

We completed several analyses in order to identify factors potentially contributing to and/or driving changes in the produced water microbial ecology. First, we constructed three non-metric multidimensional scaling (NMDS) plots using Bray-Curtis dissimilarity metrics. These ordination plots contain: (1) all samples and blanks from our current study (Supplementary Figure 6A), (2) samples and blanks from our current study as well as samples previously collected by Lipus et al. (2018) from the same sample sites (Supplementary Figure 6B), and (3) samples and blanks from our current study and all samples previously collected by Lipus et al. (2018) from the Bakken region (Supplementary Figure 6C). Next, we used analysis of similarity (ANOSIM) calculations to confirm any visible clustering present in each of the three plots. We specifically tested for the clustering by the following categories, when possible: sample origin (tank, separator, or blank), well number, formation, well age, and sampling event. We found no statistically significant clustering based on sample origin, formation, or well age (ANOSIM: p > 0.05) in our first plot (Supplementary Figure 6A). However, incorporation of previously collected data at the sample sample sites revealed clustering by sample origin (ANOSIM: R = 0.1769; p < 0.05), well age (ANOSIM: R = 0.4350; p < 0.05), and sampling event (ANOSIM: R = 0.5622; p < 0.05). Samples did not cluster by well number or formation (ANOSIM: p > 0.05 for both). Interestingly, when we completed an ANOSIM with our experimental data and the paired 2018 data (but without our sequenced blanks), there was no separation by sample origin (ANOSIM: p > 0.05). These trends remained consistent when we analyzed our sequenced blanks, experimental data, and all previously published sequences by Lipus et al. (2018) (Supplementary Figure 6C). We found no significant clustering by sample origin, well number, or formation although there was significant clustering by well age (ANOSIM: R = 0.4533; p = 0.001) and sampling event (ANOSIM: R = 0.4874; p = 0.001) (Supplementary Figure 6C).

In order to examine the relationship between the microbial community within each well and the associated geochemical profile, we completed a Mantel test. We found a statistically significant, but weak relationship (R = 0.2336; p = 0.0430). Next, we calculated and fit vectors onto a fourth ordination plot, which contained only our experimental data, in order to identify environmental factors associated with produced water microbial structure (Figure 3B and Supplementary Figure 7). We determined that 11 of the 29 geochemical measurements have a significant correlation with the microbial data visualized on the NMDS plot (Figure 3B and Supplementary Figure 7). The R² of the significant correlations was strong, ranging from 0.56 to 0.84 with an average R² of 0.65 (Supplementary Figure 7). TDS was significantly correlated with the microbial data, with an R² of 0.62. TDS is a measurement of the total amount of minerals, salts, metals, cations, and/or anions dissolved in water, therefore it is unsurprising that many of the TDS constituents also had a significant correlation with the microbial data. However, there were several notable trends. First, the Bakken Shale is known to have high NaCl concentrations compared to other shale regions (Table 2 and Supplementary Table 2) (Lipus et al., 2018). Consistent with this, sodium and chloride were significantly correlated with the microbial data (R² = 0.62 and R² = 0.59, respectively), although NaCl% was not. Second, the R² for boron, calcium, copper, potassium, and manganese were higher than the R2 for TDS, suggesting that these ions may be especially linked to microbial metabolism within the Bakken Shale Formation. Other geochemical measurements with a significant correlation to the microbial data include bromide, sulfate, and TIC (Figure 3B and Supplementary Figure 7).