Assessing recovery of alpine spoil heaps by vascular plant, bryophyte, and lichen functional traits

Functional traits are linked to ecosystem processes and services and therefore relevant in recovery assessment. However, traits of bryophytes and lichens, important components ofmany ecosystems, have received less attention than those of vascular plants.We explored the use of functional traits of multiple important organism groups in recovery assessment. We combined data on traits and species composition for vascular plants, bryophytes, and lichens from four alpine spoil heaps and their undisturbed surroundings in western Norway, collected at three time-points spanning more than two decades. We studied changes in community-weighted mean (CWM) trait values and distribution of trait-category optima over time. We analyzed temporal variation in joint functional trait composition using the ordination regression-based approach (ORBA) to predict time to recovery. We observed functional shifts along the successional gradient for all organism groups, e.g. from wind-dispersed propagules shortly after disturbance to vegetative reproduction at later successional stages. Over time, the similarity between dispersalrelated traits of vascular plants and bryophytes on the spoil heaps and in their surroundings increased, indicating that propagule influx is important in alpine restoration. The joint functional trait composition of all spoil heaps converged towards that of their surroundings: one spoil heap had recovered 34 years after construction, while the predicted time to recovery for the other three was 59–74 years. Our results indicate that inclusion of multiple organism groups improves trait-based recovery assessments and time-to-recovery predictions. Further development of trait databases is essential for future use of joint functional trait composition in recovery assessment.


Introduction
Ecosystem degradation caused by human population growth and infrastructure development poses a major threat to future generations, and its mitigation is therefore a global priority (Díaz et al. 2019). Active use of ecological restoration should play a central role in reducing negative anthropogenic impacts on the environment (Palmer et al. 2004). Alpine ecosystems are particularly vulnerable to degradation because the short growing season and slow biological processes prolong the recovery period (Hagen & Evju 2013;Rydgren et al. 2013). Development of hydropower, an energy source with substantial potential in mountainous areas (Kumar et al. 2011), causes major disturbance in alpine ecosystems. Surplus rock material from hydropower development is often deposited on site in spoil heaps, which have a considerable ecological and visual impact (Skjerdal & Odland 1995;Rydgren et al. 2011). Spoil heaps have attracted research interest as model ecosystems for alpine restoration (Auestad et al. 2018), and long-term data suitable for assessing progress towards recovery are therefore available (Rydgren et al. 2020). The slow recovery of alpine ecosystems means that the consequences of restoration failure are long lasting, and recovery should therefore be monitored regularly by comparing measurable properties of degraded sites with restoration targets (Nilsson et al. 2016).
Multiple ecosystem properties have been proposed as suitable for use in recovery assessments (Rydgren et al. 2020). They include both univariate properties, e.g. species richness (Ruiz-Jaen & Aide 2005), and multivariate properties, e.g. species composition (Shackelford et al. 2013). Functional traits, i.e. morphological, physiological, or phenological features of individuals that influence their fitness through their effects on growth, reproduction, and survival (Violle et al. 2007), have recently been recognized as useful indicators of ecosystem functioning (de Bello et al. 2010;Funk et al. 2017). Functional traits are directly associated with key ecosystem processes, e.g. primary production (Lavorel & Garnier 2002;Funk et al. 2017), and ecosystem services, e.g. carbon sequestration (de Bello et al. 2010). There is therefore growing interest in functional trait targets in the context of ecological restoration (Laughlin 2014). Functional trait composition is a measurable multivariate ecosystem property that is particularly suitable for recovery assessment because of its direct links with changes in species composition and ecosystem functioning during the restoration (Zirbel et al. 2017).
Over the past two decades, the development of trait databases has facilitated analysis of multiple traits, which is an essential basis for assessment of functional trait composition. Evaluation of restoration success by means of functional trait composition has been explored in several studies (e.g. Pywell et al. 2003;Hedberg et al. 2013;Engst et al. 2016), though generally only vascular plants have been considered. Non-vascular cryptogams (i.e. bryophytes and lichens, cf. Cornelissen et al. 2007) have received considerably less attention in trait research (but see Lang et al. 2009;Hedberg et al. 2013;Roos et al. 2019), despite their important functional roles as essential components of many ecosystems. Bryophytes and lichens constitute a substantial fraction of the biomass at high latitudes and elevations, have complex interactions with vascular plants and other biota, and influence soil temperature, hydrology, and biogeochemistry, e. g. through nitrogen input from symbiotic bacteria (Cornelissen et al. 2007;Asplund & Wardle 2017).
To the best of our knowledge, no previous study has assessed ecosystem recovery through changes in joint functional trait composition of vascular plants, bryophytes, and lichens. All of these organism groups are relevant in a comprehensive assessment of the recovery of alpine ecosystems (Rydgren et al. 2011). However, joint analysis of the functional traits of multiple organism groups is methodologically challenging due to group-specific sets of traits and ecological and morphological differences between organisms (Hedberg et al. 2013). This challenge can be addressed through ordination methods, which are tools for the extraction of gradient structure (Økland 1990;van Son & Halvorsen 2014) applicable to any kind of compositional data, including data on functional trait composition (e.g. Fukami et al. 2005;van Son et al. 2013). If a clear successional gradient can be identified in the compositional data, the distances between restored and reference sites in the ordination space can be analyzed using the novel ordination regression-based approach (ORBA; Rydgren et al. 2019). Predictions of time to recovery provided by ORBA can be used to assess the restoration process (Rydgren et al. 2019).
In this study, we combined data on species composition for vascular plants, bryophytes, and lichens collected during the restoration of four alpine spoil heaps with database and literature data on functional traits to explore changes in joint functional trait composition. The main aim of our study was to explore the potential of using the joint functional trait composition of multiple organism groups in assessing ecosystem recovery.
We sought to answer the following questions: (1) How does the functional composition of vascular plants, bryophytes, and lichens differ between spoil heaps and their surroundings, and how does it change over time? (2) Are there any parallels in temporal development in related traits across the three organism groups? (3) Is ORBA useful for analyzing data on the joint functional trait composition of multiple organism groups for use in assessments of ecosystem recovery?
By answering these questions, we aim to advance the understanding of functional changes occurring during restoration of ecosystems with multiple important organism groups, and thus contribute to the improvement of functional trait-based recovery assessments.

