Methanogens and Their Syntrophic Partners Dominate Zones of Enhanced Magnetic Susceptibility at a Petroleum Contaminated Site

Geophysical investigations documenting enhanced magnetic susceptibility (MS) within the water table fluctuation zone at hydrocarbon contaminated sites suggest that MS can be used as a proxy for investigating microbial mediated iron reduction during intrinsic bioremediation. Here, we investigated the microbial community composition over a 5-year period at a hydrocarbon-contaminated site that exhibited transient elevated MS responses. Our objective was to determine the key microbial populations in zones of elevated MS. We retrieved sediment cores from the petroleum-contaminated site near Bemidji, MN, United States, and performed MS measurements on these cores. We also characterized the microbial community composition by high-throughput 16S rRNA gene amplicon sequencing from samples collected along the complete core length. Our spatial and temporal analysis revealed that the microbial community composition was generally stable throughout the period of investigation. In addition, we observed distinct vertical redox zonations extending from the upper vadose zone into the saturated zone. These distinct redox zonations were concomitant with the dominant microbial metabolic processes as follows: (1) the upper vadose zone was dominated by aerobic microbial populations; (2) the lower vadose zone was dominated by methanotrophic populations, iron reducers and iron oxidizers; (3) the smear zone was dominated by iron reducers; and (4) the free product zone was dominated by syntrophic and methanogenic populations. Although the common notion is that high MS values are caused by high magnetite concentrations that can be biotically formed through the activities of iron-reducing bacteria, here we show that the highest magnetic susceptibilities were measured in the free-phase petroleum zone, where a methanogenic community was predominant. This field study may contribute to the emerging knowledge that methanogens can switch their metabolism from methanogenesis to iron reduction with associated magnetite precipitation in hydrocarbon contaminated sediments. Thus, geophysical methods such as MS may help to identify zones where iron cycling/reduction by methanogens is occurring.

Geophysical investigations documenting enhanced magnetic susceptibility (MS) within the water table fluctuation zone at hydrocarbon contaminated sites suggest that MS can be used as a proxy for investigating microbial mediated iron reduction during intrinsic bioremediation. Here, we investigated the microbial community composition over a 5year period at a hydrocarbon-contaminated site that exhibited transient elevated MS responses. Our objective was to determine the key microbial populations in zones of elevated MS. We retrieved sediment cores from the petroleum-contaminated site near Bemidji, MN, United States, and performed MS measurements on these cores. We also characterized the microbial community composition by high-throughput 16S rRNA gene amplicon sequencing from samples collected along the complete core length. Our spatial and temporal analysis revealed that the microbial community composition was generally stable throughout the period of investigation. In addition, we observed distinct vertical redox zonations extending from the upper vadose zone into the saturated zone. These distinct redox zonations were concomitant with the dominant microbial metabolic processes as follows: (1) the upper vadose zone was dominated by aerobic microbial populations; (2) the lower vadose zone was dominated by methanotrophic populations, iron reducers and iron oxidizers; (3) the smear zone was dominated by iron reducers; and (4) the free product zone was dominated by syntrophic and methanogenic populations. Although the common notion is that high MS values are caused by high magnetite concentrations that can be biotically formed through the activities of iron-reducing bacteria, here we show that the highest magnetic susceptibilities were measured in the free-phase petroleum zone, where a methanogenic community was predominant. This field study may contribute to the emerging knowledge that methanogens can switch their metabolism from methanogenesis to iron reduction with associated magnetite precipitation in hydrocarbon contaminated sediments. Thus, geophysical methods such as MS may help to identify zones where iron cycling/reduction by methanogens is occurring.

