Long-term dynamics of plant communities after biological remediation of oil-contaminated soils in far north

We studied the long-term dynamics of plant communities after bio and phytoremediation of oil-polluted soils. Nine plots located in European Northeast and treated using various bioremediation methods were monitored from 2002 to 2014. Geobotanical descriptions (relevés) of each plot were performed in 2006 and 2014, and Grime’s theoretical CSR (competition–stress–ruderality) framework was used to assess the vegetation state and dynamics. We observed a clear shift of communities from pioneer (where ruderal species were prevalent) to stable (where competitor species were dominant) states. However, the remediation type did not significantly impact the vegetation recovery rate. After 12 years, all methods led to a 55–90% decrease in the oil content of the soil and a recovery of the vegetation cover. The plant communities contained mainly cereals and sedges which significantly differed from the original tundra communities before the oil spill. The control plot, treated only by mechanical cleaning, had minimum oil degradation rate (50%) and vegetation recovery rates, although, in CSR terms, its vegetation assemblage resembled the background community. Cereals (Agrostis gigantea, Deschampsia cespitosa, Phalaris arundinacea, and Poa pratensis), sedges (Carex canescens, Carex limosa, and Eriophorum vaginatum), and shrubs (Salix) were found to be the most effective species for phytoremediation, exhibiting high community productivity under the harsh northern conditions.


Materials and methods
The experiments were conducted at the Verkhnevozeyskoye oil field (Usinsky district, Komi Republic: 66°37′40″ N, 57°07′56″ E), in the Lukoil Usinskneftegaz area, where oil spills occurred several times between 1989 and 1996 39 . This territory belongs to the forest-tundra zone and is characterized by harsh climatic conditions: low temperatures, strong winds, and a short growing season 49 .
The contaminated area of the experimental site was approximately 2 ha and located on a high peat bog with a peat deposit thickness of up to 1.5-2.0 m, in which the oil penetrated to a depth of 1.0 to 1.5 m. After technical recultivation (removing of surface oil) and land melioration (removal of excess water and oil from the soil by drainage), the site was plowed to a depth of 25-30 cm and then divided into test plots of 0.2 ha each (40 m × 50 m).
The experiments for assessing the efficiency of the various bioremediation methods (biopreparations, sorbents, and agricultural techniques) started in June 2002 (Table 1). In the beginning, the soil moisture ranged from 65 to 88%, while the oil concentration ranged from 150 to 450 mg g −1 39 . The oil concentration of the soil was determined by fluorimetric method using a Fluorat-02 fluid analyzer (Russia) 50 .
During the next few years, the soil moisture changed unevenly due to the fragmented siltation of the drainage channels. This caused variations in the distribution pattern of the cultivated and invading native plants. Mineral fertilizers, introduced in 2002, ensured the start of growth of sown plants and the efficient oxidation of the oil by microorganisms. Further feeding was not carried out. The development of the sown and invading plants was determined by the soil moisture level and the nature of the residual pollution which correlated with the amount of the oil-destructing bacteria. Soil contamination was uneven and varied greatly, even within the same plot (Table 2).
Two relevés were performed for the testing plots, four and twelve years into the experiment (in 2006 and 2014, respectively), using the standard geobotanical methods. We focused on the vascular plants which were the predominant species in the plots. The species nomenclature used hereafter is in accordance with www.thepl antli st.org.
We selected a willow-dwarf birch sedge-horsetail swamp as the reference plot. This was located on a site with a similar to experimental plots landscape, away from the oil-contaminated area. The most abundant species in the swamp were Salix myrtilloides, S. lapponum, Betula nana, Equisetum fluviatile, Menyanthes trifoliata, Comarum palustre, Andromeda polifolia, Eriophorum vaginatum, and Calamagrostis lapponica. The moss-lichen layer consisted mainly of sphagnum mosses (Sphagnum magellanicum and S. fuscum), with small inclusions of Polytrichum commune. The water content of the reference plot was 80-89%.
Full relevés for the experimental and reference plots are provided in the supplementary materials (Suppl. 1). The integral position of relevés in the CSR triangle was estimated using the community-weighted mean (CWM) method 55 . To reduce the influence of dominant species on the resulting scores, we used the logarithmic transformation ln (A + 1), where A denotes the species abundance 56 . The CSR species types were determined according to J.G. Hodgson 57 . Furthermore, each CSR type was presented in tertiary coordinates form from 0 to 100%, where 0 signifies the absence of traits of this type, and 100% is a clearly defined strategy 58 . The CWM score is defined as follows: www.nature.com/scientificreports/ where A i is the abundance of the i-th species in the relevés, B i is the score value along the C, S, R axes of the Grime's triangle for the i-th species, and N is the number of species. The species diversity of the relevés was estimated based on the number of recorded species and the Shannon species diversity index 59,60 , which is expressed as follows: where H is the Shannon index, p i is the contribution of the i-th species to the plant community (the ratio of a species abundance to the sum abundances of all species of a given community), and N is the number of species.

