Soil initial bacterial diversity and nutrient availability determine the rate of xenobiotic biodegradation.

Understanding the relative importance of soil microbial diversity, plants and nutrient management is crucial to implement an effective bioremediation approach to xenobiotics-contaminated soils. To date, knowledge on the interactive effects of soil microbiome, plant and nutrient supply on influencing biodegradation potential of soils remains limited. In this study, we evaluated the individual and interactive effects of soil initial bacterial diversity, nutrient amendments (organic and inorganic) and plant presence on the biodegradation rate of pyrene, a polycyclic aromatic hydrocarbon. Initial bacterial diversity had a strong positive impact on soil biodegradation potential, with soil harbouring higher bacterial diversity showing ~ 2 times higher degradation rates than soils with lower bacterial diversity. Both organic and inorganic nutrient amendments consistently improved the degradation rate in lower diversity soils and had negative (inorganic) to neutral (organic) effect in higher diversity soils. Interestingly, plant presence/type did not show any significant effect on the degradation rate in most of the treatments. Structural equation modelling demonstrated that initial bacterial diversity had a prominent role in driving pyrene biodegradation rates. We provide novel evidence that suggests that soil initial microbial diversity, and nutrient amendments should be explicitly considered in the design and employment of bioremediation management strategies for restoring natural habitats disturbed by organic pollutants.

Introduction

Polycyclic aromatic hydrocarbons (PAHs) are among the most harmful xenobiotics, belonging to a class of ubiquitous hazardous organic pollutants. Excessive bioaccumulation of PAHs due to anthropogenic activities affects the health of aquatic, terrestrial and human life (Smets and Pritchard, 2003; Megharaj et al., 2011; Bisht et al., 2015; Ghosal et al., 2016). Amongst the high‐molecular‐weight PAHs, pyrene is one of the most recalcitrant and predominant in soils from industrial sites (Haeseler et al., 1999) and is quite prevalent and persistent in the environment (Nzila et al., 2018). Because of its low biodegradability and carcinogenic effects, pyrene is considered a priority pollutant by US EPA (US EPA, 1990; IARC, 2010), attracting significant research efforts to find efficient and environmentally friendly approaches to speed up removal from contaminated sites (Megharaj et al., 2011).

Microbially mediated biotransformation and biodegradation are considered as the main natural forms of PAHs attenuation in polluted sites (i.e. bio‐attenuation) (Cébron et al., 2011; Xu and Zhou, 2017; Ławniczak et al., 2020) and, in the plethora of highly diverse and specialized microbial communities capable of PAH degradation (Fernández‐Luqueño et al., 2011; Patel et al., 2020), bacteria are considered as the most efficient and favourable (Premnath et al., 2021). However, bio‐attenuation relies on the natural ability of microorganisms such as bacteria to carry out the degradation of PAHs, and several engineered bioremediation strategies are commonly employed to improve the microbially mediated removal of these recalcitrant contaminants in soils, including biostimulation and phytoremediation (Tyagi and da Fonseca, 2011; Chen and Zhong, 2019). Biostimulation provides for the adjustment of abiotic environmental conditions (e.g. nutrient availability, soil pH) to stimulate the degradation capacity of the indigenous soil microbial community and is commonly implemented through the application of amendments such as organic or inorganic fertilizers (Nikolopoulou and Kalogerakis, 2008). By contrast, phytoremediation involves the removal of xenobiotic organic contaminants through mutual interaction of plant roots and the indigenous microbial population, whereby plant addition favours microbial growth and activity by supplying readily available nutrient sources for microorganisms in the rhizosphere (rhizoremediation), and thus can indirectly enhance biodegradation by shaping an environment conducive of increased microbial biological activity (Chaillan et al., 2006; Cheema et al., 2009; Kalantary et al., 2014; Leff et al., 2015; Xue et al., 2018; Gabriele et al., 2021).

Despite their potential, the efficacy of these bioremediation enhancement techniques remains inconsistent in field conditions (Megharaj et al., 2011; Rayu and Karpouzas, 2012). For instance, commonly employed inorganic and organic amendments differ in the content and nature of the nutrients they provide (e.g. high‐nutrient, water‐soluble compounds versus low‐nutrient, C‐rich compounds respectively). Yet, fertilization efficacy can depend on the initial carbon and nutrient requirements of the resident microbial community (e.g. carbon, nitrogen or phosphorus limited), making it difficult to generalize the most appropriate management approach to improve natural attenuation (Nervo et al., 2017; Ite and Ibok, 2019). Similarly, recent studies have shown that soil biotic and abiotic characteristics can influence the key soil processes linked to phytoremediation of hydrocarbons (Yergeau et al., 2015; Bell et al., 2016; Correa‐García et al., 2020), pointing at a pivotal role of soil nutrient status and initial levels of microbial diversity in determining phytoremediation outcomes.

Recently, the combination of these two approaches (i.e. biostimulation and phytoremediation) has been proposed as potentially effective in improving bioremediation due to their synergistic effect in promoting natural attenuation (McIntosh et al., 2017; García‐Sánchez et al., 2018). Indeed, fertilization during phytoremediation can reduce competition between plants and microorganisms for nutrients in soils, while also enhancing microbial activity (Hoang, Sarkar et al., 2021). However, very few studies have evaluated the efficacy of combined biostimulation and phytoremediation for organic pollutant degradation (Hutchinson and Banks, 2001; McIntosh et al., 2017). Additionally, whether and to what extent inherent soil biological characteristics, such as the initial level of microbial diversity and composition, can affect the outcome of these combined bioremediation approaches is unknown. Yet, differences in microbial diversity and structure can alter important ecosystem functions (Bardgett and van der Putten, 2014; Delgado‐Baquerizo et al., 2016), including specialized processes such as the biodegradation of xenobiotics (Díaz, 2008; Singh et al., 2014; Kong et al., 2018). As such, understanding the relative importance and mechanisms that underpin interactions between plants, soil amendments and microbial characteristics in the biodegradation context is essential to design effective strategies and thus achieve successful outcomes under field condition.

In this study, we used pyrene as a model xenobiotic compound for high‐molecular‐weight PAHs remediation (Kanaly and Harayama, 2010), and evaluated the individual and interactive impacts of differing soil microbial diversities (lower and higher bacterial diversity) and composition, biostimulation via nutrient amendment (inorganic and organic fertilization) and rhizoremediation via the presence of C3 and C4 grasses, on the pyrene biodegradation rate. We hypothesized that (i) nutrient amendment and plant presence would positively affect degradation rate because of direct (e.g. provision and exchange of nutrients) and indirect (e.g. improvement of soil characteristics) enrichment effects on the microbial community (Chaîneau et al., 2005; Kalantary et al., 2014; Leff et al., 2015; Xue et al., 2018); and (ii) the nutrient and plant effects on degradation rates would be influenced by the initial levels of soil microbial diversity and composition. To test these hypotheses, we conducted a full factorial microcosm experiment involving two soils, similar in pedogenic attributes but harbouring contrasting bacterial diversity (lower and higher, referred as initial microbial diversity hereafter), treated with pyrene and with or without nutrient amendment and plants, for a total of 24 unique treatments plus controls (Fig. 1). We monitored pyrene biodegradation rates, measured changes in soil nutrient availability, and assessed diversity and community structure of soil bacteria by 16S rRNA gene amplicon sequencing across three sampling time points. To further characterize the soil biodegradation potential, we investigated the abundance of the pahAc subunit of the bacterial PAH ring‐hydroxylating dioxygenases, a widespread enzyme involved in bacterial‐mediated biodegradation of PAHs.

