Toluene biodegradation in the vadose zone of a poplar phytoremediation system identified using metagenomics and toluene-specific stable carbon isotope analysis

Abstract Biodegradation is an important mechanism of action of phytoremediation systems, but performance evaluation is challenging. We applied metagenomic molecular approaches and compound-specific stable carbon isotope analysis to assess biodegradation of toluene in the vadose zone at an urban pilot field system where hybrid poplars were planted to remediate legacy impacts to an underlying shallow fractured bedrock aquifer. Carbon isotope ratios were compared spatio-temporally between toluene dissolved in groundwater and in the vapor phase. Enrichment of 13C from toluene in the vapor phase compared to groundwater provided evidence for biodegradation in the vadose zone. Total bacterial abundance (16S rRNA) and abundance and expression of degradation genes were determined in rhizosphere soil (DNA and RNA) and roots (DNA) using quantitative PCR. Relative abundances of degraders in the rhizosphere were on average higher at greater depths, except for enrichment of PHE-encoding communities that more strongly followed patterns of toluene concentrations detected. Quantification of RMO and PHE gene transcripts supported observations of active aerobic toluene degradation. Finally, spatially-variable numbers of toluene degraders were detected in poplar roots. We present multiple lines of evidence for biodegradation in the vadose zone at this site, contributing to our understanding of mechanisms of action of the phytoremediation system.


Introduction
Non-halogenated petroleum hydrocarbons, such as benzene, toluene, ethylbenzene, and xylene (BTEX), are common subsurface contaminants at numerous commercial gas stations and industrial sites worldwide (Holliger et al. 1997), and are often released into the subsurface as light, non-aqueous phase liquids (LNAPLs). If release volumes are sufficiently large, LNAPLs can reach the water table, where they float and spread laterally due to their density being less than water, and can result in persistent contaminant plumes (ITRC 2009a). Due to their high vapor pressure, LNAPLs also volatilize and migrate upward into the vadose zone forming vapor phase contaminant plumes (Kim and Corapcioglu 2003;Christophersen et al. 2005). The shallow occurrence of these contaminants in the saturated and vadose zones makes them amenable to phytoremediation (ITRC 2009b). Hybrid poplar trees have been demonstrated to be an effective application of this phytotechnology, being able to contain and even fully remediate petroleum hydrocarbon-contaminated sites (Ferro et al. 2001; Barac et al. 2009;El-Gendy et al. 2009;Landmeyer and Effinger 2016). Plant-associated bacteria play an integral role in remedial activity, especially in the vadose zone, where transpiration draws contaminants and oxygen into the rhizosphere and into contact with degrader populations (Weishaar et al. 2009). However, evaluating activity and overall phytoremediation efficacy is challenging and requires well-informed, effective monitoring strategies.
Metagenomic molecular techniques are powerful emerging tools that can be used to delineate the distribution and activity of marker degradation genes across dynamic field conditions in order to understand the ecology, capacity, and activity of degrader organisms in remediation systems (Larentis et al. 2013). Metabolism of petroleum hydrocarbons by degraders as carbon and energy sources can produce shifts in community composition; thus, differences in proportions of toluene degraders relative to total bacterial abundance across contamination gradients can serve as a line of evidence for biodegradation (Bell et al. 2013).
Tracking of degradation marker genes provides a means for evaluating remedial capability (Winderl et al. 2008;Brow et al. 2013;Key et al. 2014;) and elucidating complex ecological patterns of communities in applied remediation systems (Nebe et al. 2009;Baldwin et al. 2010;Lunsmann et al. 2016). However, to reliably overcome limitations in assay specificity and environmental variability, molecular techniques should be supported by additional lines of evidence (Smith and Osborn 2009;Iwai et al. 2011). Compoundspecific isotope analysis (CSIA) is as an effective tool that can be combined with molecular techniques to offer significant insight into in situ biodegradation processes affecting organic contaminants (Illman and Alvarez 2009). To detect biodegradation, the CSIA method makes use of the preferential cleavage of bonds between light compared to heavy isotopes that leads to a continuous enrichment of heavy isotopes in the parent compared to the daughter compound (Braeckevelt et al. 2012). Initially, CSIA was mainly applied in groundwater systems; however, recent studies have demonstrated its applicability also in the vadose zone (Bouchard et al. 2008a, b;Goli et al. 2017). Molecular techniques have been employed previously in combination with CSIA in both field and lab-scale investigations (Beller et al. 2008;Rakoczy et al. 2011), but have not yet been employed in tandem to study in situ biodegradation in the vadose zone of phytoremediation systems.
To-date, examination of plant-associated biodegradation within poplar-based petroleum hydrocarbon phytoremediation systems has often been conducted on the lab-scale, limited to culture-dependent methods (Jordahl et al. 1997;Barac et al. 2009;Weishaar et al. 2009;Timm et al. 2015), or involved uncontaminated field conditions (Ulrich et al. 2008;Gottel et al. 2011;Beckers et al. 2017). The root microbiome of poplars grown under field conditions is known to harbor a distinct assemblage of organisms compared to associated rhizosphere as a result of environmental differences, rather than opportunistic colonization of roots by dominant rhizosphere organisms (Gottel et al. 2011); however, environmental selective pressures are poorly defined. There is lack of information on the ecological differences between degraders associated with the rhizosphere, poplar root surfaces (rhizoplane), and within the root microbiome in situ. These critical knowledge gaps are significant for evaluating biodegradative processes and are relevant also for evaluating other phytoremediation activities, such as mass removal through phytoextraction (Wilson et al. 2013;Limmer et al. 2018).
In the present study, we assessed the efficacy of a smallscale (220 m 2 planted area) poplar phytoremediation system accessing toluene mass from soil and a shallow, transient, fractured bedrock flow system. Quantitative real-time PCR (qPCR), reverse transcription-qPCR, and CSIA were employed to investigate biodegradation processes affecting toluene in the vadose zone, which operate in conjunction with plant mass-removal processes. Herein, we provide results relating to plant-microbe interactions and their effects on abundance and expression of degradation genes, and on isotopic signatures of biodegradation within the vadose zone. The present study provides a greater understanding, supported through multiple lines of evidence, of the ecology and activities of bacterial toluene degrader populations in phytoremediation systems that is paramount for effective corrective action planning, implementation, and site monitoring.