Results
The oil content of the soil in the treated plots decreased from 200-360 mg g −1 in 2002 to 30-125 mg g −1 in 2014. This decrease occurred unevenly over time. Over the first period (2002)(2003)(2004)(2005)(2006), statistically significant changes were observed in all the experimental plots, except plot 5 and the control. During the second period (2006-2014), the oil content significantly decreased in the control and only for two of the treated plots (Table 2). By the fourth year of the experiment (2006), the cereals seeded during the remediation (Agrostis gigantea, Deschampsia cespitosa, Phalaris arundinacea, and Phleum pratense) became dominant. Also, the presence of a single species, which was not used in remediation, was noticed. These are intermediate CSR species (Rumex acetosella and Stellaria graminea) and ruderal types (Lapsana communis, Rumex crispus, Plantago major, and Tripleurospermum perforatum) (Fig. 1). In 2006, the total projective cover at the experimental plots was 70-80% (Table 3).
By 2014, four of the six species used for remediation had retained their presence. The Phleum pratense (CSR type according to Grime's classification) was found to be less abundant; Agrostis gigantea (CR type) was more abundant. Other seeded cereals: Deschampsia cespitosa (CS/CSR) and Phalaris arundinacea (C) were as abundant as at the beginning of the experiment.
Two species, Avena sativa (−) and Trifolium pratense (CSR) had the lowest survival rate. Avena sativa was sown on six out of nine experimental plots. In 2006, this species was observed in seedling form only on 3 plots Compared to 2006, by 2014 the total projective cover of vegetation had slightly decreased to 60-70%. By then, the number of species and the Shannon biodiversity index had also decreased for all plots, except for plot 9 ( Table 3). The number of species increased from 19 (2006) to 22 (2014). The control plot (plot 2), where vegetation was completely absent in 2006, also showed an increase in the number of species, having up to 17 species of vascular plants in 2014. However, the 2014 projective cover of vegetation at the control plot was only 25-30% smaller than that at the other experimental plots. In the reference plant community, the number of species remained stable.
The relevés analysis by NMDS ordination revealed the presence of two distinct groups (Fig. 4). All the experimental plots of 2006 were in the first group, while those of 2014 were in the second cluster. Reference and control plots lie separately on the NMDS diagram.
C, S, and R vectors indicate the correlation between NMDS axes and community-weighted mean (CWM) scores calculated for relevés.
Changes in the species composition and abundance over time led to significant changes in the average CWM scores along the C, S, and R axes ( Table 4). The competitive (C) score increased from 37-45% in 2006 to 49-55% in 2014, and the corresponding ruderality (R) score decreased from 30-40% to 10-30% (Fig. 5). The average score along the stress tolerance (S) axis did not significantly change and was 18-25% for all experiments. It should be noted that all the values for the experimental plots are substantially less than those of the control (C-39%, S-51%, R-10%) and reference (C-29-36%, S-60-62%, R-4-9%) plots.
For all plots, the individual analysis of the CWM changes revealed an approximately equal increase in the competition scores (10-20%) and a decrease in the ruderal scores (10-30%), as shown in Fig. 6. Along the stress tolerance axis, these changes are differently directed. Due to the increase in Deschampsia cespitosa (CS/CSR strategy) abundance and the appearance of Eriophorum vaginatum (S/CS strategy), the stress tolerance scores increased for most of the experimental plots. However, we observed a slight decrease in the stress tolerance scores for plots 7 and 8. This change was caused by a decrease in the Phleum pratense (CSR strategy) abundance.
The calculated Euclidean distance revealed the number of changes for all experimental plots, between 2006 and 2014, in CSR coordinates. The minimum difference, 13.4, was observed in plot 7 (agro-stimulation without biopreparation). The maximum changes, approximately 30, were observed in plots 1 and 4. For the reference plot, the difference was 8.8 (Table 3; Fig. 6).
An analysis of the trajectories in the CSR triangle reveals a slight converging of the experimental plots results with those of the control and reference territories (Fig. 6). It should be noted that the control plot has the closest position to the reference territory.

Discussion
During this study, visible oil pollution manifestations disappeared at all test plots, except for plot 8, which contained the largest initial amount of oil. In the first period (2002)(2003)(2004)(2005)(2006), almost all experimental plots (except control and plot 5) showed a statistically significant decrease in the oil content. However, in the second period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), only the control plot and two out of the eight plots treated with biopreparation showed significant changes ( Table 2).
A single use of biopreparation in 2002 encouraged the rapid development of bacteria which process oil degradation 39,63 . In the control plot, there was no such boost and hence, the oil decomposition rate was lower. Later, the bacterial count in the control plot increased naturally, which in turn increased the oil decomposition rate. www.nature.com/scientificreports/ By 2014, the total oil reduction in the control plot was 51%, which was only slightly lower than some biopreparations. For example, the minimum decomposition rate of 54-55% was observed in plots 3 and 5. In contrast, the www.nature.com/scientificreports/   www.nature.com/scientificreports/ maximum oil decomposition rate of 90-92% was observed in plots 6 and 9 ( Table 2). These plots were treated with the "Universal" (Institute of biology Komi Science Center) and "Roder" biopreparations (the Chemical Faculty of Moscow State University). However, it should be noted that by 2014, the residual oil contamination of all the plots did not return to the initial values and was higher than that of the reference plots. The vegetation recovery rate varied widely and depended on the type of remediation technology used and the soil water content. By 2014, the projective cover of the vegetation of the biopreparation-treated variants was 60-70%. The control plot treated only by mechanical cleaning had only a 25-30% projective cover. In 2006, the average projective coverage at the experimental plots was slightly higher than that of 2014 and was 70-75% (Table 3). This species-abundance decrease was associated mainly with the water regime restoration outcome of the siltation of drainage channels, and an increase in the groundwater level.
Along with the decrease in projective cover, we observed a decline in the plant species number from 15-22 in 2006 to 6-17 in 2014. The exception was plot 9 (Table 3), where the total projective cover did not change, and the number of plant species increased.
By 2014, Avena sativa and Trifolium pratense, which were sown during phytoremediation, had disappeared. Oats (Avena sativa) is an annual species with exceptional seed regeneration. Under the condition of cold and short summer typical of forest-tundra, oats do not have the time to form seeds. Therefore, all observed seedlings  www.nature.com/scientificreports/ were from seeds that did not germinate in the previous years. The depletion of the seed bank led to the complete disappearance of this species. Red clover (Trifolium pratense) is a perennial species which reproduces by seeds and vegetatively. Clover grows in this geographical area but under milder conditions, usually found in relief depressions or floodplains of large rivers. On the plateau, this species does not survive.
In addition, we observed that the ruderal species had almost completely disappeared. These species include the pioneer plants (annual and biennial), which first appeared on open ground due to rapid seed regeneration, and the perennial weeds, which reproduce mainly vegetatively. In 2006, there were 11 species with prevalent ruderal strategy (R, R/CR, or R/CSR). By 2014, there were only three of such species (Suppl. material). The presence of the ruderal species on the plots in 2006 could be attributed to the soil milling performed before the experiment.   www.nature.com/scientificreports/ The vegetation cover was mechanically destroyed, and the restoration succession began anew, leading to the appearance of these species. The ruderals were gradually replaced by species exhibiting competitive, stress-tolerant or intermediate strategies. This is the classic behavior of plant communities predicted by Grime's theory 43 . Under productive ecological conditions, ruderal species are replaced by competitors, due to their ability to capture most of the available resources (water, light, mineral nutrients). In chronically unproductive habitats, ruderal species are replaced by stress-tolerant ones, which have a high survival rate in the long term.
In 2014, on the experimental plots, we observed the presence of the species which were used in remediation (Agrostis gigantean, Deschampsia cespitosa, Phalaris arundinacea, Phleum pretense) and the appearance of several new species: sedges (Carex canescens, Carex limosa, Carex rostrata), and cereals (Calamagrostis purpurea, Poa pratensis). These species were typical native flora, having competitor or stress-tolerant strategy, and showing high resistance to oil contamination 41,64,65 . Among the newly observed species, the most widespread (6 plots out of 9) was Eriophorum vaginatum. This species is reportedly tolerant of oil contamination 7,66,67 , at least in the northern regions.
A possible reason for the increase in the diversity of sedges family species was the waterlogging of experimental plots which reduced the total project cover and decreased the competition between species. Another possible advantage of sedges is their resistance to oil spills following their natural sustainability despite oxygen deficiency in their root layers 68 .
In 2014, the first trees and shrubs of the Salix genus appeared on the experimental plots (Fig. 3). Willows have good resistance to the negative effects of oil contamination. They grow quickly, have deep-rooting characteristics and are often used in phytoremediation 4,40,69,70 . The appearance of willows and trees indicated the restoration of the vegetation cover.
Another sign of restoration is the detection of the first mosses (Calliergon giganteum and Polytrichum commune). They covered only 5-10% of the study plots, which is less than their coverage in background communities. Mosses and lichens are known for being the most vulnerable to oil pollution 5,66 and, therefore, their appearance indirectly indicated a restoration of the contaminated areas 71 .
The NMDS ordination (Fig. 4) revealed a pronounced gradient of the plant communities studied. The upper right corner shows the relevés of 2006. They are characterized by the prevalence of R-species. The relevés of 2014 display high C-scores and an intermediate position in the ordination space. The reference and control communities are dominated by S-species; they are located in the lower-left corner of the NMDS diagram. In other words, for all the experimental plots, we observed a clear shift from the prevalence of R-species towards the dominance of C and CS species. These changes are characteristic of abandoned areas 48 . The overgrowing occurs primarily due to competitive cereals (Agrostis gigantea, Phalaris arundinacea, Deschampsia cespitosa, and Poa pratensis). These species have good survival rates, produce relatively large biomass, and form a dense turf, impeding the introduction of typical bog and tundra species with S-strategy prevalent in the native vegetation.
At the control plot, where no biopreparation and plant seedling were used, we observed the typical tundrabog S-species such as Empetrum hermaphroditum, Eriophorum angustifolium, Eriophorum vaginatum, Vaccinium myrtillus, and Lycopodium clavatum, as well as young birch, spruce, and willow species. The vegetation at this plot displayed a low projective cover and the latest overgrowing onset 63,72 . However, it had the closest position to the reference plots in terms of ordination (Fig. 4) and CSR scores (Fig. 6).
Phytoremediation employs plant species that display the best resistance to oil contamination, that have a high growth rate and high productivity, and can participate in the biological cleansing processes of the soil. Thus, sedges and cereals are the best candidates 29,64,65,71 for heavily oil contaminated areas. Shrubs 40 , trees 24,73,74 , and crop plants 26,75 are usually cultivated in areas with low oil concentration. Matching the plant species used for remediation with the natural local flora is not usually considered. For example, in this experiment, we used C-species cereals, which did not allow the native S-plants to populate the recovery area. This led to the replacement of natural vegetation with an intrazonal one. In addition, to restore soils, land melioration (partial drainage of wetlands) is often carried out. This leads to a change in the hydrological regime of the soils and, as a result, stimulates changes in the plant communities.
In our opinion, the leading indicator of the eligibility of different recovery methods should be the expected long-term impact on oil-contaminated areas. The main factors to take into account are the biological features of the plants, such as resistance to oil contamination and potential biomass productivity, as well as their affinity for the natural flora 42 . The restored community should resemble, as close as possible, the original communities located on the affected territory before the oil spills. It is preferable to use plant species from the local flora, such as mire sedges or willows. They have a fast growth rate and do not significantly change the plant communities' appearance and physico-chemical properties.