Fig. 1

Experimental overview presenting fully factorial design of the growth‐chamber microcosm experiment, consisting of two soils with contrasting microbial diversity (high and low), treated with or without pyrene (200 ppm), with three nutrient amendment treatments (untreated, inorganic and organic amendments), and four plant combinations (unplanted, C3, C4 monocultures and C3+C4 combination). Triplicate pots were established for each treatment combination in a randomized design.

Results

Overall impact of treatments on pyrene biodegradation rates, microbial diversity and composition and soil physico‐chemical properties

Repeated measures ANOVA (Fig. 2, Table S1) across the four time points revealed that pyrene degradation rate was significantly impacted by the initial level of microbial diversity (F 1,48 = 70.9, P < 0.05) and its interactive effect with nutrient levels (F 2,48 = 45.5, P < 0.05) and plant presence (F 3,48 = 43.5, P < 0.05). Conversely, nutrient levels (F 2,48 = 0.90, P > 0.05) and plant presence (F 3,48 = 2.01, P > 0.05) alone or in combination (F 6,48 = 0.85, P > 0.05) had no significant effect on degradation rates. Soil physico‐chemical parameters were significantly impacted by the initial levels of microbial diversity (NH4: F 1,48 = 73.2; NO3: F 1,48 = 8.65; PO4: F 1,48 = 30.66 and pH: F 1,48 = 790, P < 0.05 for all), the application of nutrient amendments (NH4: F 2,48 = 42.6; NO3: F 2,48 = 35.2; PO4: F 2,48 = 84.9 and pH: F 2,48 = 38.2, P < 0.05 for all), and their interaction (NH4: F 2,48 = 44.01; NO3: F 2,48 = 1.1; PO4: F 2,48 = 14.7 and pH: F 2,48 = 56.4, P < 0.05 for all), while plant presence had significant impact on the levels of phosphorus (F 3,48 = 3.6, P < 0.05) and pH (F3,48 = 11.42, P < 0.05 ) exclusively. Similarly, soil microbial diversity was significantly impacted by nutrient availability (F 2,48 = 31, P < 0.05), and their interaction (F 2,48 = 13.6, P < 0.05), but it was not affected by the presence of plants (F 3,48 = 0.908, P > 0.05). Consistently, our global PERMANOVA indicated a significant effect of nutrient availability (F 2,48 = 31, P < 0.05) on the microbial community structure.

Fig. 2

Effect size (Eta squared, η2) of initial microbial diversity (IMD), nutrient amendments application, plant presence and time (weeks) on pyrene degradation rate, microbial diversity, soil parameters and their interactions. Significant (P < 0.05) effects are coloured in blue. Error bars represent ± 1 SE.

Most treatments and their interactions showed significant differences over time (Fig. 2, Tables S1 and S2). Consequently, to assess the change in treatment effect on pyrene biodegradation, soil properties and microbial diversity during the course of the experiment, we also performed one‐ and three‐way analyses of variance for the intermediate and final time points (week 4, week 10 and week 15). Similarly, to assess the impact of treatments on microbial diversity, we performed a PERMANOVA analysis at week 4 and week 15.

Impact of treatments on microbial diversity and composition

Initial microbial analysis of soils from irrigated and unmanaged plots prior to treatment manipulation confirmed that the irrigated plots (HDS) had higher microbial diversity compared to the soils from unmanaged plots (LDS) (Shannon, one‐way ANOVA, F 1,4 = 864.2, Pr(>F) <0.001) (Fig. S1A). Species richness (Chao index) was also higher in irrigated soil (3887 ± 76) compared to unmanaged (2809 ± 31) soil (one‐way ANOVA, F 1,4 = 168.8, Pr(>F) <0.001; Fig. S1B), confirming the substantial difference in microbial diversity of soils we used at the commencement of this study.

We also monitored the changes in microbial diversity and composition during the course of this experiment (weeks 4 and 15) to evaluate whether initial differences in soil diversity and structure were maintained throughout the experiment and how treatments affected these properties. Significantly, different levels of microbial diversity between the two soils were still present after 4 weeks (Shannon index F 1,60 = 59.17, P < 0.0001 and observed OTUs F 1,60 = 203.0, P < 0.0001), where nutrient amendment treatments showed significant effects (Shannon index F 2,60 = 20.49, P < 0.0001 and observed OTUs F 2,60 = 41.48, P < 0.0001), but the presence of plants did not have any effect (Shannon index F 1,60 = 0.08, P = 0.7737 and observed OTUs F 1,60 = 0.16, P = 0.6863) (Fig. S2). Similarly, at week 15, HDS had baseline higher levels of microbial diversity (Shannon index: F1,60 = 13.43, P = 0.0005; Fig. S3), and richness (Chao index: F 1,60 = 30.38, P < 0.0001; Fig. 2B, and observed OTUs: F1,60 = 25.28, P < 0.0001; Fig. S3). Also, nutrient amendment treatments showed no significant effect on microbial diversity (Shannon index: F 2,60 = 2.66, P = 0.0780; Fig. 2A) but had significant effect on richness (Chao index: F 2,60 = 4.41, P = 0.0162; Fig. 2B, and observed OTUs: F 2,60 = 4.14, P = 0.0206; Fig. S3).

When considering the treatment effects on microbial community structure across the individual time points, the two soils consistently grouped separately in the ordination space as a function of nutrient treatments (Fig. S4A,B). The PERMANOVA conducted for the individual dates corroborated these patterns, whereby community composition was impacted by nutrient amendments, initial microbial diversity and the interaction between nutrient amendments and initial microbial diversity across all time points (P < 0.05) (Fig. S4C,D).

Role of initial levels of soil microbial diversity on pyrene biodegradation

Our results of pyrene degradation suggest that initial levels of microbial diversity and microbial composition had a strong impact on the biodegradation potential of soils. Although pyrene degradation rate was found similar in the early stages of the experiment (week 4), by week 10, degradation was ~ 2 times higher (ANOVA, F 1,60 = 58.84, P < 0.0001) in HDS (95% ± 1.4) compared to LDS (51% ± 35.7; Fig. 3). Even after 15 weeks, pyrene degradation remained significantly lower in the LDS (83% ± 4.3) compared to the HDS (ANOVA, F 1,60 = 92.89, P < 0.0001), where we observed complete (100%) degradation (Fig. 3).

Fig. 3

Role of soil initial microbial diversity, nutrient amendments and plant presence on rate of pyrene biodegradation. Bar graphs comparing the pyrene degradation rates between soils with contrasting initial microbial diversity (IMD) and across the treatments at different time points. Three‐way ANOVA was performed followed by Sidak’s multiple comparisons test compares significant differences between soil initial microbial diversity mentioned as ****P < 0.0001, ***P < 0.001, **P < 0.01, *P < 0.05, and Tukey's multiple comparisons test compares across the nutrient amendment and plant presence treatments in high and low diversity soil separately, indicated as letters significant at P < 0.05.

Consistent with our observation of increased degradation rates in HDS compared to the LDS soil at week 10 (Fig. 3), we detected a significant difference in abundance of pahAc gene (a functional marker gene of PAH‐degrading bacteria) between the two soils (three‐way ANOVA, F 1,60 = 42.88, P < 0.0001; Fig. 4A), with high diversity soil harbouring an overall higher degradation gene copy number in controls and across all treatments (Fig. 4A). This trend was mirrored by a positive correlation between degradation gene abundance and pyrene degradation rate (Spearman rho = 0.27, P = 0.025) (Fig. 4B). This indicates that the number of genes associated with pyrene degradation increases with increasing degradation rate, and that this change is correlated with the initial levels of soil microbial diversity.