Site description
The investigated site was a manufacturing facility in a mixed residential/industrial neighborhood in urban South-Western Ontario, Canada. Toluene was used as a solvent in manufacturing processes and was stored in an outdoor storage tank. Toluene impacts to soil, bedrock, and groundwater were discovered in 1989 during decommissioning of the partially buried tank and a buried product distribution line ( Figure 1A).
Site geology includes a 2-m-thick overburden layer of sandy, cobbled soil with remnants of former concrete infrastructure, overlying a 21-m-thick sequence of Silurian-era, fractured dolostone bedrock. A suite of high-resolution field methodologies, collectively known as the Discrete Fracture Network (DFN) approach (Parker et al. 2012), was recently applied to characterize site hydrogeology and the distribution and phase of toluene in the fractured rock and groundwater (Fernandes 2017). Eleven multi-level wells (six to eight monitoring ports) were installed to a maximum depth of 22 m below ground surface (mbgs). Three hydrogeological units (HGUs) were delineated based on vertical hydraulic head profiles. The shallowest HGU, spanning 2-6.5 mbgs, shows a strong vertical hydraulic gradient, and low vertical connectivity to underlying HGUs. The water table fluctuates between 2.1 and 3.0 mbgs, coincident with the bedrock-overburden contact, and groundwater flows in a northern direction showing a velocity of <0.01-3.0 m/day ( Figure 1A).
A primary groundwater source zone of residual toluene was identified in the north-eastern corner of the site ( Figure  1A), wherein toluene occurs as residual LNAPL, and in the dissolved and sorbed phases within the shallow bedrock. Dissolution of LNAPL into groundwater and back-diffusion from the low-permeability dolostone bedrock matrix is sustaining a persistent plume ( Figure 1A). Approximately 95% of dissolved-phase toluene is found within the upper 2 m of the saturated zone, corresponding to the water table depth in deep overburden during high water table conditions, and uppermost bedrock. As such, the main residual toluene mass is located favorably for phyto-remedial uptake and/or biodegradation. Residual impacts to soil have not previously been characterized at this site.
A phytoremediation tree stand of 51 hybrid poplars (Populus deltoides × nigra OP-367) was installed in 2008 in a dense formation overlying the former tank and supply line ( Figure 1A). Mean tree diameter at breast height was approximately 15.97 cm in September 2016 prior to the time of this study. Goals of the phytoremediation application included plume containment and uptake and remediation of toluene in impacted soil, groundwater, and bedrock. Trees were planted in nutrient-amended boreholes to promote rapid establishment and deep rooting; boreholes were advanced to refusal up to a maximal depth of approximately 2 mbgs at the bedrock-overburden contact. This approach has been shown to increase interaction with impacted media and speed up remedial and containment activities (Ferro et al. 2003).

