Centennial recovery of recent human-disturbed forests

International commitments to restore degraded forests require global assessments of recovery timescales and trajectories of different forest attributes to inform restoration strategies. We use a meta-chronosequence approach including 125 forest chronosequences to reconstruct the past (*c. *300 years) and model future recovery trajectories of forests recovering from agriculture and logging impacts. We found recovering forests significantly differed from undisturbed ones after at least 150 years for ecosystem attributes like nitrogen stocks or species similarity and projected that difference to remain for up to 218 (38-745) or 494 (92-2,039) years, respectively. These conservative recovery metrics, however, still fail to capture the complexity of forest ecosystems, suggesting longer recovery timescales. Global restoration strategies have now the opportunity to engage in planning for a restored world that incorporates ecologically meaningful centennial implementation timescales and monitoring frameworks.

Methods
Database construction

We collected data from 16,873 plots from 125 chronosequences of forest ecosystems recovering for 50 to 295 years in 110 published primary studies. From these 125 chronosequences, we extracted 634 recovery trajectories of quantitative measures of ecosystem attributes along time, related to the six most widely included recovery metrics with enough representation to be statistically meaningful. These included biodiversity metrics (organism abundance, species diversity, and species similarity) and biogeochemical functioning metrics (carbon cycling, nitrogen stock, and phosphorus stock). We also extracted factors related to the context of where recovery and restoration happened and included the restoration strategy (passive and active), the disturbance type [agriculture (including land recovering from cultivation, grazing or combinations of both), logging and mining], the latitude, and the climatic condition (i.e., aridity index).

The trajectories related to organism abundance mainly contained biomass and density measurements. Species diversity included measurements of species richness and diversity indices. Species similarity trajectories contained information about species composition along the recovery trajectory, which were used to calculate pairwise compositional similarity at specific recovery times compared to a reference value. We used the Morisita-Horn similarity index, which accounts for species relative abundance. Abundance, diversity, and composition trajectories included five life forms: plants (including trajectories of woody plants, non-woody plants, and all plants combined), invertebrates, microorganisms other than fungi, fungi, and birds. Carbon cycling included pools and fluxes in soil, plants, litter, and microorganisms, whereas nitrogen and phosphorus stocks included bioavailable pools in soil, plants, and litter. From each trajectory, we collected all available recovery measures and compared them with a reference value.

Statistical analysis

To estimate the trajectory of forest recovery overtime, we fitted a separate linear mixed model (LMM) for the RR of each recovery metric. We included the recovery time as a fixed factor and as a random slope, and the trajectory identity as a random intercept, enabling a different slope and intercept for each trajectory (Eq. 1). As the recovery process over time may result in a wide range of trajectories from linear to more saturating shapes (65), we consider three functions to include the recovery time variable: one linear and two decelerating trends [ln(recovery time + 1) and √recovery time]. We then selected among the three options the one that best fit to the data of each recovery metric according to the minimum AICc. Their absolute values were square root transformed to meet the assumptions of general linear models and then multiplied by -1 to facilitate interpretation.

Using the resulting LMMs, we predicted the RR after 73, 146, and 219 years of recovery [i.e., one, two, and three times the global life expectancy in 2019, 73 years]. We then predicted the time needed for forest ecosystems to recover to 90% of reference values for each trajectory and recovery metric and calculated the median by metric.

Also using the resulting LMMs, we predicted the RR after 50 and 100 years of recovery for each metric and trajectory (1) to know if the recovery completeness is dependent on the metric and (2) to understand the main explanatory variables underlying the recovery process for each metric. Predictions were performed by using function predict() from stats package. We fitted linear models (LM) to analyse the difference in the RR after 50 years and after 100 years of recovery among recovery metrics.

We then fitted a separate LM for the effect of each explanatory variable studied (i.e., aridity, disturbance type or life form) on the RR predictions after 50 and 100 years of all recovery metrics together, and then for each recovery metric individually. In all the cases, the intercept of the LMMs for each trajectory was also included as a fixed factor to account for the effect of the initial state of degradation when recovery started. For the models fitted for the disturbance type and the life form, we excluded the categories with <1% of the values (i.e., “mining” for disturbance and “bird” for life form) or those including data with mixing information from other categories (i.e., “cultivation and grazing and logging” for disturbance and “woody and non-woody” for life form). We could not test the effects of the restoration strategy on the RR predictions as 88% of the recovery trajectories belong to passively restored forests.

To check the effect of the uncertainty in the recovery estimations, we first predicted the level of recovery after 50 and 100 years by using 999 random coefficients within the CI given by the models from the first stage and assuming a normal distribution. We then compared the original coefficients from the second-stage models with the resulting average coefficient calculated from 999 models coming from the randomizations. We did the comparison to check the effect of the recovery metric and of the recovery predictors (i.e. disturbance, aridity, and life form) after 50 and 100 years.