Fig. 4

(A) Difference in abundance of pahAc genes in high and low diversity soils at week 10, assessed via three‐way ANOVA followed by Sidak’s multiple comparisons test significant at **P < 0.01; *P < 0.05 and Tukey's multiple comparisons test across nutrient amendment and plant treatments.

B. Correlation (Spearman) between biodegradation rate and pahAc gene abundance. IMD, Initial Microbial Diversity.

Impact of soil nutrient availability on soil physico‐chemical properties and pyrene biodegradation rates

Extractable ammonium (), nitrate () and phosphorus (PO4) were measured to monitor the nutrient availability in soils (Fig. S5A–C respectively), to reflect the changes induced from nutrient amendment treatments (untreated, inorganic fertilizer and organic fertilizer) and to evaluate the consequences for biodegradation potential of soils. Both high and low diversity soils showed relatively low nitrogen levels, with inorganic fertilizer treatment showing high initial (week 0) ammonium levels (P < 0.0001, Fig. S5A). High diversity soils showed immediate conversion of ammonium into nitrate (Fig. S5B), in contrast to low diversity soils where this conversion was significantly slower (P < 0.001, at week 4). Untreated and organic fertilizer treatments showed similar patterns over time with relatively low ammonium (ANOVA, F 3,88 = 32.9, P < 0.001 and F 3,88 = 2.4, P = 0.076) and nitrate (ANOVA, F3,88 = 4.3, P = 0.006 and F 3,88 = 4.1, P = 0.009) levels respectively. Phosphate levels were also significantly different among nutrient amendment treatments (ANOVA, F 2,264 = 832, P < 0.001) where availability of phosphate increased after inorganic fertilizer addition compared to untreated soils (P < 0.001, Fig. S5C), while no significant changes were observed in organic fertilizer treatment (P = 0.938, Fig. S5C). Soil pH was also affected by soil nutrient amendments (ANOVA, F 2,264 = 6.46, P = 0.002; Fig. S5D), where inorganic fertilizer consistently increased pH both in high diversity (P < 0.001) and low diversity soil (P < 0.001). Organic fertilizer had no effect on soil pH in high diversity soil (P = 0.806) but had a negative effect on pH in low diversity soil (P = 0.003). Untreated high diversity soil showed no significant change in pH over 15 weeks, whereas in low diversity soil, pH fell by more than one pH unit (from 7.5 to 6.1, P < 0.001). A similar trend was observed in both inorganic and organic fertilizer treatments, where soil pH in high diversity soils was maintained at pH ~7 (inorganic ‐ pH 8.9 to 7.0; P < 0.001, and Organic ‐ pH 7.7 to 7.0; P = 0.009), while in low diversity soils pH turned acidic by dipping below 5.8 (Inorganic ‐ pH 8.3 to 5.8; P < 0.001, and Organic ‐ pH 6.8 to 5.6; P < 0.001). Overall, the trend is evident (Fig. S5D) that soil pH declined over time in all treatments, and it was stronger in inorganic amendment than organic because of higher starting pH values.

Soil nutrient amendment treatments (organic and inorganic) showed mixed effects on the pyrene biodegradation rates. The addition of inorganic fertilizer significantly increased the degradation rates in the low diversity soil, with 97% ± 1.5 of pyrene degraded on average after 10 weeks, versus 51% ± 35.7 (p < 0.001) in the untreated soil. At the same time, it showed the opposite effect on the high diversity soil, with 71% ± 21.8 of pyrene dissipated on average after 10 weeks, versus 95% ± 1.4 (P = 0.19) in the untreated soil (Fig. 3). After 15 weeks, complete degradation (100%) was observed in low diversity soil treated with inorganic fertilizer compared to untreated at 83% ± 4.3 (P < 0.001). Conversely, there was no difference in the degradation rate between untreated and inorganic amendment in high diversity soil. Degradation rates were stimulated by organic fertilizers in the low diversity soil, where the degradation rate almost doubled from 16.9% in non‐amended soils to 33.6% in organic amendment soils by week 4 (P < 0.05). However, by week 15, organic amendment had no effect, and the degradation rate remained below 76%. Similarly, organic fertilizer had a negligible effect on the degradation rates in high diversity soils throughout the experiment (week 4; P = 0.10, week 10; P = 0.99 and week 15; P = 0.81) (Fig. 3). Generally, inorganic amendments improved the biodegradation of pyrene in the low diversity soil, while having negative (inorganic) effect in the more diverse soils; conversely, organic amendment’s positive effect was limited to early stage in low diversity soil and had an overall neutral effect on pyrene degradation rates in high diversity soil (Fig. 3 and Table 1).

Table 1

Three‐way ANOVA comparing the role of interactions among soil initial microbial diversity, nutrient amendment, and plant presence on rate of pyrene biodegradation. Numbers in bold represent statistically significant interactions.

ANOVA table	Week 4		Week 10		Week‐15
Source of variation	F (DFn, DFd)	P value	F (DFn, DFd)	P value	F (DFn, DFd)	P value
α Initial microbial diversity (IMD)	F (1,60) = 0.80	0.3746	F (1,60) = 58.84	< 0.0001	F (1,60) = 92.89	< 0.0001
β Nutrient amendments	F (2,60) = 8.96	0.0004	F (2,60) = 3.58	0.0339	F (2,60) = 14.31	< 0.0001
γ Plant presence	F (1,60) = 0.31	0.5797	F (1,60) = 22.88	< 0.0001	F (1,60) = 0.72	0.3970
α IMD × β Nutrient	F (2,60) = 3.21	0.0471	F (2,60) = 68.07	< 0.0001	F (2,60) = 27.24	< 0.0001
β Nutrient × Plant presence	F (2,60) = 6.37	0.0031	F (2,60) = 4.56	0.0143	F (2,60) = 10.47	0.0001
α IMD × γ Plant presence	F (1,60) = 0.16	0.6835	F (1,60) = 15.70	0.0002	F (1,60) = 0.05	0.8289
α IMD × β Nutrient × γ Plant presence	F (2,60) = 1.77	0.1790	F (2,60) = 2.49	0.0910	F (2,60) = 1.34	0.2697

Initial microbial diversity: soils with higher and lower bacterial diversity.

Nutrient amendment treatments: Untreated, inorganic fertilizer and organic fertilizer.

Plant presence treatments: Plant (average of C3, C4 and C3C4 combination), and no plant.

Influence of the plant presence on soil physico‐chemical properties and degradation rates