Soil vapor and groundwater sampling
Soil vapor for toluene concentration analysis and CSIA was sampled from three soil vapor probes at the site (SVP1, SVP2, and SVP3; Figure 1A). For concentration analysis, soil vapor was collected monthly between July 11 and September 18, 2017 with 1-L Tedlar ® bags using an SKC XR500 air pump at a pump rate of 450 mL/min with dedicated 0.64 OD Teflon ® tubing. For CSIA, soil vapor was sampled passively with Waterloo Membrane Samplers (WMSs) (Seethapathy and Gorecki 2011). The WMSs were installed for 69 days between July 11 and September 18, 2017. This period was necessary to accumulate enough contaminant mass for CSIA. Groundwater for toluene analysis and CSIA was sampled monthly between July 11 and September 18, 2017, coinciding with soil vapor sampling for concentration analysis, from adjacent multilevel wells (M22, M24, and M29; Figure 1A). Groundwater was sampled from the shallowest port of the multilevel wells adjacent to the soil vapor probes using Geopump peristaltic pumps with dedicated 0.32 cm OD diameter Teflon ® tubing. It was not possible to sample the shallowest port of multilevel M29 on August 13 and September 18, 2017 due to the low water table, therefore the second-shallowest port was sampled instead.

Toluene concentration analysis and compound-specific carbon isotope analysis (CSIA)
Detailed descriptions of toluene concentration analysis and CSIA can be found in the Supporting Information. Briefly, toluene concentrations in the vapor phase and in groundwater were analyzed using a gas chromatograph coupled to a mass spectrometer (GC-MS). Compound-specific stable carbon isotope ratios of toluene were determined by a gas chromatograph coupled to an isotope mass spectrometer (GC-IRMS). Carbon isotope signatures were analyzed relative to the Vienna Pee Dee Belemnite (VPDB) standard and results were expressed in the delta notation: δ VPDB (‰) = (R/R std −1) × 1000, where R and R std are the isotope ratios of the sample and the standard, respectively. The measured δ 13 C values in the vadose zone were corrected by adding 1‰ to account for isotope fractionation during sampling, caused by diffusion through the WMS membrane and sorption in the WMS (Goli et al. 2017).
Rhizosphere soil, rhizoplane soil, and root samples Poplar roots and soil samples were collected over three days in November 2016 from sampling trenches excavated at three locations, upgradient of (ST1 and ST2) and above (ST3), the toluene groundwater source zone ( Figure 1A and B). Weather conditions were reflective of very mild autumn with ambient temperatures of approximately 13-18°C, the ground was not frozen, and trees had not entered senescence. Mechanical advancement of single trenches, versus multiple trenches in a given area, was necessitated by large cobbles and/or buried infrastructure in the overburden material. Trenches were advanced using a CAT ® 420 excavator with an 18-inch bucket in three depth intervals until refusal at final depths of 1.70, 1.93, and 1.37 mbgs (Supplementry data Table S1).
Root samples with adhering soil (ca. 18 g) were collected from each depth (n = 3), then separated into roots and rhizosphere soil by shaking off and collecting adhering soil. Subsamples of soil and roots were flash-frozen in the field by submerging in liquid nitrogen for 3 min, then transported on dry ice and stored at −80°C until further processing. Additional rhizosphere soil subsamples were analyzed for toluene at Maxxam Analytics (Mississauga, Canada) according to EPA 8260 C protocols, and gravimetric moisture content was determined within 48 h.
Rhizoplane samples were obtained as follows in the lab: ca. 6 g of root material was added to 250 mL sterile 1 × phosphate-buffered saline (PBS), vortexed briefly, shaken at 220 rpm at 10°C for 20 min, then rinsed in 100 mL sterile 1 × PBS. Wet soil was recovered from the solution by centrifugation at 13,700 × g for 20 min. Wet soil and surfacecleaned roots were immediately used for DNA extraction to analyze rhizoplane soil and root microbiome bacterial communities, respectively.

