Philippa Kaur

Experimental protocol
Twenty healthy male adults participated in the study. To limit the effects of differences in the participants’ anthropometry, the volunteers were matched according to their height and weight.
The leg and hand dominances of each participant was noted before running the experiment. The subjects were placed at random on the left ((#1—right hand holding the load) or the right (#2—left hand holding the load) side of the carried object. It happened by chance that three lefthanders were affected to #2. The individuals had an average height of (mean ± SD) 1.77 ± 0.07 m (#1), and 1.77 ± 0.05 m (#2), and an average weight of 74.78 ± 9.00 kg (#1), and 74.54 ± 12.38 kg (#2). The load was symmetrical in shape, and its weight was evenly balanced between the participants who were positioned randomly with respect to the object transported to counteract the effect of a dominant side (Fig. 1A).
The study was carried out with healthy individuals who wrote their informed consent to participate in the experiment and to be filmed and photographed. The experiment was non-interventional, and the movements performed by the volunteers were no more risky than those they perform in daily activities. The study was approved by the Research Ethics Committee of the University of Toulouse, France (Number IRB00011835-2019-11-26-172, Université Fédérale de Toulouse IRB #1).
The instructions given to the volunteers were: “Move the load together from point A to point B” and “Any communication between you is forbidden during the experiment”. Point A and point B laid 20 m apart. No explicit instruction was given as to how fast the volunteers should perform the task. The volunteers were tested with three conditions called CT20, CT30, and CT40 corresponding to a load representing on average 20%, 30%, and 40% of the sum of their body masses respectively. The three conditions were tested in random order for each pair. To avoid adaptations due to familiarization or learning, only one trial per pairs and condition was recorded.
Kinematics and kinetics
Thirteen MX3, and TS40 Vicon cameras (Vicon©, Oxford) were used to capture the positions of ninety-one retro-reflective markers taped on the system formed by the paired individuals and the load they carry (hereafter called Poly-Articulated Collective System—PACS): 42 markers on each individual31,32, and seven on the load (Fig. 1A). The acquisition frequency was set to 200 Hz. In order to record the walking patterns of the individuals at a stable speed, and thus to exclude the acceleration and deceleration phases at the beginning and end of each trial, the calibrated volume (30m3) corresponded to the central part of the walkway. This covered about two steps. Concerning kinematic analysis, the PACS was reconstructed with the Vicon Nexus™ 1.8.5 software. Reconstruction was impossible for one pair of individuals who had lost one reflective marker. The two lateral handles on each side of the load were equipped with a 6-axis force sensor (Sensix®, France) (Fig. 1B), allowing to record the reaction forces and moments at a sampling frequency of 2000 Hz. The kinematic and kinetic measurement errors were 1 mm for 1 m for the positions (Vicon system) and ± 0.01 N for the forces (Sensix sensors), respectively. The sensors frames were located with the help of screwed reflective markers. The data were filtered with 4th order Butterworth filters with a cut-off frequency of 5 Hz for the kinematic data, and of 10 Hz for the kinetic data. To ensure at least one complete walking cycle for each subject of a pair, the gait cycle of the PACS was defined from the first heel strike of individual #1 to the third heel strike of individual #2.
COM determination and related parameters
The carried object, which constituted the 33th segment of the PACS, was built in aluminum and was therefore extremely rigid. It was completely symmetrical about its sagittal plane (Fig. 1B) and therefore its weight was evenly balanced between the participants. The object was equipped with a rod at its center where standard cast iron discs could be slid to increase its weight. The CoM of the object was determined at the intersection point of the vertical lines obtained by hanging the object without discs with a thread fixed at different positions. When the object was loaded, the position of its CoM was then adjusted by taking into account the added cast iron discs and by considering a homogeneous mass distribution inside the discs.
The De Leva Anthropometric table33 allowed us to estimate the mass mi as well as the CoM of each segment i (CoMi) of the PACS, and thus to compute its global CoM (CoMPACS) as follows:

$${varvec{G}}{text{PACS}} = frac{1}{{m_{PACS} }}mathop sum limits_{i = 1}^{n = 33} m_{i} {varvec{G}}_{i}$$
(1)

with GPACS corresponding to the 3D position of the CoMPACS in the frame R (the global coordinate system), mPACS to the mass of the PACS, n the number of PACS segments (i.e. 16 segments per volunteer, plus one segment for the box), and Gi corresponding to the 3D position of the CoMi in R.
The vertical amplitude (Az = Zmax − Zmin, in meters) of the CoMPACS trajectory along two consecutive steps, and the length of two consecutive steps by each individual were also computed.
Assessment of energetic exchanges
To assess energetic exchanges, forward kinetic work, vertical work, and external work of the forces applied to the CoMPACS were computed25.
Forward kinetic work (Wkf) was defined as the positive work needed to move the CoMPACS forward, and it was calculated as the sum of the increments of forward kinetic energy (Ekf) along the time curve:

$$E_{{{text{kf}}}} = { }frac{1}{2} m overrightarrow{V_{f}}left( {text{t}} right)^{{2}} _{/R}$$
(2)

with m being the mass of the individual, and (overrightarrow {{V_{f} }}) (t)/R the linear forward velocity of the CoMPACS in the frame R. The x-, y- and z-axis of the frame R, corresponding to the medio-lateral, antero-posterior, and vertical directions respectively, are illustrated in Fig. 1A.
Vertical work (Wv) was defined as the positive work needed to move the CoMPACS against gravity, and it was calculated as the sum of the increments of the vertical kinetic energy (Ekv) plus the potential energy (Epot) along the time curve with:

$$E_{{{text{kv}}}} = { }frac{1}{2} m overrightarrow{V_{v}}left( {text{t}} right)^{{2}} _{/R}$$
(3)

and

$$E_{{{text{pot}}}} = mgh_{{/{text{R}}}}$$
(4)

where (overrightarrow {{V_{v} }}) (t)/R is the linear vertical velocity of the CoMPACS in R, g = 9.81 m s−2 is the acceleration due to gravity, and h/R is the height of the CoMPACS in R.
The external work (Wext), corresponding to the positive external work needed to raise and accelerate the CoMPACS, was computed as the sum of the increments of the external mechanical energy (Eext) along the time curve with:

$$E_{{{text{ext}}}} = E_{{{text{pot}}}} + E_{{{text{kv}}}} + E_{{{text{kf}}}}$$
(5)

The energy recovered (called recovery rate (RR)10) by the CoMPACS in the sagittal plane was computed with the following formula17:

$$RR = { 1}00frac{{W{text{kf}} + W{text{v}} – W{text{ext}}}}{{W{text{kf}} + W{text{v}}}}$$
(6)

RR is the percentage of kinetic energy converted into potential energy7,24,25,34,35 and vice versa.
In the present study, internal work was also considered in order to encompass the coordination between all body segments. Based on the assumption of a conservative Poly-Articulated Model (PAM), internal work (Wint) was computed as the sum of the increments of the Eint,k along the time curve with:

$$E_{{{text{int}},{text{k}}}} = frac{1}{2}~mathop sum limits_{{i = 1}}^{{33}} (m_{i} overrightarrow {{V_{{~i}} }} left( {text{t}} right)^{{text{2}}} _{{/{text{R}}*}} + m_{{text{i}}} K_{{i^{2} }} {text{ }} times vec{omega }^{2} _{i} /_{{{text{R}}*}} )$$
(7)

where mi is the mass of the ith segment, (overrightarrow {{V_{i} }})(t)/R* the linear velocity of its CoM in the sagittal plane of the barycentric coordinate system (R*), Ki its radius of gyration around its CoM, and (vec{omega }_{i})2/R* its angular velocity in R* 36.
The total mechanical energy of the PACS (Etot) was computed as the sum of the internal kinetic energy (Eint,k) of each segment, plus the potential energy (Epot), and the forward (Ekf ) and vertical (Ekv ) kinetic energy of the CoMPACS in the sagittal plane21,25,37,38:

$$E_{{{text{tot}}}} = E_{{{text{int}},{text{k}}}} + E_{{{text{pot}}}} + E_{{{text{kf}}}} + E_{{{text{kv}}}}$$
(8)

Finally, the total mechanical power (PmecaTot) was used to assess the amount of energy spent or gained by the CoMPACS per unit of time (Δt):

$$P_{{{text{mecaTot}}}} = frac{{W{text{ext}}}}{Delta t} + frac{{{ }W{text{int}}}}{Delta t} = P_{{{text{ext}}}} + P_{{{text{int}}}}$$
(9)

Calculation of internal efforts
The resultant joint moments at the wrist, elbow, shoulder, neck, and back joints were calculated using a bottom-up Newton–Euler recursive algorithm39. Cardanic angles were used to represent the rotation of the segments coordinate system relative to the global coordinate system40. The segment masses, inertia tensors, and center of mass locations were estimated for each subject according to the scaling equations proposed in Dumas et al. (2007)41. In order to estimate the muscular torque produced at all the joints of the upper-limbs, shoulders, neck, and back, the Moment Cost Function (MCF in kg m2 s−2, 42) was computed as follows:

$${text{MCF}} = sqrt {M_{L_wt}^{2} } + sqrt {M_{R_wt}^{2} } + sqrt {M_{L_el}^{2} } + sqrt {M_{R_el}^{2} } + sqrt {M_{L_sh}^{2} } + sqrt {M_{R_sh}^{2} } + sqrt {M_{back}^{2} } + sqrt {M_{neck}^{2} }$$
(10)

where ML_wt, MR_wt, ML_el, MR_el, ML_sh, MR_sh, Mback, and Mneck are the mean values over a PACS gait cycle of the three-dimensional left and right wrist, left and right elbow, left and right shoulder, top of the back and neck moments, respectively. (sqrt {{text{M}}^{2} }) represents the Euclidian norm of M, i.e. (sqrt {sumnolimits_{i = 1}^{3} {left( {M_{i} } right)^{2} } }), with Mi the i-th component of the vector M.
We summed the MCF values of the two individuals of each pair to obtain the total moment cost function (TotMCF). The TotMCF allows to quantify the global muscular effort developed at the upper-limbs of the PACS during one gait cycle of the carrying of the load. The MCF difference (∆MCF) between the two individuals was also computed to investigate whether the volunteers developed the same efforts while carrying the object.
Kinetic synergy analysis
We extracted the synergies by using a principal component analysis (PCA) applied to the wrist, elbow, shoulder, back, and neck joint moment on the right and left sides of the body. The PCA was used to reduce data dimensionality13,35,43. It consisted in the eigen-decomposition of the co-variance matrix of the joint moment data (Matlab eig function). The joint moments data were arranged in time × joint moment matrices. We called the eigenvectors extracted from the PCA synergy vectors13. The number of synergies was determined from the VAF (Variance Accounted For), which corresponds to the cumulative sum of the eigenvalues, ordered from the greatest to the lowest value, normalized by the total variance computed as the sum of all eigenvalues. We defined the number of synergies as the first number for which the VAF was greater than 0.9. The synergy vectors retained were then rotated using a Varimax rotation method to improve interpretability44.
We extracted the synergy vectors for each experimental condition separately. We first performed an analysis on each individual separately. In this analysis the initial data matrices were constituted of all available time frames in line, concatenated, and of eight columns corresponding to each joint moment, namely the right wrist, left wrist, right elbow, left elbow, right shoulder, left shoulder, back, and neck. The values in the matrix corresponded to the norm of the joint moment vector at a given time frame. We then performed a second analysis to identify possible co-variations between the joint moments of the two participants in each pair. The columns of the initial matrices were thus constituted of the joint moments of the two loaded arms, i.e., the right wrist, elbow, and shoulder joint moments of participant #1, plus the left wrist, elbow and shoulder joint moments of participant #2. The synergy vectors were compared across conditions by computing Pearson’s r correlations on their PCA weightings, after being matched together, also using Pearson’s r to identify the best matches.
Statistical analysis
We used generalized linear mixed models (GLMM)45 to compare the velocity and the vertical amplitude of the CoMPACS, the length and duration of the gait cycles, the recovery rate, the external, internal, and total mechanical power produced by the PACS, the TotMCF and ∆MCF, the number of synergies, as well as the Pearson’s r-values across conditions.
The experimental condition was entered as a fixed factor in the model, and individuals as a random variable. We used a Gaussian GLMM for all variables, except for the comparison of the number of synergies across conditions, for which a Poisson GLMM was used. For Gaussian GLMMs we systematically inspected the normality of the model residuals with Q-Q plots. We used the functions lmer() and glmer() of the R package lme4 46 to run the Gaussian and the Poisson mixed models, respectively. The effect of experimental conditions was tested by comparing the deviance of the model with and without the fixed factor with a χ2 test. Multiple comparisons across experimental conditions were performed with the function glht() of the multcomp R package47 using the default Tukey test as post-hoc test. Pearson’s r were Fisher Z-transformed before running the analyses. The significance threshold was set to 0.05. All data in the text are given as mean ± SD. Since our sample size was low, which could lead to inflate the Type II error (not rejecting H0 when H0 is false), we followed the recommendations of Nakagawa & Foster (2004)49 and provide in the Supplemental Table S1 the value of Cohen d standardized effect size50, along with its 95% confidence interval51, for each studied parameter and each paired comparison between conditions. A confidence interval that largely extends on both sides of zero indicates an absence of effect that would probably not change with increasing the sample size.
Ethics statement
All methods used in this study were carried out in accordance with relevant guidelines and regulations. More

1.
Tilman, D., Balzer, C., Hill, J. & Befort, B. L. Global food demand and the sustainable intensification of agriculture. Proc. Natl Acad. Sci. USA 108, 20260–20264 (2011).
ADS  CAS  PubMed  Article  Google Scholar 
2.
Alexandratos, N. & Bruinsma, J. World Agriculture Towards 2030/2050: The 2012 Revision (FAO, ESA Working paper, Rome, 2012).

3.
Carlson, K. M. et al. Greenhouse gas emissions intensity of global croplands. Nat. Clim. Change 7, 63–68 (2017).
ADS  CAS  Article  Google Scholar 

4.
Houghton, R. A. et al. Carbon emissions from land use and land-cover change. Biogeosciences 9, 5125–5142 (2012).
ADS  CAS  Article  Google Scholar 

5.
Cassman, K. G., Dobermann, A., Walters, D. T. & Yang, H. Meeting cereal demand while protecting natural resources and improving environmental quality. Annu. Rev. Environ. Resour. 28, 315–358 (2003).
Article  Google Scholar 

6.
Laurance, W. F., Sayer, J. & Cassman, K. G. Agricultural expansion and its impacts on tropical nature. Trends Ecol. Evol. 29, 107–116 (2014).
PubMed  Article  Google Scholar 

7.
Lobell, D. B., Schlenker, W. & Costa-Roberts, J. Climate trends and global crop production since 1980. Science 333, 616–620 (2011).
ADS  CAS  Article  Google Scholar 

8.
Asseng, S. et al. Rising temperatures reduce global wheat production. Nat. Clim. Change 5, 143–147 (2015).
ADS  Article  Google Scholar 

9.
Rosenzweig, C. et al. Assessing agricultural risks of climate change in the 21st century in a global gridded crop model intercomparison. Proc. Natl Acad. Sci. USA 111, 3268–3273 (2014).
ADS  CAS  Article  Google Scholar 

10.
Mueller, N. D. et al. Closing yield gaps through nutrient and water management. Nature 490, 254–257 (2012).
ADS  CAS  PubMed  Article  Google Scholar 

11.
Jägermeyr, J., Pastor, A., Biemans, H. & Gerten, D. Reconciling irrigated food production with environmental flows for Sustainable Development Goals implementation. Nat. Commun. 8, 15900 (2017).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

12.
Schauberger, B. et al. Consistent negative response of US crops to high temperatures in observations and crop models. Nat. Commun. 8, 13931 (2017).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

13.
Tack, J., Barkley, A. & Hendricks, N. Irrigation offsets wheat yield reductions from warming temperatures. Environ. Res. Lett. 12, 114027 (2017).
ADS  Article  Google Scholar 

14.
Troy, T. J., Kipgen, C. & Pal, I. The impact of climate extremes and irrigation on US crop yields. Environ. Res. Lett. 10, 054013 (2015).
ADS  Article  Google Scholar 

15.
Elliott, J. et al. Constraints and potentials of future irrigation water availability on agricultural production under climate change. Proc. Natl Acad. Sci. USA 111, 3239–3244 (2014).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

16.
Jägermeyr, J. et al. Integrated crop water management might sustainably halve the global food gap. Environ. Res. Lett. 11, 025002 (2016).
ADS  Article  Google Scholar 

17.
Rosegrant, M. W., Ringler, C. & Zhu, T. J. Water for agriculture: maintaining food security under growing scarcity. Annu. Rev. Environ. Resour. 34, 205–222 (2009).
Article  Google Scholar 

18.
Siebert, S. & Döll, P. Quantifying blue and green virtual water contents in global crop production as well as potential production losses without irrigation. J. Hydrol. 384, 198–217 (2010).
ADS  Article  Google Scholar 

19.
Li, X. & Troy, T. J. Changes in rainfed and irrigated crop yield response to climate in the western US. Environ. Res. Lett. 13, 064031 (2018).
ADS  Article  Google Scholar 

20.
Neverre, N., Dumas, P., Nassopoulos, H. Large-scale water scarcity assessment under global changes: insights from a hydroeconomic framework. Hydrol. Earth Syst. Sci. Discuss. Preprint at https://doi.org/10.5194/hess-2015-502 (2016).

21.
Lesk, C., Rowhani, P. & Ramankutty, N. Influence of extreme weather disasters on global crop production. Nature 529, 84–87 (2016).
ADS  CAS  PubMed  Article  Google Scholar 

22.
Lobell, D. B., Sibley, A. & Ortiz-Monasterio, J. I. Extreme heat effects on wheat senescence in India. Nat. Clim. Change 2, 186–189 (2012).
ADS  Article  Google Scholar 

23.
Gourdji, S. M., Sibley, A. M. & Lobell, D. B. Global crop exposure to critical high temperatures in the reproductive period: historical trends and future projections. Environ. Res. Lett. 8, 024041 (2013).
ADS  Article  Google Scholar 

24.
Espe, M. B. et al. Point stresses during reproductive stage rather than warming seasonal temperature determine yield in temperate rice. Glob. Change Biol. 23, 4386–4395 (2017).
ADS  Article  Google Scholar 

25.
Martre, P. et al. Multimodel ensembles of wheat growth: many models are better than one. Glob. Change Biol. 21, 911–925 (2015).
ADS  Article  Google Scholar 

26.
Li, T. et al. Uncertainties in predicting rice yield by current crop models under a wide range of climatic conditions. Glob. Change Biol. 21, 1328–1341 (2015).
ADS  CAS  Article  Google Scholar 

27.
Müller, C. et al. Global gridded crop model evaluation: benchmarking, skills, deficiencies and implications. Geosci. Model Dev. 10, 1403–1422 (2017).
ADS  Article  Google Scholar 

28.
Wang, X. et al. Emergent constraint on crop yield response to warmer temperature from field experiments. Nat. Sustain. https://doi.org/10.1038/s41893-020-0569-7 (2020).

29.
Folberth, C. et al. Uncertainty in soil data can outweigh climate impact signals in global crop yield simulations. Nat. Commun. 7 https://doi.org/10.1038/ncomms11872 (2016).

30.
Elliott, J. et al. The Global Gridded Crop Model Intercomparison: data and modeling protocols for Phase 1 (v1.0). Geosci. Model Dev. 8, 261–277 (2015).
ADS  Article  Google Scholar 

31.
Roberts, M. J. et al. Comparing and combining process-based crop models and statistical models with some implications for climate change. Environ. Res. Lett. 12, 095010 (2017).
ADS  Article  Google Scholar 

32.
Raftery, A. E., Gneiting, T., Balabdaoui, F. & Polakowski, M. Using Bayesian model averaging to calibrate forecast ensembles. Mon. Weather Rev. 133, 1155–1174 (2005).
ADS  Article  Google Scholar 

33.
Monfreda, C., Ramankutty, N. & Foley, J. A. Farming the planet:2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000. Glob. Biogeochem. Cycles 22, GB1022 (2008).
ADS  Article  CAS  Google Scholar 

34.
Fekete, B. M., Vorosmarty, C. J., & Grabs, W. Global, composite runoff fields based on observed river discharge and simulated water balances. Report No. 22. (World Meteorological Organization–Global Runoff Data Center, Koblenz, Germany, 1999).

35.
Wang, X. X. et al. Taking account of governance: Implications for land-use dynamics, food prices, and trade patterns. Ecol. Econ. 122, 12–24 (2016).
Article  Google Scholar 

36.
Neumann, K. et al. Exploring global irrigation patterns: a multilevel modelling approach. Agric. Syst. 104, 703–713 (2011).
Article  Google Scholar 

37.
Fekete, B. M. & Vörösmarty, C. J. The current status of global river discharge monitoring and potential new technologies complementing traditional discharge measurements. IAHS Publ. 309, 129–136 (2007).
Google Scholar 

38.
Hanasaki, N., Yoshikawa, S., Pokhrel, Y. & Kanae, S. A global hydrological simulation to specify the sources of water used by humans. Hydrol. Earth Syst. Sci. 22, 789–817 (2018).
ADS  Article  Google Scholar 

39.
Siebert, S., Henrich, V., Frenken, K., & Burke, J. Update of the digital global map of irrigation areas to version 5 (Rheinische Friedrich-Wilhelms-Universität, Bonn, Germany, Food and Agriculture Organization of the United Nations, Rome, Italy, 2013).

40.
Rosa, L. et al. Potential for sustainable irrigation expansion in a 3 °C warmer climate. Proc. Natl Acad. Sci. USA 117, 29526–29534 (2020).
ADS  CAS  PubMed  Article  Google Scholar 

41.
Rosa L., Chiarelli D. D., Rulli M. C., Dell’Angelo J. & D’Odorico P. Global agricultural economic water scarcity. Sci. Adv. 6, eaaz6031 (2020).

42.
Rosegrant, M. W. & Cai, X. Global water demand and supply projections. Water Int. 27, 170–182 (2002).
Article  Google Scholar 

43.
Berkoff, J. China: the South–North Water Transfer Project—is it justified? Water Policy 5, 1–28 (2003).
Article  Google Scholar 

44.
Schmitz, C. et al. Blue water scarcity and the economic impacts of future agricultural trade and demand. Water Resour. Res. 49, 3601–3617 (2013).
ADS  Article  Google Scholar 

45.
Biewald, A., Rolinski, S., Lotze-Campen, H., Schmitz, C. & Dietrich, J. P. Valuing the impact of trade on local blue water. Ecol. Econ. 101, 43–53 (2014).
Article  Google Scholar 

46.
Portmann, F. T., Siebert, S. & Döll, P. MIRCA2000—Global monthly irrigated and rainfed crop areas around the year 2000: a new high-resolution data set for agricultural and hydrological modeling. Glob. Biogeochem. Cycles 24, GB1011 (2010).
ADS  Article  CAS  Google Scholar 

47.
Muller, C. et al. The Global Gridded Crop Model Intercomparison phase 1 simulation dataset. Sci. Data 6, 50 (2019).
PubMed  PubMed Central  Article  Google Scholar 

48.
Kobayashi, K. & Salam, M. U. Comparing simulated and measured values using mean squared deviation and its components. Agron. J. 92, 345–352 (2000).
Article  Google Scholar 

49.
Deutscher Wetterdienst Frankfurt M. Mitteilungen des Deutschen Wetterdienstes: Windschutzanlagen auf der hohen Rhön (Dt. Wetterdienst P.11, 1954).

50.
Gornott, C. & Wechsung, F. Statistical regression models for assessing climate impacts on crop yields: a validation study for winter wheat and silage maize in Germany. Agric. Meteorol. 217, 89–100 (2016).
Article  Google Scholar 

51.
Harris, I., Jones, P. D., Osborn, T. J. & Lister, D. H. Updated high-resolution grids of monthly climatic observations – the CRU TS3.10 Dataset. Int J. Climatol. 34, 623–642 (2014).
Article  Google Scholar 

52.
Jägermeyr, J. et al. Water savings potentials of irrigation systems: global simulation of processes and linkages. Hydrol. Earth Syst. Sci. 19, 3073–3091 (2015).
ADS  Article  Google Scholar 

53.
Pastor, A. V., Ludwig, F., Biemans, H., Hoff, H. & Kabat, P. Accounting for environmental flow requirements in global water assessments. Hydrol. Earth Syst. Sci. 18, 5041–5059 (2014).
ADS  Article  Google Scholar 

54.
Poff, N. L. et al. The natural flow regime: a paradigm for river conservation and restoration. BioScience 47, 769–784 (1997).
Article  Google Scholar  More

1.
Tinius, A. & Russell, A. P. Points on the curve: An analysis of methods for assessing the shape of vertebrate claws. J. Morphol. 278, 150–169 (2017).
Article  Google Scholar 
2.
Ribas, S. C. et al. Structure of claws and toes of two tropidurid lizard species of Restinga from Southeastern Brazil: Adaptations to the vertical use of the habitat. Rev. Chilena Hist. Nat. 77, 599 (2004).
Google Scholar 

3.
Burnham, D. A., Feduccia, A., Martin, L. D. & Falk, A. R. Tree climbing: A fundamental avian adaptation. J. Syst. Paleontol. 9, 103–107 (2011).
Article  Google Scholar 

4.
Crandell, K. E., Herrel, A., Sasa, M., Losos, J. B. & Autumn, K. Stick or grip? Co-evolution of adhesive toepads and claws in Anolis lizards. Zoology 117, 363–369 (2014).
Article  Google Scholar 

5.
Zani, P. A. The comparative evolution of lizard claw and toe morphology and clinging performance. J. Evol. Biol. 13, 316–325 (2000).
Article  Google Scholar 

6.
Feduccia, A. Evidence from claw geometry indicating arboreal habits of Archaeopteryx. Science 259, 790–793 (1993).
ADS  CAS  Article  Google Scholar 

7.
Manning, P. L., Payne, D., Pennicott, J., Barrett, P. M. & Ennos, R. A. Dinosaur killer claws or climbing crampons? Biol. Lett. 2, 110–112 (2006).
Article  Google Scholar 

8.
Pattrick, J. G., Labonte, D. & Federle, W. Scaling of claw sharpness: mechanical constraints reduce attachment performance in larger insects. J. Exp. Biol. 221, jeb188391 (2018).
Article  Google Scholar 

9.
Dai, Z., Gorb, S. N. & Schwarz, U. Roughness-dependent friction force of the tarsal claw system in the beetle Pachnoda marginata (Coleoptera, Scarabaeidae). J. Exp. Biol. 205, 2479–2488 (2002).
PubMed  Google Scholar 

10.
Heethoff, M. & Koerner, L. Small but powerful: The oribatid mite Archegozetes longisetosus Aoki (Acari, Oribatida) produces disproportionately high forces. J. Exp. Biol. 210, 3036–3042 (2007).
Article  Google Scholar 

11.
Büscher, T. H. & Gorb, S. N. Complementary effect of attachment devices in stick insects (Phasmatodea). J. Exp. Biol. 222, 209833 (2019).
Article  Google Scholar 

12.
Pfingstl, T., Kerschbaumer, M. & Shimano, S. Get a grip-evolution of claw shape in relation to microhabitat use in intertidal arthropods (Acari, Oribatida). PeerJ 8, e8488 (2020).
Article  Google Scholar 

13.
Herrel, A. Correlations between habitat use and body shape in a phrynosomatid lizard (Urosaurus ornatus): A population-level analysis. Biol. J. Lin. Soc. 74, 305–314 (2001).
Article  Google Scholar 

14.
Tulli, M. J., Cruz, F. B., Herrel, A., Vanhooydonck, B. & Abdala, V. The interplay between claw morphology and microhabitat use in neotropical iguanian lizards. Zoology 112, 379–392 (2009).
CAS  Article  Google Scholar 

15.
Klingenberg, C. P. & Gidaszewski, N. A. Testing and quantifying phylogenetic signals and homoplasy in morphometric data. Syst. Biol. 59, 245–261 (2010).
CAS  Article  Google Scholar 

16.
Muñoz-Muñoz, F., Talavera, S. & Pagès, N. Geometric morphometrics of the wing in the subgenus Culicoides (Diptera: Ceratopogonidae): From practical implications to evolutionary interpretations. J. Med. Entomol. 48, 129–139 (2011).
Article  Google Scholar 

17.
Khang, T. F., Soo, O. Y. M., Tan, W. B. & Lim, L. H. S. Monogenean anchor morphometry: Systematic value, phylogenetic signal, and evolution. PeerJ 4, e1668 (2016).
Article  Google Scholar 

18.
Blomberg, S. P. & Garland, T. Tempo and mode in evolution: Phylogenetic inertia, adaptation and comparative methods. J. Evol. Biol. 15, 899–910 (2002).
Article  Google Scholar 

19.
Pfingstl, T., Baumann, J. & Lienhard, A. The Caribbean enigma: the presence of unusual cryptic diversity in intertidal mites (Arachnida, Acari, Oribatida). Organ. Divers. Evol. 19, 609–623 (2019).
Article  Google Scholar 

20.
Blomberg, S. P., Garland, T. & Ives, A. R. Testing for phylogenetic signal in comparative data: Behavioral traits are more labile. Evolution 57, 717–745 (2003).
Article  Google Scholar 

21.
Caumul, R. & Polly, P. D. Phylogenetic and environmental components of morphological variation: Skull, mandible, and molar shape in marmots (Marmota, Rodentia). Evolution 59, 2460–2472 (2005).
Article  Google Scholar 

22.
Klingenberg, C. P. There’s something afoot in the evolution of ontogenies. BMC Evol. Biol. 10, 221 (2010).
Article  Google Scholar 

23.
Edwards, E. J. & Donoghue, M. J. Pereskia and the origin of the cactus life-form. Am. Nat. 167, 777–793 (2006).
Article  Google Scholar 

24.
Ge, D., Yao, L., Xia, L., Zhang, Z. & Yang, Q. Geometric morphometric analysis of skull morphology reveals loss of phylogenetic signal at the generic level in extant lagomorphs (Mammalia: Lagomorpha). Contrib. Zool. 84, 267–284 (2015).
Article  Google Scholar 

25.
Pfingstl, T. Habitat use, feeding and reproductive traits of rocky-shore intertidal mites from Bermuda (Oribatida: Fortuyniidae and Selenoribatidae). Acarologia 53, 369–382 (2013).
Article  Google Scholar 

26.
Pfingstl, T., Lienhard, A. & Baumann, J. New and cryptic species of intertidal mites (Acari, Oribatida) from the Western Caribbean: An integrative approach. Int. J. Acarology 45, 10–25 (2019).
Article  Google Scholar 

27.
Rohlf, F. J. TpsDig2. Version 2.31 (2017).

28.
Klingenberg, C. P. MorphoJ: An integrated software package for geometric morphometrics. Mol. Ecol. Resour. 11, 353–357 (2011).
Article  Google Scholar 

29.
Hammer, Ø., Harper, D. A. T. & Ryan, P. D. PAST: Paleontological statistics software package for education and data analysis. Palaeontol. Electron. 4, 1 (2001).
Google Scholar 

30.
Bookstein, F. L. Morphometric Tools for Landmark Data: Geometry and Biology (Cambridge University Press, Cambridge, 1991).
Google Scholar 

31.
Adams, D. C., Rohlf, F. J. & Slice, D. E. Geometric morphometrics: Ten years of progress following the ‘revolution’. Italian J. Zool. 71, 5–16 (2004).
Article  Google Scholar 

32.
Cavalcanti, M. J. Geometric morphometric analysis of head shape variation in four species of hammerhead sharks (Carcharhiniformes: Sphyrnidae). Morphometrics 14, 97 (2004).
Article  Google Scholar 

33.
Schäffer, S., Stabentheiner, E., Shimano, S. & Pfingstl, T. Leaving the tropics: The successful colonization of cold temperate regions by Dolicheremaeus dorni (Acari, Oribatida). J. Zool. Syst. Evol. Res. 56, 505–518 (2018).
Article  Google Scholar 

34.
Otto, J. C. & Wilson, K. Assessment of the usefulness of ribosomal 18S and mitochondrial COI sequences in Prostigmata phylogeny. Acarology: Proceedings of the 10th International Congress (2001).

35.
Schäffer, S., Krisper, G., Pfingstl, T. & Sturmbauer, C. Description of Scutovertex pileatus sp. nov. (Acari, Oribatida, Scutoverticidae) and molecular phylogenetic investigation of congeneric species in Austria. Zoologischer Anzeiger 247, 249–258 (2008).
Article  Google Scholar 

36.
Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35, 1547–1549 (2018).
CAS  Article  Google Scholar 

37.
Ronquist, F. & Huelsenbeck, J. P. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574 (2003).
CAS  Article  Google Scholar 

38.
Guindon, S. et al. New algorithms and methods to estimate maximum-likelihood phylogenies: Assessing the performance of PhyML 3.0. Syst. Biol. 59, 307–321 (2010).
CAS  Article  Google Scholar 

39.
Rambaut, A. FigTree: tree figure drawing tool. version 1.3.1. Software available at http://tree.bio.ed.ac.uk/software/figtree/ (2009).

40.
Rambaut, A. & Drummond, A. Tracer: MCMC trace analysis tool. version 1.6. 0. http://tree.bio.ed.ac.uk/software/tracer/ (2007).

41.
Adams, D. C. A generalized K statistic for estimating phylogenetic signal from shape and other high-dimensional multivariate data. Syst. Biol. 63, 685–697 (2014).
Article  Google Scholar 

42.
Adams, D. C. & Collyer, M. L. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution 73, 2352–2367 (2019).
Article  Google Scholar 

43.
Adams, D. C., Collyer, M. L., Kaliontzopoulou, A. & Balken E. geomorph: Software for geometric morphometric analyses. R package version 3.3.2. https://cran.r-project.org/package=geomorph (2021).

44.
Collyer, M. L. & Adams, D. C. RRPP: Linear model evaluation with randomized residuals in a permutation procedure. https://cran.r-project.org/web/packages/RRPP (2021).

45.
Collyer, M. L. & Adams, D. C. RRPP: An R package for fitting linear models to high-dimensional data using residual randomization. Methods Ecol. Evol. 9(2), 1772–1779 (2018).
Article  Google Scholar 

46.
Pennell, M. W. et al. geiger v2.0: An expanded suite of methods for fitting macroevolutionary models to phylogenetic trees. Bioinformatics 30, 2216–2218 (2014).
CAS  Article  Google Scholar  More

1.
Acuto, M., Parnell, S. & Seto, K. C. Building a global urban science. Nat. Sustain. 1, 2–4 (2018).
Article  Google Scholar 
2.
Anderies, J. M., Folke, C., Walker, B. & Ostrom, E. Aligning key concepts for global change policy: robustness, resilience, and sustainability. Ecol. Soc. 18, 8 (2013).
Article  Google Scholar 

3.
Brenner, N. & Schmid, C. Planetary urbanization. In Implosions/explosions: Towards A Study of Planetary Urbanization (ed. Brenner, N.) 142–463 (Jovis Verlag, 2014).

4.
Elmqvist, T. et al. Sustainability and resilience for transformation in the urban century. Nat. Sustain. 2, 267 (2019).
Article  Google Scholar 

5.
McPhearson, T. Scientists must have a say in the future of cities. Nature 538, 165–166 (2016).
CAS  Article  Google Scholar 

6.
Groffman, P. M. et al. Moving towards a new Urban Systems Science. Ecosystems https://doi.org/10.1007/s10021-016-0053-4 (2016).

7.
Pataki, D. E. Grand challenges in urban ecology. Front. Ecol. Evol. 3, 57 (2015).
Article  Google Scholar 

8.
National Science Foundation. What is Convergence?|NSF-National Science Foundation. https://www.nsf.gov/od/oia/convergence/index.jsp (Accessed April 23, 2019).

9.
Ramaswami, A. et al. Sustainable Urban Systems: Articulating a Long-Term Convergence Research Agenda. Vol. 31 (National Science Foundation, 2018).

10.
Grimm, N. B., Pickett, S. T. A., Hale, R. L. & Cadenasso, M. L. Does the ecological concept of disturbance have utility in urban social-ecological-technological systems? Ecosyst. Health Sustain 3, e01255 (2017).
Article  Google Scholar 

11.
McPhearson, T. et al. Advancing urban ecology towards a science of cities. BioScience 66, 198–212 (2016).
Article  Google Scholar 

12.
United Nations. Transforming Our World: the 2030 Agenda for Sustainable Development. https://sustainabledevelopment.un.org/post2015/transformingourworld (Accessed March 6, 2020).

13.
Seto, K. C. et al. Urban land teleconnections and sustainability. Proc. Natl Acad. Sci. USA 109, 7687–7692 (2012).
CAS  Article  Google Scholar 

14.
Folke, C., Biggs, R., Norstrom, A. V., Reyers, B. & Rockstrom, J. Social-ecological resilience and biosphere-based sustainability science. Ecol. Soc. 21, 41 (2016).
Article  Google Scholar 

15.
Grimm, N. B., Cook, E. M., Hale, R. L. & Iwaniec, D. M. A broader framing of ecosystem services in cities: benefits and challenges of built, natural, or hybrid system function. In The Routledge Handbook of Urbanization and Global Environmental Change (eds. Seto, K. C., Solecki, W. D. & Griffith, C. A.) 203–212 (Routledge, 2016).

16.
Pelling, M. & Manuel-Navarrete, D. From resilience to transformation: the adaptive cycle in two Mexican urban centers. Ecol. Soc. 16, 11 (2011).
Article  Google Scholar 

17.
Meerow, S., Newell, J. P. & Stults, M. Defining urban resilience: a review. Landsc. Urban Plan. 147, 38–49 (2016).
Article  Google Scholar 

18.
Shane, D. G. Recombinant Urbanism: Conceptual Modeling in Architecture. (John Wiley & Sons, 2005).

19.
McHale, M. R. et al. The new global urban realm: complex, connected, diffuse, and diverse social-ecological systems. Sustainability 7, 5211–5240 (2015).
Article  Google Scholar 

20.
Ellin, N. Integral urbanism: a context for urban design. In Resilience in ecology and urban design: linking theory and practice for sustainable cities (eds. Pickett, S. T. A., Cadenasso, M. L. & McGrath, B.) 63–78 (Springer, 2013).

21.
Marcotullio, P. J. & Solecki, W. What is a city? an essential definition for sustainability. In Urbanization and Sustainability: Linking Urban Ecology, Environmental Justice, and Environmental Change (eds. Boone, C. G. & Fragkias, M.) 11–25 (Springer, 2013).

22.
Burch, W. R., Jr., Machlis, G. E. & Force, J. E. The Structure and Dynamics of Human Ecosystems: toward A Model for Understanding and Action. (Yale University Press, 2017).

23.
Redman, C., Grove, J. M. & Kuby, L. Toward a Unified Understanding of Human Ecosystems: Integrating Social Sciences Into Long-term Ecological Research. Vol. 13 (LTER Network, 2000).

24.
Barnett, R. & Margetts, J. Disturbanism in the South Pacific: disturbance ecology as a basis for urban resilience in small island states. In Resilience in Ecology and Urban Design: Linking Theory and Practice for Sustainable Cities (eds. Pickett, S. T. A., Cadenasso, M. L. & McGrath, B.) 443–459 (Springer, 2013).

25.
Folke, C. et al. Resileince and Sustainable Development: Building Adaptive Capacity in A World of Transformations. (Ministry of the Environment, 2002).

26.
Scheffer, M., Westley, F., Brock, W. A. & Holmgren, M. Dynamic interaction of societies and ecosystems–linking theories from ecology, economy, and sociology. In Panarchy: Understanding Transformations in Human and Natural Systems (eds. Gunderson, L. H. & Holling, C. S.) 195–239 (Island Press, 2002).

27.
Pickett, S. T. A. et al. Dynamic heterogeneity: a framework to promote ecological integration and hypothesis generation in urban systems. Urban Ecosyst. 20, 1–14 (2017).
Article  Google Scholar 

28.
Wu, J. G. & Loucks, O. L. From balance of nature to hierarchical patch dynamics: a paradigm shift in ecology. Q. Rev. Biol. 70, 439–466 (1995).
Article  Google Scholar 

29.
Boone, C. G. et al. Reconceptualizing land for sustainable urbanity. In Rethinking Urban Land Use in A Global Era (eds. Seto, K. C. & Reenberg, A.) 313–330 (MIT Press, 2014).

30.
Machlis, G. E., Force, J. E. & Burch, W. R. The human ecosystem 1. The human ecosystem as an organizing concept in ecosystem management. Soc. Nat. Resour. 10, 347–367 (1997).
Article  Google Scholar 

31.
Cadenasso, M. L. & Pickett, S. T. A. Three tides: the development and state of the art of urban ecological science. In Resilience in Ecology and Urban Design: Linking Theory and Practice for Sustainable Cities (eds. Pickett, S. T. A., Cadenasso, M. L. & McGrath, B.) 29–46 (Springer, 2013).

32.
Collins, S. L. et al. An integrated conceptual framework for long-term social-ecological research. Front. Ecol. Environ. 9, 351–357 (2011).
Article  Google Scholar 

33.
Naveh, Z. The total human ecosystem: integrating ecology and economics. BioScience 50, 357–361 (2000).
Article  Google Scholar 

34.
Pickett, S. T. A. & Cadenasso, M. L. Ecosystem as a multidimensional concept: meaning, model and metaphor. Ecosystems 5, 1–10 (2002).
Article  Google Scholar 

35.
Alberti, M. Advances in Urban Ecology: Integrating Humans and Ecological Processes in Urban Ecosystems. (Springer, 2008).

36.
Pickett, S. T. A. & Grove, J. M. Urban ecosystems: what would Tansley do? Urban Ecosyst. 12, 1–8 (2009).
Article  Google Scholar 

37.
Lachmund, J. Greening Berlin. (MIT Press, 2013).

38.
Rademacher, A., Cadenasso, M. L. & Pickett, S. T. A. From feedbacks to coproduction: toward an integrated conceptual framework for urban ecosystems. Urban Ecosyst. https://doi.org/10.1007/s11252-018-0751-0 (2018).

39.
Johnson, E. A. & Miyanishi, K. (eds.) Plant Disturbance Ecology: the Process and the Response. (Academic Press, Burlington, 2007).

40.
Pickett, S. T. A. & White, P. S. (eds.) The Ecology of Natural Disturbance and Patch Dynamics. (Academic Press, Orlando, 1985).

41.
Schumpeter, J. A. The Theory of Economic Development: an Inquiry Into Profits, Capital, Credit, Interest, and the Business Cycle. (Transaction Books, 1983).

42.
Peters, D. P. C. et al. Cross-system comparisons elucidate distrubance complexities and generalities. Ecosphere 2, art 81 (2011).
Article  Google Scholar 

43.
Holling, C. S. Engineering resilience versus ecological resilience. In Engineering within Ecological Constraints (ed. Schulze, P. C.) 31–44 (National Academies of Engineering, 1996).

44.
Gunderson, L. H. & Holling, C. S. (eds.) Panarchy: understanding transformations in human and natural systems. (Island Press, Washington DC, 2002).

45.
Walker, B., Holling, C. S., Carpenter, S. R. & Kinzig, A. Resilience, adaptability and transformability in social-ecological systems. Ecol. Soc. 9, Article 5 (2004).
Article  Google Scholar 

46.
McPhearson, T., Andersson, E., Elmqvist, T. & Frantzeskaki, N. Resilience of and through urban ecosystem services. Ecosyst. Services 12, 152–156 (2015).
Article  Google Scholar 

47.
Tidball, K., Frantzeskaki, N. & Elmqvist, T. Traps! An introduction to expanding thinking on persistent maladaptive states in pursuit of resilience. Sustain. Sci. 11, 861–866 (2016).
Article  Google Scholar 

48.
Biggs, R., Westley, F. R. & Carpenter, S. R. Navigating the back loop: fostering social innovation and transformation in ecosystem management. Ecol. Soc. 15, 9 (2010).
Article  Google Scholar 

49.
Changnon, S. A., Kunkel, K. E. & Reinke, B. C. Impacts and responses to the 1995 heat wave: a call to action. Bullet. Am. Meteorol. Soc. 77, 1497–1506 (1996).
Article  Google Scholar 

50.
Borden, K. A. & Cutter, S. L. Spatial patterns of natural hazards mortality in the United States. Int. J. Health Geogr. 7, 64 (2008).
Article  Google Scholar 

51.
Park, R. E. & Burgess, E. W. The City (University of Chicago Press, 1925).

52.
Jacobs, J. The Death and Life of Great American Cities (Random House, 1961).

53.
Lynch, K. Good City Form (MIT Press, 1981).

54.
Shane, D. G. Urban Design Since 1945–A Global Perspective (John Wiley & Sons, Ltd, 2011).

55.
Hamstead, Z., Farmer, C. & McPhearson, T. Landscape-based extreme heat vulnerability assessment. J. Extreme Event. 5, 1–23 (2018).
Google Scholar 

56.
Uejio, C. K. et al. Intra-urban societal vulnerability to extreme heat: the role of heat exposure and the built environment, socioeconomics and neighborhood stability. Health Place 17, 498–507 (2011).
Article  Google Scholar 

57.
Rosenthal, J. K., Kinney, P. L. & Metzger, K. B. Intra-urban vulnerability to heat-related mortality in New York City. 1997–2006. Health Place 30, 45–60 (2014).
Google Scholar 

58.
Madrigano, J., Ito, K., Johnson, S., Kinney, P. L. & Matte, T. A case-only study of vulnerability to heat wave–related mortality in New York City (2000–2011). Environ. Health Perspect. 123, 672–678 (2015).
Article  Google Scholar 

59.
Allen, T. F. H. & Starr, T. B. Hierarchy: Perspectives for Ecological Complexity (2nd edn.) (University of Chicago Press, Chicago, 2017).

60.
McGrath, B. & Shane, G. Introduction: metropolis, megalopolis, and metacity. In The SAGE Handbook of Architectural Theory (eds. Crysler, C. G., Cairns, S. & Heynen, H.) (SAGE, 2012).

61.
Mihaljevic, J. R. (2012). Linking metacommunity theory and symbiont evolutionary ecology. Trends Ecol. Evol. 27, 323–329 (2012).
Article  Google Scholar 

62.
McGrath, B., Sangawongse, S., Thaikatoo, D. & Corte, M. B. The architecture of the metacity: land use change, patch dynamics and urban form in Chiang Mai, Thailand. Urban Plan. 2, 53–71 (2017).
Article  Google Scholar 

63.
Leibold, M. A. The metacommunity concept and its theoretical underpinnings. In The Theory of Ecology (eds. Scheiner, S. M. & Willig, M. R.) 163–183 (University of Chicago Press, 2011).

64.
McGrath, B. & Pickett, S. T. A. The metacity: a conceptual framework for integrating ecology and urban design. Challenges 2011, 55–72 (2011).
Article  Google Scholar 

65.
Batty, M. The New Science of Cities. (MIT Press, 2013).

66.
Gandy, M. Where does the city end? In Implosions/explosions: Towards A Study of Planetary Urbanization (ed. Brenner, N.) 86–89 (jovis Verlag, 2014).

67.
McPhearson, T., Kremer, P. & Hamstead, Z. Mapping ecosystem services in new york city: applying a social-ecological approach in urban vacant land. Ecosyst. Service 11–26, https://doi.org/10.1016/j.ecoser.2013.06.005 (2013).

68.
Kremer, P., Hamstead, Z. & McPhearson, T. A social-ecological assessment of vacant lots in New York City. Landsc. Urban Plann. 218–233, https://doi.org/10.1016/j.landurbplan.2013.05.003 (2013).

69.
Burkholder, S. The new ecology of vacancy: rethinking land use in shrinking cities. Sustainability 4, 1154–1172 (2012).
Article  Google Scholar 

70.
Bowman, A. O. M. & Pagano, M. A. Transforming America’s cities: policies and conditions of vacant land. Urban Affairs Rev. 35, 559–581 (2000).
Article  Google Scholar 

71.
Kabisch N., et al. Nature-Based Solutions to Climate Change Adaptation in Urban Areas—Linkages Between Science, Policy and Practice. 91–109 (Springer, 2017).

72.
Schwarz, K., Berland, A. & Herrmann, D. L. Green, but not just? Rethinking environmental justice indicators in shrinking cities. Sustain. Cities Soc. 41, 816–821 (2018).
Article  Google Scholar 

73.
McDonnell, M. J. & Hahs, A. K. The future of urban biodiversity research: moving beyond the ‘low-hanging fruit’. Urban Ecosyst. 16, 397–409 (2013).
Article  Google Scholar 

74.
Pickett, S. T. A., Kolasa, J. & Jones, C. G. Ecological Understanding: The Nature of Theory and the Theory of Nature (Academic Press, 2007).

75.
Depietri Y. & McPhearson T. Integrating the grey, green, and blue in cities: Nature-based solutions for climate change adaptation and risk reduction. In Nature-based solutions to climate change adaptation in urban areas. Theory and practice of urban sustainability transitions (eds. Kabisch N., Korn H., Stadler J. & Bonn A.) (Springer, Cham, 2017). More

Aquarium experiment and sampling
To examine the effect of changes in species composition on the behaviour of eDNA, we conducted aquarium experiments using two mock fish communities comprising H. neglectus, C. temminckii, O. latipes, R. flumineus, and M. anguillicaudatus. Mock community 1 (MC1) consisted of one individual of each of the five fish species, whereas mock community 2 (MC2) consisted of three H. neglectus individuals and one individual of each of the other four fish species (Fig. 2). We used two aquaria (A and B). Each aquarium was used four times, twice for each mock community, giving two replicates (R1 and R2). This resulted in eight experimental units (2 mock fish communities × 2 aquaria × 2 replicates). Figure 2 shows the experimental setup used in this study.
Figure 2

Experimental setup of the aquarium experiments.

Full size image

To set up the aquaria, 20 L of tap water was added into each aquarium (GEX Co. Ltd., Osaka, Japan) and heated with a heater (Spectrum Brands, Wisconsin, US) until the water temperature reached 25 °C. Water in the two aquaria was maintained at 25 °C and constantly circulated with an aeration device. Before adding fish to the aquaria, the water was sampled for the negative control. The first experimental samples (day 0) were taken 1 h after adding the fish and subsequent samples were taken each day until day 4. At each sampling, two 1-L samples of surface water were collected from each aquarium and then 2 L of tap water was added to each aquarium to maintain the volume of water. The weight of individual fish species was measured using an electronic balance immediately after the final sampling. After each experiment, the two aquaria were bleached before being reused.
In Japan, experiments on fish do not require any legal procedures or permission. However, in order to avoid causing pain to the specimens, the experiments in this study were conducted in accordance with the ARRIVE guidelines, Japanese laws and guidelines for mammals, birds, and reptiles as below; Act on Welfare and Management of Animals (Notice of the Ministry of the Environment No. 105 of October 1, 1973), Standards relating to the Care and Keeping and Reducing Pain of Laboratory Animals (Notice of the Ministry of the Environment No. 88 of 2006), Fundamental Guidelines for Proper Conduct of Animal Experiment and Related Activities in Academic Research Institutions under the jurisdiction of the Ministry of Education (Notice of Ministry of Education No. 71, 2006), and Guidelines for Proper Conduct of Animal Experiments (established by the Science Council of Japan on June 1, 2006).
DNA extraction
Each 1-L water sample was filtered immediately through a GF/F glass fibre filter (nominal pore size = 0.7 μm, diameter = 47 mm; GE Healthcare Japan Corporation, Tokyo, Japan). Filter funnels and measuring cups were bleached after filtration to prevent cross-contamination among the water samples. All filters were stored separately at − 20 °C until DNA extraction. Total eDNA was extracted from each filter using a DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany) and Salivette tubes (Sarstedt AG & Co. KG, Nümbrecht, Germany). Extraction methods were as previously described18 with modifications. A filter sample was placed in the upper part of the Salivette tube and 220 μL of solution containing Buffer AL (200 μL) and Proteinase K (20 μL) was added. The tube containing the filter was incubated at 56 °C for 30 min, then centrifuged at 5000 × g for 3 min, and the solution was collected in the base of the tube. To increase eDNA yield, 220 μL Tris-EDTA (TE) buffer was added to the filter sample and centrifuged at 5000 × g for 1 min. Then, ethanol (200 μL) was added to the collected solution, and the mixture was transferred to a spin column. Total eDNA was eluted in buffer AE (100 μL), following the manufacturer’s instructions. All eDNA samples were stored at − 20 °C prior to qSeq and dPCR.
Quantitative sequencing
Simultaneous quantification and sequencing of the extracted eDNA were performed by qSeq as previously described15,16. First, SPE was performed. The SPE reaction mixture (20 µL) consisted of 1 × PrimeSTAR Max premix (Takara Bio Inc., Kusatsu, Japan), 300 nM of the primer qSeq-MiFish-U-F (Table 1), and extracted DNA (2 µL). The SPE primer qSeq-MiFish-U-F contains an eight-base length random sequence tag, which creates 65,536 different variations, enabling the quantification of up to approximately 1.0 × 105 copies of DNA15. This amount of variation was sufficient to quantify the abundance of eDNA in this study. SPE was initiated by denaturation at 94 °C for 1 min, followed by cooling to 60 °C at 0.3 °C/s, incubation at 60 °C for 1 min, and final extension at 70 °C for 10 min. Subsequently, the excess primer was completely digested by adding exonuclease I (4 µL, 5 U/µL; Takara Bio Inc.) to the SPE mixture. The digestion was performed at 37 °C for 120 min, followed by inactivation of the exonuclease I at 80 °C for 30 min. The first-round PCR mixture (25 µL) contained PrimeSTAR Max premix (12.5 µL), primers qSeq-MiFish-U-R and F2 (300 nM each; Table 1), and the SPE product (2 µL). Following 40 cycles of amplification at 98 °C for 10 s, 55 °C for 5 s, and 72 °C for 5 s, the amplification product was subjected to agarose gel electrophoresis, and the band of the expected size was removed and purified using Nucleospin Gel and PCR Clean-up column (Takara Bio Inc.). The qSeq-MiFish-U-R primer also contains eight N bases to increase the complexity, which improves the sequencing quality, and thus PhiX was not added in this study. Finally, a 2nd-round PCR was performed to add an index for Illumina sequencing as described elsewhere15. The indexed PCR amplicon was purified using AMPure XP beads (Beckman Coulter, Indianapolis, IN) followed by sequencing using a MiSeq platform with MiSeq Reagent Kit v3 for 600 cycles (Illumina). The sequence data obtained in this study were deposited in the DDBJ database under accession numbers SAMD00219124–SAMD00219214.
Table 1 Oligonucleotide sequences used in this study.
Full size table

Data analysis
First, all sequences were assembled and screened by length and quality of reads using the mothur software package (v1.39.5)22. The processed sequence reads were classified using the MiFish pipeline (http://mitofish.aori.u-tokyo.ac.jp/mifish/), with the parameters as previously described23. Subsequently, the representative sequences of individual operational taxonomic units (OTUs) were extracted using the Usearch program (http://www.drive5.com/usearch/). The random sequence tags (RST) at the end of sequences in the OTUs were counted to quantify the environmental DNA from each fish species as described elsewhere16. For comparison, the relative proportion of eDNA from individual species in each sample was calculated from the composition of the sequences of the fish species obtained by qSeq.
Microfluidic digital PCR
Quantification of eDNA was also performed by microfluidic dPCR using the BioMark Real-time System and 12.765 Digital Array (Fluidigm Corporation, South San Francisco, CA, United States) as previously described13. For each sample, the PCR mixture (6 µL) contained 2 × Probe qPCR mix (3.0 µL; Takara Bio Inc.), 20 × binding dye sample loading reagent (0.6 µL; Fluidigm Corporation), forward and reverse primers (900 nM), TaqMan probe (125 nM), ROX solution (0.015 µL), and sample DNA (1.0 µL). We used three sets of primers and probes to quantify the eDNA of H. neglectus, O. latipes, and M. anguillicaudatus (Table 1). PCR was initiated at 98 °C for 2 min, followed by 50 cycles of 98 °C for 10 s and 60 °C for 1 min. The amplification curves obtained from individual reaction chambers of the microfluidic chip were analysed using Fluidigm Digital PCR analysis software (Fluidigm Corporation) to obtain abundance of DNA molecules.
Statistical analysis
We employed Gaussian Type II regression models with the standardised major axis method to determine the relationship between the log10 eDNA abundances obtained from qSeq and dPCR analyses with the “sma” function of the “smatr” ver. 3.4.8 package in R ver. 3.6.024. Zero values were disregarded for the modelling. We employed the Gaussian Type II model because our preliminary evaluation showed higher R2 values for Type II regression models with a Gaussian distribution than for those with a logarithmic distribution in all cases. We compared the differences in the coefficient values by overlapping the 95% confidence interval (CI) ranges. More

1.
IPBES. Global Assessment Report on Biodiversity and Ecosystem Services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (eds. Brondizio, E.S., Settele, j., Díaz, s. & Ngo, H.T.) (IPBES Secretariat, Bonn, Germany, 2019).
2.
Champness, B. S., Palmer, G. C. & Fitzsimons, J. A. Bringing the city to the country: relationships between streetscape vegetation type and bird assemblages in a major regional centre. J. Urban Ecol. 5, juz018 (2019).
Article  Google Scholar 

3.
Frantzeskaki, N. et al. Nature-based solutions for urban climate change adaptation: linking the science, policy and practice communities for evidence based decision-making. BioScience 69, 455–566 (2019).
Article  Google Scholar 

4.
Willeme, A. Rotterdam drops €233 million on green spaces—and they look INCREDIBLE. DutchReview. https://dutchreview.com/cities/rotterdam-drops-233-million-on-green-spaces-and-they-look-incredible/ (18 June 2020).

5.
Díaz, S. et al. The IPBES Conceptual Framework – connecting nature and people. Curr. Opin. Environ. Sustain. 14, 1–16 (2015).
Article  Google Scholar 

6.
Elmqvist, T. et al. (eds). Urbanization, Biodiversity and Ecosystem Services: Challenges and Opportunities. A Global Assessment (Springer, Dordrecht, 2013).

7.
Ives, C. et al. Cities are hotspots for threatened species. Glob. Ecol. Biogeograph. 25, 117–126 (2016).
Article  Google Scholar 

8.
Parris, K. M. & Hazell, D. L. Biotic effects of climate change in urban environments: the case of the grey-headed flying-fox (Pteropus poliocephalus) in Melbourne, Australia. Biol. Conserv. 124, 267–276 (2005).
Article  Google Scholar 

9.
Prévot, A. C., Cheval, H., Raymond, R. & Cosquer, A. Routine experiences of nature in cities can increase personal commitment toward biodiversity conservation. Biol. Conserv. 226, 1–8 (2018).
Article  Google Scholar 

10.
ICLEI CBC. Edinburgh Process for Subnational and Local Governments on the Development of the Post 2020 Global Biodiversity Framework (ICLEI, 2020). https://cbc.iclei.org/edinburgh-process-for-subnational-and-local-governments-on-the-development-of-the-post-2020-global-biodiversity-framework/.

11.
Nilon, C. et al. Planning for the future of urban biodiversity: a global review of city-scale initiatives. BioScience 67, 332–342 (2017).
Article  Google Scholar 

12.
Garrard, G. E., Williams, N. S. G., Mata, L., Thomas, J. & Bekessy, S. A. Biodiversity sensitive urban design. Conserv. Lett. 11, e12411 (2018).
Article  Google Scholar 

13.
Bush, J. & Doyon, A. Building urban resilience with nature-based solutions: how can urban planning contribute? Cities 95, 102483 (2019).
Article  Google Scholar 

14.
Prober, S. M., Doerr, V. A. J., Broadhurst, L. M., Williams, K. J. & Dickson, F. Shifting the conservation paradigm: a synthesis of options for renovating nature under climate change. Ecol. Monog. 89, e01333 (2019).
Article  Google Scholar 

15.
Canaway, J. Unveiling the misunderstood magical mistletoes of Australia. ABC Life. https://www.abc.net.au/life/the-misunderstood-magical-mistletoes-of-australia/11505510 (20 December 2019).

16.
Burgin, S. What about biodiversity? Redefining urban sustainable management to incorporate endemic fauna with particular reference to Australia. Urban Ecosys. 19, 669–678 (2016).
Article  Google Scholar 

17.
Rigolon, A. A complex landscape of inequity in access to urban parks: a literature review. Landsc. Urban Plann. 153, 160–169 (2016).
Article  Google Scholar 

18.
Fuller, R. A. et al. Psychological benefits of greenspace increase with biodiversity. Biol. Lett. 3, 390–394 (2007).
Article  Google Scholar 

19.
Sugiyama, T., Carver, A., Koohsari, M. J. & Veitch, J. Advantages of public green spaces in enhancing population health. Landsc. Urban Plann. 178, 12–17 (2018).
Article  Google Scholar 

20.
Hockings, M. et al. COVID-19 and protected and conserved areas. Parks 26.1, 7–24 (2020).
Article  Google Scholar 

21.
Frantzeskaki, N. Seven lessons for planning nature-based solutions in cities. Environ.Sci. Pol. 93, 101–111 (2019).
Article  Google Scholar 

22.
Ahern, J. From fail-safe to safe-to-fail: sustainability and resilience in the new urban world. Landsc. Urban Plann. 100, 341–343 (2011).
Article  Google Scholar 

23.
Kabisch, N., van den Bosch, M. & Lafortezza, R. The health benefits of nature-based solutions to urbanization challenges for children and the elderly—a systematic review. Environ. Res. 159, 362–373 (2017).
CAS  Article  Google Scholar 

24.
Tengö, M. et al. Weaving knowledge systems in IPBES, CBD and beyond—lessons learned for sustainability. Curr. Opin. Environ. Sustain. 26–27, 17–25 (2017).
Article  Google Scholar 

25.
Davies, C. & Lafortezza, R. Transitional path to the adoption of nature-based solutions. Land Use Pol. 80, 406–409 (2019).
Article  Google Scholar 

26.
Mumaw, L. M. & Bekessy, S. A. Wildlife gardening for collaborative public–private biodiversity conservation. Austral. J. Environ. Manage. 24, 242–246 (2017).
Article  Google Scholar 

27.
Eilam, E. & Garrard, G. E. Perception of space among children studying their local grasslands: examining attitudes and behavioural intentions. Sustainability 9, 1660 (2017).
Article  Google Scholar 

28.
ICLEI. 6th Global Biodiversity Summit of Local and Subnational Governments. Event Report (2018) https://cbc.iclei.org/wp-content/uploads/2019/10/Egypt-Summit-EVENT-REPORT-FINAL-digital_compressed.pdf. More

Observed forest disturbances
We focused on the vulnerability of European forests to three major natural disturbances: forest fires, windthrows and insect outbreaks (bark beetles, defoliators and sucking insects). In order to identify/calibrate/validate vulnerability models (details on model development in the following sections) we used a large number of records of forest disturbances collected over the 2000–2017 period (Supplementary Fig. 1, step1). Fires were retrieved from the European Forest Fire Information System (EFFIS, https://effis.jrc.ec.europa.eu/) and count 15,818 records. Windthrows were acquired from the European Forest Windthrow dataset62 (FORWIND, https://doi.org/10.6084/m9.figshare.9555008) with 89,743 records. Insect outbreaks were retrieved from the National Insect and Disease Survey (IDS, http://foresthealth.fs.usda.gov) database of the United States Department of Agriculture (USDA) which includes 50,777 records. Each disturbance record is represented by a vector feature describing the spatial delineation of the damaged forest patch obtained by visual photointerpretation of aerial and satellite imagery or field surveys.
Even if the study focuses on Europe, for insect diseases we used the IDS-USDA database due to the lack of an analogous monitoring system and related dataset for Europe. Therefore, the models of vulnerability to insect outbreaks were identified/calibrated/validated on US data and then applied in predictive mode to Europe (see following sections for details). To assure the transferability of such models, we developed models for functional groups instead of working on species-specific models. For this purpose, we classified records based on functional groups of the pest (bark beetles, defoliators and sucking insects) and on the PFT of the host tree species. Records were considered if the host plant belonged to the following PFTs: broadleaved deciduous, broadleaved evergreen, needle leaf deciduous and needle leaf evergreen.
Reconstruction of annual biomass time series
In order to evaluate the biomass loss expected given a disturbance event occurs, multi-temporal information of biomass is required. However, there is still no single technology for direct and continuous monitoring of such variable in time. In order to reconstruct the temporal variations in biomass over the 2000–2017 period we integrated a static 100-m above ground biomass map acquired for the year 2010 from multiple Earth Observation systems63 with forest cover changes derived from the Global Forest Change (GFC) maps recorded at 30-m spatial resolution from Landsat imagery21. The GFC maps include three major layers: “2000 Tree Cover”, “Forest Cover Loss” and “Forest Cover Gain”. “2000 Tree Cover” (TC2000) is a global map of tree canopy cover (expressed in percentage) for the year 2000. “Forest Cover Loss” is defined as the complete removal of tree-cover canopy at the Landsat pixel scale (natural or human-driven) and is reported annually. “Forest Cover Gain” reflects a non-forest to forest change and refers to the period 2000–2012 as unique feature without reporting the timing of the gain.
The data integration approach built a on the assumption that changes in biomass are fully conditioned by the changes in tree cover. First, we quantified the percentage of tree cover in 2010 (TC2010) by masking out all pixels where forest loss occurred over the 2000–2010 period from the TC2000 map.
Then, in order to characterize to what extent an increase or decrease in tree cover may affect biomass, we quantified the density of biomass per percentage of tree cover lost (ρloss) and gained (ρgain) as follows:

$$rho _{mathrm{{loss}}} = frac{{B_{2010}}}{{{mathrm{{TC}}}_{2010,{mathrm{{loss}}}}}},$$
(1)

$$rho _{mathrm{{gain}}} = frac{{B_{2010}}}{{{mathrm{{TC}}}_{2010,{mathrm{{gain}}}}}},$$
(2)

where (B_{2010}) is the static biomass map available for the year 2010 (ref. 63). ({mathrm{{TC}}}_{2010,{mathrm{{loss}}}}) is the ({mathrm{{TC}}}_{2010}) masked over the pixels where there has been a forest loss during the 2011–2017 period. This filtering provides a picture of forests that were intact in 2010 but removed since then. Similarly, ({mathrm{{TC}}}_{2010,{mathrm{{gain}}}}) is the ({mathrm{{TC}}}_{2010}) masked over the pixels where there has been a forest gain and identifies the reforested and afforested areas. Since the map of forest gain is a binary map referring to the year 2012, forest gain pixels lack any information on their tree cover as their value in 2000 is zero. We therefore associated to forest gain pixels the maximum of tree cover percentage computed in a moving window with a radius of 2.5 km. This value represents the maximum potential tree cover in the local environmental conditions and refers to the whole 2000–2012 period (({mathrm{{TC}}}_{2012,{mathrm{{gain}}}})). Then, we assumed that forest gain proceeds at a constant rate over time and that the associated tree cover thus grows linearly:

$$frac{{{mathrm{{TC}}}_{2010,{mathrm{{gain}}}}}}{{left( {2010 – 2000} right)}} = frac{{{mathrm{{TC}}}_{2012,{mathrm{{gain}}}}}}{{left( {2012 – 2000} right)}} to {mathrm{{TC}}}_{2010,{mathrm{{gain}}}} = 0.83 cdot {mathrm{{TC}}}_{2012,{mathrm{{gain}}}},$$
(3)

Both ({mathrm{{TC}}}_{2010,{mathrm{{loss}}}}) and ({mathrm{{TC}}}_{2010,{mathrm{{gain}}}}) were resampled to the (B_{2010}) spatial resolution (100 m). Supplementary Figure 13 shows the frequency distribution of (rho _{mathrm{{loss}}}) and (rho _{{mathrm{{gain}}}}) over a test area in Southern Finland. As expected, the density of biomass associated with forest losses is higher than that one associated to forest gain. Indeed, biomass of new forest plantations is generally lower than the biomass of an old one (e.g. a forest that is typically harvested).
The obtained maps of (rho _{{mathrm{{loss}}}}) and (rho _{{mathrm{{gain}}}}) in Eqs. (1) and (2) refer to sparse and isolated pixels where there have been forest gain or loss. To obtain continuous fields, such density values were spatialized by computing their median over a 0.1° grid. Annual maps of biomass were finally obtained at 100 m spatial resolution as follows:

$$B_t = B_{2010} + alpha cdot rho _{{mathrm{{loss}}}} cdot {mathrm{{TC}}}_{t,{mathrm{{loss}}}} – rho _{{mathrm{{gain}}}} cdot {mathrm{{TC}}}_{t,{mathrm{{gain}}}} cdot frac{{left( {2010 – t} right)}}{{10}},$$
(4)

where t is the year (over the 2000–2017 period) and α takes the value of +1 for t  5% were selected (Supplementary Fig. 1, step 3). In the case of windthrows, we noted that maximum wind speeds retrieved from 0.5° spatial resolution of reanalysis data may largely underestimate effective maximum winds. This was particularly evident for tornado events, given their limited spatial extents compared to the grid cell, and the storm event Klaus that occurred in 2009 and for which we noticed an underestimation of the effective wind speed of the 78% (retrieved ~12 ms−1 instead of observed maximum wind speed of 55 ms−1 (ref. 67)). Therefore, such events were excluded from our analysis.
Possible missing data in the environmental variables were corrected by the median value of the variable-specific distributions (Supplementary Fig. 1, step 4). Potential effects of spatial dependence structure in the observational datasets were reduced by resampling ({mathrm{{BL}}}_{{mathrm{{rel}}}}), F, C and L along the gradients of the three principal components (PC) derived from the initial set of predictors. To this aim, we used 20 bins of equal intervals for each PC dimension spanning the full range of values. The resampling procedure was stratified by splitting the records in training and testing sets. For each year between 2000 and 2017, we randomly extracted 60% of the records. The extracted subset (({mathrm{{BL}}}_{{mathrm{{rel}}}}), F, C and L) was then binned in the PC space using the average as aggregation metric weighted by the areal extents of each disturbance record. The remaining 40% of records were similarly processed and used as a separate validation set (Supplementary Fig. 1, steps 5–7). The cover fraction of each PFT was resampled using the same approach and renormalized within each bin. Only bins with at least three records were retained for model development.
The resampled training and testing sets were used to calibrate and validate an “approximate” RF model using the full set of variables (A) as predictors initially identified based on literature review (Supplementary Fig. 1, step 8 and Supplementary Table 1). With the RF algorithm importance scores for each environmental variable can be calculated31. These scores reflect how important each covariate is in determining the fitted values of relative biomass loss. The RF implemented here uses 500 regression trees, whose depth and number of predictors to sample at each node were identified using Bayesian optimization. To reduce potential redundancy effects across predictors and facilitate the interpretability of results, we implemented a feature selection procedure. Based on the “approximate” RF model the importance of each predictor was quantified. We then computed the Spearman correlation between each pair of predictors and when it exceeded 0.8, the predictor with the lower variable importance was excluded (Supplementary Fig. 1, step 9 and Supplementary Table 1). The remaining predictors (I) were then used for a second set of RF runs, in which we iteratively evaluated RF performance on a reduced set of predictors, excluding in each new run the less important variable computed on the new reduced set of features. The set of predictors which maximizes the R2 was finally selected (Q hereafter for short) (Supplementary Fig. 1, step 10 and Supplementary Table 1). The implemented iterative feature selection procedure identifies a reasonable compromise between computing cost and model performance. The general equation describing the vulnerability is as follows:

$${mathrm{{BL}}}_{{mathrm{{rel}}}} = vleft( {{Q}} right),$$
(6)

where v is the vulnerability model implemented in the RF regression algorithm, and describes the relative biomass losses as a function of a selected Q set of environmental variables.
Such automatic feature selection process was complemented with visual interpretation of the PDPs68 based on the RF algorithm. PDP is used to visualize the relationship between explanatory covariates (environmental predictors) and ({mathrm{{BL}}}_{{mathrm{{rel}}}}), independent of other covariates (Supplementary Figs. 2–4). PDP results were analysed in combination with a detailed study of the literature and allowed us to understand and interpret the response functions to natural disturbances (see details in the main text and Fig. 2). Consistency of PDPs at the boundaries of the observational ranges was carefully checked to reduce possible artefacts generated when the models are used to extrapolate outside the range of training conditions.
Vulnerability models were further refined by retrieving v functions separately for each PFT. For PFT-specific vulnerability models, only resampled records in the PC space with a cover fraction >5% were retained and used for the model development (Supplementary Fig. 1, step 11). Model performances were ultimately evaluated on the testing set in terms of coefficient of determination (R2), root mean square error (RMSE), percent bias (PBIAS)69 and RE.
Regarding the insect-related disturbance, we initially implemented specific RF models for different insect groups (bark beetles, sucking insect and defoliators). However, due to the limited sample size of the first two groups, RF was not able to represent their effects on biomass losses reliably. We therefore opted to merge all three groups in a unique insect disturbance class (hereafter referred as insect outbreaks). We recognize that different ecological processes may characterize each insect group and therefore the use of a unique insect class may potentially mask some distinctive features. The resulting vulnerability models can therefore identify only drivers and patterns common to all groups (e.g., susceptibility to temperature anomalies70,71).
Interacting processes
The co-occurrence of multi-dimensional environmental factors resulting from the combination of interacting physical processes (compound events) may amplify or dampen ecosystem responses29. Tree-based models consider all variables together in the model and account for nonlinear feature interactions in the final model31,68. The inherent ability of RF models to detect interacting variables allows avoiding the prescription of specific relations between variables based on “a priori” knowledge—as for instance required in parametric regression frameworks—by letting the model learn automatically these relations from data.
In order to detect feature interactions and assess their strength in the developed RF-based vulnerability models we computed the Friedman’s H-statistic50. Here, we derived the H-statistic to assess second-order interactions by quantifying how much of the variation of the prediction depends on two-way interactions. To speed up the computation, we sampled 50 equally spaced data points over the environmental gradients.
We complemented this analysis by estimating the amplification or dampening effect (({Delta}{{P}})) associated to each feature interaction. To this aim, we quantified the difference in the peak values between the response function which incorporates interacting processes (two-way partial dependences) and those ones decomposed without interactions (one-dimensional partial dependences) and expressed in terms of relative variations.
The H and ({Delta}{{P}}) metrics were computed for each pair of features, and averaged for different combinations of predictor categories (forest, climate, landscape).
Spatial and temporal patterns of vulnerability and its key drivers
The RF models were used to evaluate the vulnerability of forests annually between 1979 and 2018 for each grid cell (0.25°) of the spatial domain covering the geographic Europe (including Turkey and European Russia). To this aim, vulnerability models were used in predictive mode using as input spatial maps of predictors, preliminary resampled to the common resolution, and with results expressed in terms of potential relative biomass loss (({mathrm{{PBL}}}_{{mathrm{{rel}}}})). Estimates of ({mathrm{{PBL}}}_{{mathrm{{rel}}}}) are obtained as the average from all trees in the RF ensemble. The ongoing changes in climate features were also accounted for in our framework. Climate predictors were kept dynamic for backward RF runs, while the remaining forest and landscape features were fixed to their current values averaged over the 2009–2018 period. Doing so, we implicitly assume that the sampling of response variables and predictors is representative for the whole temporal period. However, over longer time periods (from decades to century) additional ecosystem processes may play a role, such as adaptation phenomena driven by species change and shifting biomes, which could also affect vulnerability trends. The lack of multi-temporal monitoring of most of the forest and landscape predictors hampered the integration of their dynamics in the backward RF runs.
Results of PFT-specific vulnerability models were averaged at grid-cell level with weighting based on the cover fractions of PFTs (Supplementary Fig. 1, steps 12–13). This resulted in annual maps of vulnerability to each natural disturbance. Spatial and temporal variations in vulnerability were both expressed in relative and absolute terms. Absolute biomass losses were retrieved by multiplying estimates of potential relative biomass loss by the available biomass. Therefore, vulnerability values in a given grid cell reflect the biomass (relative or absolute) that would be affected if exposed to a disturbance under its specific local and temporal environmental conditions.
Grid-cell uncertainty of predicted vulnerability values were quantified in terms of standard error (SE) derived by dividing standard deviations of the computed responses over the ensemble of the grown trees of the model by the square root of the ensemble size (Supplementary Fig. 7).
We then calculated the “current” vulnerability as the average vulnerability over the 2009–2018 period. To factor out the local dependence of the current vulnerability on each predictor we retrieved the Individual Conditional Expectation72 (ICE) for each grid cell. ICE plots show the relationship between the predicted target variable (({mathrm{{PBL}}}_{{mathrm{{rel}}}})) and one predictor variable for individual cases of the predictor dataset. In our application, an individual case is a specific combination of F, C and L data for a given grid cell. To summarize and map the ICE of each grid cell in a single number, we fitted by linear regression the partial dependence of ({mathrm{{PBL}}}_{{mathrm{{rel}}}}) versus the corresponding predictor variable and mapped the slope of this regression, hereafter referred as “local sensitivity” (Supplementary Figs. 5–7), similarly to the approach presented in ref. 30. The marginal contribution ((Z_{mathrm{{marg}}})) of each environmental category of predictors (F, C and L, hereafter referred as X for short) on the current vulnerability was derived as follows:

$$Z_{{mathrm{{marg}}},X} = 100 times frac{{mathop {sum }nolimits_{i in X} left| {s_i} right|}}{{mathop {sum }nolimits_{j in Q} left| {s_j} right|}},$$
(7)

where s represents the slope of ICE, i runs over all predictors of X, whereas j runs over all available predictors Q. Therefore (Z_{mathrm{{marg}},X}) values range between 0 (no dependence of current vulnerability on X predictors) and 100% (full dependence of current vulnerability on X predictors).
Long-term linear trends in vulnerability ((delta {mathrm{{PBL}}}_{{mathrm{{rel}}}})) were quantified over the 1979–2018 period for each grid cell and their significance evaluated by the two-sided Mann–Kendall test. In order to isolate the key determinants of the emerging trends in vulnerability, a set of factorial simulations was performed. To this aim, we estimated the vulnerability due to the temporal variations in a given k climate predictor (({mathrm{{PBL}}}_{{mathrm{{rel}}}}^k)), by applying the RF models to a data array in which the k climate variable is dynamic while all the remaining features are kept fixed to their “current” value (average value over 2009–2018). The resulting trends in vulnerability associated to the k factor (({mathrm{{PBL}}}_{{mathrm{{rel}}}}^k)) are then calculated by linear regression and subject to the Mann–Kendall test.
Spatial and temporal patterns were visualized at grid-point scale and averaged over geographic macro-regions (Supplementary Fig. 14 and Supplementary Tables 2 and 3). Zonal statistics were obtained by averaging grid-cell results weighted by their forest areal extent. Forests with cover fraction lower than 0.1 were excluded from the analyses. Uncertainty in spatial averages were based on the 95% bootstrap confidence interval computed with 100 bootstrap samples.
In order to derive statistics minimally affected by potential extrapolation errors of the RF models, we replicated the aforementioned analyses by excluding areas outside the observational ranges of climatological temperature and precipitation (Supplementary Fig. 8).
Combining forest vulnerability to multiple natural disturbances
To quantify the total vulnerability to multiple disturbances we defined the OVI, similarly to the multi-hazard index developed in ref. 73. We assumed that the considered disturbances are independent and mutually non-exclusive and the potential biomass loss of single disturbances is spread homogeneously within each grid cell. From the inclusion-exclusion principle of combinatorics the potential biomass loss associated to the OVI can be expressed for a given year as follows:

$${mathrm{{PBL}}}_{{mathrm{{rel}}}}left( {{mathrm{{OVI}}}} right) = mathop {bigcup}nolimits_{p = 1}^D {{mathrm{{PBL}}}_{{mathrm{{rel}}},p}} = mathop {sum }limits_{q = 1}^D left( {left( { – 1} right)^{q – 1} cdot mathop {sum }limits_{{G subset left{ {1, ldots ,D} right}} atop {left| G right| = q} } {mathrm{{PBL}}}_{{mathrm{{rel}}},G}} right),$$
(8)

where p refers to the disturbance-specific ({mathrm{{PBL}}}_{{mathrm{{rel}}}}), D is the number of disturbances considered, the last sum runs over all subsets G of the indices {1, …, D} containing exactly q elements, and

$${mathrm{{PBL}}}_{{mathrm{{rel}}},G}: = mathop {bigcap}nolimits_{p in I} {{mathrm{{PBL}}}_{{mathrm{{rel}}},p}} ,$$
(9)

expresses the intersection of all those ({mathrm{{PBL}}}_{{mathrm{{rel}}},p}) with index in G. Maps of current overall vulnerability and trends were ultimately analysed following the approach adopted for single disturbances.
This approach does not account for the potential reduction in exposed biomass following the occurrence of a given disturbance. Furthermore, possible amplification/dampening effects due to interacting disturbances could also occur3,74. A strong interaction effect has been documented for instance between windthrows and bark beetle disturbances. Uprooted trees are virtually defenseless breeding material supporting the build-up of beetle populations and the consequent increase in vulnerability to insect outbreaks3,59. Insect outbreaks, in turn, may potentially affect the severity of subsequent forest fires by altering the abundance of available fuel60. The magnitude of these effects varies with insect type and outbreak timing. Despite the relevance of these interactions, the lack of reference observational data of compound events hampered the integration of their effects in our modelling framework. Therefore, estimates of OVI can only partially capture the overall vulnerability resulting from multiple disturbances and should be viewed in light of these limitations.
Spatial maps of current overall vulnerability and trends in OVI were then normalized separately based on the min–max method and combined by simple multiplication into a single index, hereafter referred as space-time integrated OVI. High values of space-time integrated OVI depict forest areas that are currently susceptible to multiple disturbances and their vulnerability have experienced a substantial increase over the 1979–2018 period. The space-time integrated OVI is used to identify currently fragile ecosystems that might in the future become even more susceptible to natural disturbances.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.
Taberlet, P., Bonin, A., Zinger, L, & Coissac, E. Environmental DNA, for Biodiversity Research and Monitoring (Oxford Univ. Press, 2018).
2.
Jo, T., Arimoto, M., Murakami, H., Masuda, R. & Minamoto, T. Particle size distribution of environmental DNA from the nuclei of marine fish. Environ. Sci. Technol. 53, 9947–9956 (2019).
CAS  PubMed  Article  Google Scholar 

3.
Wilcox, T. M., McKelvey, K. S., Young, M. K., Lowe, W. H. & Schwartz, M. K. Environmental DNA particle size distribution from Brook Trout (Salvelinus fontinalis). Conserv. Genet. Resour. 7, 639–641 (2015).
Article  Google Scholar 

4.
Thomsen, P. F. & Willerslev, E. Environmental DNA – an emerging tool in conservation for monitoring past and present biodiversity. Biol. Conserv. 183, 4–18 (2015).
Article  Google Scholar 

5.
Seymour, M. et al. Executing multi-taxa eDNA ecological assessment via traditional metrics and interactive networks. Sci. Total Environ. 729, 138801 (2020).
CAS  PubMed  Article  Google Scholar 

6.
Jarman, S. N., Berry, O. & Bunce, M. The value of environmental DNA biobanking for long-term biomonitoring. Nat. Ecol. Evol. 2, 1192–1193 (2018).
PubMed  Article  Google Scholar 

7.
Jeunen, G.-J. et al. Species-level biodiversity assessment using marine environmental DNA metabarcoding requires protocol optimization and standardization. Ecol. Evol. 9, 1323–1335 (2019).
PubMed  PubMed Central  Article  Google Scholar 

8.
Turner, C. R. et al. Particle size distribution and optimal capture of aqueous microbial eDNA. Methods Ecol. Evol. 5, 676–684 (2014).
Article  Google Scholar 

9.
Koziol, A. et al. Environmental DNA metabarcoding studies are critically affected by substrate selection. Mol. Ecol. Resour. 19, 366–376 (2019).
CAS  PubMed  Article  Google Scholar 

10.
Tsuji, S., Takahara, T., Doi, H., Shibata, N. & Yamanaka, H. The detection of aquatic macroorganisms using environmental DNA analysis – a review of methods for collection, extraction, and detection. Environ. DNA 1, 99–108 (2019).
Article  Google Scholar 

11.
Shu, L., Ludwig, A. & Peng, Z. Standards for methods utilizing environmental DNA for detection of fish species. Genes 11, 296 (2020).
CAS  PubMed Central  Article  PubMed  Google Scholar 

12.
Deiner, K., Walser, J.-C., Mächler, E. & Altermatt, F. Choice of capture and extraction methods affect detection of freshwater biodiversity from environmental DNA. Biol. Conserv. 183, 53–63 (2015).
Article  Google Scholar 

13.
Jeunen, G.-J. et al. Environmental DNA (eDNA) metabarcoding reveals strong discrimination among diverse marine habitats connected by water movement. Mol. Ecol. Resour. 19, 426–438 (2019).
CAS  PubMed  Article  Google Scholar 

14.
Thomas, A. C., Howard, J., Nguyen, P. L., Seimon, T. A. & Goldberg, C. S. ANDeTM: a fully integrated environmental DNA sampling system. Methods Ecol. Evol. 9, 1379–1385 (2018).
Article  Google Scholar 

15.
Schumer, G. et al. Utilizing environmental DNA for fish eradication effectiveness monitoring in streams. Biol. Invasions 21, 3415–3426 (2019).
Article  Google Scholar 

16.
Zinger, L. et al. DNA metabarcoding – need for robust experimental designs to draw sound ecological conclusions. Mol. Ecol. 28, 1857–1862 (2019).
PubMed  Article  Google Scholar 

17.
Bessey, C. et al. Maximizing fish detection with eDNA metabarcoding. Environ. DNA 2, 493–504, https://doi.org/10.1002/edn3.74 (2020).
Article  Google Scholar 

18.
Harrison, J. B., Sunday, J. M. & Rogers, S. M. Predicting the fate of eDNA in the environment and implications for studying biodiversity. Proc. R. Soc. Ser. B 286, 20191409 (2019).
CAS  Article  Google Scholar 

19.
Seymour, M. et al. Acidity promotes degradation of multi-species environmental DNA in lotic mesocosms. Commun. Biol. 1, https://doi.org/10.1038/s42003-017-0005-3 (2018).

20.
Deiner, K. & Altermatt, F. Transport distance of invertebrate environmental DNA in a natural river. PLoS ONE 9, e88786 (2014).
PubMed  PubMed Central  Article  CAS  Google Scholar 

21.
Mächler, E., Deiner, K., Spahn, F. & Altermatt, F. Fishing in the water: effect of sampled water volume on environmental DNA-based detection of macroinvertebrates. Environ. Sci. Technol. 50, 305–312 (2016).
PubMed  Article  CAS  Google Scholar 

22.
Hanfling, B. et al. Environmental DNA metabarcoding of lake fish communities reflects long-term data from established survey methods. Mol. Ecol. 25, 3101–3119 (2016).
PubMed  Article  CAS  Google Scholar 

23.
Cantera, I. et al. Optimizing environmental DNA sampling effort for fish inventories in tropical streams and rivers. Sci. Rep. 9, 3085 (2019).
PubMed  PubMed Central  Article  CAS  Google Scholar 

24.
McQuillan, J. S. & Robidart, J. C. Molecular-biological sensing in aquatic environments: recent developments and emerging capabilities. Curr. Opin. Biotechnol. 45, 43–50 (2017).
CAS  PubMed  Article  Google Scholar 

25.
Schabacker, J. C. et al. Increased eDNA detection sensitivity using a novel high-volume water sampling method. Environ. DNA 2, 244–251 (2020).
Article  Google Scholar 

26.
Mariani, S., Baillie, C., Colosimo, G. & Riesgo, A. Sponges as natural environmental DNA samples. Curr. Biol. 29, R395–R402 (2019).
Article  CAS  Google Scholar 

27.
Keesing, J., Webber, B.L. & Hardiman, L. Ashmore Reef Marine Park Environmental Assessment. Final report to director of National Park (2020).

28.
Kirtane, A., Atkinson, J. D. & Sassoubre, L. Design and validation of passive environmental DNA samplers using granular activated carbon and montmorillonite clay. Environ. Sci. Technol. https://doi.org/10.1021/acs.est.0c01863 (2020).
Article  PubMed  Google Scholar 

29.
Taberlet, P., Coissac, E., Hajibabaei, M. & Rieseberg, L. H. Environmental DNA. Mol. Ecol. Resour. 21, 1789–1793 (2012).
CAS  Article  Google Scholar 

30.
Fonseca, V. G. Pitfalls in relative abundance estimation using eDNA metabarcoding. Mol. Ecol. Resour. 18, 923–926 (2018).
CAS  Article  Google Scholar 

31.
Lamb, P. D. et al. How quantitative is metabarcoding: a meta-analytical approach. Mol. Ecol. 28, 420–430 (2019).
PubMed  Article  Google Scholar 

32.
Derocles, S. A. P. et al. Biomonitoring for the 21st century: integrating next-generation sequencing into ecological network analysis. Adv. Ecol. Res. 58, 1–62 (2018).
Article  Google Scholar 

33.
Prosser, J. I. Replicate or lie. Environ. Microbiol. 12, 1806–1810 (2010).
CAS  PubMed  Article  Google Scholar 

34.
MacKenzie, D. I. What are the issues with presence-absence data for wildlife managers? J. Wildl. Manag. 69, 849–860 (2005).
Article  Google Scholar 

35.
Liang, Z. & Keeley, A. Filtration recovery of extracellular DNA from environmental water samples. Environ. Sci. Technol. 47, 9324–9331 (2013).
CAS  PubMed  Article  Google Scholar 

36.
Renshaw, M. A., Olds, B. P., Jerde, C. L., McVeigh, M. M. & Lodge, D. M. The room temperature preservation of filtered environmental DNA samples and assimilation into a phenol-chloroform-isoamyl alcohol DNA extraction. Mol. Ecol. Resour. 15, 168–176 (2015).
CAS  PubMed  Article  Google Scholar 

37.
Eichmiller, J. J., Miller, L. M. & Sorensen, P. W. Optimizing techniques to capture and extract environmental DNA for detection and quantification of fish. Mol. Ecol. Resour. 16, 56–68 (2016).
CAS  PubMed  Article  Google Scholar 

38.
Majaneva, M. et al. Environmental DNA filtration techniques affect recovered biodiversity. Sci. Rep. 8, 4682 (2018).
PubMed  PubMed Central  Article  CAS  Google Scholar 

39.
Stier, A. C., Bolker, B. M. & Osenberg, C. W. Using rarefaction to isolate the effects of patch size and sampling effort on beta diversity. Ecosphere 7, e01612 (2016).
Article  Google Scholar 

40.
Yates, M. C., Fraser, D. J. & Derry, A. M. Meta-analysis supports further refinement of eDNA for monitoring aquatic species-specific abundance in nature. Environ. DNA 1, 5–13 (2019).
Article  Google Scholar 

41.
Strickland, G. J. & Roberts, J. H. Utility of eDNA and occupancy models for monitoring an endangered fish across diverse riverine habitats. Hydrobiologia 826, 129–144 (2019).
CAS  Article  Google Scholar 

42.
Deagle, B. E. et al. Counting with DNA metabarcoding studies: how should we convert sequence reads to dietary data? Mol. Ecol. 28, 391–406 (2019).

43.
Shogren, A. J. et al. Controls on eDNA movement in streams: transport, retention, and resuspension. Sci. Rep. 7, 5065 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

44.
Berry, T. E. et al. DNA metabarcoding for diet analysis and biodiversity: a case study using the endangered Australian sea lion (Neophoca cinerea). Ecol. Evol. 7, 5435–5453 (2017).
PubMed  PubMed Central  Article  Google Scholar 

45.
Deagle, B. E. et al. Studying seabird diet through genetic analysis of faeces: a case study on Macaroni penguins (Eudyptes chrysolophus). PLoS ONE 2, e831 (2007).
PubMed  PubMed Central  Article  CAS  Google Scholar 

46.
Murray, D. C., Coghlan, M. L. & Bunce, M. From benchtop to desktop: important considerations when designing amplicon sequencing workflows. PLoS ONE 10, e0124671 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

47.
Benson, D. A. et al. GenBank. Nucleic Acids Res. 42, D32–D37 (2014).
CAS  PubMed  Article  Google Scholar 

48.
Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).
CAS  PubMed  PubMed Central  Article  Google Scholar 

49.
Paradis, E. APE 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528 (2019).
CAS  Article  PubMed  Google Scholar 

50.
Baselga, A. & Orme, C. D. L. Betapart: an R package for the study of beta diversity. Methods Ecol. Evol. 3, 808–812 (2012).
Article  Google Scholar 

51.
Dixon, P. VEGAN, a package of R functions for community ecology. J. Veg. Sci. 14, 927–930 (2003).
Article  Google Scholar 

52.
Herve, M. RVAideMemoire, testing and plotting procedures for biostatistics. https://cran.r-project.org/web/packages/RVAideMemoire/index.html (2018). More