Plant presence had a modest impact on soil physico‐chemical properties (Fig. S6), whereby significant (P < 0.05) effects were limited to pH in untreated soils and soils treated with organic fertilizers at week 4, phosphate at week 4 in untreated and inorganic fertilizer treatment and nitrate at week 15 in untreated soils. Interestingly, and at odds with our predictions, plant treatments (both monocultures and C3‐C4 combinations) did not substantially improve pyrene degradation rate in any of the treatment combinations compared to no plant. There were a few exceptions, where C3 plants significantly improved biodegradation only in untreated and organic fertilizer treatments (P < 0.05) in low diversity soil (Fig. S7). Therefore, the C3, C4 and C3‐C4 combinations were combined as ‘Plant’ treatment for further comparisons with no plant and nutrient amendment treatments (Fig. 3). Plant presence in combination with inorganic fertilizer treatment showed the least degradation rates in high diversity soils at week 10 (58%, P < 0.0001) and week 15 (92%, P = 0.232), compared to ‘untreated+plant’ treatment (90% and 97% respectively). By contrast, ‘organic fertilizer + plant’ treatment slightly improved the degradation rates to 95% (P = 0.962) and 98% (P = 0.892) respectively at weeks 10 and 15, compared to ‘organic fertilizer’ alone, which improved degradation from 85% and 94%. Otherwise, there was no significant effect of plant growth on degradation rates in high diversity soil. In case of low diversity soil, ‘plant presence’ treatment negatively affected the pyrene degradation consistently across all the nutrient amendment treatments at week 10 by 52% to 31% (P = 0.15) in untreated, 97% to 79%, (P = 0.32) in inorganic and, 52% to 38% (P = 0.662) in organic amendment treatments with no significant differences at week 15. Largely, plant presence showed no effect at week 4 (F 1,60 = 0.31, P = 0.579) and week 15 (F 1,60 = 0.72, P = 0.397) (Fig. S7), while showing moderate effects at week 10 (F 1,60 = 57.38, P = 0.0029) (Fig. 3).

Interactive effects of microbial diversity, plant and nutrient availability on the degradation potential of soils

The above results highlighted the importance of management practices in enhancing pyrene biodegradation rates, which is directly linked to soil microbial diversity and nutrient availability in the low nutrient soils we studied. Consistently, we found a strong interaction between soil microbial diversity and nutrient availability (F 2,60 = 68.07, P < 0.0001; week 10) (Table 1). While initial microbial diversity seems to be the most prevailing factor to determine biodegradation potential of a soil (F 2,60 = 92.89, P < 0.0001; week 10), nutrient availability is one of the main rate‐limiting factor, which can significantly influence the degradation rate (week 10; F 2,60 = 3.823, P = 0.0274) (Table 1). The other factor that may contribute to enhance degradation rate is plant presence, by direct or indirect interaction with microbial diversity and nutrient availability. Although influence of plant presence was negligible at week 4 (F 1,60 = 0.31, P = 0.579) during their initial establishment, we saw some impact during the plant vegetative growth phase at week 10 (F 1,60 = 22.88, P < 0.0001) (Table 1). Additionally, there was a significant interaction between nutrient amendment and plant during this phase (F1,60 = 4.56, P = 0.0146), although the interaction of nutrient amendment with soil microbial diversity was stronger (F 2,60 = 68.07, P < 0.0001) (Table 1). In the final stage of biodegradation process (by week 15), initial microbial diversity was the only factor (F 1,60 = 92.89, P < 0.0001) that facilitated the complete degradation of pyrene in untreated high diversity soils (Fig. 3) (Table 1). Conversely, a similar degradation rate was achieved in low diversity soil with the aid of inorganic amendment by improving the nutrient availability of soil (F 1,60 = 27.24, P < 0.0001) (Table 1). These interactions seemingly articulated our initial hypothesis that soils harbouring higher microbial diversity will show higher PAHs degradation potential, while nutrient amendment and plant’s growth will also positively affect degradation rate via direct and indirect changes in microbial activities.

Effect of nutrient amendments on microbial community composition and indicator taxa

Seven major phyla, namely, Firmicutes, Proteobacteria, Actinobacteria, Acidobacteria, Planctomycetes, Gemmatimonadetes and Cyanobacteria, accounted for more than 98% of the total microbial abundance in the pyrene‐treated samples, as well as between the two reference soils (i.e. the original pyrene‐untreated soil). Particularly, the high diversity reference soil were characterized by significantly (P < 0.05) higher initial abundance of Proteobacteria, Acidobacteria, Cyanobacteria, Bacteroidetes, Planctomycetes and Gemmatimonadetes, while the low diversity soil had higher initial proportion of Actinobacteria and Firmicutes (Fig. S8). A total of 189 and 64 OTUs were identified as significant indicators of the pyrene‐treated high‐ and low diversity soils, respectively, and were assigned to the phyla Proteobacteria, Planctomycetes, Actinobacteria, Acidobacteria, Chloroflexi, Bacteroidetes and Firmicutes (Fig. 5). The composition of the microbial indicator communities changed mainly in relationship to the inorganic fertilizer application in the low diversity soil at week 4. In this soil, we observed a number of Firmicutes (genera Bacillus, Clostridiales, Cohnella, Sporosarcina) and Proteobacteria (genera Cytophaga, Pedobacter, Cohnella, Sporosarcina, Paludisphaera, Delftia, Herbaspirillum, Noviherbaspirillum, Variovorax, Vulgatibacter, Devosia, Kaistia, Microvirga, Pseudorhodoplanes, Sphingopyxis) uniquely associated with inorganic amendments compared to other treatments (Table S3), indicating a possible increase in microbial activity related to fertilizer application at early stages of the experiment. However, this difference in indicator taxa composition was dissipated at the end of the experiment (week 15), whereby inorganic fertilized low‐diversity samples shared most of the indicator OTUs with other treatments (Fig. 5, Table S3). By contrast, in the high diversity soils, the composition of the microbial indicator community remained overall similar across treatments at both time points, except for a sharp decrease in the Chloroflexi in the organic‐treated soils at week 15 (Fig. 5, Table S4).

Fig. 5

Bar plot showing the phylum‐level composition of the bacterial indicator community significantly (P < 0.01) and strongly (IndVal > 0.7) associated with different nutrient amendments at weeks 4 and 15 in low and high diversity soils. Detailed assignments are reported in Table S3 and S4. HDS, High Diversity Soil; LDS, Low Diversity Soil.

Direct and indirect impacts of soil biotic and abiotic variables on the degradation rate

We constructed a SEM for the final sampling point (week 15) to infer the potential mechanisms through which nutrient and plant community manipulation, as well as differing initial microbial community diversity, affected pyrene degradation via effects on microbial community structure and diversity (Fig. 6A and B). Our initial model also included the pahAc gene copy numbers; however, the addition of this variable resulted in non‐satisfactory fit metrics, and thus pahAc abundance was dropped in further analyses. We found that initial soil microbial diversity was the main driver of enhanced degradation rate and as it influenced both soil diversity and community structure. Microbial community effects were also influenced by fertilizer and plant + fertilizer treatments through their effects on pH, nitrogen and phosphorus concentration, while plant treatment alone indirectly affected the degradation rate through its effect on nitrogen. The total standardized effects from SEM confirmed that initial soil microbial diversity was a major determinant of degradation rates, followed by nutrient content (nitrogen, phosphorus) and soil pH (Fig. 6B). The different goodness‐of‐fit metrics (chi‐square = 6.486, df = 7, P = 0.518, CFI = 0.999, RMSEA = 0, P RMSEA = 0.630, SRMR = 0.016, Bollen‐Stine = 1) indicated that the a priori model had a satisfactory fit to the data.

Fig. 6

(A) Structural equation model for pyrene degradation rate for the final data collection point (week 15), based on the effects of soil initial microbial diversity, nutrient amendments (inorganic and organic fertilizer), plant presence, soil nutrient (NO3/NH4 ratio and PO4) content and biological attributes (bacterial community diversity and structure). Numbers on arrows are standardized path coefficients, analogous to partial regression weights and indicative of the effect size of the relationship. Dashed arrows indicate negative relationships. Only significant (P < 0.05) coefficients are reported. R 2 indicates the proportion of variance explained.

B. Total (i.e. both direct and indirect effects) standardized effects of each variable on the degradation rate derived from the SEM model.

Discussion