Extraction of nucleic acids and cDNA synthesis
The quantity and purity of all extracted DNA and RNA were assessed using a NanoDrop™ 8000 Spectrophotometer (Thermo Fisher Scientific Corp.). Total DNA and RNA were co-extracted from rhizosphere soil (ca. 1.4 g) using an RNA PowerSoil® Total RNA Isolation Kit and DNA Elution Accessory Kit (Mobio Laboratories Inc.) according to the manufacturer's instructions. Extracted DNA used in qPCR was stored at −20°C until analysis. The RNA extracts were immediately DNase-treated using RQ1 RNase-Free DNase (Promega Corp.) in quadruplicate reactions according to the manufacturer's instructions. The RNA was immediately converted to cDNA using a High Capacity cDNA Reverse Transcription Kit (Applied Biosystems Corp.) according to the manufacturer's instructions using 1 µL of SUPERase•In™ RNase Inhibitor (Life Technologies Corp.). DNase treatment efficacy was evaluated by visualizing DNase-treated and untreated RNA on agarose gels and by qPCR analysis.
Surface-cleaned roots (ca. 3 g) were homogenized using an MM400 Mixer Mill (Retsch GmbH) by agitating flashfrozen grinding jars twice at 25 Hz for 20 s. Total root DNA was extracted from root homogenate (0.3 g each in two pooled extractions) using a PowerSoil ® DNA Isolation Kit (Mobio Laboratories Inc.), and gravimetric moisture content of root homogenate was determined.
Total rhizoplane soil DNA was extracted using a PowerSoil ® DNA Isolation Kit (Mobio Laboratories Inc.) with modification to the manufacturer's protocol. DNA was extracted following kit instructions for wet soil samples; centrifugation was increased to 3 min, and the process was repeated until achieving ca. 0.4 g of pelleted material. Pellet weights were later used for standardizing rhizoplane gene quantities.

Quantitative PCR
Total bacteria and toluene-degrading bacterial populations were assessed by quantifying genes and transcripts (cDNA) in the metagenomic samples using previously published primers. The 335f/769r primer set was used to quantify the V3-V4 regions of the bacterial 16S rRNA gene to estimate total bacterial abundance (Dorn-In et al. 2015). Aerobic degradation pathways were targeted using toluene dioxygenase (TOD), ring-hydroxylating monooxygenase (RMO), and phenol hydroxylase (PHE) primers (Baldwin et al. 2003). Anaerobic pathways were targeted using the 7772f/8543r primer set that target benzyl succinate synthase alpha subunit (bssA)-related genes encoding fumarate-adding enzymes (von Netzer et al. 2013). All qPCR reactions were carried out on a CFX96™ Real-Time PCR Detection System (Bio-Rad Laboratories Inc.) and conformed to MIQE guidelines (Bustin et al. 2009). Full details on primers used and reaction conditions can be found in the Supporting Information.
Gene and transcript values were normalized to grams dry soil (rhizosphere soil), dry root, or grams of pelleted extractant (rhizoplane soil) to relay results on a biologically significant scale. Negligible differences in DNA extraction efficiency between samples were assumed. Relative toluene degrader abundances were calculated by standardizing degradation gene quantities to 16S rRNA gene quantities in order to determine if changes in toluene degradation gene abundance reflected changes in bacterial abundance or shifts in community composition.