Study Area
The study was carried out at four alpine spoil-heap sites in western Norway (Figs. S1 & S2, Table S1). The spoil heaps, each covering 2.2-4.1 ha, were established between 1974 and 1984. The bedrock in the study area comprises a mixture of granite, gneiss, phyllite, and other metamorphic rocks (NGU 2018), and all spoil heaps were made up of local rock material. Three spoil heaps consisted of blasted rocks and one (Kleådalen) was made up of fine-grained material from full-profile drilling (Skjerdal & Odland 1995). Spoil heaps differed from their surroundings in soil organic matter content and pH; soil organic matter content was most similar to the surroundings in the Kleådalen spoil heap, and pH was most similar in the Fossane spoil heap (Rydgren et al. 2020). Annual precipitation in the study area in the period 1971-2000 ranged between approximately 1,500 and 1,900 mm (NVE 2018).
All spoil heaps were seeded with a commercial grass mixture and fertilized with artificial fertilizer, without prior addition of topsoil (Skjerdal & Odland 1995 (Rydgren et al. 2011). All the study sites are in alpine areas used as summer pasture for free-ranging domestic sheep. No active measures (e.g. fencing) have been taken to exclude sheep from the spoil heaps, and there has been lowintensity grazing at all study sites.

Species Composition Dataset
Species composition was sampled at all four study sites at three time-points: 1990s (1991for Fossane, Kleådalen, and Svartavatn, 1994 for Øydalen), 2008 and 2015. In the 1990s, only spoil heaps were sampled, while in 2008 and 2015, we sampled both spoil heaps and their undisturbed surroundings; these are classified as two statistical "treatments" in this study. In the 1990s, stratified random sampling along baselines was used instead of sampling in blocks (Skjerdal & Odland 1995). In 2008, we subjectively placed eight blocks (each 5 × 10 m) at each site (Rydgren et al. 2011): five blocks on each spoil heap to cover variation in vegetation and ecological conditions, and three in the surroundings to provide a reference for recovery. Blocks in the surroundings were placed in well-drained locations where hydrological conditions were similar to those on the generally dry spoil heaps (Rydgren et al. 2011). In the 1990s, 10-20 non-permanent plots were placed on each spoil heap (total n = 64). In 2008, we placed three plots randomly within each block on the spoil heaps and in the surroundings (total n = 96), and permanently marked them with metal tubes. Permanent plots were resampled in 2015, except for three plots in one block on the Kleådalen spoil heap, which were buried under newly deposited rock material between 2008 and 2015. All plots measured 0.5 × 0.5 m and were divided into 16 equal-sized subplots. At each time-point, all species (or species aggregates for genera in which species could not be distinguished, see Supplement S1) of vascular plants (n = 107), bryophytes (n = 89), and lichens (n = 48) in each subplot were recorded. Presence in subplots (0-16) represented species abundance. The species composition dataset consisted of 253 plot-by-time combinations, which formed 20 "plot groups," i.e. combinations of 4 sites × 3 timepoints × 2 "treatments" (for the time-point 1990s, only 1 "treatment," i.e. spoil heap. was sampled).

Trait Dataset
We selected traits for the trait dataset using the following criteria: (1) relevance for growth, reproduction, and survival; (2) relevance for successional and ecosystem processes (Violle et al. 2007); and (3) data availability. We used trait data from databases and other literature sources deriving primarily from Norway, Northern Europe and European alpine regions, i.e. the Alps and the Carpathians. As most traits were categorical and we analyzed community-level averages of trait values, the potential bias due to intraspecific trait variation was considered negligible (cf. Auger & Shipley 2013). Missing species-by-trait combinations (n = 35) were replaced by the median value for the given trait for species from the same genus and alpine vegetation zone present in the trait dataset (n = 15) or present in a database/ other literature source (n = 19). For one missing species-by-trait combination without a suitable replacement (n = 1), we used the median value for the trait for all species present in the trait dataset (cf. van Son et al. 2013). We aggregated several originally distinct traits/trait categories into compound traits (for details about trait data, see Table 1).
Vascular plant traits included maximum height, seed mass, specific leaf area, main dispersal mode, lateral spread, and growth form. At our study sites, tall shrubs and trees do not reach the potential maximum heights reported in the databases. We therefore used 1 m as an approximation for the local maximum height of all tall shrub species (Juniperus communis, Salix phyllicifolia, and S. pentandra), and 2 m for tree species (Betula pubescens).
Bryophyte traits included maximum size, mean spore size (i.e. the mean of maximum and minimum spore size values), frequency of sporophytes, frequency of vegetative propagules, sexuality, and growth form.
Lichen traits included maximum thallus size, main reproduction type, photobiont association and growth form. No trait database was available for lichens (Branquinho et al. 2015), and we used trait data from multiple sources (Table 1). Maximum thallus size was included because of associations with specific thallus mass, water-holding capacity, and production of secondary compounds (Asplund & Wardle 2017) and because similar (but not fully analogous) traits were used for vascular plants and bryophytes. The other three traits were included because of their relevance for growth, survival, and reproduction (Ellis & Coppins 2006;Koch et al. 2013;Nelson et al. 2015b). The trait dataset consisted of 1,368 species-by-trait combinations.

Statistical Analyses of Changes in Joint Functional Trait Composition
The joint functional trait composition dataset was prepared in two steps. In the first step, we transformed continuous trait data into categorical data to ensure comparability with a priori categorical traits. Values for each continuous trait for all species were sorted in ascending order, and divided into five categories based on sample quintiles (cf. van Son et al. 2013): small, medium-small, medium, medium-large, and large. Each species was subsequently assigned to one category per trait. A priori categorical traits were left untransformed. In the second step, we obtained the joint functional trait composition of each plot-by-time combination by calculating trait-category abundances as the weighted averages of trait categories with species abundances as weights (for details, see Supplement S2). The joint functional trait composition dataset consisted of 70 trait categories × 253 plot-by-time combinations.
Recovery assessment using the ordination regression-based approach (ORBA; Rydgren et al. 2019) required a successional gradient quantified by ordination methods. To quantify the successional gradient, we subjected the joint functional trait composition dataset to multiple parallel ordination (MPO; van Son & Halvorsen 2014), using detrended correspondence analysis (DCA; Hill & Gauch 1980) and global (GNMDS; Minchin 1987) and local (LNMDS; Sibson 1972) non-metric Restoration Ecology multidimensional scaling (NMDS). The MPO procedure, in which more than one ordination method is applied to the same dataset, facilitated artifact detection and established a consolidated gradient structure (van Son & Halvorsen 2014). We conducted all ordinations using package "vegan" (Oksanen et al. 2017). For details of ordination methods, see Supplement S3. Based on the selection process described in Supplement S3, we chose the two-dimensional LNMDS for interpretation and further analyses.
To confirm that the LNMDS ordination adequately quantified the successional gradient, we used linear mixed-effects models (LMM) as implemented in package "lme4" (Bates et al. 2015). Separate models were obtained for LNMDS axes 1 and 2, using plot scores as response variables. As the fixed effect, we used time-point × "treatment" combinations, i.e. a factor variable with five levels: 1990s spoil heap (1990s surroundings not sampled), 2008 spoil heap, 2008 surroundings, 2015 spoil heap, and 2015 surroundings. To address repeated sampling, i.e. using non-permanent plots in the 1990s and permanent plots in 2008 and 2015, we randomly assigned permanent plot IDs to the non-permanent plots from the same site in the models. Repeated surveys of plots and the hierarchical sampling design (plots nested in blocks nested in study sites) were accounted for by including appropriate random effects in the models. We assessed differences between fixed-effect factor levels by multiple comparisons between all pairs of means using Tukey contrasts, as implemented in the function "glht" of package "multcomp" (Hothorn et al. 2008).
In the multiple comparisons, the LNMDS-axis 1 scores differed significantly between all fixed-effect factor levels, except Table 1. Overview of used traits and data sources. Abbreviations: (trait) C, compound trait (i.e. several originally distinct traits/trait categories from data sources were aggregated), F, fuzzy coding (i.e. assignment of species/species aggregates to multiple trait categories; see Supplement S2 for details); (relevance) G, growth, R, reproduction, S, survival; description (scale) cat., categorical, cont., continuous. Data sources (

Trait
Relevance Description (scale) Units or Categories Data Sources Vascular plants Restoration Ecology for the comparison between the surroundings in 2008 and 2015 (Fig. S3A). Plot-score means increased along LNMDS-axis 1 from spoil heap plots sampled in the 1990s, via 2008 and 2015 to plots in the surroundings (Fig. S3A). We found no significant differences in plot scores along LNMDS-axis 2 between the spoil heaps and the surroundings (Fig. S3B). We therefore concluded that plot scores along LNMDS-axis 1 adequately represent the successional gradient, which is required for recovery assessment using ORBA. We used the distance along LNMDS-axis 1 between each spoil-heap plot-by-time combination and the centroid of plots from the corresponding surroundings, referred to as successional distance (Rydgren et al. 2019), to assess the recovery status of a given spoil-heap plot at a given time-point. As the surroundings were not sampled in the 1990s, scores for the surroundings in 2008 were used to estimate successional distances in the 1990s (for justification, see Supplement S4). We defined recovery as reached when the successional distance between the centroid of spoil-heap plots and the centroid of plots from the corresponding surroundings (at a given site and time-point) was lower than or equal to a threshold. As the threshold, we used the mean absolute deviation (MAD), that is, the mean distance of the plots from the surroundings to their centroid in 2015.
We used ORBA with a dynamic reference (Rydgren et al. 2019) to model the recovery process (see Fig. S4). Time to recovery was modeled separately for each spoil-heap site using an asymptotic LMM. We used log-transformed successional distances as the response and spoil-heap age as the fixed effect. Repeated and spatially nested sampling was accounted for by specifying plot nested within block within site as random effects (Rydgren et al. 2020). We calculated 95% confidence intervals (95% CI) around model predictions as ±1.96 × SE of fitted values.

Variation in Single Traits
To explore variation in continuous traits (Table 1) over time and between spoil heaps and their surroundings, we first calculated the community-weighted mean (CWM; Garnier et al. 2004) for each trait in each plot-by-time combination. Continuous trait values from the trait dataset were centered and standardized by subtracting the mean and dividing by the SD before calculating the CWM, using species abundances as weights. We applied LMM to each trait and each site with plot-wise CWM values as the response, "plot group" affiliation (i.e. combination of site × time-point × "treatment") as fixed-effect factor, and plots nested in blocks as random effects. We subsequently assessed differences between "plot groups" within each site by multiple comparisons of means using Tukey contrasts (Hothorn et al. 2008).
Changes in a priori categorical traits (Table 1) along the successional gradient were studied by comparing positions of traitcategory optima. These optima were estimated as weighted averages (WA) of LNMDS-axis 1 scores using trait-category abundances as weights, as implemented in function "wascores" in package "vegan" (Oksanen et al. 2017). We divided the WA values, sorted in ascending order, by sample quintiles into five successional stages (early, early-intermediate, intermediate, intermediate-late, and late). This reflected the successional gradient from "young" spoil heaps to the vegetation of the surroundings. However, as the successional gradient is continuous, any division into discrete stages is arbitrary, serving illustrational purposes only.
All analyses were conducted in R version 3.2.0 (R Core Team 2015). For flowchart of datasets and analyses, see Figure S5.

Variation in Single Traits of Vascular Plants, Bryophytes, and Lichens
In general, trait community-weighted means (CWMs) of spoilheap plots became more similar to the surroundings over time (Fig. 1). However, CWMs of several continuous traits differed between spoil heaps, with respect to both change over time and difference from their surroundings.
The specific leaf area (SLA) of vascular plants decreased over time on the spoil heaps towards the level in the surroundings (Fig. 1A). Seed mass on the spoil heaps gradually converged towards that of the surroundings at all sites except Svartavatn, where, however, there was no significant difference between the spoil heap and the surroundings at any time-point. Seed mass increased at Kleådalen and Øydalen, and decreased at Fossane and Svartvatn spoil heaps (Fig. 1B). Bryophyte maximum size increased with time on all spoil heaps, becoming more similar to that in the surroundings. However, at Fossane and Øydalen, there were still significant differences between lichen and bryophyte maximum size on the spoil heaps and in their surroundings in 2015 (Fig. 1C, E). Mean bryophyte spore size on the spoil heaps became gradually more similar to that in the surroundings, but the convergence pattern was weaker than for vascular plant seed mass, and differed between sites. In 2015, there was still a significant difference in bryophyte spore size between the spoil heaps and their surroundings at Fossane and Kleådalen (Fig. 1D).
The distribution of optima of a priori categorical traits along the successional gradient revealed functional shifts ( Fig. 2; for details, see Table S2). The early-successional stage was associated with epizoochorous dispersal of vascular plants and abundant production of sporophytes of bryophytes, while anemochorous dispersal and frequent production of sporophytes prevailed in the early-intermediate stage. Hemicryptophytes (forbs, e.g. Alchemilla alpina) and turfforming bryophytes (e.g. Polytrichum alpinum) were common at the early-successional stage, while foliose lichens (e.g. Peltigera spp.) characterized the early-intermediate stage.

Recovery Assessment Based on Joint Functional Trait Composition
Although the spoil-heap plot scores along LNMDS-axis 1 gradually became more similar to the scores of plots in the surroundings (Fig. 3), they were still significantly different in 2015 (Fig. S3A). Only the Kleådalen spoil heap had reached the recovery threshold in 2015 (Fig. 4B). According to ORBA estimates, the threshold was reached in 2010, 29 years after establishment (95% CI: 24-34 years). For the other three spoil heaps (Fig. 4A,C,D), the predicted time from establishment to recovery varied from 74 years (95% CI: 69-80 years) at Fossane and 65 years (95% CI: 60-72) at Øydalen to 59 years (95% CI: 52-69) at Svartavatn. The uncertainty of predictions expressed as the 95% CI varied from 10 years for Kleådalen to 17 years for Svartavatn (Fig. 4, Table S3).   Table 1) on the spoil heaps and in the surroundings at three time-points, with site-wise multiple comparisons of means of "plot groups" with Tukey contrasts (i.e. all-pair comparisons). "Plot groups" within each site with no common letter are significantly different at the α = 0.05 level. Note that the Y-axes in (A-E) differ in scale, and that the different study sites (names indicated at the top) are separated by vertical lines. X-axis labels: "Plot group," i.e. time-point and "treatment" -Spoil heap (H; gray boxes) or surroundings (S; green boxes). No lichen species were recorded in the 1990s at Fossane and Svartavatn  -point (vertical -1990s, diagonal -2008, and horizontal -2015). "Plot-group" centroids are indicated by "+" for spoil-heap "plot groups" and "×" for those of the surroundings.

Variation in Single Traits of Vascular Plants, Bryophytes, and Lichens
Patterns of functional traits of vascular plants observed in this study generally agreed with those reported in other studies of succession (Garnier et al. 2016). There was a gradual decrease in community-weighted mean (CWM) values of specific leaf area (SLA) on the spoil heaps, reflecting a shift from rapid acquisition to increased conservation of resources (Navas et al. 2010). At all sites, the CWM values for seed mass on the spoil heaps either became more similar to those of their surroundings, or were not significantly different from the surroundings to start with (Svartavatn). In the former case, the increasing similarity resulted from a combination of increasing (Kleådalen, Øydalen) and decreasing (Fossane) spoil-heap CWM values over time. A plausible explanation for this paradoxical result is that recovery of the vascular plant community on the spoil heaps is primarily mediated by local seed influx from the surroundings, possibly with dispersal mechanisms varying between sites. The importance of anemochory and epizoochory in early-successional stages highlights the role of wind and animals (e.g. domestic sheep) in local seed dispersal. Plant establishment on the spoil heaps may thus be limited by the availability of suitable microsites rather than by seed availability (Rydgren et al. 2011), as might be expected in an open, wind-exposed alpine landscape (Erschbamer et al. 2008). The sequence of optima for vascular plant growth forms along the successional gradient, and the increasing abundance of tree and (dwarf) shrub species that are dominant in the surroundings, correspond to the general successional pattern (Garnier et al. 2016). The significant differences between bryophyte and lichen maximum size on the spoil heaps and in the surroundings at Fossane and Øydalen suggest that the bryophyte and lichen communities on these spoil heaps are not yet functionally restored. The abundant sporophyte production and smaller spores observed early in early-successional stages are consistent with the inverse relationship between spore mass and dispersal potential reported by During (1992). However, the pattern of convergence is weaker for spore mass than for seed mass, indicating that sexual reproduction plays a smaller role in bryophytes than in vascular plants (During 1992). The difference between mean spore size on the Kleådalen spoil heap, which is close to functional recovery, and in its surroundings, may reflect the increasing importance of vegetative reproduction as succession proceeds. This may result from a trade-off between the production of sporophytes and of vegetative propagules (Austrheim et al. 2005b). The sequence of bryophyte growth form optima along the successional gradient, ending in wefts, that is, layers of intertwining branched shoots, fits with patterns typical of the establishment of closed vegetation during primary succession (Gimingham & Birse 1957).
Like Nelson et al. (2015a), we found that lichen species reproducing by soredia prevailed in the early-successional stage ("young" spoil heaps in our case). Soredia are highly mobile vegetative propagules particularly suitable for rapid colonization of new areas (Nelson et al. 2015a) because they disperse the mycobiont and photobiont simultaneously. In lichens, only the mycobiont reproduces sexually and must subsequently find a new photobiont partner. This may explain why sexual reproduction of vascular plants and bryophytes prevailed at an earlier successional stage than sexual reproduction of lichens. Thallus fragmentation, on the other hand, became more important towards later successional stages, as did vegetative reproduction in vascular plants and bryophytes. The sequence of optima we observed for lichen growth forms, i.e. foliose-crustose-fruticose, deviated from the sequence typical of succession, i.e. crustose-foliosefruticose (During 1992). This supports the view of During (1992) that the latter sequence does not occur in all successions. However, the aggregation of lichen growth forms into three main categories may also have influenced the result, as the "crustose" category contained squamulose lichens only.

Recovery Assessment Based on Joint Functional Trait Composition
The joint functional trait composition of vascular plants, bryophytes, and lichens on all spoil heaps converged towards that of their surroundings, although at different speeds at different sites. Similar patterns were found for the temporal dynamics of total cover, species richness, and physical environment on the spoil heaps (Rydgren et al. 2020). The fact that the ordination regression-based approach (ORBA) provides ecologically plausible results for joint functional trait compositional data suggests that it is an appropriate analytical method for such data.
Only one of the spoil heaps studied, Kleådalen, had reached the recovery threshold for functional trait composition by 2015, the last time-point. A likely explanation for the rapid recovery at Kleådalen is that local soil characteristics, especially the more fine-grained topsoil, facilitated plant establishment and accumulation of organic matter (Rydgren et al. 2013;Rydgren et al. 2020). In addition, Kleådalen is the study site at the lowest altitude, and the longer growing season may encourage the establishment of vegetation. Joint functional trait composition at Kleådalen shows an interesting late-successional pattern, with successional distances apparently stabilizing around the recovery threshold. However, more data will be needed to follow post-recovery developments, and to identify the causes of any deviations from the current ORBA models, e.g. persistent differences in environmental conditions between the spoil heaps and their surroundings. Nevertheless, the decrease over time in successional rates observed at all sites shows that this general property of succession (Rydgren et al. 2019) also applies to successional gradients in functional trait composition.

Joint Functional Trait Composition: Relationship With Species Composition and Use in Recovery Assessment
We observed faster recovery in functional trait composition than in species composition at the same study sites (Rydgren et al. 2020), in accordance with Engst et al. (2016). Furthermore, ordering the sites by time to recovery gives different results for the two metrics.
For example, the Svartavatn spoil heap was found to show the second fastest recovery of functional trait composition, but the slowest recovery of species composition (Rydgren et al. 2020).
Thus, a comparison between our results for functional traits composition and those of Rydgren et al. (2020) for species composition reveals that the two metrics provide somewhat different signals about ecosystem recovery at the study sites. This is interesting, as the two metrics are closely related. For example, based on the mass ratio hypothesis (Grime 1998;Garnier et al. 2004), species abundances are used in calculating functional trait CWMs, which represent functional trait composition. There are examples indicating that recovery assessment based on functional trait composition alone may be misleading, e.g. if invasive alien species are present in the restored area ). However, recovery assessment based on species composition alone may fail to reveal situations where a restored ecosystem has achieved functional similarity with the reference, but with a different species composition (Cadotte et al. 2011;Engst et al. 2016). We therefore recommend the use of functional trait composition and species composition in parallel, and suggest that the two approaches should be considered as complementary (cf. Engst et al. 2016). In alpine and other ecosystems where bryophytes and lichens are fundamental elements (Cornelissen et al. 2007), we recommend using joint functional trait composition rather than vascular plant functional trait composition alone for recovery assessments.
Time-to-recovery predictions by ORBA for our study sites showed similar degrees of uncertainty, expressed as 95% CI of model estimates, for functional trait composition (this study) and species composition (Rydgren et al. 2020). Laughlin et al. (2017) found functional trait-based metrics to be less variable and more predictable than species composition-based metrics for a variety of restoration treatments, while Abella et al. (2018) questioned the existence of a general hierarchy of predictability. Thus, the two types of metrics should be considered as complementary for predictions of time to recovery as well as for recovery assessments. Predictions based on both functional trait composition and species composition provide a more reliable basis for making decisions on restoration than predictions based on either metric alone. If there is substantial disagreement between recovery assessments based on the two metrics, the results should be discussed in the context of specific restoration goals, such as ecosystem functioning or conservation of particular species or species groups.