Despite the mounting evidence on microbially derived biodegradation of xenobiotics, most studies have mainly focused the individual role of microbial diversity, nutrient availability and plant presence in performing bioremediation. To the best of our knowledge, no previous study has empirically and statistically examined the relative importance of microbial (bacterial) diversity, nutrient availability and plant presence in enhancing biodegradation potential of soils. This study provides evidence that the initial levels of soil microbial diversity are an important predictor of biodegradation rates and biostimulants efficiency in nutrient‐deficient soils. Our results demonstrated that soils with originally higher bacterial diversity can accomplish complete pyrene degradation within 10 weeks of incubation without any inputs, while in soils with lower bacterial diversity, pyrene degradation was still incomplete after 15 weeks. These results on pyrene biodegradation, a microbially driven specialized function (Delgado‐Baquerizo et al., 2016), are consistent with previous observations, where a change in diversity was reported to reduce the rate of key specialized soil functions such as total petroleum biodegradation (Dell’Anno et al., 2012), pesticide mineralization (Singh et al., 2014), nitrate availability, methane and nitrous oxide flux (Trivedi et al., 2019). Soils with higher microbial diversity underpin a higher rate of biological processes (Delgado‐Baquerizo et al., 2016; Maron et al., 2018), which results in a high rate of biodegradation (Venail and Vives, 2013; Crampon and Bodilis, 2018) and higher availability of N and P (Nannipieri et al., 2003), as observed in our study.

In addition, our results also suggest that initial bacterial structure and diversity levels influence the efficiency of bioremediation‐enhancement approaches such as biostimulation (Wu et al., 2020). Indeed, the impact of two typical inorganic and organic amendments on pyrene degradation rate was markedly distinct between the two soils used in this study, whereby organic amendments had an initial impact on pyrene degradation rates, while inorganic amendments enhanced pyrene degradation but only in the lower diversity soil. By contrast, in the higher diversity soils, the effect of inorganic nutrient amendments was negative compared to untreated soils, despite the fact that the fertilizer‐treated soils had initially more bioavailable P and N. This inconsistent effect from inorganic fertilizer addition could be ascribed to the rapid nutrient metabolization that occurred in the higher diverse soil, in which the pool of added nutrients was rapidly depleted and thus possibly made unavailable to the pyrene‐degrading microbial subcommunity in the critical first few weeks of exposure to pyrene. Conversely, in less diverse soils, inorganic amendments took longer to be metabolized, increasing nutrient availability for pyrene degraders (Sarkar et al., 2005; Kalantary et al., 2014; Lukić et al., 2016), which in turn had a boosting effect on the biodegradation rate (Leff et al., 2015; Zhu et al., 2016; Xue et al., 2018; Sivaram et al., 2019). Consistent with these findings, inorganic nutrient addition in the low diversity soil supported a unique assemblage of bacterial indicator taxa with known bioremediation abilities, such as members of Proteobacteria and Firmicutes (Venail and Vives, 2013; Ren et al., 2016; Zada et al., 2021), while the differentiation of the indicator community in nutrient‐amended samples was relatively weaker in the high diversity soil. Overall, biostimulation of the indigenous microflora is a strategy of choice to achieve effective remediation of PAHs in soils, which are chronically contaminated with diesel (Ruberto et al., 2009), phenanthrene (Kalantary et al., 2014), and other petroleum hydrocarbons (Jiang et al., 2016). However, given the inconsistent effect of nutrient amendments in the two soils used in this study, we suggest that the use of fertilizers to improve bioremediation should be tailored to the specific nutrient requirements of the resident microbial community.

The potential of the C3 and C4 plant species used in this study for remediation of PAHs‐contaminated soils has previously been explored (Sivaram et al., 2018, 2020), and a positive effect of plant addition on PAHs biodegradation was reported. Yet, contrary to our expectation, we observed a very limited effect of plant presence on the rate of pyrene biodegradation in both soils. The phytoremediation potential of plants, however, varies significantly among species and also depends on the edaphic factors such as PAHs concentration in soil, aging and soil properties (Cheema et al., 2010). Importantly, previous studies reported shifts in plant bioremediation abilities due to changes in soil edaphic factors and microbial activities (Afzal et al., 2011), suggesting that discrepancies in plant effect could be ascribed to differences in soil physico‐chemical and biological features between soil types (i.e. sandy, loamy sand and loam). Additionally, in this study, we measured the plant effect on biodegradation on the bulk soil only. However, in the case of phytoremediation processes, the synergetic interactions of plant and rhizospheric microbial communities are crucial (Tang et al., 2010; Mohanram and Kumar, 2019; Kotoky and Pandey, 2020), as they alter the mobility and bioavailability of PAHs while also contributing to their removal by enhancing uptake by plants (Pretorius et al., 2018). As such, it is possible that, rather than in the bulk soil, the plant effect might be more apparent in the rhizosphere space, where PAHs are more bioavailable to plants (Wu et al., 2020). Nonetheless, our observations on the limited impact of plants with known phytoremediation abilities for pyrene bioremediation indicate that soil biotic and abiotic properties might override plants in driving bioremediation outcomes, especially in soil not in direct contact with the plant roots. As such, the importance of such parameters should be considered in the design of bioremediation applications that include the use of plants (Kong et al., 2018; Haider et al., 2021). Overall, to bring the proven remediation technologies into practical management practice, it is critical to acquire more knowledge about the spatial influence and mechanisms of microbially mediated rhizoremediation for the selection of the most appropriate phytoremediation approaches (Yang et al., 2020; Haider et al., 2021; Hoang, Lamb et al., 2021).

Our general findings were further validated by a SEM model, which identified initial soil microbial diversity and composition as the major drivers of biodegradation rates through direct and indirect effects on soil biotic (community diversity and composition) and abiotic‐driven (N and P contents) characteristics. High diversity soils had a higher initial diversity and abundance of phyla harbouring bacteria with known biodegradation abilities, such as Gemmatimonadetes, Chloroflexi, Bacteroidetes and Actinobacteria, which likely supported a higher baseline rate of biodegradation (Bell et al., 2013). To achieve a similar level of biodegradation rate, soils with lower microbial diversity required external inputs of nutrients that increased growth and activity of the resident microbial community (Chaîneau et al., 2005; Leewis et al., 2016). Interestingly, the qPCR results also indicated that the high diversity soil possessed a considerable higher abundance of the functional gene pahAc involved directly in PAH degradation. As expected, (Liang and Huang, 2019), pahAc gene abundance had also a moderate positive correlation with PAH degradation. However, we observed a poor association of degrading gene number and biodegradation rates in our SEM models, indicating that degrading gene abundance is a modest predictor of bioremediation potential in our system. Such decoupling between biodegradation gene abundance and biodegradation rates could be attributed to bias in primer selectivity of PCR amplification, or to the contribution of other enzymes not measured here to the degradation of pyrene, as suggested in recent metagenomics analyses (Suenaga et al., 2009; Cébron et al., 2011; Gao et al., 2021). Additionally, we did not assess gene expression levels, although it is possible that pahAc gene expression, rather than abundance per se, might be a more reliable and accurate measure of the biodegradation dynamics in some cases. Overall, the systems‐based approach undertaken in our study emphasizes the importance of validating management practices to enhance PAHs biodegradation rates, which are directly linked to soil microbial diversity and nutrient availability (Tyagi et al., 2011) in low nutrient soils.