Statistical analyses
Statistical analyses (ANOVA) were performed to evaluate fixed effects and to compare means of gene and transcript quantities and relative gene abundances using SigmaPlot (v. 12.5) software. Significant differences in means were determined using the Holm-Šídák multiple comparisons test. Where necessary, Box-Cox power transformations (λ = −2 or 0) were performed prior to analysis to ensure normal distribution and equal variance of data. Significance was determine using a threshold of α = 0.05.

Results and discussion
Toluene concentrations in soil, groundwater, and soil vapor Toluene measurements in rhizosphere soil showed a spatial distribution of mass between depths and locations across the site ( Figure 1B). Toluene was detected in soil from all three sampled depths of ST1 (0.025-0.26 µg/g soil ). In ST2, toluene was only detected in the deepest interval (0.045 µg/g soil ). No toluene was detected in soil obtained from ST3, near the groundwater source zone. During trench excavation at ST3, excavation depth was limited to 33-56 cm shallower than the other locations due to buried infrastructure, possibly limiting our ability to sample soil at the groundwater interface, where we expect toluene levels to be highest.
Toluene concentrations in groundwater in the shallowest port of the multilevel wells showed minor temporal variation over the study period ( Figure 2). The lowest concentrations in groundwater (0.79-15.52 µg/L water ) were measured in the inferred upgradient well (M24; Figure 2A). In contrast, concentrations were elevated (94,850-323,200 µg/L water ) in the inferred downgradient direction (M22; Figure 2B) and in the source zone (M29; Figure 2C). When the uppermost port at M29 was dry, toluene concentrations measured in the deeper port were somewhat lower (2683-13,800 µg/ L water ), consistent with the bulk of toluene mass residing at the approximate depth of the water table ( Figure 2C).
Toluene was also detected in vapor collected from all three soil vapor probes installed in the vadose zone, at concentrations significantly lower than in groundwater (0.00279-0.111 µg/L vapor ; Figure 2). The lowest toluene concentrations in vapor over the study period were typically detected under seasonally high water table conditions ( Figure 2). This suggests a potential relationship between exposure of the high toluene concentration zone in the shallow bedrock and vapor concentrations in the vadose zone that merits further investigation.

Stable carbon isotope ratios of toluene
Carbon isotope signatures of toluene in groundwater exhibited small temporal variations during the sampling period ( Figure 3). The δ 13 C values adjacent to the source zone (M29) were in a range between −29.82‰ and −29.26‰ and can be used as an approximation for the carbon isotope composition of the toluene groundwater source ( Figure 3C). The δ 13 C signatures in the upgradient and downgradient multilevel wells (M22 and M24) were consistent with signatures at M29, despite variation in toluene concentration between those locations.
In the vadose zone at the upgradient location (SVP1), a depletion of 13 C in toluene, i.e. a more negative δ 13 C value, was observed just above the water table compared to the groundwater (Δδ 13 C = 1.6‰), whereas at shallower depths the δ 13 C values were similar to the groundwater ( Figure   Figure 2. Measured toluene concentrations in soil vapor probes 1 (A), 2 (B), and 3 (C) in µg/L vapor (upper three data points), and in the shallowest port of the adjacent multilevel wells (M24, M22, and M29) in µg/L water (lowest data point) between July 11 and September 18, 2017. Horizontal lines indicate groundwater levels on sampling dates. The shallowest port in M29 was dry on August 14 and September 18, 2017 due to the low water table, therefore toluene concentrations from the second-shallowest port are presented instead (C). Note that toluene concentrations are plotted on a logarithmic scale to facilitate comparison between groundwater and vapor phase. 3A). In contrast, a less negative δ 13 C value and thus, an enrichment of 13 C in toluene was observed in the vadose zone compared to the groundwater at the downgradient location (Δδ 13 C = 2.1‰) and at the location close to the groundwater source zone (Δδ 13 C = 1.9‰), where high toluene concentrations were detected ( Figure 3B and C). For a system without biodegradation, a depletion of 13 C equaling a more negative δ 13 C value is expected in the vadose compared to the saturated zone, as isotopically light molecules diffuse faster than heavy molecules (Wanner and Hunkeler 2015). Thus, observed enrichment of 13 C in the vadose zone downgradient and adjacent to the source zone in comparison to groundwater provides one line of evidence for microbial toluene degradation in the vadose zone (Bouchard et al. 2008a;Bouchard et al. 2008b). Molecular tools targeting bacterial toluene degradation genes and transcripts were applied to soil samples to provide additional lines of evidence for toluene biodegradation within the vadose zone at this site.

