Legacy Effects of Phytoremediation on Plant-Associated Prokaryotic Communities in Remediated Subarctic Soil Historically Contaminated with Petroleum Hydrocarbons

Subarctic ecosystems provide key services to local communities, yet they are threatened by pollution caused by spills and disposal of petroleum waste. Finding solutions for the remediation and restoration of subarctic soils is valuable for reasons related to human and ecosystem health, as well as environmental justice. ABSTRACT Phytoremediation of petroleum hydrocarbons in subarctic regions relies on the successful establishment of plants that stimulate petroleum-degrading microorganisms, which can be challenging due to the extreme climate, limited nutrients, and difficulties in maintaining sites in remote locations. A long-term phytoremediation experiment was initiated in Alaska in 1995 with the introduction of grasses and/or fertilizer to petroleum hydrocarbon (PHC)-contaminated soils that were subsequently left unmanaged. In 2011, the PHC concentrations were below detection limits in all soils tested and the originally planted grasses had been replaced by volunteer plant species that had colonized the site. Here, we sought to understand how the original treatments influenced the structure of prokaryotic communities associated with plant species that colonized the soils and to assess the interactions between the rhizospheric and endophytic communities of the colonizing vegetation 20 years after the experiment was established. Metataxonomic analysis performed using 16S rRNA gene sequencing revealed that the original type of contaminated soil and phytoremediation strategy influenced the structure of both rhizospheric and endophytic communities of colonizing plants, even 20 years after the treatments were applied and following the disappearance of the originally planted grasses. Our findings demonstrate that the choice of initial phytoremediation strategy drove the succession of microorganisms associated with the colonizing vegetation. The outcome of this study provides new insight into the establishment of plant-associated microbial communities during secondary succession of subarctic areas previously contaminated by PHCs and indicates that the strategies for restoring these ecosystems influence the plant-associated microbiota in the long term. IMPORTANCE Subarctic ecosystems provide key services to local communities, yet they are threatened by pollution caused by spills and disposal of petroleum waste. Finding solutions for the remediation and restoration of subarctic soils is valuable for reasons related to human and ecosystem health, as well as environmental justice. This study provides novel insight into the long-term succession of soil and plant-associated microbiota in subarctic soils that had been historically contaminated with different sources of PHCs and subjected to distinct phytoremediation strategies. We provide evidence that even after the successful removal of PHCs and the occurrence of secondary succession, the fingerprint of the original source of contamination and the initial choice of remediation strategy can be detected as a microbial legacy in the rhizosphere, roots, and shoots of volunteer vegetation even 2 decades after the contamination had occurred. Such information needs to be borne in mind when designing and applying restoration approaches for PHC-contaminated soils in subarctic ecosystems.