In conclusion, based on comprehensive system‐based approach, our findings provide novel evidence that initial soil microbial diversity and nutrient requirements interact to affect the rate of biodegradation and should be considered explicitly when designing bioremediation strategies including the use of plant and nutrient inputs. Taken together, our study provides a holistic view and a systematic understanding of the importance of soil initial bacterial diversity, nutrient amendments and plant presence when considering an effective approach to enhance the bioremediation potential in natural soils. Our results demonstrated that a systematic approach that explicitly considers soil microbial diversity and nutritional status can help in designing effective bioremediation strategies. Such an approach could facilitate eco‐friendly and economic remediation strategies that lead to rapid remediation in situ without the need for excavation and transport of polluted soil.

Experimental procedures

Soil preparation and spiking

Two soil types were used in this study with contrasting initial microbial diversity: high diversity soil (HDS; soils harbouring higher bacterial diversity) and low diversity soil (LDS; soils with lower bacterial diversity compared to HDS) (Figs 1 and 2). Soil samples were collected from the Hawkesbury Forest Experiment (HFE) site (33°36'41.0"S, 150°44'25.0"E), in Richmond (NSW, Australia). The soil is a sandy loam with relatively low nutrient concentrations (Table 2), low organic matter content (0.7%) and low water holding capacity (Barton et al., 2010; Hu et al., 2015). A complete description of HFE plots setup and soil characteristics can be found in Hu et al. (2015). Briefly, HFE is a forest plantation experiment established in 2007 that comprises four soil management treatments (control, irrigated, fertilizer‐treated, irrigated with fertilizer), replicated in four plots (38.5 × 41.6 m) planted with Eucalyptus saligna monocultures in a randomized complete block design. For this study, soils were collected from control (unmanaged) and irrigated plots based on the observations from previous work, which demonstrated contrasting bacterial diversity and structures between irrigated (HDS) and unmanaged (LDS) soils (Colombo et al., 2016). These two soils have similar pH and carbon levels, but differ in , and availability (Table 2).

Table 2

Soil physico‐chemical properties among unmanaged (LDS) and irrigated (HDS) soils presented as mean values ± standard error (n = 12).

Soil properties	Low diversity soil (unmanaged)	High diversity soil (irrigated)
pH	7.4 ± 0.11	7.6 ± 0.09
Total C a	10.92 ± 2.09	10.34 ± 1.59
Total N a	0.78 ± 0.14	0.76 ± 0.04
PO43‐ b	9.3 ± 0.65	2.2 ± 0.5
NH4+ b	4.17 ± 0.17	1.83 ± 0.22
NO3‐ b	5.79 ± 0.21	9.0 ± 0.35

Colombo et al. (2016).

Presented in mg per kg‐1 of dry soil.

For each low (unmanaged) and high diversity (irrigated) soil, 150 kg of soil (0–10 cm) was collected from three replicated plots in September 2018. Soils were subdivided into two lots, and one lot was spiked with pyrene, dissolved in hexane to yield a final concentration of 200 mg kg‐1 (Sivaram et al., 2018). After adding pyrene, soils were mixed thoroughly to achieve homogeneity, and they were left uncovered overnight until the solvent was evaporated completely. The remaining soil was used for untreated controls (Fig. 1). We did not measure the pyrene immobilization in each type of soil. However, as the soils considered in this study had nearly identical levels of nutrient and carbon content, we expected minimal effects of soil physico‐chemical parameters on pyrene mobilization.

Microcosm setup with nutrient amendment and plant treatments

The microcosm experiment was set up in a climate‐controlled growth chamber (Biochambers, Manitoba, Canada). The experimental design involved the use of two common soil nutrient amendments (organic and inorganic fertilizers) to assess the impact of nutrient availability, and two native grass species (Rhytidosperma sp., and Dichanthium sp.) to assess the influence of plant presence and type (photosynthetic pathway) on enhancing the rate of xenobiotic biodegradation in two soils with contrasting initial microbial diversity (Fig. 1). For the preparation of microcosms, we selected two common commercial inorganic and organic amendments, Thrive® (Yates, NSW, Australia) and Converte™ (Converte Pty Ltd, Griffith, NSW, Australia) respectively. Thrive® is an all‐purpose liquid plant food, which mainly constitutes potassium dihydrogen phosphate (H2KO4P), potassium and urea (CH4N2O) and Converte™ is an organic liquid plant food derived from natural resources, containing organic carbon, potassium and nitrogen in addition to mineral trace elements, phytoproteins, vitamins and plant hormones. Both inorganic and organic fertilizers were applied to the corresponding soil treatments following the manufacturers’ recommendations and mixed well. The inorganic fertilizer treatment received N, P and K at a rate of 95, 20 and 49 kg ha‐1, respectively, while the organic (biostimulant) treatment received C, N, P and K at a rate of 0.58, 0.006, 8.25 × 10‐5 and 0.007 kg ha‐1 respectively. Microcosms were prepared in pots (inner diameter 13 cm, height 15 cm; 1.8 L) filled with approximately 1.7 kg soil. To explore the phytoremediation potential of C3 and C4 grasses and the interactive effect of rhizoremediation on the rate of PAHs biodegradation, we used two Australian widespread native grass species (C3 – Rhytidosperma sp., and C4 – Dichanthium sp.) with known phytoremediation abilities (Sivaram et al., 2018, 2020; Sivaram et al., 2018). Plant treatments were established by sowing ~ 12 seeds directly into each pot. Seedlings were thinned after three weeks to achieve 6 plants per pot. The first two plant treatments were set up using monocultures of C3 and C4, with a combination treatment including both these plant species, in addition to an equal number of no‐plant control pots. In total, we created 48 different treatments for a total of 144 microcosms (two soil types [low and high diversity] × two xenobiotic treatments [pyrene, no pyrene] × three nutrient amendment treatments [no addition, inorganic, organic] × four plant treatments [C3, C4, C3C4, no plant] × three biological replicates) (Fig. 1). The experimental design allowed us to quantify the relative importance of initial soil microbial diversity, nutrient amendment and plant presence/type on pyrene biodegradation rate and to identify which combination (i.e. higher/lower microbial diversity, +/− nutrient amendment, +/− plant presence) (Fig. 1) expedited pyrene degradation.

Plant growth conditions and sampling

Daily average temperatures of growth chamber were maintained at 24°C (day) and 16°C (night), and daylight settings were set at 16/8 h. All the pots were maintained in well‐watered conditions (~ 60% of water holding capacity) throughout the experiment. Soils were sampled at four time points: week 0; (i.e. before sowing the seeds), week 4 and week 10 (intermediate sampling points), and the final harvest at week 15 of the experiment. During the intermediate sampling, ~ 15 g of bulk soil (i.e. soil not in contact with the roots) was collected by a punch hole of 8 cm with a 15‐ml centrifuge tube, and at the final harvest, sufficient soils were collected after separating plant roots by sieving. At each time of sampling, soils were subdivided into three portions. The first portion was immediately processed for pyrene quantification, the second was stored at 4°C before soil physico‐chemical properties analysis, and the third was stored at −20°C for DNA extractions.

Extraction and analysis of pyrene from soil to monitor pyrene biodegradation rate

Extraction of pyrene from soils was performed by modified acetone: hexane extraction method (Haleyur et al., 2016). Briefly, one gram of soil was taken in a glass vial and spiked with 100 µg of phenanthrene (Sigma‐Aldrich, Castle Hill, NSW, Australia) as an internal standard. Extracted twice with 4 ml of acetone: hexane (1:1 v/v) for 20 min by vortex and intermediate shaking. The supernatant was collected in a different tube and evaporated using a CentriVap vacuum concentrator and the residues were dissolved in 1.5 ml of acetonitrile. Finally, samples were filtered through a 0.45 µm, 13‐mm diameter Millex® nylon syringe filter (Merck Millipore, Darmstadt, Germany). Extractions were done in duplicate and were analysed by Gas Chromatography‐Mass Spectrometry (GC‐MS) on an Agilent 7890A gas chromatography system coupled to the Agilent 5975C inert XL MSD Triple Axis Detector (Agilent Technologies, CA, USA), using a DB‐17MS capillary column (60 m × 0.25 mm i.d. × 0.25 μm film, (50%‐phenyl‐methyl polysiloxane), Agilent J&W GC Columns). Helium was used as the carrier gas and maintained at a flow rate of 1 ml min‐1. The column temperature was set to initial hold at 60 °C for 1 min and then raised to 150 °C at 6 °C min‐1, 170 °C at 2 °C min‐1, and to 300 °C at 15 °C min‐1 with a final hold for 2 min. Injector and detector temperatures were 230 and 280 °C, respectively in a split mode (10:1) with a split flow of 12 ml min‐1. Pyrene was identified by comparing the retention time and mass fragmentation pattern of standard (Sigma) and quantified by comparing with that of the internal standard. The pyrene concentration µg g‐1 of soil was calculated by dry weights and the percentage of pyrene degradation was calculated using the following equation: where CI was the initial concentration of pyrene in soils at the beginning of the experiment, and CT was the final pyrene concentration after the number of respective weeks at each time point. Pyrene degradation was monitored over 15 weeks by measuring the pyrene removal (%) at week 4, week 10, and the final harvest at week 15 by comparing the respective degradation levels with initial levels.

Soil physico‐chemical properties

Soil pH was measured using 2.5 ± 0.01 g of each fresh sample mixed with 12.5 ml of Mili‐Q water and kept on a shaker for 1 h at 180 r.p.m. and allowed to settle before reading with a Delta pH meter (Mettler‐Toledo Instruments, Columbus, OH, USA). Soil moisture content was calculated and expressed as the percentage of soil weight by measuring the difference in weight of soil samples before and after oven drying at 105°C for 24 h.

Extractable NO3, NH4, and PO4

Extractable nitrate () and ammonium () were determined on fresh soils extracted with 2 M KCl (Keeney and Nelson, 1982). Extractable phosphate () was determined from air‐dried soil (40°C for 48 h) by the Colwell‐P method (Colwell, 1963). Filtered extracts (Whatman No. 42) were analysed for NO3 ‐, NH4 +and PO4 3‐ on a SEAL AQ2 Analyzer (SEAL Analytical, Maquon, WI, USA).

Microbial community structure and diversity

Total genomic DNA was isolated from soils at week 0 (reference pyrene‐untreated soils only), week 4 (pyrene‐treated soils only) and week 15 (pyrene‐treated soils only) using DNeasy PowerSoil Pro® DNA Isolation kit (Qiagen Inc., Valencia, CA, USA) according to the manufacturer’s instructions with a slight modification, whereby the initial cell lysis step was performed by bead beating (speed of 5.5 m s‐1 for 30 s) using a FastPrep bead beating system (Bio‐101, Vista, CA, USA). The quality and quantity of DNA extracts were assessed using a NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). For samples collected at weeks 4 and 15, the bacterial diversity and community composition were determined using 16s rRNA gene amplicon sequencing with primers 341F (5’ CCT ACG GGN GGC WGC AG 3’)/805R (5’GAC TAC HVG GGT ATC TAA TCC 3’) (Klindworth et al., 2013). Library preparation and sequencing on an Illumina MiSeq platform (2x300 PE) were performed by the Next‐Generation Sequencing Facility of the Western Sydney University (Richmond, NSW, Australia).