Outlook for the Functional Trait-Based Approach to Recovery Assessment
It is only possible to use joint functional trait composition in recovery assessments if reliable, relevant trait data can be obtained from trait databases. While the number of vascular plant species and traits covered by available databases has increased strongly in recent years, information about intraspecific (including regional) variation in traits remains sparse (Moran et al. 2016). This variation may be substantial (de Bello et al. 2011;Asplund & Wardle 2017). However, interspecific trait variation generally dominates in the presence of a strong gradient (Auger & Shipley 2013), such as a successional Restoration Ecology gradient (Rydgren et al. 2019). The issue of intraspecific trait variation is also less relevant for categorical traits than for continuous traits. Judicious use of trait-database data is therefore justified in recovery assessments of ecosystems where the interspecific trait variation is likely to dominate over the intraspecific (Hedberg et al. 2013), and when the assessed traits are mostly categorical. However, further development of trait databases for bryophytes and lichens, e.g. LIAS light (Rambold et al. 2014), is essential to enable the use of joint functional trait composition in future recovery assessments.
As far as we know, this study is the first to use joint functional trait composition of multiple organism groups for assessment of ecosystem recovery. We have demonstrated that ORBA (Rydgren et al. 2019) can be used to analyze functional traits of multiple organism groups simultaneously, solving the methodological challenge of combining traits of morphologically different organisms in one analysis (cf. Hedberg et al. 2013). ORBA thus has the potential to become a powerful tool for assessing recovery on the basis of functional traits, and we recommend that more experience of this approach should be gained by testing it on datasets from other ecosystems.
Analysis of functional traits of multiple functionally important groups will make it possible to develop two approaches needed in restoration ecology in the 21st century: using multiple variables for monitoring and assessing restoration projects (Ruiz-Jaen & Aide 2005; Shackelford et al. 2013), and focusing on ecosystem function in restoration (Choi 2007;Shackelford et al. 2013). Our results indicate that using appropriate quantitative methods to analyze joint functional trait composition of multiple organism groups may be an important step towards the goal of making restoration ecology a predictive science (Brudvig 2017).

Supporting Information
The following information may be found in the online version of this article: Table S1. Overview of studied spoil heaps. Table S2. Optima for categorical traits as WA-scores at LNMDS-axis 1. Table S3. Predicted time to recovery (model estimate and 95% CI) in years since establishment (i.e. spoil heap age). Figure S1. Map of the study area and locations of the four study sites in western Norway Figure S2. The four studied spoil heap sites. Figure S3. Multiple comparisons of means of LNMDS-axis scores Figure S4. Illustration of the data input for regression step in ordination regressionbased approach (ORBA) Figure S5. Flowchart of datasets and analyses Supplement S1. Species aggregates Supplement S2. Preparation of the joint functional trait composition dataset Supplement S3. Ordinationsjustification and settings Supplement S4. Successional distances in the 1990sjustification