Cold Lake Winter Blend dilbit, sourced from pipeline stocks provided by the Canadian Association of Petroleum Producers (CAPP), was weathered on June 19–20, 2019 for 36h by exposing 9 kg of oil to sunlight and air movement in a 1.1 m diameter stainless steel pan over ∼220 L of lake water, as described in Palace et al. (2021). Weathering was performed to simulate oil that would typically reach the shoreline environment after a spill event. During the weathering process, oil undergoes evaporative loss of lighter compounds (e.g., BTEX), altering its chemical and physical characteristics, including an increase in viscosity and density that make recovery of an aquatic spill more difficult (Lee et al., 2015). After weathering, the oil was collected from the surface of the water in the weathering pan using slotted stainless steel spoons and 1,444 g was applied in the EE at the near shore (within 50 cm of the bank) environment on 21 June 2019. Shorelines were treated 4 days after oil addition (25 June 2019) to allow floating oil to strand naturally on the shoreline and to simulate conservative response times in remote locations. The EE received primary recovery of shoreline washing for 12 minutes through pumping enclosure water (1,600 L) over the shoreline with a low pressure manifold and collecting oil on the water surface with polypropylene sorbent media (Spill Ninja, MEP Brothers, Winnipeg, MB, Canada). Cleanup treatment followed response industry standards for wetland (vegetated) shorelines. Primary recovery removed 49.4 g of dilbit, leaving 1,394 g (97%) of residual oil on the shoreline and in the enclosure. The selected shoreline treatment method was non-invasive and did not include any attempts to scrub oil from surfaces or activities that would alter the shoreline and native vegetation. Residual oil was then treated in the EE on 26 June 2019, with secondary remediation of an EFW located in the middle to far shore environment of the enclosure. The RE also received shoreline washing on 27 June 2019, followed by EFW addition. In order to conduct such research, contingency measures were implemented to minimize further impact to the lake ecosystem beyond the experimental regions (Supplementary Material, Section 1.1).

EFWs were BioHaven Floating Islands sourced and produced by Martin Ecosystems (LA, United States), under license by Floating Islands International (MT, United States). The platform was made of recycled polyethylene terephthalate plastic and coated in polyurea for protection, with holes to transplant emergent vegetation (Floating Islands International, 2024). EFWs (1.14 m²; 21 holes of 7.62 cm diameter) were originally established in 2017 with local emergent wetland plants, Gold Label coconut mulch, and soil, and re-transplanted in June 2019 to replace dead plant material with coconut mulch and Typha sp. (cattail) and three Carex spp. (sedge); C. atherodes/utriculata, C. lasiocarpa, and C. pseudocyperus/hystericina (identified by Chris Penner, personal communication with Dr. Richard Grosshans, IISD, n.d.). Noting, Carex are particularly difficult to identify due to the high number of species, classified by minor changes in plant features (Naczi, 1992, as cited in Starr et al., 2009). EFWs plants were measured for heights and number of live plants, however, due to insect activity, measurements are not discussed herein.

