1. Human-induced disturbance has substantially influenced the structure and function of terrestrial ecosystems globally. However, the extent to which multiple ecosystem functions (multifunctionality) recover following anthropogenic disturbance (ecosystem recovery) remains poorly understood.

2. We report on the first study examining the temporal dynamics in recovery of multifunctionality from 3 to 12 years after the commencement of rehabilitation following mining-induced disturbance, and related this information to changes in biota. We examined changes in 57 biotic (plants, microbial) and functional (soil) attributes associated with biodiversity and ecosystem services at four open-cut coal mines in eastern Australia. 

3. Increasing time since commencement of rehabilitation was associated with increases in overall multifunctionality, soil microbial abundance, plant productivity, plant structure and soil stability, but not nutrient cycling, soil carbon sequestration nor soil nutrients. However, the temporal responses of individual ecosystem properties varied widely, from strongly positive (e.g., litter cover, fine and coarse frass, seed biomass, microbial and fungal biomass) to strongly negative (groundstorey foliage cover). We also show that sites with more developed biota tended to have greater ecosystem multifunctionality. Moreover, recovery of plant litter was closely associated with recovery of most microbial components, soil integrity and soil respiration. Overall, however, rehabilitated sites still differed from reference ecosystems a decade after commencement of rehabilitation.

4. Synthesis and applications. The dominant role of plant and soil biota and litter cover in relation to functions associated with soil respiration, microbial function, soil integrity and C and N pools suggests that recovering biodiversity is a critically important priority in rehabilitation programs. Nonetheless, the slow recovery of most functions after a decade indicates that rehabilitation after open-cut mining is likely to protracted.

Plant and soil sampling

Within a 20 m by 20 m Biodiversity Assessment Method (BAM; DPIE, 2020) floristics plot we recorded the foliage cover of all native and exotic vascular plant species, the length of logs (> 0.1 m diameter), the cover of litter (detached leaves, sticks, bark) and characteristics of the soil surface. Within each BAM plot we also established eight 25 x 25 cm quadrats; four beneath trees and four in the open. All litter was collected from each quadrat and pooled to obtain one tree and one open sample per plot. Litter was air dried and sorted in the laboratory into the following categories: sticks and bark, tree and shrub leaves, grass and forb leaves, coarse frass (> 4 mm diameter and < 40 mm long), fine frass (litter fragments < 4 mm but > 2 mm diameter), reproductive material (seeds/seedpods, fruits and flowers), woody material (woodchips and mulch), and vertebrate dung. Following litter collection, a single soil sample (0-5 cm depth) was collected from the centre of each of the eight quadrats and pooled by quadrat type (beneath trees or open) to obtain one composite sample of each type. For each composite sample, a small subsample was immediately separated and frozen (-18^oC) and later freeze-dried for DNA- and phospholipid- based microbial analyses. 

We used a field-based protocol (Landscape Function Analysis, hereafter “LFA” Tongway 1995) to assess the characteristics of the soil surface within the BAM plot. Within the plot we measured 12 surface attributes: surface roughness, crust resistance, crust brokenness, surface stability, surface integrity (cover of uneroded surface), cover of deposited materials, biocrust cover, plant foliage cover, plant basal cover, litter cover, litter origin, and the degree of litter incorporation (Appendix S2). These 12 surface features provide a measure of the health of the soil surface. Indices derived from these measures have been shown to be highly correlated with ecosystem functions related to soil stability, nutrient cycling and infiltration (Maestre & Puche, 2009; Eldridge et al., 2019; Appendix S2). 

Soil chemical and microbial analyses

Extractable phosphorus (Colwell-P) and pH (1:5 soil water extract) were determined following standard procedures (Rayment & Lyons, 2011). Total nitrogen and soil carbon fractions (total organic carbon, particulate organic carbon, humic organic carbon, resistant organic carbon) were estimated using mid-infrared (MIR) spectroscopy techniques (Baldock et al., 2013; Baldock et al., 2014; Appendix S3).

DNA was extracted from 0.25 g of freeze-dried soil using the DNeasy PowerSoil Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Amplicons targeting the bacterial 16S rRNA gene (341F-805R, Herlemann et al. 2011) and the fungal ITS region (FITS7-ITS4R, Ihrmark et al., 2012) were sequenced at the Ramaciotti Centre for Genomics, University of New South Wales (Sydney, Australia), on the Illumina MiSeq platform (Appendix S4). Raw sequencing data can be accessed at the NCBI Sequence Read Archive under BioProject <accession to be added prior to publication>. DNA sequence processing methods, including assignment of DNA sequencing reads to operational taxonomic units (OTUs) are described in Appendix S4. The OTU abundance tables were rarefied to an even number of sequences per sample (6677 and 9073 sequences for bacteria and fungi, respectively; the minimum number of sequences observed). Alpha diversity metrics were then calculated using the ‘vegan’package (Oksanen et al., 2019) in R statistical software. 