Abundance of toluene degradation genes in poplarassociated rhizosphere soil
Toluene degraders equipped with the genetic capacity for aerobic and anaerobic toluene biodegradation were detected in all rhizosphere soil sampled across the planted area and gradient of detected toluene concentrations in rhizosphere soil (Supplementry data Table S2). Biodegradation genes were detected in the range of 10 5 -10 6 (RMO and PHE), 10 4 -10 5 (TOD), and 10 6 (bssA) copies per gram dry soil. Detected quantities of 16S rRNA genes (10 7 -10 9 copies per gram dry soil) revealed that total bacterial abundance decreased significantly with increasing depth in ST1 and ST3 (Supplementry data Table S2). Relative abundances of degradation genes, as proportions of total bacteria, reflected differences in bacterial community composition across the site, and indicated enrichment of bacterial populations capable of both aerobic and anaerobic toluene degradation at greater depths and toluene concentrations (Figure 4). The effects of sample location and depth on enrichment were evaluated; a significant increase in relative abundances of all toluene degradation genes was observed with increasing depth. Average relative abundances increased between depths 1 and 3 by 0.18% (RMO), 0.13% (PHE), 0.01% (TOD), and 0.82% (bssA) (Figure 4). Relative abundances were not significantly different between locations for TOD and RMO targets ( Figures 4B and C) but changed significantly with depth for bssA and PHE targets (Figures 4A and D). Most notably, PHE gene enrichment was significantly greater in depth 3 of ST1 and ST2 (where toluene was highest; Figure 1B) compared to ST3.
Degrader enrichment is a known marker for contaminant biodegradation (Lovley 2003). Therefore, the overall enrichment of toluene-degrading populations relative to total bacterial numbers and the patterns described above provide a line of evidence for toluene biodegradation in the vadose zone. Toluene is found in vapor phase across the vadose Aerobic transcripts as ratios to gene copies for RMO (E) and PHE (F). Error bars are standard error of replicates. For relative gene abundances, n = 3, except ST3 depth 3, where n = 2. For transcript to gene ratios, n = 3, except RMO (ST1 depth 2, ST2 depth 1, and ST3 depth 1) and PHE (ST3 depth 1), where n = 2. zone ( Figure 2); thus, evidence of biodegradation in rhizosphere soil in which toluene was not detected was not unexpected (Hohener et al. 2006). Biodegradation is likely influencing community composition more strongly at depth, where nutrients are more limited (Jobbágy and Jackson 2001) and where soil toluene concentrations were typically observed to be higher, thus, giving organisms with the capability to biodegrade toluene a competitive advantage. Increased abundance of degrader organisms in poplar rhizospheres has previously been shown to be a result of a general rhizosphere effect of the trees on contaminant-selected soil communities (Jordahl et al. 1997), and enriched degraders in the rhizosphere have been found to drop below detection once contaminant sources are fully removed or remediated (Barac et al. 2009). Therefore, we do not attribute degrader enrichment in rhizosphere soil at this site to specific biochemical and ecological plant-microbe interactions that are involved in other phytoremediation systems (Siciliano et al. 2002;Thijs et al. 2016).
Enrichment of PHE genes most strongly followed patterns of toluene concentrations detected in rhizosphere soil ( Figure 4D). This suggests that organisms equipped with toluene biodegradation pathways that use PHE-encoded enzymes are dominating in locations within the vadose zone where toluene concentrations in rhizosphere soil are highest. Toluene monooxygenase genes phylogenetically related to the tbmD gene of Burkholderia sp. strain JS150, which are amplified by PHE but not RMO primers, may be responsible for observed enrichment of these populations (Nebe et al. 2009). We observed enrichment of organisms having anaerobic bssA genes ( Figure 4A). These results do not, however, necessarily suggest that biodegradation was substantially occurring in the vadose zone under anoxic conditions. Increased oxygen infiltration into the vadose zone is typical in poplar phytoremediation systems due to transpiration activity and resultant depression of the water table (Weishaar et al. 2009). We observed a seasonal flux in the height of the water table in this manner below the poplar stand (Figures 1 and 2). Moreover, subsequent preliminary sampling of the three soil vapor probes has demonstrated that the vadose zone exists under oxic conditions (data not shown). It is therefore possible that facultative anaerobic toluene degraders possessing both anaerobic and aerobic toluene degradation pathways, such as Thauera sp. strain DNT-1 (Shinoda et al. 2004), were being enriched.
Toluene degradation gene expression in poplarassociated rhizosphere soil Active expression of degradation gene targets in rhizosphere soil was evaluated by quantifying mRNA transcripts of RMO and PHE gene targets to determine whether these genes were simply present, or indeed actively being utilized. Aerobic toluene biodegradation genes were actively being transcribed at the time of sampling, reinforcing that biodegradation is occurring in the vadose zone of the phytoremediation system (Supplementry data Table S2). Transcript quantities ranged from non-detect to 10 4 (RMO) and 10 2 -10 5 (PHE) copies per gram dry soil. Transcripts were detected at all locations, except for RMO in depth 3 of ST3, where genes were potentially being transcribed at levels below our limit of detection. On average, RMO transcripts differed significantly between locations, and were more abundant in ST1 and ST2, where rhizosphere soil toluene concentrations were highest. However, PHE transcripts did not follow the same trend, despite PHE gene enrichment at these locations ( Figure 4D). Expression of BTEX degradation genes has previously been shown in contaminated soils to be highest at locations with the greatest availability of substrate and/or electron acceptor (Key et al. 2014). However, degradation gene expression is reflective of metabolic processes that are transient in nature, with numerous fluctuating environmental factors, such as nutrient and oxygen status, pH, moisture, and temperature affecting associated transcript quantities (Saleh-Lakha et al. 2005). Therefore, expression patterns may be reflecting conditions at the time of sampling where metabolic activity and therefore gene induction were low. Relative transcript quantities (as proportions of corresponding gene quantities) were evaluated in an attempt to identify zones where high-level expression might be occurring (Brow et al. 2013), but did not show distinct patterns ( Figure 4E and F). Induction of degradation gene expression occurs in the presence of toluene; but terminal electron acceptors are occasionally sufficient inducers, even in the absence of toluene (Brow et al. 2013). Moreover, both aerobic and anaerobic genes can be expressed at low levels in preparation for changing environmental conditions that favor biodegradation, such as temperature, oxygen status, and substrate concentration (Shinoda et al. 2004). Oxygen concentrations were not measured at the time of sampling, but oxygen has been detected throughout the vadose zone in preliminary sampling of the three soil vapor probes (data not shown). Thus, PHE gene targets may have had sufficient inducers of expression available across the site, while induction of RMO targets may be more dependent on access to toluene. Together, these results provide evidence of a dynamic response of bacterial toluene degrader populations in the vadose zone to environmental conditions and toluene concentrations and support the finding that enrichment of degraders in the rhizosphere is a result of biodegradation.
Abundance of toluene degradation genes in poplar roots and root-associated rhizoplane soil Toluene degradation gene abundances (RMO, TOD, and PHE) in rhizoplane soil and root microbiomes indicated that degrader communities changed across the soil-root interface at this site (Supplementry data Table S3). Quantities of 16S rRNA genes ranged from 10 8 to 10 9 copies per gram wet rhizoplane soil and 10 9 -10 10 copies per gram dry root. Degradation gene quantities per gram wet rhizoplane soil ranged from 10 5 to 10 6 (RMO), 10 4 -10 5 (TOD), and 10 5 -10 7 (PHE). Gene copies per gram dry root ranged from 10 6 to 10 7 (RMO) and 10 5 -10 6 (TOD). PHE genes were present in all root samples (10 6 -10 7 ); however, quantification was limited by a secondary PCR product (∼100 bp) that amplified in some samples. Off-target plant DNA can sometimes interfere with metagenomic study of bacteria in the root microbiome (Arenz et al. 2015); therefore, specificity of PHE primers merits further examination when looking at root-associated communities. Degrader abundances did not differ significantly between roots where toluene was detected in associated rhizosphere soil compared to roots with rhizosphere soil having no detected toluene (data are presented in the Supporting Information). All roots are considered toluene-exposed due to vapor-phase toluene in the vadose zone (Figure 2), and the small-scale nature of the site precluded access to non-exposed control roots.
The abundance of RMO and TOD genes relative to 16S rRNA genes ( Figure 5) showed changes in the composition of toluene degrader populations associated with poplar roots. Relative degradation gene abundance changed significantly with averages of 0.304% (RMO) and 0.0165% (TOD) in rhizoplane soil, and 0.471% (RMO) and 0.0217% (TOD) in roots. Relative TOD abundance differed by sampling location, illustrating significantly different enrichment patterns between and within rhizoplane soil and roots ( Figure 5B and E). A greater relative abundance of TOD genes in roots in depth 3 was observed in ST3 compared to ST1 and ST2 ( Figure 5E). These results suggest that distinct degrader assemblages are populating poplar root microbiomes across the site, and that degraders are more active and are being enriched in deeper roots of trees overlying more heavilyimpacted groundwater in the toluene groundwater source zone. Given that poplars at this site are of the same genotype and age, access to toluene is likely the dominant selective pressure on community composition in their root microbiomes (Ulrich et al. 2008). Efforts to measure in planta concentrations of toluene and degrader abundances in above-ground tissues along transpiration streams of trees at this site are currently underway. Wilson et al. (2013) previously showed that biodegradation prior to plant uptake directly impacts the concentrations of BTEX in poplar transpiration streams; however, the influence of biodegradative processes along the transpiration stream has remained relatively unexplored and merits further investigation.