Water was collected 3 days pre-oil addition and on days 1 through 6, 8, 12, 20, 38, 66, 87, 209, 263, 356, and 411 post-oil addition for PAC chemistry (where day 0 was the date of oil addition). Total PACs (n = 44) included 16 United States Environmental Protection Agency (2014) priority parent PACs and 28 alkylated PACs (which includes parent and alkylated heterocyclic PACs) (Supplementary Table 1). Water was collected in 1 L amber bottles using a peristaltic pump from the far shore environment in the enclosures (∼0.5 m from enclosure end and collecting approximately halfway through the water column) and transported in a cooler to the Centre for Oil and Gas Research and Development (COGRAD) laboratory at the University of Manitoba (Winnipeg, MB, Canada). Within 24 h of collection, 250 mL of each sample was vacuum filtered (1.2 μm, Whatman GF/C filter) and liquid-liquid extraction was completed following methods previously described in Dearnley (2022). Briefly, filtrates were transferred into separatory funnels which received 20 μL of 5 ng/μL recovery internal standard (suite of d8-naphthalene, d8-acenaphthylene, d10-acenaphthene, d10-fluorene, d10-phenanthrene, d10-pyrene, d12-benz(a)anthracene, d12-chrysene, d12-benzo(b)fluoranthene, d12-benzo(k)fluoranthene, d12-benzo(a)pyrene, d12-indeno(1,2,3-c,d)pyrene, d14-dibenzo(a,h)anthracene, and d14-benzo(g,h,i)perylene) and 10 g of sodium chloride (NaCl). Samples were extracted twice with 50 mL dichloromethane (DCM) by gently swirling liquids for 1 min. DCM was collected into round bottom flasks (total 100 mL DCM) and was reduced to approximately 2 mL using rotary evaporation. Sodium sulphate (Na₂SO₄) was added to remove any residual water from the extraction process (Dearnley, 2022) and samples were transferred to 12 mL round bottom borosilicate glass tubes with a Teflon lined screw caps. The remaining sodium sulphate was rinsed three times with Optima ™ grade hexanes (ThermoFisher Scientific), which was then transferred into the glass test tubes. Solvent exchange from DCM to Optima ™ hexanes was conducted using a nitrogen gas evaporator (N-EVAP™111, Organomation Associates Inc., MA, USA and OA-SYS Heating System), until 1 mL of sample remained. Each sample was spiked with 20 μL of an instrument performance internal standard (5 ng/μL d10-anthracene) and was transferred into an amber GC vial with a Teflon cap. PACs were detected and quantified on an Agilent 7890 gas chromatograph coupled to a 7000C triple quadrupole mass spectrometer (GC-MS/MS) following methods and GC-MS/MS conditions outlined in Idowu et al. (2018).

During each sampling period, a field blank was collected by opening an amber bottle containing 200 mL Optima ™ water (ThermoFisher Scientific) for the duration of sample collection. PAC extraction of field blanks followed the same protocols described above for experimental samples. Final concentrations were calculated using internal standard recoveries and blank correction. Of all samples analyzed, GC-MS/MS runs for two samples (EE and RE) on day 209 failed and are not included.

Water was also regularly monitored for basic quality parameters and was sampled bi-weekly for nutrient chemistry, discussed in Supplementary Material, Sections 1.2 and 1.3.

Sediment cores were collected with a sediment corer (4.6 cm diameter x 30 cm long clear polycarbonate tube) in three random locations (n = 3) from each site 3 days pre-oil addition and on days 38, 66, and 91 post-oil addition. Single sediment samples were collected on days 263 and 642. Sediment samples were collected in triplicate for PAC chemistry, unlike the single water samples, to capture spatial heterogeneity of PAC concentrations resulting from variation in oil submergence and/or sediment characteristics. The top 5 cm of each sediment core was transferred into a 125 mL amber jar and refrigerated overnight to remove overlying water once particles had settled and then frozen at −20°C. Sediment PACs were extracted at the COGRAD laboratory using accelerated solvent extraction methods described in Dearnley (2022) and quantified on GC-MS/MS as described above (Idowu et al., 2018).