S ubarctic regions are often impacted by fuel and oil spills associated with petroleum extraction, transportation, storage, and use (1)(2)(3)(4)(5). This release of petroleum hydrocarbons (PHCs) into soils and/or water can pose exposure risk to humans, impair environmental health, and disrupt diverse microbial communities and their ecological functions in the environment (6,7). Since natural attenuation of contaminants proceeds more slowly in cold regions than in warmer areas, it is not always a suitable approach for their restoration (8). Therefore, the development of appropriate remediation strategies is of great environmental significance in these regions.
Phytoremediation, the use of plants and plant-associated microorganisms to remediate contaminated areas, can be a cost-effective and sustainable strategy for the ecological restoration of PHC-polluted sites (9,10). Unlike commonly used physical and chemical remediation strategies, phytoremediation does not necessarily cause further disturbance of the landscape and can help re-establish vegetation in the ecosystem (11). To date, there have been a limited number of studies examining the phytoremediation of PHC-contaminated sites at high latitudes (12)(13)(14)(15)(16). The short growing season, low precipitation, and nutrient limitation due to slow chemical weathering and biological decomposition dramatically limit the productivity of plants and their associated microbes in some cold regions (8,17). Moreover, high concentrations of PHCs can impair the water movement in soil and gas exchange between soil and air. Such processes can restrict plant growth and lower the activity of soil microorganisms, thereby decreasing soil health (18)(19)(20). Therefore, identifying plant species that can thrive under these extreme conditions and successfully remediate PHC-contaminated soils is important to the development of successful phytoremediation strategies for subarctic regions.
(Re)introduction of plants to contaminated areas has been shown to promote ecosystem recovery by increasing soil health through the improvement of soil structure and the accumulation of organic carbon and limiting nutrients, such as nitrogen or phosphorus (21,22). Microorganisms living in the rhizosphere and endosphere play important roles in these processes (23)(24)(25). Both the rhizosphere, the narrow zone of soil directly influenced by the roots, and the endosphere, the plant interior, have been shown to harbor bacteria that often stimulate plant growth through fixation of atmospheric nitrogen, solubilization of inorganic phosphate, or production of phytohormones and siderophores or help to protect their host plants by warding off phytopathogens and herbivores (25,26). In addition, some endophytes have been found to mitigate the impacts of contaminant-associated stress by directly degrading pollutants (27,28). Overall, plant-associated bacteria enable their hosts to adapt to changing environmental conditions, and their beneficial properties could potentially be even more important under the harsh conditions that prevail in subarctic regions (29,30). While the composition of plant-associated bacterial communities is generally influenced by a variety of factors, such as plant species, soil characteristics, or climate (25), the exact mechanisms that drive the assembly of plant microbiota are likely to be a combination of factors and may vary across different environments.
In this study, we built upon a long-term phytoremediation experiment with PHC-contaminated soils established in Fairbanks, Alaska, in 1995 (3,(31)(32)(33). The original study determined how the introduction of different grass species, fertilization, or their combination would affect the remediation of PHC-contaminated soils over time. We re-examined the field site after it had experienced no active site management for almost 20 years. By that time, the PHC concentrations were below detection limits in all soils and the originally planted grasses had already disappeared and had been replaced by volunteer vegetation depending on the originally applied phytoremediation strategy (12). Here, we hypothesized that the legacy of previously applied bioremediation strategies would still be detectable in microbial communities associated with volunteer vegetation (Fig. 1). Together, our investigations of this research site have helped to uncover how different factors contribute to both plant and microbial secondary succession in PHC-contaminated subarctic environments.

RESULTS
Diversity and structure of plant-associated prokaryotic communities. The prokaryotic diversities in the rhizosphere and both the root and shoot endospheres of colonizing plants ( Fig. 2A). Both Shannon and Simpson indices differed significantly among the plant-associated microbial habitats studied (P # 0.001, Kruskal-Wallis test); the prokaryotic diversity decreased in the order of rhizosphere, root endosphere, shoot endosphere. Nonmetric multidimensional scaling (NMDS) using weighted UniFrac distances revealed that the communities in the rhizosphere and shoot endosphere clustered separately (Fig. 2B), indicating that rhizospheric communities differed from those found in the shoots, while the communities in the roots were more similar to both rhizospheric and shoot endophytic communities. There was no clear separation based on the colonizing plant species, indicating that communities hosted by different species were not necessarily dissimilar (Fig. 2B).
Association of prokaryotic communities with the original remediation strategies. The assessment of prokaryotic diversity in the rhizosphere, root endosphere, and shoot endosphere across different treatments is shown in Fig. S1 in the supplemental material. Out of the studied factors, the original soil type (crude oil contaminated or diesel contaminated) was found to be significantly associated with the prokaryotic diversity (P # 0.05, Kruskal-Wallis test) in the rhizosphere, while the original phytoremediation strategy significantly influenced the prokaryotic diversity in the root endosphere (P # 0.05, Kruskal-Wallis test). In the shoot endosphere, both the original soil type and original phytoremediation strategy were significantly associated with prokaryotic diversity (P # 0.01 and P # 0.001, respectively, Kruskal-Wallis test). Fertilization was not found to be significantly associated with prokaryotic diversity in any studied habitat (P . 0.05, Kruskal-Wallis test). Additionally, the interaction between the original type of soil and the phytoremediation strategy was found to be significantly associated with the structure of prokaryotic communities in all three habitats studied: the rhizosphere and root and shoot endospheres (P # 0.05, permutational multivariate analysis of variance [PERMANOVA]) ( Table 1).
The dominant effect of the original soil type on prokaryotic communities in all studied habitats was further demonstrated by redundancy analysis (RDA). Prokaryotic communities in both soil types formed distinct clusters in the ordination space and were separated by the first RDA axis (Fig. 3). Additionally, the rhizospheric soil communities were separated based on the original phytoremediation strategy applied as well (Fig. 3). When included in the PERMANOVA model, the colonizing plant species were not found to be significantly associated with the microbial community structure (tested separately on rhizospheric, root endophytic, and shoot endophytic data sets) (P . 0.05). It should be noted that since the experimental plots were not uniformly vegetated due to distinct patterns of secondary succession (Fig. 1), we were not able to sample the same number of individuals for all plant species. Thus, the statistical power was likely not the same for all plant species sampled. For that reason, the variable of colonizing plant species was excluded from the final PERMANOVA analysis (Table 1) and the association of colonizing plant species with community structure was investigated directly by analyzing the relative abundances of the top 25 prokaryotic genera in the rhizosphere, roots, and shoots of volunteer vegetation (Fig. 4). The most abundant genera were found to be present in similar relative abundances across plant species (Fig. 4). The relative abundances of the top 25 prokaryotic genera in the rhizosphere, root, and shoot samples across different treatments are shown in Fig. S2.
Significantly enriched prokaryotic genera in the rhizosphere and endospheres. To investigate the legacy effects of phytoremediation in more detail, differential abundance analysis was conducted to identify prokaryotic amplicon sequence variants (ASVs) that were significantly enriched (adjusted P value [P adj ] of ,0.01 [Wald test with Benjamini-Hochberg multiple testing correction]). The samples of crude oil-contaminated soil and  diesel-contaminated soil were analyzed separately, as there was a significant interaction between original soil type and phytoremediation strategy (Table 1). For each soil type, three pairwise comparisons were conducted among treatments: (i) planting with one grass species versus the no-remediation control (1P versus 2P), (ii) planting with two grass species versus the control (12P versus 2P), and (iii) 1P versus 12P. The use of different strategies led to the enrichment of distinct ASVs in both the rhizosphere and endospheres of the plants that succeeded in crude oil-contaminated soil (Fig. 5) and diesel-contaminated soil (Fig. 6).

DISCUSSION
In this study, we investigated the prokaryotic communities associated with plants that had colonized a phytoremediation study site previously contaminated with PHCs (diesel or crude oil) in subarctic Alaska. The studied soils were originally subjected, in 1995, to different phytoremediation strategies; specifically, they were sown with one or two grass species or were left unplanted (1P, 12P, and 2P, respectively), with or without fertilization (2F and 1F) (3,12,31,33). The site was subsequently left unmanaged for almost 20 years. In the meantime, both types of soils successfully reached remediation targets, most likely with the help of both the original treatments and colonizing vegetation, which replaced the initial grasses planted on the experimental site and followed different successional trajectories based on the original treatment (12). Here, we aimed to decipher whether and how the original treatments shaped the prokaryotic communities, not only in the rhizosphere but also in the root and shoot endospheres of the colonizing vegetation, 20 years after the original study was established. We found that of the factors studied, the interaction between the original type of contaminated soil and the initial phytoremediation strategy was primarily associated with the structures of the rhizospheric and endophytic prokaryotic communities of the colonizing plants.
Overall, the rhizospheres, root endospheres, and shoot endospheres of the colonizing plants differed from each other in both the structure and diversity of prokaryotic communities, demonstrating an effect of plant physiology on all microbial habitats studied (Fig. 2). Such a finding is in line with previous studies that showed that below-ground and above-ground plant-associated habitats host different microbial communities-for  (35). More importantly, even after nearly 20 years and the disappearance of the originally planted grasses, not only did the original soil type (diesel-or crude oil-contaminated soil) significantly influence the structures of prokaryotic communities in all studied habitats of colonizing plants, but so did the phytoremediation strategy (2P, 1P or 12P) ( Table 1; Fig. 3). In fact, there was an interactive effect of these two variables. The observed influence of the original soil type (crude oil-or diesel-contaminated soil) on community structure is not very surprising, as the two experimental soils differed in texture as well as the original source of PHCs. This finding is in accordance with previous studies (36)(37)(38)(39)(40)(41), which demonstrate that soil characteristics, together with plant species, are generally among the major factors driving plant-associated communities. To our surprise, we did not find the identity of the host plant to be significantly associated with the structure of the prokaryotic community in the rhizosphere, root endosphere, or shoot endosphere, based on PERMANOVA analysis (Table 1). Moreover, there was no clear separation of samples according to plant species either in NMDS ordination (Fig. 2B) or in a heat map investigating the relative abundances of prokaryotic genera across sampled plant species (Fig. 4). It is possible that, rather than employing host-specific mechanisms of microbial selection, which can benefit individual plant species during competition (42), the plants colonizing the disturbed site generally selected microbes that alleviated the contaminant-associated stress induced by the original presence of PHCs. For instance, Oliveira et al. (27) demonstrated such a phenomenon in plants growing in PHC-contaminated soils, which selected endophytes that were able to degrade PHCs. In addition, several studies showed that distantly related plant species growing at different sites in cold regions harbor similar microbial communities: a so-called "cold-adapted plant microbiome" characterized by psychrotolerance and low host-species specialization (43)(44)(45). Taking that into consideration, it could also be that a combination of a high concentration of PHCs and the cold climate created a selective pressure that reduced the pool of microorganisms in the soil that were able to associate with the colonizing plants, and hence, little to no effect of colonizing plant species on community composition occurred. The significant influence of the initial phytoremediation strategy (2P, 1P, or 12P) on the diversity and structure of prokaryotic communities associated with colonizing host plants (Table 1; Fig. 3, 5, and 6; Fig. S1) may be indirectly attributed to plant-driven processes that alter physical and chemical attributes of soil, i.e., rhizodeposition (46). Rhizodeposition has been found to play an important role in shaping soil microorganisms that utilize plant-derived compounds as sources of carbon and/or energy (23,47,48). The rate of microbial degradation of such compounds depends on their chemical structure and availability in soil and can be influenced by the presence of pollutants. Alternatively, plant compounds were shown to affect the degradation of pollutants (49). Moreover, certain pools of soil organic matter can take years or decades to be transformed by microorganisms (50). Therefore, even after the disappearance of the originally planted grasses (1P and 12P), the initial input of organic carbon to PHCcontaminated soils through rhizodeposition and subsequent plant litter deposition likely shaped prokaryotic communities associated with the next generations of plant successors. In addition, microbial communities associated with the planted grasses (Lolium multiflorum and Festuca rubra) could have persisted in PHC-contaminated soils to later inhabit the colonizing plants, which could explain the significant association of the prokaryotic communities with the phytoremediation strategies found herein (Table 1; Fig. 5 and 6). Several recent studies indicated that soil legacy, including the We thus hypothesize that the transfer of potentially beneficial microorganisms among different generations of plants could be important in contaminated soils, including those in extreme environments. Such a process could aid in establishing and maintaining biodiversity at previously disturbed sites, especially in the case of cold-adapted plant microbiota with a wide spectrum of host specificity. Here, for instance, ASVs belonging to the plant growth-promoting diazotrophic genera Bradyrhizobium, Mesorhizobium, and Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium were found to be differentially abundant among phytoremediation treatments. While rhizobia may have contributed to the distinct patterns of plant and microbial succession in these nutrientpoor subarctic soils previously contaminated with crude oil and diesel, testing such hypotheses is beyond the scope of this study and requires further research.
Our analyses indicate that each bioremediation strategy (2P, 1P, and 12P) led to the significant enrichment of specific bacterial ASVs associated with colonizing plants. Accordingly, as there was a significant interaction between the original soil type and the phytoremediation strategy (Table 1), the enriched ASVs were distinct for crudeand diesel-contaminated soils (Fig. 5 and 6). Several of these enriched ASVs belonged to genera that were previously reported to contain members that degrade PHCs; for   (55)(56)(57)(58). It should be noted, though, that the soil texture also varied between the two contaminated soils (12), which may have contributed to different patterns of microbial succession as well. Overall, our results suggest that the originally planted grasses had either served as a source of some of the facultative endophytes that had persisted in PHC-contaminated soils and/or provided an environment favorable to these specific ASVs that later became associated with the colonizing vegetation that succeeded in crude oiland/or diesel-contaminated soils.
Finally, fertilization has previously been shown to stimulate the bioremediation of polluted soils through the enhancement of plant growth and microbial degradation of pollutants (59)(60)(61)(62). Such effects have been observed in PHC-contaminated subarctic soils upon the addition of a fertilizer (12,63). In this study, the interaction of two factors implemented at the beginning of the original study, the application of mineral fertilizer and the original soil type, was still significantly associated with the prokaryotic communities in the rhizosphere and roots of plants growing at the site almost 20 years later. Taken together, the significant association of fertilization, original soil type, and phytoremediation strategy with prokaryotic communities of colonizing plants further points toward the long-lasting impacts that a single treatment can have on plant-associated microbiota in subarctic environments previously disturbed by PHC contamination.
To conclude, we present novel information that extends the understanding of microbial succession in soil and plants in cold regions. By using high-throughput amplicon sequencing, we were able to characterize the structure and diversity of plant-associated prokaryotic communities following initial phytoremediation treatment of PHC-contaminated soils, using different grass species (Lolium multiflorum and Festuca rubra) as phytoremediation agents. Even after almost 20 years and the disappearance of the initially planted grasses, the interaction of the original soil type and phytoremediation strategy, rather than the current host plant, were among the drivers of prokaryotic communities associated with plants that had successfully colonized the site. Furthermore, each phytoremediation strategy led to the significant enrichment of distinct bacterial ASVs in the rhizosphere and endospheres of colonizing plants, which may have initially helped the colonizing plants to survive and adapt to PHC-contaminated soils. Nevertheless, the reader should be once again reminded that the original treatments influenced patterns of secondary succession at the research site (12). Thereby, it is possible that, even though we did not find a strong association between the colonizing plant species and plant-associated microbial community composition, the legacy of the original treatments detected in the rhizosphere, roots, and shoots of volunteer vegetation might be, at least partially, a secondary effect of distinct secondary succession strategies that occurred at the site. Overall, the outcome of this study provides evidence that in previously disturbed subarctic ecosystems, the non-host-specific soil-to-root transfer of microorganisms seems to be of great importance for the succession of plant-associated microbiota. Future research might therefore aim to decipher how different environmental factors promote either horizontal or vertical transmission of microorganisms among different generations and species of plants in subarctic regions. Finally, the results of our study highlight the role of phytoremediation not only in the removal of PHCs but also as a legacy that determines the long-term microbial succession of soils previously contaminated by PHCs. Understanding the microbial assemblages associated with vegetation colonizing restored soils is valuable for the management of contaminated subarctic soils with the goals of sustainably maintaining ecosystem function and resilience. Microbiology Spectrum the experiment included crude oil-contaminated soil collected from a gravel pad at a pump station on the Trans-Alaska pipeline and diesel-contaminated soil collected during the removal of an underground storage tank. In addition to the difference in the original source of PHCs (crude oil or diesel), the soils differed in structure: the crude oil-contaminated soil was gravel with a large grain size, while the dieselcontaminated soil was finer in texture with more organic matter (12). For that reason, the crude oil-and diesel-contaminated soils are referred to herein as the original soil types. The total initial PHC concentrations were approximately 3,420 and 800 mg/kg in crude oil-contaminated soil and diesel-contaminated soil, respectively. Each soil type was separately placed in a lined and bermed area approximately 21 by 3 m and 60 cm deep. Each area was later subdivided into seven treatment plots. The treatments included three levels of vegetation and two levels of nutrient addition for each soil type, as follows: the plots were sown with annual ryegrass (Lolium multiflorum) (1P) or a 1:1 mixture of annual ryegrass and Arctared fescue (Festuca rubra) (12P) or left unplanted (2P) and were additionally treated with commercially available mineral fertilizer, which was surface applied at approximately 620 g/m 2 of N, P, and K (granular 20-20-10) (1F), or left unfertilized (2F) (3,12,31,33). Soil and plant sampling for this study was carried out in September of 2014, by which point the PHC concentrations were below detection limits (,0.5 ppm; EPA method 8015M) in all treatment plots (64). The plots had been successfully colonized by a variety of volunteer herbaceous and woody plants depending on the original phytoremediation strategy (12). The plants (referred to herein as colonizing plants) that were sampled belonged to the following species: Picea glauca (white spruce; WS), Epilobium angustifolium (fireweed; FW), Achillea millefolium (yarrow; Y), Salix spp. (willow; W), Poa spp. (blue grass; L), Populus balsamifera (poplar; P), Shepherdia canadensis (buffalo berry; BB), Betula neoalaskana (birch; B), and Trifolium hybridum (clover; C). Since the site gradually turned from experimental to observational and the experimental plots were not colonized uniformly (12), not all species were harvested from all plots. Sampling was carried out according to the scheme in Fig. 1, with five colonizing plants of five species (i.e., 5 replicas) per plot being harvested.

MATERIALS AND METHODS
Processing of plant and rhizosphere soil samples. After the removal of bulk soil, the rhizosphere soil samples were collected by shaking off the soil directly adhering to the roots of harvested plants. The roots and shoots were subsequently surface sterilized by immersion in 2% sodium hypochlorite for 15 min. Immediately after surface sterilization, plant samples were rinsed with sterile deionized water three times and finally washed with 10 mM MgSO 4 to remove sodium hypochlorite. Twenty-microliter amounts of the final wash solutions were spread on plate count agar (PCA) plates and incubated at 25°C for 2 days to confirm sterility. The plant and soil samples were stored at 280°C prior to further analyses.
DNA extraction and 16S rRNA gene sequencing. To isolate metagenomic DNA, 0.5 g of the rhizosphere soil and 0.1 to 0.7 g (depending on the quantity harvested) of surface-sterilized roots and shoots were used. Prior to DNA isolation, plant samples were subjected to repeated freezing at 280°C and thawing at 25°C. All plant and rhizosphere soil samples were then homogenized by using the FastDNA spin kit for soil (MP Biomedicals, OH, USA) and the FastPrep instrument for 40 s at 30 m/s followed by 3 min at 15 m/s or for 80 s at 12 m/s, respectively. The DNA was extracted using the same kit following the manufacturer's instructions and then purified and concentrated using the Genomic DNA Clean and Concentrator-10 kit (Zymo Research Irving, CA, USA).
The metataxonomic analysis of plant-associated prokaryotic communities using 16S rRNA gene amplicon sequencing was done as described previously (14). Briefly, for the soil DNA samples, the 515 forward primer (59-GTGYCAGCMGCNGCGG-39; Sigma-Aldrich, USA) and 926 reverse primer (59-CCGYCAATTYMTTTRAGTTT-39; Sigma-Aldrich, USA) were used to target the V4-V5 hypervariable region of the 16S rRNA gene (49). The first 15-mL PCR mixture contained 0.02 U/mL Kapa HiFi hot start ready mix (Kapa Biosystems, USA), 0.3 mM each primer (Sigma-Aldrich, USA), 1 mL of template DNA (;15 ng/mL), and PCR-grade water (Sigma-Aldrich, USA). The temperature cycling conditions were as follows: initial denaturation at 95°C for 5 min, 25 to 35 cycles of 20 s at 98°C, 15 s at 56°C, 15 s at 72°C, and final extension at 72°C for 5 min (65). A volume of 0.5 mL of the PCR product was used as a template in a second PCR with the same primers modified with internal barcodes and sequencing adapters of various lengths (49). This round of PCR was performed under the same conditions, but the final reaction mixture volume was 25 mL, the concentration of each primer was 1 mM, and the number of cycles and annealing temperature were decreased to 8 to 14 and 50°C, respectively.
For DNA samples extracted from plant biomass, peptide nucleic acids were used to prevent the amplification of plant organellar DNA (66). The first 15-mL PCR mixture contained 0.3 mM each of peptide nucleic acid probes targeting mitochondrial genes (mPNAs; 59-GGCAAGTGTTCTTCGGA-39) and plastid genes (pPNAs; 59-GGCTCAACCCTGGACAG-39) (PNA Bio, Thousand Oaks, CA), 0.02 U mL 21 of Kapa HiFi hot start ready mix (Kapa Biosystems, USA), 0.3 mM 515 forward primer, 0.3 mM 1068 reverse primer (59-CTGRCGRCRRCCATGCA-39; Sigma-Aldrich, USA), 1 mM template DNA (;15 ng), and PCR-grade water (Sigma-Aldrich, USA) (14). The thermal cycling conditions started with a 5-min denaturation at 95°C, followed by 35 cycles at 98°C, 15 s at 72°C (annealing of PNAs), 15 s at 56°C, 15 s at 72°C, and a final extension for 5 min at 72°C. Each sample was prepared in 8 copies that were pooled and separated by electrophoresis on a 1.5% agarose gel. Bands 550 bp in size were excised from the gel and purified using a Zymoclean gel DNA recovery kit (Zymo Research, USA). Then, 0.5 mL of the purified PCR product was used as a template in the same 2-step PCR process as described above for soil DNA samples, except that each round of PCR included the PNAs and an additional step of PNA annealing (15 s at 72°C).
The resulting PCR products were purified using SPRI magnetic beads (Beckman Coulter, USA), and further library preparation and sequencing analysis on an Illumina MiSeq instrument were performed at the Institute of Arctic Biology Genomics Core Laboratory at the University of Alaska Fairbanks, USA. The concentrations of amplicons were normalized to 1 to 2 ng/mL using a SequalPrep kit (Thermo Fisher Scientific, USA), and the samples corresponding to each plate were pooled and either not diluted or diluted 1.5-fold and 3-fold. These three technical replicates of the amplicon libraries were then subjected to 8-cycle PCR to add specific Illumina adapters and sequenced.
Data analyses. The sequence data were analyzed using the DADA2 pipeline (67) and R software (version 4.1.0). Several modifications were made to the DADA2 standard operating procedure (SOP) based on the analysis of a mock community, which consisted of 12 bacterial strains (68) that were amplified in parallel with our DNA samples. Briefly, the primer sequences were trimmed off when found present, otherwise the whole read was discarded. To manage the lower sequence quality toward read ends, forward and reverse reads were truncated to lengths of 247 and 170 bp, respectively, and filtered according to their quality using the following parameters: maxN=0, maxEE=1, truncQ=2. After dereplication and the application of DADA2-based removal of sequencing errors, denoised forward and reverse reads were merged and chimeric sequences were removed using the method="pooled". Based on the mock community analysis, the sequences differing by one base were clustered and the most abundant sequence was considered the valid one. The technical replicates were merged, and only the sequences that occurred in all technical replicates were kept, resulting in a total number of 6,008,156 reads. Finally, taxonomy was assigned using silva_nr_v132_train_set.fa.gz (69) to create a database of amplicon sequence variants (ASVs). ASVs of plant origin, including mitochondria and chloroplasts, were discarded, which accounted for 12.5% of all sequences. The remaining data set was rarefied to the smallest sample size, 2,300 reads per sample, to ensure the comparability of rhizospheric and endophytic data sets.
Further multivariate statistical analyses of sequence data were conducted using the packages phyloseq (70), vegan (71), DESeq2 (72), ggplot2 (73), and ampvis2 (74). A maximum-likelihood phylogenetic tree (GTR1G1I) was constructed using the packages DECIPHER and phangorn as described by Callahan et al. (75). Prokaryotic diversity was assessed by calculating the Shannon and Simpson alpha diversity indices (75), and the nonparametric Kruskal-Wallis test was used to test whether the prokaryotic diversity differed significantly among plant-associated habitats. To test the statistical significance of the effects of the original soil type (diesel-or crude oil-contaminated soil), fertilization, and phytoremediation strategy on prokaryotic diversity, the Kruskal-Wallis test was performed using Shannon diversity index values calculated separately for the rhizospheric, root endophytic, and shoot endophytic data sets. To investigate the prokaryotic community structures in the rhizosphere, roots, and shoots of colonizing plants, nonmetric multidimensional scaling (NMDS) with weighted UniFrac distances was used. Heat maps were constructed using the package ampvis2 (74). The statistical significance of the effects of the original soil type, fertilization, phytoremediation strategy, and colonizing host plant species on the prokaryotic communities of colonizing plants was determined separately on the rhizospheric, root endophytic, and shoot endophytic data sets of ASVs by permutational multivariate analysis of variance (PERMANOVA) (76) implemented in the adonis2 function of the vegan package (71) and based on Bray-Curtis dissimilarity. For each data set, R 2 values were calculated for all individual variables using the adonis function, and subsequently, the order of the variables in the final PERMANOVA model was sorted in descending order based on their respective R 2 values. Redundancy analysis (RDA) based on centered log ratio-transformed data (77) was then conducted to investigate the association of plant-associated prokaryotic communities with the factors studied (the original soil type, phytoremediation strategy, fertilization, and colonizing host plant species) at the ASV level. To identify prokaryotic ASVs that were significantly enriched among treatments, differential abundance analysis was conducted on the unrarefied data set of ASVs using the Wald test with Benjamini-Hochberg multiple testing correction (P adj , 0.01) as implemented in the DESeq2 package (72). The function lfcShrink was used to shrink the log 2 fold changes.
Data availability. The unprocessed FASTQ files for all samples were deposited in SRA under BioProject accession number PRJNA771088.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 0.3 MB.