Conclusions
The remediation methods applied on the experimental plots reduced the oil content of the soil by 55-90% and allowed the development of productive plant communities. Despite the use of various methods of soil remediation and the amount of residual oil pollution, the vegetation recovery showed the same patterns over the twelve-year period. The vegetative projective cover in 2014 was approximately 50-70% and vegetation assemblages were also similar.
For all the plots, a significant increase in the number of competitive species and decrease in ruderal species was noticed. This indicated the clear shift of the communities from pioneer to stable state. The relatively small amount of stress-tolerant species and the resulting low scores along the S-axis indicate significant differences between the background vegetation (S-species prevalence) and the recovered plant communities (C-species dominance) at all experimental plots. Only the control plot (subjected only to mechanical cleaning; no plant www.nature.com/scientificreports/ seedlings and biopreparation) exhibited a close position to the reference territory in the CSR term. However, this plot displayed the lowest vegetation restoration rate; its projective cover of vegetation in 2014 was 2 to 2.5 times lower than those at the experimental plots. Based on the theoretical principles of CSR theory, it can be assumed that the communities of the experimental plots will be relatively stable. The presence of a large number of competitor species (primarily cereals) hampers the growth of typical tundra and mire plant species at the remediated plots and prevents the restoration of background ecosystems. Employing only mechanical cleaning methods allowed the communities to recover a vegetation similar to the background vegetation; however, the recovery rate of this community was minimal.
In future studies, we plan to investigate whether the observed changes are stable in the long-term or whether the vegetation will return to its original state in a few years. Furthermore, it will be imperative to evaluate the differences in the rate of vegetation restoration of the treated and control plots.