For raw pair‐end reads, primers at the beginning and end of each sequence were trimmed off and the quality control of raw reads was performed using the USEARCH pipeline (Edgar, 2010). Briefly, the maximum of expected error (ee) was set as 1.0 for the merged reads filtering as default. zOTUs (denoised sequences, or amplicon sequence variants; ASVs) were gained by denoising (error correction) the amplicon reads using unoise3 at 100% sequence identity (Edgar, 2016). Chimeras were examined and discarded using the UCHIME algorithm. zOTUs’ representative sequences were annotated against the Silva database (Quast et al., 2013) in QIIME (Caporaso et al., 2010) using UCLUST (Edgar, 2010). Approximately 8.95 M high‐quality merged sequences were mapped for all the samples. A normalization procedure was then performed at a sequencing depth of 10,000 reads per sample before downstream analyses; 5 samples were discarded due to the low yield. Microbial alpha diversity indices (Shannon Index, Richness, and Chao1 index) were derived from the ‘diversity’ function in the R phyloseq package (McMurdie and Holmes, 2013).

All the data and associated metadata are publicly available at 10.6084/m9.figshare.15145086.

PAH‐RHD α subunit quantification

To assess the impact of functional genes directly involved in PAH degradation, we investigated whether there was any correlation between the abundance of microbial pyrene degradation genes and biodegradation rate. The abundance of functional genes was obtained using a quantitative‐PCR approach. Following Ding et al. (2010), we targeted the PAH‐RHD α subunit (pahAc gene) of the bacterial PAH ring‐hydroxylating dioxygenases, a widespread enzyme involved in bacterial‐mediated biodegradation of PAHs. Notably, the PAH‐RHD α subunit drives the initial oxidation of PAHs, a rate‐limiting step in PAH degradation. As such, this subunit is generally used as a functional biomarker of PAH degradation (Debruyn and Chewning, 2007; Cébron et al., 2008; Liang et al., 2019). For this analysis, we focused on samples at week 10 as those samples showed the most pronounced differences in degradation rates among the pyrene‐treated samples.

Briefly, the qPCR assay was conducted in 96‐well plates on a CFX96 Touch™ Real‐Time PCR Detection System (BioRad, Australia). Each 10 μl‐reaction contained 5 μl of SensiFAST™ SYBR no‐ROX Master mix (Bioline, Australia), 0.9 μl nuclease‐free water (Ambion, Life Technologies, Australia), bacterial primers AH‐RHDα‐396F (5′‐ATT GCG CTT AYC AYG GBT GG‐3′) and PAH‐RHDα‐696R (5′‐ATA GGT GTC TCC AAC RAA RTT‐3′) (Sigma‐Aldrich) at a final concentration of 250 nM, and 5 ng of DNA template. The PCR conditions were 95°C for 3 min followed by 40 cycles of 95°C for 15 s, 53°C for 15 s and 72°C for 30 s. Fluorescent signals were measured after this last step at each cycle. Each plate included duplicate reactions per DNA sample, set of standards and negative control samples (without DNA). Standard curves were generated using a tenfold serial dilution of an amplicon obtained from a composite soil sample composed by pooling 10 soils randomly selected among the high diversity and low diversity soils. Melt curve analysis was conducted after each assay (ramp from 65 to 95°C, rising by 0.5°C each step) to confirm the specificity of amplification products.

Statistical analysis