Three methods were used to derive proxies of microbial activity: 1) potential soil enzyme activity, 2) soil respiration, and 3) phospholipid derived fatty acids (which is a surrogate of active microbial biomass). We measured the potential activity of six enzymes (Appendix S5, Table S5) as proxies for carbon, nitrogen and phosphorus degradation using fluorometry with, as described in Bell et al. (2013) with some modifications (Appendix S5). The MicroResp^TM technique was used to determine substate-induced respiration rates using seven substrates: fructose, glucose, maltose, raffinose, sucrose, threonine and xylose. Basal respiration was assessed with filtered sterile deionized water. Thirty microlitres each of filtered sterile substrate and deionized water were added to the pre-incubated soil in deep-well plates. Absorbance was measured at 570 nm immediately before and after the 24 h incubation period using a microplate reader (SpectraMax M2^e) (Campbell et al., 2003; Burton et al., 2007). The rate of CO₂ respiration per gram of dry soil was calculated according to the formula as described in MicroResp™ manual (MicroResp™, James Hutton Ltd, Aberdeen, UK). Microbial phospholipid fatty acids (PLFA) were assessed as indicators of total microbial biomass in soil and the relative abundance of general microbial functional groups such as bacteria (gram negative, gram positive), Actinobacteria, fungi, mycorrhizae, as well as of protists, which are soil microfauna (Appendix S6).

Statistical analyses

Here we report on the recovery of 57 attributes described above for 90 sites (Appendix S7). Of the 90 sites, 42 ranged in time since the commencement of rehabilitation treatments from 3 to 12 years (hereafter ‘Rehabilitated’), and 48 were largely intact communities that may have been altered slightly but represent highly functional (hereafter ‘Reference’) sites that showed various levels of disturbance or change typical of the remaining remnants for these vegetation communities (Table S1.2). Our attributes were classified as either biotic structural attributes (n=22), i.e., they represent biological components such as the plant cover or microbes, or functional (n =35), i.e., they contributed to resource and/or energy flows, e.g., soil respiration (Table S7.1). We also had information on pH, which we used to explore downstream relationships with microbial biomass and diversity. Prior to statistical analyses, we scaled all quadrat-level data to the plot level using the relative cover of trees at each site.

First we standardised (z-transformed) the value of each attribute, for each of the 42 rehabilitation sites, and arranged the 57 biotic/functional attributes into seven categories (soil carbon sequestration, decomposition, microbes, plant productivity, plant structure, soil stability, soil nutrients, Appendix S7). The multifunctionality values for the seven categories were calculated as the mean standardised value for those attributes within a given category. The multifunctionality approach allows us to compare attributes that might vary markedly in their mean and range by bringing the value of all attributes to a common scale with a mean = 0 and SD = 1. For example, the soil stability multifunctionality index was calculated as the mean standardised values of six attributes (soil brokenness, deposited materials, soil integrity, surface resistance, soil surface roughness, crust stability). This is possible because increasing values of each attribute represent increasing function. We then explored whether the multifunctionality index of the seven categories changed in relation to time since rehabilitation. 

Next, we allocated the seven categories to either a biotic or functional category and averaged the functional attributes (soil carbon, decomposition, nutrients, plant productivity, soil stability) to form an average value of ecosystem multifunctionality and did the same for the microbes and plant structure to create an average measure of biotic structure. This averaging is possible because the values are standardised and therefore bounded. The relationship between biotic structure and multifunctionality allowed us to explore whether rehabilitation of the average functional effect was associated with changes in the average biotic effect.

We then used a Relative Interaction Intensity (RII) index (Armas et al., 2004) to explore how the raw values of the 57 individual attributes changed with time since restoration after accounting for values at their relevant reference (control) site. The RII index is calculated as (X_R-X_C)/(X_R+X_C), where X is the value of the attribute, R = Rehabilitated and C = Control (Reference). This relativisation process results in an index that is bounded by −1 and 1, with RII values > 0 indicating relatively greater values in rehabilitated sites. The control (reference) sites used in this index were those from the same mine. The significance of any relationships between biota structure and multifunctionality, time since commencement of rehabilitation or RII and any individual attributes was tested using linear and non-linear (quadratic) regressions, with mine identity as a random effect. Akaike’s Information Criterion (AIC) was used to decide the model that provided the best fit in each case. We then visualised Pearson’s correlation, among the 22 relativised biotic attributes and the 35 relativised functional attributes in a heat map.

Relevant methodological considerations are given in the Supplementary material accompanying the manuscript.