INTRODUCTION
Magnetic susceptibility (MS) is a measure of the ability of a material to become magnetized when placed in a magnetic field and is dependent on the concentration and type of the magnetic minerals present. It generally reflects the presence and volume of the iron oxide magnetite. Because of the abundance of iron in terrestrial environments, MS is commonly used to investigate different geological processes. Example applications of MS include: (1) determining the magnetic stratigraphy in loess deposits, lake sediments and marine sediments to reconstruct the paleoenvironmental evolution and erosional history (e.g., Zhang et al., 2012), (2) examining variations in the detrital input of magnetic minerals related to climate change (e.g., Da Silva et al., 2013), (3) aiding archeological investigations (Walach et al., 2011), (4) serving as a proxy for extreme hydrologic events and rainfall (Li et al., 2019), (5) monitoring soil pollution (Hoffmann et al., 1999;D'Emilio et al., 2007;Zhang et al., 2012), and (6) investigating the movement of tectonic plates. Further, the MS technique has also been used for oil exploration (e.g., Guzmán et al., 2011), for evaluation of oil bearing units (Perez-Perez et al., 2011;Emmerton et al., 2013) and reservoir properties (Potter, 2007;Potter and Ivakhnenko, 2008).
There is also a growing interest in relating environmental magnetic studies to microbial mediated iron cycling (Porsch et al., 2010;Mewafy et al., 2011;Byrne et al., 2015), and oil biodegradation after spills (Rijal et al., , 2012Atekwana et al., 2014;Beaver et al., 2016;Lund et al., 2017). In particular, Rijal et al. (2010Rijal et al. ( , 2012 and Atekwana et al. (2014) observed large MS anomalies within the water table fluctuation zone (WTFZ) associated with the presence of hydrocarbons relative to background measurements and related these changes to the precipitation of magnetite resulting from microbial degradation of the hydrocarbons by iron-reducing bacteria.
Biodegradation is a natural process performed by microbes that is dependent on the terminal electron acceptors (TEAs) available in the environment (Meckenstock et al., 2015). Oil degradation using aerobic respiration provides microbes with the most energy, but oxygen is eventually exhausted and sediments turn anoxic. After oxygen is depleted, microbial communities switch to different TEAs. Nitrate reduction is the next most energetically favorable process. When nitrate is used up, iron reduction occurs. After Fe(III) is depleted, sulfate reduction and finally fermentation and methanogenesis will take over as predominant microbial processes. To produce methane, methanogens can use methylated compounds such as acetate or methanol as substrates or just carbon dioxide (CO 2 ) and hydrogen gas (H 2 ). In fact, the methanogens' use of H 2 leads to a locally very low hydrogen partial pressure, making hydrocarbon fermentation energetically favorable. This interspecies hydrogen transfer represents a syntrophic relationship between fermenters and methanogens. In other syntrophic associations between respiratory bacteria and methanogens, electrons are transferred via certain bacterial appendages such as nanowires (direct interspecies electron transfer -DIET) (Cruz Viggi et al., 2014;Rotaru et al., 2014;Zhuang et al., 2015).
Magnetite is often found in sediments where iron reduction has occurred. When the activity of iron-reducing bacteria increases the local concentration of Fe(II), Fe(II) will sorb to ferrihydrite, causing the precipitation of the mixed valent magnetite (Hansel et al., 2005).
Therefore an increase in MS can indicate if iron-reducing microorganisms are active while degrading hydrocarbons (Porsch et al., 2010). Nonetheless the relationship between MS and the activity of iron-reducing bacteria is not completely understood and the use of MS as a proxy for biodegradation is tenuous at best. This is underscored by the observation that the MS signature can be transient due to the creation and destruction of magnetic minerals by biotic and abiotic processes (e.g., Ameen et al., 2014;Lund et al., 2017). Some of these magnetic minerals may be used by microbes to transfer electrons, causing changes in mineral structure. Redox cycling using magnetite by Rhodopseudomonas, an iron oxidizer, and Geobacter, an iron reducer, has been demonstrated under laboratory conditions, and MS was shown to decrease during the oxidation phase from Fe dissolution and increase during the reduction phase from magnetite precipitation (Byrne et al., 2015). In the field, Rijal et al. (2012) and Ameen et al. (2014) documented temporal decreases in MS at a contaminated former military air force base in Hradčany, Czech Republic, and speculated this to be due to the dissolution of magnetite. Lund et al. (2017) also documented a >90% loss in MS at the National Crude Oil Spill Fate and Natural Attenuation Research Site administered by the United States Geological Survey (USGS) near Bemidji, MN, United States, between 2011 and 2015 when compared to initial MS observations from this aquifer (Mewafy et al., 2011;Atekwana et al., 2014). It was hypothesized that the reduction of MS was due to the transformation of magnetite to siderite coupled with the exhaustion of ferrihydrite.
Currently, we only have a poor understanding of the microbial communities driving the changes in MS and this uncertainty limits the effective application of MS as a proxy for biodegradation and iron cycling studies. The relationship between oil biodegradation and MS is complex as bacterial oil biodegradation coupled to iron reduction could result in the precipitation of magnetite and cause an increase in MS response. Once the bioavailable iron is depleted, the microbes may switch to other metabolic processes that use magnetite. This, in turn, can cause a decrease in the MS response through the oxidation of magnetite to maghemite and then hematite. Interestingly, initial characterization of the microbial community composition at the Bemidji site revealed that not iron-reducing bacteria were prevalent in zones where high MS was measured, but a microbial community dominated by the syntrophic bacterium Smithella and methanogenic Archaea (Beaver et al., 2016).
In this publication, we extend the Beaver et al. (2016) study, which was limited in scope because it showed data based on 1 year only, from just four samples along a 10 m long core, and it used clone libraries with the 16S rRNA gene, which is a limited molecular approach for the characterization of microbial communities. Here, we present a multi-year analyses of the microbial community structure with closely spaced samples along complete cores using high-throughput DNA sequencing at the same hydrocarbon contaminated site exhibiting transient MS signatures. Our objective was to determine the key microbial populations in zones of elevated magnetic susceptibility and to provide insights into possible microbial-driven MS changes.

Site Description
This study was conducted at the National Crude Oil Spill Fate and Natural Attenuation Research Site administered by the USGS. Here, a pipeline transporting crude oil broke in 1979 spilling about 1,700,000 L of crude oil. After cleanup efforts were completed in 1980, about 400,000 liters of oil remained in the unsaturated zone on the water table. This oil continues to be a source of groundwater contamination.
The spill resulted in two plumes in two low topographic areas designated as the north and south pools. The north oil pool is the location of the rupture site whereas the south oil pool resulted from the sprayed oil within a small wetland region ( Figure 1A). The hydrocarbon contamination at the site is characterized by different hydrocarbon phases including a freefluid phase approximately 0.75 m thick floating on the water table, a dissolved phase in the ground water, and a vapor phase in the unsaturated zone. Seasonal fluctuations of the water table have resulted in a smear zone characterized by residual phase smearing aquifer sediments (Essaid et al., 2011). Our study focuses on cores obtained from contaminated sediments in the free-phase petroleum in the north oil pool ( Figure 1A).
The geology has been described as ∼20-25 m of moderately calcareous silty sand and outwash glacial deposits, overlying an unknown thickness of clayey till (Bennett et al., 1993). The upper 3 to 6 m are characterized by mostly poorly sorted, interbedded coarse sand, gravel, and silt. Beneath the upper layer, a 15-20 m laminated silt, very fine sand, and clay make up the majority of the lithology and this is underlain by the lower till layer.

Core Retrieval
All five cores from the free-phase oil zone were drilled within 1.9 m of each other and, according to their retrieval year, named C1110, C1215, C1308, C1407, and C1507 ( Figure 1B). According to the fluctuating water table, these five cores were 8-10 m long and their sections are described in Supplementary  Table 1. Core 1411 was retrieved from a location in the dissolved phase zone as indicated in Figure 1A. Cores were retrieved via percussion drilling using a core barrel containing either a 1 or 2 polycarbonate liner pushed through the sediment in a hollow stem auger. For the lowest core section crossing the water table, a sample-freezing drive shoe was used to retain the water from the aquifer by pumping liquid CO 2 to the bottom of the core which allowed the retention of the final section (Trost et al., 2018). A detailed description of the core sections can be found in Supplementary Table 1. Core C1110 was shipped on ice to Oklahoma State University (OSU, Stillwater, OK, United States), and stored frozen until MS analysis and sediment sample extraction for microbial analysis. Cores collected in years 2012, 2013, 2014, and 2015 were transported on ice to Western Michigan University, where they were stored at 4 • C while the sampling of sediments for DNA isolation occurred in less than 10 days after arrival at the laboratory. Sediment samples were stored at −20 • C until DNA isolation. The cores were subsequently shipped on ice to OSU, where they were frozen until MS analysis.

Magnetic Susceptibility (χ) Measurements
Magnetic susceptibility is the degree to which a sediment or mineral can be magnetized when a magnetic field is applied. MS is the ratio between M and H, where M is the magnetization of the material in the magnetic field and H is the field intensity (Blakely, 1995 .65 kHz ± 1%; maximum resolution of 2 × 10 −6 SI) coupled to a MS2B dual frequency sensor. For cores retrieved in 2014 and 2015, the volume MS values were acquired using the same meter, but it was coupled to a MS2C core scanning sensor, allowing the measurements to be taken from the complete cores that were slid through the sensor.

Gas Concentrations
Concentrations of O 2 + Ar, CO 2 , and CH 4 were obtained from vapor well 601G ( Figure 1A) as described by Mewafy et al. (2011) from the USGS National Water Information System database accessed on 11/7/2020 (U. S. Geological Survey, 2020). Based on limitations of the gas-chromatographic method used, oxygen and argon could not be separated and are plotted together.

Sampling and DNA Extraction
Sediment samples were obtained from holes cut through the side of the polycarbonate core liner. The sampling occurred in an anaerobic chamber (Shel Lab Model BactronEZ, Sheldon Manufacturing, Cornelius, OR, United States) for all years except for 2011. The anaerobic chamber was continuously fed with 90% nitrogen, 5% CO 2 and 5% H 2 . Sediment samples were taken along the complete cores from top to bottom in ∼15 cm distances. Sterile spatulas were used to extract the sediments, which were placed aseptically into plastic centrifuge tubes. The sediment samples were then stored at −20 • C in order to preserve the DNA until it was extracted. The PowerSoil DNA isolation kit (MO BIO, now Qiagen, Germantown, MD, United States) was used to isolate total DNA as directed in the protocol, but 0.5 g of sediment was substituted for 0.25 g in order to increase the DNA yield for cores obtained in 2011, 2012, and 2014. This method still did not produce enough DNA for analysis in some The black star indicates the location of the five cores that were collected from the free-phase petroleum zone. The yellow star represents the location of the core that was collected in the dissolved-phase plume. The white dot represents the location of well 601G, which is ∼5 m away from the location of the cores. The blue arrow indicates the direction of the ground water (GW) flow. BTEX, benzene, toluene, ethylbenzene, xylene. The image was retrieved from the USGS Bemidji Research Site (https://mn.water.usgs.gov/projects/bemidji/spatial/ index.html) and modified. (B) Photographic image of the location of the five wells from which the free-phase cores were drilled in five consecutive years. The five wells are not more than 1.9 m apart.
Frontiers in Earth Science | www.frontiersin.org of the samples for those years, so a different method that can extract more DNA from the sediments was used for 2013 and 2015. The method used for the samples obtained in 2013 and 2015 was the DNA extraction method described by Zhou et al. (1996), except that 4 g of sediment was used instead of 5 g due to the size of the centrifuge tubes. Although the use of different extraction methods has been shown to cause some variation in microbial community structure analysis (Guo and Zhang, 2013), To determine if any contaminating DNA was present in the DNA extraction kits or other reagents, control samples were obtained by following the extraction procedures using all the components without the addition of any sediment.

Sequencing and Statistical Analyses
Total DNA of the samples was sent to the RTSF Genomics Core at Michigan State University (East Lansing, MI, United States) for the construction of 16S rRNA gene libraries and high-throughput DNA sequencing. The libraries were constructed with the primers F515 and R806, targeting the V3/V4 variable regions of the 16S rRNA gene. The Illumina MiSeq platform (Illumina, San Diego, CA, United States) was used for sequencing the libraries by 250 base paired-end sequencing. Sequences were binned and annotated with QIIME. Forward and reverse reads were merged and quality controls were performed with Pandaseq (Masella et al., 2012). If the merged reads consisted of sequences shorter than 247 or longer than 275 base pairs, less than 45 base pairs as the minimum overlap of the forward and reverse read, or had an alignment score of less than 0.9, the read was discarded. Chimeras were identified with identify chimeric_seqs.py and eliminated, then operational taxonomic units (OTUs) were picked with the pick_open_reference_otus.py command using VSEARCH (Rognes et al., 2016) in QIIME (Caporaso et al., 2010). Only OTUs with 2 or more sequences were kept as part of the analysis. The reference set used with QIIME was the "Silva_119_rep_set97.fna" set which used a 97% similarity cutoff and taxonomic identification was completed with the "majority_taxonomy_7_levels.txt" file found in SILVA QIIME release 128. Contaminant OTUs found in the control samples without soil were removed from the core sample data manually in the Excel program (Microsoft Office Suite, Version 15.32).

Generation of Graphs and Statistical Analysis
Phyla and UniFrac graphs were generated and statistical analyses were performed with R (R Core Team, 2016), specifically the following R packages:  (Auguie, 2017). The nmds graph was generated with vegan and ggplot2. Predicted functional genes were determined by PICRUSt (Langille et al., 2013). The LEfSe analysis was generated with dada2 v. 1.0.3 (Callahan et al., 2016), microbial (Guo and Gao, 2020) and phyloseq.

Description of Cores
Sediment cores C1110, C1215, C1308, C1407, and C1507 were retrieved from the center of the free-phase oil at the north pool (Figure 1). Complete cores were retrieved from the surface at the elevation of ∼433 m down to below the water table, which was at ∼424 m. The units for the elevations are meters above sea level (m.a.s.l.) and pertain to all elevations mentioned in here, but are omitted for readability. The cores in the free-phase oil zone were taken no more than 1.9 m apart from each other to minimize spatial variability ( Figure 1B). The cores were characterized by five distinct zones based on sediment characteristics from visual inspection and gas data (Figures 2A,E Table 1). The first zone was the upper vadose (UV) zone, which was not contaminated by any free-phase or residual oil and consisted of layers of sand or clay. This zone extended down to ∼ 428.5 m, approximately where the carbon dioxide and oxygen concentration graphs cross over in Figures 2C,G,K. The second zone, the lower vadose (LV) zone, was characterized by decreasing oxygen and increasing CO 2 and methane concentrations and was part of the unsaturated zone between the upper vadose zone and the smear zone, ending at ∼425.5 m. In the smear zone (SZ), residual phase hydrocarbons adhered to the sand particles above the free-phase oil. The length of the smear zone was dependent on the depth of the water table. The fourth zone, the free-phase oil (FP), consisted of a ∼0.75 m thick layer of free-phase petroleum on top of the water table, which was at approximately 424 m. The saturated (S) zone, with dissolved phase hydrocarbons, was below the free-phase oil zone. The geological description reveals mostly fine-grained sand with clay and silt across the water table.

Magnetic Susceptibility (MS)
Sediments were sampled from the cores, and subjected to MS measurements in the laboratory, augmenting our previous in situ measurements that were taken at the Bemidji site   Atekwana et al., 2014;Lund et al., 2017). For core 1110, volume MS measurements ranged from a low of 29 × 10 −5 SI at 427.0 m to a peak of 287 × 10 −5 SI at 424.8 m at the water table ( Figure 2B). In addition, there were several intermediate peaks found with one more pronounced peak in the vadose zone at 427.3 with 209 × 10 −5 SI and another one in the saturated zone at 423.6 with a MS value of 188 × 10 −5 SI. The MS measurements for core C1215 collected in 2012 varied from a low of 36 × 10 −5 SI at 422.6 m to a peak of 959 × 10 −5 SI at 430.1 m in the vadose zone ( Figure 2F). The second largest MS peak occurred in the smear zone at 425.2 m, with a measurement of 549 × 10 −5 SI. Several other peaks were found around the water table and in the lower vadose zone with values between 350 and 500 × 10 −5 SI. No MS measurements were taken in 2013.
Core C1407 displayed several higher peaks around the water table with values above 350 × 10 −5 SI, including the highest peak of 561 × 10 −5 SI at 423.4 m ( Figure 2J). Interestingly, core 1507 collected in 2015 showed an extremely high peak with MS values up to 2,000 × 10 −5 SI in the upper vadose zone (Figure 2N), at a similar depth compared to highest MS peak found in the core retrieved in 2012. Several other peaks with MS values above 500 × 10 −5 SI occurred in the lower vadose zone. Unfortunately, a core section from around the water table was destroyed during transport and MS could not be measured. The elevated MS values found around the water table and in the free-phase petroleum in the cores from 4 different years are at similar depths that exhibited elevated MS values during the in situ measurements at the Bemidji site (Atekwana et al., 2014).

Gas Concentrations
Gas concentrations were measured in years 2011, 2012, and 2014 in well 601G, which is closest (∼5 m) to the core locations in the free-phase petroleum. Based on the limitations of the gaschromatographic method used, oxygen and argon could not be separated and are therefore plotted together (Earth's atmosphere consists of 20.95% and 0.934% argon). Figures 2C,G,K show that oxygen concentrations decreased with increasing depths, while methane concentrations increased in the lower depths in all 3 years. Maximum methane concentrations were measured in 2011 and 2014 with 17% at the lowest depths. The CO 2 concentrations varied between 9 and 16% below 428.8 m. The O 2 -CO 2 transition zone was at ∼429.4 m, whereas the O 2 -CH 4 transition zone was at ∼428.8 m. The changes in concentration of these gasses indicate that sediments were anoxic below 428.8 m.
No gas data were available from well 601G for 2015.

Analysis of Replicates
Sediment samples were collected from the center of the retrieved cores with the spacing of 15 cm along the complete core length and DNA was isolated from these samples. First, we wanted to determine the reliability of our DNA isolation method and performed the DNA isolation and subsequent high-throughput DNA sequencing for core C1407 in duplicates. A UniFrac analysis was performed, which is a phylogenetic method for comparing microbial communities. It shows the relatedness of the microbial communities from each sample in a 2-dimensional space (Figure 3). This figure illustrates that (i) the replicate samples clustered very close together, and (ii) that there are three larger groupings which represent the free-phase oil zone (bottom left), the lower (bottom right) and the upper vadose zone (top). The two replicate smear zone samples were spaced somewhat away from the free-phase samples in the direction toward the lower vadose samples (Figure 3).

Comparison of Two Different Extraction Methods
During the 5 years of this study, we optimized the DNA extraction techniques. For samples collected in 2011-2013, we used a commercially available kit, while for the soil samples collected in 2014 and 2015, we used a more laborious method, not involving a kit, which was published by Zhou et al. (1996). Some variations are expected through the use of different extraction methods (Guo and Zhang, 2013). Nevertheless, Supplementary Figure 1 shows that the two different methods applied to the same sediment sample collected from the free-phase oil zone displayed a similar microbial community composition. Deltaproteobacteria were found in similar high abundance with both methods (∼40%), while there seems to be a slight bias toward Actinobacteria in the DNA extracted with the commercial kit (25.1% vs. 14.2%) and a bias toward Euryarchaeota (3.9% vs. 8.5%) and Betaproteobacteria (5.5 % vs. 10.9 %) in the DNA extracted with the more laborious DNA extraction method. Overall, despite these biases, both methods revealed a quite similar microbial community composition (Supplementary Figure 1).

Analysis of Microbial Phyla
The abundance of microbial phyla according to depth is shown in Figure 4 for the samples collected during the five different years. Each zone has a characteristic phyla distribution with little year-to-year variations. In the free-phase oil zone, shown at the bottom of the graphs, Deltaproteobacteria were found in high relative abundance. They were accompanied by methanogenic Euryarchaeota. The smear and lower vadose zone were dominated by Alpha-and Beta-Proteobacteria, representing methanotrophs and iron-reducers, respectively. Highly abundant phyla in the upper vadose zone were common soil bacteria from the Actinobacteria and Acidobacteria phyla (Figure 4).

Comparison of Free-Phase Samples With Samples From the Dissolved Phase
A comparison between the microbial communities present in the free-phase oil zone versus those present in a location in the dissolved-phase plume was warranted. Lower MS values were found around the water table in the dissolved-phase plume than in the free-phase oil (Atekwana et al., 2014). It was no problem isolating sufficient microbial DNA quantities for downstream applications from the cores retrieved in the free-phase petroleum zone. It seems that the overabundance of hydrocarbons at the water table resulted in active microbial communities producing sufficient biomass for making our DNA isolations successful throughout the length of the core. Unfortunately, and probably because of the much lower hydrocarbon concentration ∼100 m away from the free-phase oil (Figure 1A), there seems to be much less biomass present in the dissolved plume. We were only able to isolate sufficient DNA quantities for reliable DNA sequencing from samples around the water table from the core retrieved in the dissolved phase. A UniFrac analysis shows that the samples from around the water table of the dissolved phase core clustered together on the one side of the diagram, whereas the samples taken from around the water table from free-phase cores clustered on the other side, far away from each other (Figure 5). This is an indication that the microbial community composition differs greatly between these two groups of samples. Also, the microbial phyla composition and the top 10 most abundant genera at this site were completely different from the ones found in the free-phase oil zone (Supplementary Figure 2 and Supplementary Table 2). Table 1 shows the results using the linear discriminant analysis (LDA) effect size (LEfSe) method, which identifies the most abundant and characteristic bacterial taxa in different sample groups (Segata et al., 2011). The most abundant and characteristic bacterial genera for a specific sample group receive a high LDA score and are considered to be a "biomarker" for that specific group of samples and are discussed below according to the zonation in the free-phase oil cores.

Free-phase oil zone
The genus found in the free-phase oil zone with the highest LDA score was Smithella. Laboratory-grown Smithella is capable of propionate oxidation, and it was found to be the most abundant long-chain alkane degrader in methanogenic enrichment cultures (Liu et al., 1999;de Bok et al., 2001;Gray et al., 2011;Cheng et al., 2013;Tan et al., 2014;Wawrik et al., 2016;Qin et al., 2017;Xia et al., 2019;Ji et al., 2020). To metabolize the longchain alkanes, Smithella activates them by fumarate addition and then excretes CO 2 and H 2 , which, in turn, are being used by methanogens. The genus with the second highest score found in all samples of the free-phase oil zone was a hydrogenotrophic methanogen, Methanoregula (Brauer et al., 2011). Hydrogenotrophic methanogens use CO 2 and H 2 to synthesize methane. Methanoregula species have been isolated from acidic peat bogs and from methanogenic sludge used for  the treatment of brewery effluents (Brauer et al., 2011;Yashiro et al., 2011). In addition, an unclassified genus of the Candidatus Yanofskybacteria was found to be a biomarker in the freephase oil samples. This candidate genus belongs to the newly described superphylum Patescibacteria and is a member of the Parcubacteria (formerly OD1), which are prevalent in anoxic groundwater, sediments and aquifers. Parcubacteria are of very small size, have reduced genomes, and seem to depend on other members of the microbial community for obtaining essential metabolites Tian et al., 2020). Candidatus Yanofskybacteria have frequently been found to be associated with anoxic carbon transformation and hydrogen generation Solden et al., 2016;Jaffe et al., 2020). Another biomarker in the free-phase oil zone was Desulfosporosinus. This spore-forming bacterium has been isolated from gasolinecontaminated aquifers (Robertson et al., 2000) and has been identified as a syntrophic toluene degrader in a methanogenic microbial community (Abu Laban et al., 2015). Another highscoring microorganism in the free-phase oil zone was Pelolinea from the Chloroflexi phylum, which has been shown to be involved in the degradation of phenolic compounds (Chen et al., 2020). All the above-mentioned genera are known to be strict anaerobes.

Smear zone
Uncultured Coriobacteriaceae, which belong to the Actinobacteria phylum, had the highest LDA score in the smear zone. These obligate anaerobes are capable of both, fermentation and anaerobic respiration. Lentini et al. (2012) found representatives of uncultured Coriobacteriaceae in enrichment cultures with ferrihydrite, goethite and hematite as electron acceptor. Coriobacteriaceae have also recently been identified in methanic environments such as terrestrial serpentinizing fluid seeps (Woycheese et al., 2015), deep groundwater (Hernsdorf and Amano, 2017), marine enrichment cultures undergoing benzoate degradation with concurrent iron reduction and methanogenesis (Aromokeye et al., 2020) and at a thermophilic electromethanogenic biocathode (Kobayashi et al., 2017). Thermincola, a Gram-positive bacterium of the Peptococcaceae family, also exhibited a high score for the smear zone. In an enrichment culture, Thermincola ferriacetica was found to be dominant under ferrihydrite-reducing conditions (Kato et al., 2019), and the formation of magnetite was observed when amorphous Fe(III) hydroxide was used as electron acceptor for T. ferriacetica (Zavarzina et al., 2007). Magnetite could serve as an alternative electron acceptor for T. ferriacetica in these sediments (Zavarzina et al., 2007). Also having a high score in the smear zone was Paenibacillus, a facultative anaerobic endosporeformer of the Firmicutes phylum. Paenibacillus strains have been frequently isolated from soil and they can degrade a variety of aliphatic and aromatic hydrocarbons (Grady et al., 2016). Despite considered being aerobic, the Alphaproteobacterium Beijerinckia was also a biomarker of the smear zone, and in addition to its ability to fix atmospheric nitrogen, it is able to oxidize a wide variety of organic compounds, including methane and methanol (Wegner et al., 2020).

Lower vadose zone
The highest LDA scores for the lower vadose zone were assigned to two members of the Betaproteobacteria, Polaromonas and an uncultured genus of the Comamonadaceae. Strains of Polaromonas were initially isolated in Antarctica, but have now been found in many environments. Polaromonas has been reported to be instrumental in the degradation of pollutants, including aliphatic, aromatic, polycyclic and chlorinated hydrocarbons (Mattes et al., 2008). The uncultured member of the Comamonadaceae seems to be related to Rhodoferax, an anaerobic bacterium that respires short chain fatty acids by iron reduction (Finneran et al., 2003). A well-studied iron-reducing and hydrocarbon-degrading bacterium, Geobacter, was also a marker for the lower vadose zone (Lovley et al., 1989). Geobacter is known for synthesizing magnetite during the reduction of Fe(III)-oxyhydroxide (Byrne et al., 2015) and magnetite has been shown to stimulate direct electron transfer (DIET) during syntrophic methanogenesis (Kato et al., 2012;Lovley, 2017). Another biomarker for the lower vadose zone was Methylocystis. The aerobic, methane-oxidizer Methylocystis has been found to be widely distributed in wetlands, peatlands and aquifers. This wide distribution is explained by Methylocystis' extended metabolic capability to use acetate or ethanol in the absence of methane and to fix dinitrogen (Han et al., 2018). Methylocystis has also been detected in groundwater contaminated with heavy metals including arsenic, displaying genes for As(III) oxidation and As(V) reduction (Shi et al., 2018;Danczak et al., 2019). High LDA scores for the lower vadose zone were also assigned to the Alphaproteobacteria Devosia and Rhodopseudomonas. Devosia species have been isolated from diesel-contaminated soil (Ryu et al., 2008), while Rhodopseudomonas is known to be an iron-oxidizer that is able to take up electrons from solid iron and cathodes (Bose et al., 2014).

Upper vadose zone
The two most important biomarkers in the upper vadose zone were members of the Actinobacteria phylum. The first, carrying the designation MB-A2-108, was identified in methane hydrate-bearing, deep marine sediments off the Japanese coast (Reed et al., 2002;Curd et al., 2018). The second was an uncultured member of the order Gaiellales.
Not much is known about this deep-branching phylogenetic lineage, although one Gaiella species has been isolated from a deep mineral water aquifer in Portugal (Albuquerque et al., 2011). Also having a high LDA score in the upper vadose zone was the Chthoniobacterales DA101 soil group of the Verrucomicrobia phylum, which has been found to be extremely common in many soil samples analyzed with molecular methods (Brewer et al., 2016), as well as uncultured members of the Gemmatimonadaceae of the Gemmatimonadetes phylum that are widespread in nature and make up around 2% of soil microbial communities (Zhang et al., 2003). The genus Pseudarthrobacter was also indicative of the upper vadose zone. Pseudarthrobacter species have been isolated from sewage and were found to be able to degrade chlorophenol (Kim et al., 2008).

Analysis of Microbial Metabolism With PICRUSt
The high-throughput DNA sequencing of the 16S rRNA gene resulted in an overview of the microbial community structure at the different depths of the free-phase cores. From this phylogenetic information, we can estimate the metabolic capabilities of the indigenous microorganisms present with a program called PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States; Langille et al., 2013). This program predicts the predominant metabolic pathways based on marker gene survey data (such as the 16S rRNA gene). Although this program was developed for the analysis of the intestinal microbiome, we extracted three pertinent metabolic pathways: methanogenesis, iron reduction and methanotrophy, which are plotted according to depth in Figures 2D,H,L,O. Confirming earlier reports (Bekins et al., 1999;Amos et al., 2005), it is obvious that the free-phase zone is methanic with a majority of methanogenesis genes found below 424.5 m. The majority of the iron reduction genes are present above the methanic zone, from ∼424.5 to 428 m. The majority of the genes that enable microorganisms to oxidize methane were found to overlap with the upper iron reduction zone at ∼426-429 m.

Iron Mineral Analysis
The iron mineral phases were determined for samples from selected depths of core C1407 and compared to the MS (Supplementary Figure 3). The highest MS was measured at 423.2 m, which is inside the free-phase oil zone. Coinciding, the highest concentrations of Fe(II) and magnetite were measured at the same depth. In addition, the FeII/FeIII ratio was the highest at this depth. Reactive and reducible Fe(III) in the form of ferrihydrite and hematite/goethite, respectively, was depleted in the lower depths, but showed higher concentrations between 424.4 and 424.8 m and close to the surface.

Relationship Between Microbial Populations and Magnetic Susceptibility
Our multi-year analysis of the microbial community structure at the aged subsurface hydrocarbon spill site at Bemidji, MN, United States, revealed that a methanogenic microbial community is present in the free-phase petroleum zone. This confirms earlier studies that reported methane emissions at this site (Amos et al., 2005) and a preliminary microbial community analysis based on clone libraries, which showed the presence of methanogenic Archaea (Beaver et al., 2016).
The main methane-generating microorganisms belonged to the genus Methanoregula. These are hydrogenotrophic methanogens, using hydrogen and carbon dioxide to produce methane. The prominent hydrogen provider at this site seems to be the syntrophic Deltaproteobacterium Smithella that was found in all years with an abundance of up to 46% in the free-phase samples.
A comprehensive analysis revealed that Smithella is widespread at different hydrocarbon spill sites all over the world (Gray et al., 2011) and it has recently been shown that Smithella is able to attack long-chain hydrocarbons by fumarate addition and subsequent beta-oxidation (Ji et al., 2020). At the Bemidji site, the free-phase petroleum floats on the water table 9 m below the surface. Above the methanic zone, we found microorganisms that are able to oxidize methane aerobically, mainly belonging to the alphaproteobacterial genus Methylocystis. The same zone was also characterized by the presence of iron-reducing bacteria, specifically Geobacter. A conceptual model is shown in Figure 7. How can we relate these results with the MS values that were measured at Bemidji? First, we have to mention that the high MS values that were measured in year 2011 around the water table in the free-phase petroleum decreased substantially (∼90%) in the following years when measured down the borehole (Lund et al., 2017). On the other hand, here we report that the elevated MS values measured around the water table in the retrieved cores did not substantially decrease over 5 years when measured in the laboratory (Figure 2). The MS values of sediments usually depend on their magnetite content. Magnetite can form biotically, in particular through the action of ironreducing bacteria. Hansel et al. (2005) proposed a model for magnetite formation under iron reducing conditions when there is >2 mM aqueous Fe(II) in the presence of ferrihydrite. Piepenbrock et al. (2011) reported magnetite formation when the ferrihydrite concentration exceeded 5 mM. These values are comparable to the values that we measured in core 1407 (Supplementary Figure 3). It remains enigmatic, however, why the highest MS values were found in the methanic zone around 423/424 m and not in the iron-reducing zone one to two meters above (Figure 2). Possible explanations entail (1) a washing out of magnetic minerals during recharge events from the iron-reducing zone above to the methanic zone below or (2) the ability of the methanogens themselves to switch their metabolism from methanogenesis to iron reduction, resulting in the formation of magnetite. This idea is supported by growing evidence in the literature suggesting that iron FIGURE 7 | Conceptual model of the predominant bacteria and the metabolic processes in and above the free-phase oil zone.
Frontiers in Earth Science | www.frontiersin.org utilization by methanogens mediates the formation of different iron mineral phases (Sivan et al., 2016;Amiel et al., 2020;Shang et al., 2020). Nevertheless, to support our data further, whole metagenome DNA sequencing or RNA-based methods, such as a transcriptome analysis, are necessary.
The presence of magnetite can even be beneficial for methanogenesis. Magnetite was shown to increase the rate of methane production in laboratory experiments with syntrophic bacteria and methanogenic Archaea, because they are capable of direct interspecies electron transfer (DIET) using conductive minerals like magnetite (Kato et al., 2012;Cruz Viggi et al., 2014;Yamada et al., 2015;Yang et al., 2015). Contributing evidence that electron transfer may be occurring between Smithella and Methanoregula is that a self-potential anomaly resulting from a biogeobattery was observed at this site (Heenan et al., 2017). In the Heenan et al. (2017) study, self-potential measurements using a vertical array of electrodes showed a dipolar pattern across the water table fluctuation zone associated with a current source. Such biogeobatteries are postulated to occur in the capillary fringe of contaminant plumes. The large gradient in redox potential at the capillary fringe might allow microbes to transfer electrons from one to another through conductive minerals like magnetite or through biological structures such as nanowires or biofilms, operating in a similar manner to fuel cells (Revil et al., 2010). Interestingly, Methanoregula has been enriched on the cathode in microbial fuel cells when current was applied Zhang and Lu, 2016;Xia et al., 2019), which indicates that it can accept electrons from conductive materials.
Factors contributing to the dissolution of magnetite and therefore the decreasing MS response can be the activities of ironoxidizing bacteria such as Rhodopseudomonas, which was found with a high LDA score in the lower vadose zone. In laboratory experiments, Rhodopseudomonas was shown to oxidize magnetite with and without photosynthesis and it was concluded that magnetite can serve as an electron sink or source (Bose et al., 2014;Byrne et al., 2015). To evaluate the fate of magnetite at hydrocarbon contaminated sites, experiments are currently being carried out at the Bemidji site that study the long-term fate of submerged magnetite minerals at the free-phase petroleum location (L. Slater, pers. communication).
An enigma remains concerning the occurrence of large single MS peaks in the upper vadose zone which were measured at 430 m in two of the 4 years (Figure 2). These outliers could be due to technical errors; nevertheless, this is unlikely because i) at least in 2015, the high MS was represented by more than one measuring point and ii) it was found at the same depth in two different years.