Root samples were extracted using the PowerSoil Max DNA extraction kit (cat#12988–10, QIAgen). Between 5 and 10 g of roots were weighed and transferred into a PowerMax Bead Tube containing 15 mL of PowerBead Solution. The extraction was done following manufacturer’s guidelines. The nucleotide acid elution was done with 3 mL of DNAse-free Ambion water. The extracted genomic DNA (gDNA) was quantified using the Quant-iT™ PicoGreen™ dsDNA Assay Kit (cat#. P7589, ThermoFisher Scientific).

Taxonomic 16S (prokaryotic) and 18S (eukaryotic) ribosomal ribonucleic acid (rRNA) gene amplifications were done on all the gDNA samples. Bacterial communities, including eubacteria and archaea, were targeted using the primer set 515F-Y (5′- GTGYCAGCMGCCGCGGTAA-3′)/926R (5′-CCGYCAATTYMTTTRAGTTT-3′). The 18S rRNA gene region of the eukaryotic ribosomal DNA was amplified using the primer set 565F (5′- CCAGCASCYGCGGTAATTCC-3′)/948R (5′- ACTTTCGTTCTTGATYRA-3′). Primers contained the required Illumina adaptors at the 5′end of the primer sequences (5′-TCG​TCG​GCA​GCG​TCA​GAT​GTG​TAT​AAG​AGA​CAG-3′ for the forward primer and 5′-GTC​TCG​TGG​GCT​CGG​AGA​TGT​GTA​TAA​GAG​ACA​G-3′ for the reverse primer).

The V4-V5 region of the bacterial 16S rRNA gene and the V4 region of the eukaryotic rRNA gene were amplified in 25 μL volumes using 12.5 μL of KAPA HiFi HotStart ReadyMixPCR Kit (Roche), 10 μM of each primer, Ambion nuclease-free water, 20 μg of BSA and 1 μL of gDNA. Thermal cycling conditions were as follows: initial denaturation at 95°C for 3 min; 30 cycles at 95°C for 30 s, 55°C (for 16S primers) or 50°C (for 18S primers) for 30 s, 72°C for 30 s; and a final elongation at 72 °C for 7 min. PCRs products were visualized on SyBrSafe-stained 1% agarose gels. PCR amplicons were purified using 0.8X volume of magnetic beads solution (Agencourt AMPure XP, Beckman Coulter Life Science) according to the manufacturer’s protocol. Unique codes were added to each sample by amplifying 5 μL of the purified PCR product with 12.5 μL of KAPA HiFi HotStart ReadyMixPCR Kit, 150 nM of each Nextera XT Index Primer (Illumina Inc., San Diego, CA, USA) and Ambion nuclease-free water for a total volume of 25 μL. Thermal cycling conditions were as follows: 3 min at 95°C, eight cycles of 30 s at 95°C, 30 s at 55°C, 30 s at 72°C, and a final elongation step of 5 min at 72°C. Indexed amplicons were purified with the magnetic beads as previously described and quantified using the picogreen fluorescence method and combined in an equimolar ratio.

Paired-end sequencing (2 × 250 bp) of the 16S and 18S rRNA gene amplicon pools was carried out on an Illumina MiSeq sequencer at the National Research Council (Montreal, QC, Canada) using the MiSeq Reagent kit V2, 500 cycles (cat# MS-102–2003, Illumina).

Sequencing data was analyzed using AmpliconTagger (Tremblay and Yergeau, 2019). Briefly, raw reads were scanned for sequencing adapters and PhiX spike-in sequences. Remaining paired-end reads were processed to remove primer sequences (pTrimmer v1.3.3; Zhang et al., 2019) and discard reads having an average quality Phred score lower than 20. The remaining sequences were processed for generating Amplicon Sequence Variants (ASVs) (DADA2 v1.12.1; Callahan et al., 2016). Since the quality filtering step was performed in a separate upstream step, we used more lenient parameters for the DADA2 workflow which is summarized as follows: filterAndTrim (maxEE = 8, truncQ = 0, maxN = 0, minQ = 0). Errors were learned using the learnErrors (nbases = 1e8) function for both forward and reverse filtered reads. Reads were then merged using the mergePairs (minOverlap = 10, maxMismatch = 0) function. Chimeras were removed with DADA2’s internal remove BimeraDeNovo (method = ‘consensus’) method followed by UCHIME reference (Rognes et al., 2016). ASVs were assigned a taxonomic lineage with the RDP classifier (Wang et al., 2007) using the complete SILVA release 138 database (Quast et al., 2012) limited to 16S and 18S entries. Taxonomic lineages were combined with the ASV abundance matrix obtained above to generate a raw ASV table, from which the prokaryotic or eukaryotic ASV tables were generated. A total of 6,365 and 3,263 ASVs were obtained for 16S and 18S data types respectively.

Total dissolved-phase PAC concentrations before oil addition were 290.21 ng/L, 27.38 ng/L, and 57.90 ng/L in the EE, RE, and LR, respectively (Figure 1). After oil addition to the EE, concentrations on day 1 were 1,591.23 ng/L, 95.84 ng/L, and 63.33 ng/L, respectively. The maximum concentration in the EE occurred on day 20, reaching 12,440.83 ng/L, within the chronic (>96 h) toxicity range for total polycyclic aromatic hydrocarbons (PAH) causing mortality (LC₅₀) for freshwater biota (Lee et al., 2015). Total PACs on day 20 were 1,667.17 ng/L and 795.11 ng/L in the RE and LR, respectively. This peak was estimated to result from a precipitation event on day 19, a record rainfall of 145.6 mm in 24 h (IISD-ELA meteorological data), causing remobilization and reintroduction of oil from the shoreline and/or sediments. This was supported by the declining trend in total dissolved-phased PAC concentrations on days 8 and 12 after oil addition. Concentrations in the EE then declined to 2,336.03 ng/L on day 38, and later to 359.91 ng/L, nearing background conditions on day 66. Total PACs were relatively stable for the remainder of the experiment, where concentrations were last recorded on day 411 at 467.97 ng/L, 202.56 ng/L, and 119.51 ng/L in the in EE, RE, and LR, respectively. Figure 1 presents the first 87 days of exposure. The remaining data are presented in Supplementary Figure 2.

Total PACs in the EE were primarily 2-, 3- and 4-ring alkylated PACs, with fewer ≥5-ring compounds (Supplementary Figures 3, 4), and with dominance (mean >100 ng/L) of alkylated naphthalenes (C2, C3, and C4), dibenzothiophenes (C1 and C2), C2-pyrene, and C1-fluorene. Following oil addition, alkylated PACs had a mean (range; n = 16) of 90.65% (66.79%–99.28%), 70.56% (16.89%–98.49%), and 60.06% (4.28%–99.87%) of total PACs in the EE, RE and LR, respectively. Alkylated PACs are more persistent compared to parent compounds, have greater chronic toxicity to biota (Zhu et al., 2001; Lee et al., 2015), and increase with continued weathering and evaporative processes (Lee et al., 2015). Total PACs were higher than previous oil spill research at the IISD-ELA (Palace et al., 2021; Rodriguez-Gil et al., 2021), which may be a result of the oil product, PACs measured, analytical methods (Rodriguez-Gil et al., 2021), remediation efforts, site volume, and environmental conditions. For example, during this study we experienced an extreme precipitation event which may have caused redistribution and mixing of oil into the water column from the oil impacted shoreline, whereas in previous studies oil had sunk to the lake bottom (Stoyanovich et al., 2021) with potential accumulation in sediment through deposition, partitioning, and/or by direct contact (Stoyanovich et al., 2022). This illustrates the limitations of comparing to prior studies and past spill events.

Significant variability was evident among total sediment PAC concentrations within replicates and there were high concentrations observed in unoiled sites, making it difficult to discern differences between sites over the exposure period. This was in part due to significant contributions from retene, a three-ring PAC, associated with combustion of wood (i.e., forest fires) or diagenic processes (Ramdahl, 1983; Gabos et al., 2001). Forest fires are a regular and essential disturbance for boreal forest rejuvenation (Brandt et al., 2013), and the study region is often impacted by nearby fires resulting in ash deposition and smoke, and also has signs of a historic fire in the watershed (1979; personal communication with Paul Fafard, IISD-ELA, February 2023). In all sampling periods, percent retene of total PACs ranged from 0.00% to 92.29%, 22.03%to96.80%, and 28.04%to97.68% in the EE, RE, and LR, respectively. The differences in ranges between sites suggest the EE may have other compounds contributing to total PACs in some replicates and/or sample rounds, potentially those of petrogenic origin. To assess changes in sediment PACs between sites, we reanalyzed PAC data excluding retene from total PAC concentrations (n = 43) (Figure 2) (trends of the complete dataset, including retene is included in Supplementary Figure 5).

Mean total PACs (±SD; triplicate samples at each site) excluding retene (n = 43) during pre-oil conditions were 10.46 ± 1.26 μg/g, 10.86 ± 3.94 μg/g, and 5.78 ± 3.76 μg/g in the EE, RE, and LR, respectively, with large contributions from C1-pyrene (14.50%–44.05%). Concentrations increased with high variability on day 38 in the EE to 42.64 ± 54.65 μg/g, with large contributions from C3-dibenzothiophene (4.54%–20.59%), C2-dibenzothiophene (5.31%–19.13%), and C1-pyrene (5.68%–15.50%), however the largest change occurred on day 66 at 528.39 ± 909.17 μg/g, mostly of C4-phenanthrene (1.72%–87.47%) and C1-pyrene (6.54%–23.58%) in some samples. The variability among replicates indicates potential influence of sample location, similarly observed in Stoyanovich et al. (2022). Any oil globules contained in the EE sediment samples would have contributed to increased PACs and replicate variability. Sediment PACs in the RE were 18.17 ± 20.60 μg/g and 4.23 ± 3.44 μg/g on day 38 and 66, respectively, and were 10.36 ± 1.50 μg/g and 9.76 ± 3.37 μg/g in the LR, respectively, mostly C1-pyrene. Concentrations decreased in the EE on day 91, at 339.42 ± 200.17 μg/g (n = 2), almost fully contributed by C1-pyrene (98.89%–99.00%), with little contribution from alkylated dibenzothiophenes or phenanthrene, suggesting potential degradation, loss, or influence of sample location. Concentrations in the RE and LR were 34.77 ± 27.62 μg/g and 103.93 ± 167.73 μg/g, respectively. By day 263, concentrations were similar amongst all sites at 30.31 μg/g, 31.30 μg/g, and 36.61 μg/g, respectively (88.65%–93.81% C1-pyrene) (Supplementary Figure 6). On day 642, sediment PACs remained at 31.80 μg/g in the EE, while the RE and LR had 19.16 μg/g and 1.25 μg/g, respectively. It is clear there is a large contribution of C1-pyrene in the sediment, with and without oil exposure, and it is possible C1-pyrene and other sediment PACs, may have resulted from forest fires (Vergnoux et al., 2011), however the increase on day 38 and 66 in the EE, is likely contributed from dilbit.

Despite background PACs skewing interpretation of total PAC trends between oiled and unoiled sites, increases in mean trends of parent and alkylated PACs in the EE indicate potential influence from dilbit (Supplementary Figure 5). However, there was large variability captured within site replicates, which is likely a result of spatial heterogeneity due to sediment characteristics and/or oil submergence and PAC sorption. To assess whether PACs arose from petrogenic or pyrogenic sources, sediment PAC ratios are discussed in Supplementary Material, Section 1.4.

Root microbial 16S (prokaryotes) and 18S (eukaryotes) rRNA gene Bray-Curtis beta diversity of rarefied ASVs was analyzed using a PCoA to identify differences between EFWs over the experimental period. The PCoA for prokaryotic and eukaryotic beta diversity had 40.13% and 41.28% of explained variation on the first two axes, respectively (Figure 3). The communities of all three EFWs on day −3 were clustered together on both PCoAs, indicating similarity before EFW deployment to their respective sites. The greatest differences in the prokaryotic community were between the LR and enclosures on days 38 and 91, which were clustered together, suggesting there may be an enclosure effect. Differences over the exposure period were observed for the enclosure sites, with greater differences between sites on day 91, indicating a potential seasonal or condition effect. There were no changes in the eukaryotic beta diversity in the LR, while there were changes in the enclosures over the exposure period, with greatest differences between the enclosure sites on day 91, displaying a potential site condition and/or seasonal effect.

Prokaryotes had greater richness, and diversity and evenness than eukaryotes, suggesting there may be no dominating prokaryotic organisms, while certain eukaryotic organisms may be dominant or selected for. While highly variable, the enclosures appear to have greater mean richness than the LR on days 38 and 91, but diversity indices are variable, further suggesting a potential enclosure effect (Supplementary Table 3). The enclosure design may limit wind and wave aeration, nutrient cycling, and other lake processes (see Supplementary Material, Section 1.2) that may influence microbial diversity and abundance. Similar enclosure effects on the aquatic prokaryotic community were identified in the Freshwater Oil spill Remediation Study which took place at the IISD-ELA studying the use of enhanced Monitored Natural Recovery and shoreline washing agent as secondary remediation methods following dilbit spills (Kharey et al., 2024).