To investigate the impact of initial soil microbial diversity, plant treatment and amendment treatments on biodegradation rate, soil microbial diversity and soil physico‐chemical properties over the entire study, we used a repeated measures analysis ANOVA, estimating effect sizes – the relative importance of the manipulated factors – with partial eta‐squared (η2). Time (weeks) and the random effect of pot identity were also added to the analyses, including two‐way interactions between weeks and the experimental factors. When the effect of manipulation treatments differed between dates, we also conducted separate ANOVAs for each time point, followed by Sidak’s multiple comparison tests to perform pairwise comparisons and Tukey’s multiple comparison tests to determine significant interactions across treatments, using GrapPad Prism 9.1 software. Similarly, the treatment effect on bacterial community composition was assessed with PERMANOVA. Pairwise Bray‐Curtis dissimilarity metrics were calculated and differences in bacterial community composition across time and treatments were estimated using a PERMANOVA with 9999 permutations, using the ‘adonis’ function in the vegan package (Oksanen et al. 2007). Bacterial community composition was visualized using principal coordinate analysis (PCoA) based on the ASV feature tables in RStudio (Version 1.1.442) using the phyloseq package (McMurdie and Holmes, 2013).

To identify the pyrene‐associated bacterial taxa that were enriched by nutrient amendment treatments at early (week 4) and later (week 15) stages of the experiment, we used an indicator analysis combining both the abundance and occurrence of a given genus across all treatments in the pyrene‐treated soils. The indicator (IndVal) values were calculated using the multipatt function in the ‘indicspecies’ R package, with n = 999 random permutations (Cáceres and Legendre, 2009; Cáceres et al., 2012), and only significant (P < 0.01) and strong (IndVal > 0.7) associations were retained. The relationships between pyrene biodegradation rate and pahAc gene abundance was assessed by Spearman correlation using the R 4.0.2 statistical software (R Core Team, 2018).

Structural equation modelling

For the final sampling point (week 15), structural equation modelling (Grace, 2006) was used to build a system‐level understanding and evaluate the multiple effects of plant and nutrient amendments (fertilizer, plant and interaction of plant and fertilizer) and initial soil microbial diversity on the pyrene biodegradation rate, acting via effects on soil nutrient content (measured as changes in NO3/NH4 ratio and PO4 concentration), pH and microbial composition and diversity assessed via 16S rDNA amplicon sequencing. The full rationale for the model is reported in the appendices (Appendix S1). In the models, the initial soil microbial diversity was considered as a categorical variable (0 = low diversity soil and 1 = high diversity soil). Moreover, the different manipulations (plant presence, nutrient amendments and interaction of nutrient x plant treatment) were treated as categorical variables with two levels: 1 (e.g. nutrient amendment) and 0 (remaining treatments + no treatment control). This approach allowed us to compare the effect of a particular manipulation on the pyrene degradation rate with the average of the remaining treatments. Shannon index of each bacterial community (normalized to 10,000 reads per sample) was used as a measure of alpha diversity. We also performed principal coordinates analysis (PCoA) in the R package vegan (R version 4.0.2) (Dixon, 2003) based on the Bray‐Curtis dissimilarities of square‐root‐transformed bacterial community operational taxonomic units (OTUs) composition and PCoA scores of the first and second axis were used as proxies for changes in community composition in structural equation model (SEM) analysis. Before conducting SEM, soil PO4, soil NO3/NH4 ratio and pyrene degradation rate were square‐root‐transformed to improve linearity. The model was then parameterized and its overall goodness of fit tested using the Chi‐square test, the Comparative Fit Index (CFI), the Root Mean Square Error of approximation (RMSEA), and the (Standardized) Root Mean Square Residual (SRMR), as suggested in Kline, 2012. Additionally, because some variables resulted in residuals that were not normally distributed, the model fit was confirmed using the Bollen‐Stine bootstrap test. To further identify the relative importance of abiotic and biotic drivers in controlling degradation rates, the total standardized effects (i.e. both direct and indirect effects) of each variable were also calculated from SEM. All SEM analyses were conducted using the R package lavaan (Rosseel, 2012).

Conflict of interest

None declared.

Fig. S1. Estimated diversity (a, Shannon), and richness (b, Chao; and c, observed OTUs) of the bacterial communities in initial soils (week 0). One‐way ANOVA determines significant differences between initial microbial diversity (IMD).

Fig. S2. Estimated diversity (a, Shannon), and richness (b, Chao; and c, observed OTUs) of the bacterial communities at week 4, compared between initial microbial diversity (IMD) of soils, across nutrient amendment and plant treatments by three‐way ANOVA (tables d, e and f respectively). Sidak’s multiple comparisons test determines the significant differences among IMD levels at ****P < 0.0001, ***P < 0.001, **P < 0.01.

Fig. S3. Estimated diversity (a, Shannon), and richness (b, Chao; and c, observed OTUs) of the bacterial communities at week 15, compared between initial microbial diversity (IMD) of soils, across nutrient amendment, and plant presence treatments by three‐way ANOVA (tables d, e and f respectively). Sidak’s multiple comparisons test determines the significant differences among IMD levels at ****P < 0.0001, ***P < 0.001, **P < 0.01, *P < 0.05.

Fig. S4. PCOA ordination of bacterial communities at (a) week 4, and (b) week 15, separated by initial microbial diversity (IMD) and treatments. PERMANOVA results are reported for week 4 (c), and week 15 (d). ***P < 0.01

Fig. S5. Impact of initial microbial diversity (IMD) of soils and nutrient amendments on (a) Ammonium, (b) Nitrate, (c) Phosphate levels presented in mg per kg‐1 of dry soil, and (d) soil pH in (i) untreated, (ii) inorganic fertilizer and (iii) organic fertilizer nutrient amendment treatments. Two‐way ANOVA was performed in each nutrient treatment, followed by Sidak’s multiple comparisons test compares significant differences between IMD mentioned as ****P < 0.0001, ***P < 0.001, **P < 0.01 and, *P < 0.05. Different coloured letter compares significant differences across time points in high and low diversity separately, indicating mean contrasts (significant at P < 0.05) using Tukey’s multiple comparison test.

Fig. S6. Impact of plant presence on Ammonium (a), Nitrate (b), Phosphate levels (c) presented in mg per kg‐1 of dry soil, and soil pH (d) in untreated (i), inorganic fertilizer (ii) and organic fertilizer (iii) nutrient amendment treatments. Three‐way ANOVA was performed for each plant combination, followed by Tukey’s multiple comparisons test to compares significant differences by plant presence mentioned as ***P < 0.001, **P < 0.01 and, *P < 0.05, and ns – nonsignificant.

Fig. S7. Influence of plant presence and soil amendments on pyrene degradation rates. Interleaved symbols at mean with SD showing pyrene degradation rates across plant combinations at different time points in (a) untreated control (b) Inorganic fertilizer and (b) organic fertilizer. Letters indicate mean contrasts using Tukey’s multiple comparison test and differences are significant at P < 0.05.

Fig. S8. Boxplots depicting differences in relative abundance of the major bacterial phyla in the reference high and low diversity soils at week 0 (t‐test). *P < 0.05; **P < 0.01; ns = nonsignificant.

Table S1. Repeated measure ANOVA comparing the role of interactions among soil initial microbial diversity, nutrient amendment, plant presence and time (weeks) on pyrene biodegradation rate. IMD, Initial Microbial Diversity.

Table S2. PERMANOVA comparing the role of interactions among soil initial microbial diversity, nutrient amendment, plant presence and time (weeks) on bacterial community structure. IMD= Initial Microbial Diversity

Table S3. Heatmap showing the top (InVal > 0.7, green square) indicator OTUs significantly (P < 0.01) associated with no fertilization (Un), inorganic fertilization (Inor) and organic fertilization (Org) at week 4 in low diversity soil.

Table S4. Heatmap showing the top (InVal > 0.7, green square) indicator OTUs significantly (P < 0.01) associated with no fertilization (Un), inorganic fertilization (Inor) and organic fertilization (Org) at week 15 in low diversity soil.

Appendix S1. SEM Rationale.

Click here for additional data file.