MS as a Tool
Magnetic susceptibility can be used as a tool to identify hotspots of biodegradation at a single point of time as described by Beaver et al. (2016) or to track changes in microbial activity over the long term as reported by Lund et al. (2017). Peaks in MS found in cores retrieved from contaminated sediments can be used for pinpointing the most important biochemical or microbial samples to take. After core retrieval, the wells left behind can be monitored for MS changes to determine the predominant metabolic processes occurring during natural attenuation or bioremediation of hydrocarbon spills.

Conclusion
The microbiological data presented here are the first that allow a spatial and temporal analysis of this spill site because the sediment cores were retrieved in five consecutive years and the samples were collected along the complete length of the cores. The data for the years 2011-2015 confirm previous research findings that the free-phase oil is being degraded by a consortium consisting of the syntrophic Smithella and the methanogenic Methanoregula. The smear zone is a transition zone from methanogenesis to iron reduction, with Thermincola being the most common iron reducer in this zone. Iron cycling and methanotrophy seem to be the dominant processes in the lower vadose zone where iron reducers such as Geobacter, iron oxidizers such as Rhodopseudomonas, and the methanotroph Methylocystis are present (Figure 7). The upper vadose is inhabited by aerobic organisms, mainly belonging to the Actinobacteria and Acidobacteria phyla (Figure 2).
Interestingly, despite the notion that elevated MS values are indicative of higher magnetite concentrations precipitated by the actions of iron-reducing bacterial populations, higher MS values were measured in the methanic free-phase oil zone around the water table and not in the iron-reducing zone above. This finding could be one of the first field evidences that show the metabolic potential of methanogens to switch their metabolism from methanogenesis to iron reduction with concomitant magnetite precipitation. The results presented here provide additional evidence for the role of microbes in environmental magnetism and the magnetic record of terrestrial and marine environments. Further studies are needed in order to determine if changes in MS can be linked to different stages of the bioremediation process.

AUTHOR CONTRIBUTIONS
CB, SR, EA, BB, DN, and LS conceived the study. CB performed the microbiological experiments and analyses. EA, DN, and LS performed the geophysical analysis. CB, SR, and EA wrote the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
The authors acknowledge funding from the Chevron Energy Technology (Grant #CW852844) and the National Crude Oil Spill Fate and Natural Attenuation Research Site, a collaborative venture of the United States Geological Survey (USGS), Enbridge Energy, Limited Partnership, the Minnesota Pollution Control Agency and Beltrami County. Any use of trade, product, or firm names in this publication is for descriptive purposes only and does not imply endorsement by the US Government.

ACKNOWLEDGMENTS
We thank Jared Trost and Andrew Berg (USGS) for field support, and Farag Mewafy, Brittany Ford, Sundeep Sharma, and Drs. Allison Enright and Carl Rosier (Oklahoma State University) for their assistance with geophysical data collection. We also acknowledge Dr. Denise Akob and three reviewers for their helpful suggestions which have improved this manuscript.