Conclusions
In this study, we applied molecular and stable isotope methods as lines of evidence supporting biodegradation of toluene in the vadose zone of a 9-year-old poplar phytoremediation system established over a shallow, fractured bedrock aquifer contaminated with toluene. Our findings corroborate previous evidence in the literature that a general rhizodeposition effect influences toluene degradation gene abundance in the rhizosphere, and that access to contaminants strongly affects community composition. Our results indicate that organisms inhabiting poplar rhizosphere and root microbiome compartments play a role in biodegradation in the vadose zone. Finally, we provide evidence herein suggesting that tree proximity to impacted Figure 5. Relative gene abundances of toluene-degrading bacteria associated with poplar rhizoplane soil (A-C) and roots (D and E). Relative gene abundances as proportions of bacterial 16S rRNA genes (%) in rhizoplane soil and roots for RMO (A, D) and TOD (B, E). Note, PHE relative abundances in roots are not reported due to co-amplification of nonspecific targets in root samples. Error bars are standard error of replicates (n = 3, except rhizoplane soil from ST1 depth 3, where n = 2). groundwater may be shaping degrader community structure in poplar roots at greater depths, and we propose that strategies monitoring phytoextraction of petroleum hydrocarbons should also investigate the degradative capacity of endophytic communities along transpiration streams. Together, degrader enrichment in rhizosphere and roots highlight how phytoremediation systems couple contaminant uptake with enhanced biodegradation for broader remedial activity. This research provides improved knowledge of the ecological structures of poplar-associated degraders and will aid in future efforts to enhance phytoremediation application and monitoring strategies, and to quantify processes to inform remedial time-scales.