Ecology

1.
Han, J. & Zhang, Y. Land policy and land engineering. Land Use Policy 40, 64–68 (2014).
Article  Google Scholar 
2.
Liu, Y. S. Introduction to land use and rural sustainability in China. Land Use Policy. 74, 1–4 (2018).
Article  Google Scholar 

3.
Liu, Y. S., Li, J. T. & Yang, Y. Y. Strategic adjustment of land use policy under the economic transformation. Land Use Policy 74, 5–14 (2018).
Article  Google Scholar 

4.
Yi, L. et al. Spatial-temporal change of major reserve resources of cultivated land in China in recent 30 years. Trans. CSAE (in Chinese) 29, 1–12 (2013).
Google Scholar 

5.
Zhang, Y., Li, Z. B. & Dong, Q. G. Analysis of water and nutrients leakage on barren gravel land with different structure of topsoil. J. Soil Water Conserv. 32, 162–166 (2018).
Google Scholar 

6.
Wu, X. L. et al. Improvement of new farmland soil in loess area under different fertilization treatments. Sci. Soil Water Conserv. 13, 99–104 (2015).
Google Scholar 

7.
Hu, S. G., Wang, J. & Wang, Z. Q. Characterization of quantity-quality grade conversion coefficient of supplementary cultivated land in Jianghan Plain, China. Proc. SPIE Int. Soc. Opt. Eng. 7841, 1–9 (2009).
Google Scholar 

8.
Liu, Y. S. et al. Progress of research on urbanrural transformation and rural development in China in the past decade and future prospects. J. Geogr. Sci. 26, 1117–1132 (2016).
ADS  Article  Google Scholar 

9.
Sun, Z. & Han, J. Effect of soft rock amendment on soil hydraulic parameters and crop performance in Mu Us Sandy Land, China. Field Crops Res. 222, 85–93 (2018).
Article  Google Scholar 

10.
Yang, J., Liu, L., Sun, C. M. & Liu, X. F. Soil components and fertility improvement of added cultivated land. Trans. CSAE 7, 102–105 (2008).
Google Scholar 

11.
Li, Y., Liu, Y., Long, H. & Cui, W. Community-based rural residential land consolidation and allocation can help to revitalize hollowed villages in traditional agricultural areas of China: Evidence from Dancheng County, Henan Province. Land Use Policy 39, 188–198 (2014).
CAS  Article  Google Scholar 

12.
Sağlam, M. & Selvi, K. C. Effects of different tillage managements on soil physical quality in a clayey soil. Environ. Monit. Assess. 187, 1–12 (2015).
Article  CAS  Google Scholar 

13.
Márcio, R., Douglas, L. K. & Moorman, T. B. Tillage intensity effects on soil structure indicators—A US meta-analysis. Sustainability 12, 1–17 (2020).
Article  Google Scholar 

14.
Blanco-Moure, N., Moret-Fernández, D. & López, M. V. Dynamics of aggregate destabilization by water in soils under long-term conservation tillage in semiarid Spain. CATENA 99, 34–41 (2012).
CAS  Article  Google Scholar 

15.
Nunes, M. R., van Es, H. M., Schindelbeck, R., Ristow, A. J. & Ryan, M. No-till and cropping system diversification improve soil health and crop yield. Geoderma 328, 30–43 (2018).
ADS  CAS  Article  Google Scholar 

16.
Steele, M. K., Coale, F. J. & Hill, R. L. Winter annual cover crop impacts on no-till soil physical properties and organic matter. Soil Sci. Soc. Am. J. 76, 2164–2173 (2012).
ADS  CAS  Article  Google Scholar 

17.
Laura, T., Jacynthe, D. R. & Caron, J. Short-term improvement in soil physical properties of cultivated histosols through deep-rooted crop rotation and subsoiling. Agron. J. 111, 2084–2096 (2019).
Article  CAS  Google Scholar 

18.
Sun, X. F. et al. Subsoiling practices change root distribution and increase post-anthesis dry matter accumulation and yield in summer maize. PLoS ONE 12, 1–18 (2017).
Google Scholar 

19.
Huang, M., Liang, T., Wang, L. & Zhou, C. Effects of no-tillage systems on soil physical properties and carbon sequestration under long-term wheat–maize double cropping system. CATENA 128, 195–202 (2015).
CAS  Article  Google Scholar 

20.
Martínez, I. et al. Two decades of no-till in the Oberacker long-term field experiment: Part II. Soil porosity and gas transport parameters. Soil Tillage Res. 163, 130–140 (2016).
Article  Google Scholar 

21.
Suzuki, L. E. A. S., Reichert, J. M. & Reinert, D. J. Degree of compactness, soil physical properties and yield of soybean in six soils under no-tillage. Soil Res. 51, 311–321 (2013).
Article  Google Scholar 

22.
López-Fando, C. & Pardo, M. T. Changes in soil chemical characteristics with different tillage practices in a semi-arid environment. Soil Tillage Res. 104, 278–284 (2009).
Article  Google Scholar 

23.
Song, K. et al. Effects of tillage and straw return on water-stable aggregates, carbon stabilization and crop yield in an estuarine alluvial soil. Rep. U.K. 9, 4586 (2019).
Google Scholar 

24.
Zhou, H., Li, B. G. & Lu, Y. Z. Micromorphological analysis of soil structure under no tillage management in the black soil zone of northeast china. J. Mt. Sci. 6, 173–180 (2009).
Article  Google Scholar 

25.
Chang, I., Im, J., Prasidhi, A. K. & Cho, G. C. Effects of Xanthan gum biopolymer on soil strengthening. Constr. Build. Mater. 74, 65–72 (2015).
Article  Google Scholar 

26.
Jiang, M., Zhang, F., Hu, H., Cui, Y. & Peng, J. Structural characterization of natural loess and remolded loess under triaxial tests. Eng. Geol. 181, 249–260 (2014).
Article  Google Scholar 

27.
Lin, B. & Cerato, A. B. Applications of SEM and ESEM in microstructural investigation of shale-weathered expansive soils along swelling-shrinkage cycles. Ecol. Eng. 177, 66–74 (2014).
Google Scholar 

28.
Mikha, M. M., Benjamin, J. G., Vigil, M. F. & Nielsen, C. Cropping intensity impacts on soil aggregation and carbon sequestration in the central Great Plains. Soil Sci. Soc. Am. J. 74, 1712–1719 (2010).
ADS  CAS  Article  Google Scholar 

29.
Sodhi, G. P., Beri, V. & Benbi, D. K. Soil aggregation and distribution of carbon and nitrogen in different fractions under long-term application of compost in rice-wheat system. Soil Tillage Res. 103, 412–418 (2009).
Article  Google Scholar 

30.
Mikha, M. M., Vigil, M. F. & Benjamin, J. G. Long-term tillage impacts on soil aggregation and carbon dynamics under wheat-fallow in the central Great Plains. Soil Sci. Soc. Am. J. 77, 594–605 (2013).
ADS  CAS  Article  Google Scholar 

31.
Matson, P., Lohse, K. A. & Hall, S. J. The globalization of nitrogen deposition: Consequences for terrestrial ecosystems. Ambio 31, 113–119 (2002).
PubMed  Article  Google Scholar 

32.
Li, X. G. et al. Nitrogen fertilization decreases the decomposition of soil organic matter and plant residues in planted soils. Soil Biol. Biochem. 112, 47–55 (2017).
CAS  Article  Google Scholar 

33.
Hou, X. Q. et al. Effects of rotational tillage practices on soil properties, winter wheat yields and water-use efficiency in semi-arid areas of north-west China. Field Crops Res. 129, 7–13 (2012).
Article  Google Scholar 

34.
Chen, H. Q. et al. Effects of 11 years of conservation tillage on soil organic matter fractions in wheat monoculture in Loess Plateau of China. Soil Tillage Res. 106, 85–94 (2009).
Article  Google Scholar 

35.
Hernanz, J. L., Sánchez-Girón, V. & Navarrete, L. Soil carbon sequestration and stratification in a cereal/leguminous crop rotation with three tillage systems in semiarid conditions. Agric. Ecosyst. Environ. 133, 114–122 (2009).
CAS  Article  Google Scholar 

36.
Mazzoncini, M. et al. Soil carbon and nitrogen changes after 28 years of no-tillage management under Mediterranean conditions. Eur. J. Agron. 77, 156–165 (2016).
CAS  Article  Google Scholar 

37.
Sombrero, A. & Benito, A. D. Carbon accumulation in soil. Ten-year study of conservation tillage and crop rotation in a semi-arid area of castile-leon, spain. Soil Tillage Res. 107, 64–70 (2010).
Article  Google Scholar 

38.
Bottinelli, N. et al. Tillage and fertilization practices affect soil aggregate stability in a Humic Cambisol of Northwest France. Soil Tillage Res. 170, 14–17 (2017).
Article  Google Scholar 

39.
Du, Z., Ren, T., Hu, C., Zhang, Q. & Blanco-Canqui, H. Soil aggregate stability and aggregate-associated carbon under different tillage systems in the north China plain. J. Integr. Agric. 12, 2114–2123 (2013).
Article  Google Scholar 

40.
Bhattacharyya, R., Prakash, V., Kundu, S., Srivastva, A. K. & Gupta, H. S. Soil aggregation and organic matter in a sandy clay loam soil of the Indian Himalayas under different tillage and crop regimes. Agr. Ecosyst. Environ. 132, 126–134 (2009).
Article  Google Scholar 

41.
Tripathi, R. et al. Soil aggregation and distribution of carbon and nitrogen in different fractions after 41 years long-term fertilizer experiment in tropical rice-rice system. Geoderma 213, 280–286 (2014).
ADS  CAS  Article  Google Scholar 

42.
Meng, Q. et al. Distribution of carbon and nitrogen in water-stable aggregates and soil stability under long-term manure application in solonetzic soils of the Songnen plain, northeast China. J. Soil Sediment 14, 1041–1049 (2014).
CAS  Article  Google Scholar 

43.
Six, J., Bossuyt, H., Degryze, S. F. & Denef, K. A history of research on the link between (micro)aggregates, soil biota, and soil organic matter dynamics. Soil Tillage Res. 79, 7–31 (2004).
Article  Google Scholar 

44.
Gu, X., An, T. T., Li, S. Y., Li, H. & Wang, J. K. Effects of application of straw on organic carbon in brown soil aggregates by δ13c method. J. Soil Water Conserv. 28, 243–247 (2014).
Google Scholar 

45.
Cremeens, D. L. Guidelines for analysis and description of soil and regolith thin sections. Geoarchaeology 19, 507–509 (2004).
Article  Google Scholar 

46.
Wang, F., Tong, Y. A., Zhang, J. S., Gao, P. C. & Coffie, J. N. Effects of various organic materials on soil aggregate stability and soil microbiological properties on the Loess Plateau of China. Plant Soil Environ. 59, 162–168 (2013).
CAS  Article  Google Scholar 

47.
Blanco-Canqui, H. & Lai, R. No-tillage and soil profile carbon sequestration: An onfarm assessment. Soil Sci. Soc. Am. J. 141, 355–362 (2007).
CAS  Google Scholar 

48.
Pikul, J. L. & Zuzel, J. F. Soil crusting and water infiltration affected by long-term tillage and residue management. Soil Sci. Soc. Am. J. 58, 1524–1530 (1994).
ADS  Article  Google Scholar 

49.
Kinoshita, R., Schindelbeck, R. R. & Harold, M. Quantitative soil profile-scale assessment of the sustainability of long-term maize residue and tillage management. Soil Tillage Res. 174, 34–44 (2017).
Article  Google Scholar 

50.
Das, A. et al. Tillage and cropping sequence effect on physico-chemical and biological properties of soil in Eastern Himalayas, India. Soil Tillage Res. 180, 182–193 (2018).
Article  Google Scholar 

51.
Xiao, J. et al. Comment on “The whole-soil carbon flux in response to warming”. Science 359, 0218 (2018).
Article  CAS  Google Scholar 

52.
Awe, G. O., Reichert, J. M. & Wendroth, O. O. Temporal variability and covariance structures of soil temperature in a sugarcane field under different management practices in southern Brazil. Soil Tillage Res. 150, 93–106 (2015).
Article  Google Scholar 

53.
Blanco-Canqui, H. & Ruis, S. J. No-tillage and soil physical environment. Geoderma 326, 164–200 (2018).
ADS  Article  Google Scholar 

54.
Li, J. et al. Effect of mixing ratio and planting years on the micro structure of soft rock and sand compound soil: A scanning electron microscope study. Appl. Ecol. Environ. Res. 17, 10599–10610 (2019).
Google Scholar 

55.
Zhang, W. et al. Residue incorporation enhances the effect of subsoiling on soil structure and increases soc accumulation. J Soils Sediments 20, 3537–3547 (2020).
CAS  Article  Google Scholar 

56.
Kinoshita, R., Schindelbeck, R. R. & van Es, R. M. Quantitative soil profile-scale assessment of the sustainability of long-term maize residue and tillage management. Soil Tillage Res. 174, 34–44 (2017).
Article  Google Scholar 

57.
Sharma, P., Singh, G. & Singh, R. P. Conservation tillage and optimal water supply enhances microbial enzyme (glucosidase, urease and phosphatase) activities in fields under wheat cultivation during various nitrogen management practices. Arch. Agron. Soil Sci. 59, 911–928 (2013).
CAS  Article  Google Scholar 

58.
Daveiga, M., Horn, R., Reinert, D. & Reichert, J. Soil compressibility and penetrability of an Oxisol from southern Brazil, as affected by long-term tillage systems. Soil Tillage Res. 92, 104–113 (2007).
Article  Google Scholar 

59.
Karlen, D. & Rice, C. Soil degradation: Will humankind ever learn?. Sustainability 7, 12490–12501 (2015).
CAS  Article  Google Scholar 

60.
Rabot, E., Wiesmeier, M., Schlüter, S. & Vogel, H. J. Soil structure as an indicator of soil functions: A review. Geoderma 314, 122–137 (2018).
ADS  Article  Google Scholar 

61.
Xue, B. et al. Effects of organic carbon and iron oxides on soil aggregate stability under different tillage systems in a rice–rape cropping system. CATENA 177, 1–12 (2019).
CAS  Article  Google Scholar 

62.
Bremner, J. M. Determination of nitrogen in soil by the Kjeldahl method. J. Agric. Sci. 55, 11–33 (1960).
CAS  Article  Google Scholar 

63.
Tiessen, H. & Moir, J. O. Total and organic carbon. In: Carter, M.R. (Ed.), Soil sampling and methods of analysis. J. Environ. Qual. 38, 187–199 (1993).
Google Scholar 

64.
Nimmo, J. R. & Perkins KS. Aggregates Stability and Size Distribution. In Methods of Soil Analysis, Part4-Physical Methods, 317–328 (Soil Sci. Soc. Am. J., Inc., Madison, 2002).

65.
Karami, A., Homaee, M., Afzalinia, S., Ruhipour, H. & Basirat, S. Organic resource management: Impacts on soil aggregate stability and other soil physico-chemical properties. Agric. Ecosyst. Environ. 148, 22–28 (2012).
CAS  Article  Google Scholar 

66.
Liu, Z. et al. Long-term effects of different planting patterns on greenhouse soil micromorphological features in the North China Plain. Sci. Rep. U.K. 6, 31118 (2016).
ADS  Article  CAS  Google Scholar 

67.
Zhu, G., Shangguan, Z. & Deng, L. Soil aggregate stability and aggregate-associated carbon and nitrogen in natural restoration grassland and Chinese red pine plantation on the Loess Plateau. CATENA 149, 253–260 (2017).
CAS  Article  Google Scholar  More

Study system and species
This study was performed in Welgevonden Game Reserve (WGR), a privately owned game reserve in the Limpopo province, South Africa (24° 10′ S; 27° 45′ E to 24° 25′ S; 27° 56′ E). The reserve is located in the mountainous Waterberg region. WGR was established on former agricultural lands in the early 1980s and the main occurring vegetation types are Waterberg Mountain Bushveld and Sour Bushveld. The Waterberg region has a temperate climate, with two distinct seasons, characterized by the rainfall regime: a dry season ranging from April to September and a wet season ranging from October to March, with an average annual precipitation in WGR of 634 mm. Our study area is an enclosed breeding camp within WGR, with a size of approximately 1200 ha. Main predator species such as lion, cheetah and spotted hyena were excluded from this study area, as well as elephant and rhino.
WGR equipped 35 impala (Aepyceros melampus), 34 blue wildebeest (Connochaetes taurinus), 35 plains zebra (Equus burchellii) and 34 common eland (Taurotragus oryx) with a GPS and accelerometer sensor equipped custom made collar; an estimated 23% of the individual impalas present in the area, 48% of the eland, 40% of the wildebeest and 40% of the zebra. However, due to malfunctioning and errors made in the sensor development process, only 83 of the sensors yielded data at any point in time, thus lowering the effective density of sentinel animals. During the experimental intrusions (see below), the median number of data-yielding sensors was 47, and minimally 30. The animal movement data were recorded day and night and transmitted wirelessly in near real-time to five long-range low-power LoRa radiocommunication gateways in the study area, from where data packages were routed to an on-line data warehouse via a 3G/4G backhaul. The deployment of these sentinel animals were approved by the board and CEO of WGR as a management action and was performed in accordance with relevant guidelines and regulations (see Supplementary GPS Collaring letter).
Experimental intrusions
Between September 2017 and March 2018, WGR employees performed experimental intrusions (lasting ca. 2 h) on foot and by car through the study area, at varying locations and movement routes through the study area, independent from the locations of the sentinel animals. The movement of the intrusions were tracked by GPS, and the relevant metadata for each intrusion recorded (mode of transport, group size, start time, end time). The intrusions were distributed in a stratified way over the mornings, middays and afternoons (with time slots relative to specific solar positions: sunrise, solar noon and sunset). Furthermore, the intrusions were temporally spread in such a way to avoid a disturbance overflow for the sentinel animals, by performing a maximum of five experiments per week and a maximum of two experiments per day (and then only with one intrusion in the morning and one in the afternoon).
Data gathering
The animal sensors gathered location data via GPS and overall dynamic body accelerations47 (ODBA) via a tri-axial accelerometer (range ± 2 g; sampling frequency 100 Hz, down-sampled to 10 Hz prior to analysis). The GPS was scheduled to record spatial position at irregular intervals depending on the level of activity as gauged by ODBA. All sensors were scheduled to record locations every 15 min in the absence of sufficient activity (given that successive fixes were further than 5 m apart, else a geofence was applied and the new coordinate was omitted to save bandwidth and battery power, thereby assuming that the animal still was at its previous location). The GPS fix rate was increased up to 2- or 10-min intervals (depending on two different sensor settings) when ODBA indicated sufficient activity (after checking for the geofence). ODBA data were sampled continuously and summarized per 15 s window in a mean, maximum and variance value.
The experimentally intruding groups were outfitted with handheld GPS devices that recorded their location every 5 s and these groups logged and timestamped all their pre-defined activities and metadata on a tablet using CyberTracker48 during their intrusion. Most cars traveling through the study area were tracked by GPS as well to filter the animal data for disturbances by cars unrelated to the experimental intrusions.
Weather data (temperature, radiation, precipitation and wind) in the study area were recorded on a 3-min resolution with a weather station in the north of the study area. We assumed the 1200 ha study area to be sufficiently small to assume the weather station data to be representable for the prevailing weather conditions throughout the study area. GIS data of the study area (summarized in Supplementary Table S1) consisted of information on topography, infrastructure (e.g., fences, roads, powerlines, etc.) and vegetation cover (supervised classification of 25 cm resolution aerial imagery into four classes: trees, herbaceous/grass, sand/soil and other/built-up area).
All further data processing and algorithm development was done in the software R3.5.049.
Data pre-processing
To link the animal location data with the intrusion location data, as well as to correct for the substantial level of positional noise present in the animal location data, we modelled the animal location data to regular 1-min resolution trajectories using the following five steps. First, we filtered out large obvious errors (e.g., obvious outliers and irregularities such as locations far outside the study area) from the data. Second, we corrected systematic medium-scale outliers: ‘spikes’ that occurred due to positional outliers. Such spike-like outliers were visible during sensor testing while following known straight-line trajectories along an airstrip, thereby confirming that these spike-like geometries most likely resulted from positional error rather than true animal movement. Points were classified as anomalous spike points when (a) the displacement to and away from this point was high ( > 500 m), (b) when the distance between the locations before and after this point was small, and (c) when the turning angle at this point approached 180 degrees. Therefore, we corrected the locations that were classified as spike-like anomalies by shifting them closer to the straight line between the neighboring points. The extent of this shift was set relative to the degree of spikiness of the points (the spikier the pattern, the larger the shift towards the midpoint of the adjacent coordinates). Third, after filtering and correcting the original locations we smoothed the timeseries of x/y coordinates at each original timepoint with a Kalman smoother using a dynamic linear model. Fourth, we linearly interpolated the locations to a 10 s resolution based on ODBA, where we considered the animal to be stationary between multiple timepoints if the accelerometer signal suggested the animal was not moving. Fifth, we fitted an X-spline through the data, where we gave the linearly ODBA-interpolated locations a smaller weight, and sampled the fitted spline on a regular 1-min resolution. These pre-processing steps resulted in the modelled animal trajectory data, composed of spatial locations every minute, and averaged ODBA statistics per step (i.e., the segments between consecutive coordinates). These data were used as input for the next steps in the analyses. In contrast to the animal data, the raw intrusion data were of a high temporal resolution and spatial accuracy so that we only needed to subset the data in order to acquire 1-min resolution time-synchronized intrusion trajectories.
The first three parts of the data pre-processing were only needed because of firmware issues in our custom-made sensors. Without these issues, a simple denoising technique like a Kalman filter will suffice.
Feature engineering and processing
We computed a plethora of human-engineered features from the animal trajectories, ODBA data, weather data and several GIS layers with environmental data from the study area (summarized in Supplementary Table S1). All features were computed such that they could not directly be linked to specific points in space or time (by computing movement features relative to the environmental variables), so that only behavioral patterns and abnormalities therein could be linked to intrusion presence. After engineering these base features, we transformed certain features (after visual inspection of the histograms) to approximately symmetric distributions using logarithms. Then we truncated the distributions to the lower and upper 0.001 percentile to correct possible outliers. After that, we standardized all computed features to zero mean and unit variance per species. We also computed scaled versions of selected features by subtracting the mean and dividing by the variance of the selected features per reference set to capture deviations from normal behavior: (1) per area (characterized by a 30 by 30 m neighborhood around each grid cell), (2) per time of day (morning, midday, afternoon) in a period of 5 weeks around each intrusion or control, (3) per area per time of day per 5 weeks, and (4) per individual sentinel per time of day per 5 weeks (Supplementary Table S1). Furthermore, after computing and standardizing the features, we computed more features by applying moving window computations (5 min centered, 10 and 20 min lagging, and the difference between these: 5 min centered minus 10 and 20 min lagging) on the standardized features to capture (the change in) the recent history of animal movement descriptors (mean and standard deviation of all features, fitted Mean Squared Displacement exponential function parameters, net-gross distance ratio and variance of log First Passage Times). Finally, we discretized all features to ordinal values to avoid odd-, fat- and heavy-tailed distributions. In total we computed 2117 features describing different aspects of movement geometry of individual trajectories, herd topology and the interactions with landscape variation.
Subsetting and dimensionality reduction
Before analyzing the computed animal movement features, we applied some filtering on the data. We removed all periods with an experimental intrusion during which there were less than 30 active animal sensors in total. We also removed data of both animals and intrusion when they were close to the reserve’s main gate in order to avoid dilution of the data with other known disturbances. This resulted in 57 intrusions that were selected for further analyses. For every intrusion we selected control data of the same period one or 2 days earlier or later during which no intrusion took place, resulting in an approximately balanced intrusion-control dataset. Furthermore, we removed data from animals that were located within 250 m and within 20 min of a vehicle moving through the area that was not part of our experiment.
For each feature, we computed 4 importance metrics based on binary labelled data: records associated to locations within 1 km from the intrusion (subscript 1) versus an equally-sized random selection of data points during control periods (subscript 0): Mahalanobis distance, marginality (computed as (frac{{mu }_{1}-{mu }_{0}}{{sigma }_{0}}), for sample mean (mu) and sample standard deviation (sigma)), specialization (computed as (frac{{sigma }_{1}}{{sigma }_{0}})) and the Mean Decrease Accuracy of a Random Forest classifier (with default hyperparameters). We then ranked the features according to their importance and selected a feature for further analyses if it occurred in the top 125 features for any of the 4 importance measures described above (resulting in a total of 361 selected features). Subsequently, we converted the selected features per main feature class (Supplementary Table S1) to principal components, keeping those principal components that capture the most variation (in total 95%), which resulted in 99 selected components in total. Finally, we transformed these components again via a second principal component analysis, now across all the selected 99 components. In subsequent training of the animal behavior classifier, we optimized the total number of included components as a hyperparameter, which resulted in the first 8 principal components in the best performing classifier.
Labelling
We labelled the sentinel movement data through visual inspection of the animal and intruder trajectories, where we considered the animals’ behavior to be undisturbed when the animal was not near an intrusion, or when the animal was close to an intrusion yet did not visually display a change in behavior. However, when the animal was near the intrusion and displayed a sudden or gradual behavioral change in response to intrusion proximity, we labelled the data as ‘flight’ (changing the movement direction away from the intrusion, possibly with increased speed) or ‘regroup’ (when individuals clustered together). In total, only ca. 1% of the animal data were associated to either flight or regroup behavior (which we will refer to as ‘response’ behavior). A few animals also appeared to exhibit behavior we could label as ‘freeze’, i.e., halting movement in the proximity of the intrusion, yet this class was too underrepresented to be accurately predicted and hence dropped from the final dataset. Furthermore, we assigned a qualitative measure of intensity to each labelled behavioral response (‘low’, ‘medium’, ‘high’) to describe how visually pronounced this response was. Besides the supervised labelling based on visual inspection of behavioral responses via video animations of the trajectories, we also labelled data using an unsupervised k-means nearest neighbor classifier, where we clustered the feature space consisting of the 99 features selected as described above into 25 clusters per species.
Animal behavior classification
We trained an RBF kernel C-classification Support Vector Machine (SVM) with a subsequent moving window over the outputted probabilities to distinguish undisturbed versus response behavior. In the training datasets we only included the data separated by more than 1 km from the intrusion and labelled as ‘undisturbed’, and removed 90% thereof to train algorithms with a more balanced dataset. Furthermore, we only trained and validated on data with intrusions present in the area. We trained another SVM to distinguish the flight response from the regroup response. All computations were done in R 3.5.0 with the e1071 package on the Linux High Performance Cluster of Wageningen University and Research. We optimized the following hyperparameters and model settings during the training phase for the Average Precision via a grid search (with the selected values between brackets):

gamma (undisturbed-response: 10–3.2; flight-regroup: 10–2.0);

cost (undisturbed-response: 10–2.2; flight-regroup: 10–1.5);

number of principal components to include as features (undisturbed-response: 8; flight-regroup: 12);

species-specific models versus one model with species dummy variables included in the features (species-specific models);

specific models for the different times of day versus one model with time of day dummy variables (one model);

response intensities to include in the training data (only medium and high intensities);

weights to assign to the classes (equal weights);

the quantile to be computed of the SVM probabilities by the moving window (100%, i.e., maximum value);

the alignment of the moving window (centered);

the size of the moving window (15 min on both sides).

The best model was selected via a leave-one-intrusion-out cross-validation approach. We summarized the predictive performance by computing the Average Precision of the least occurring class (i.e., ‘response’ for the undisturbed-response model: 46%, Supplementary Fig. S2; and ‘regroup’ for the flight-regroup model: 80%, Supplementary Fig. S3). After having computed these probabilities with an SVM and a temporal window smoother, we tried to improve the predicted performance by including the predicted animal response probabilities of nearby animals. However, this spatial explicit approach hardly improved the predictive performance, indicating that the spatial contextualization of behavioral response was sufficiently captured by the computed features. We therefore did not include this spatial contagion effect of predicted animal response probabilities in the final analysis.
System classification—detection
Based on the predicted SVM response probabilities and feature cluster analysis, we computed summary features per 15 min of each intrusion and control period. These summary features related to the odds ratios of the probability of association of unsupervised clusters with intrusions versus controls, the SVM predicted probabilities of behavioral response, and several features describing the values (and its spatial structure, e.g., clustering or autocorrelation) of these SVM predicted response probabilities. After computing summary features per 15 min, we summarized them even further for the intrusions versus controls using the following eight statistics: mean, standard deviation, minimum, maximum, mean of the lagged differences, standard deviation of the lagged differences, minimum of the lagged differences and maximum of the lagged differences.
After computing the summary features, we build a logistic regression classifier to distinguish intrusions from controls. To create a parsimonious model, we iteratively added features to the model and evaluated its performance after each iteration. We evaluated the performance based on the model accuracy and performed validation through 25 times twofold cross-validation in a stratified way (by 25 times choosing a balanced random sample of intrusions and controls). We determined the sequence of adding features to the model by performing an independent two-sample t-test for each feature between the intrusions and controls. The feature with the largest t-value was then added to the model. After each feature addition, we removed its correlation with the remaining features using linear regressions with the added feature as independent variable and the remaining features as dependent variables, from which we extracted the residuals, standardized them to zero mean and unit variance, and applied the t-tests again. The (original) feature with the largest t-value was then added to the model again. This procedure was repeated until all features were ordered corresponding to their “importance”. We then performed logistic regressions without interactions between the features for an increasing number of features (Supplementary Fig. S4). The model already performed quite accurately with only 7 features (86.1% accuracy ± SD 3.3%, precision 82.6% ± SD 6.9%, recall 89.2% ± SD 5.1%). However, with 20 features and 2-way interactions the model achieved the maximum accuracy (90.9%).
System classification—localization
The data gathered during intrusions that were correctly predicted as such by the detection classifier were used to train the intrusion localization algorithm. The probability surface of the location of the intrusion was fitted relative to that of the sentinel animals using:

$${O}_{i,j} sim frac{{p}_{j}left({f}_{wn}left({theta }_{i,j},{mu }_{j},{rho }_{1}right) {f}_{ln}left({gamma }_{i,j},{mu }_{1},{sigma }_{1}right)right) left(1-{p}_{j}right) left({f}_{wn}left({theta }_{i,j},{mu }_{j},{sigma }_{0}right) {f}_{ln}left({gamma }_{i,j},{mu }_{0},{sigma }_{0}right)right)}{{f}_{wn}left({theta }_{i,j},{mu }_{j},{rho }_{0}right) {f}_{ln}({gamma }_{i,j},{mu }_{0},{sigma }_{0})}$$

where ({O}_{i,j}) is the odds ratio of intrusion presence at location (i) evaluated for individual (j), ({p}_{j}) is the SVM-predicted probability that individual (j) is exhibiting response behavior. The function ({f}_{wn}) is the wrapped normal probability density function, ({theta }_{i,j}) is the direction from location (i) to the location of the focal animal (j), ({mu }_{j}) is the movement direction of individual (j), ({rho }_{1}) and ({rho }_{0}) are the standard deviations of the unwrapped distributions. The function ({f}_{ln}) is the lognormal probability density function, where ({gamma }_{i,j}) is the distance of location (i) to (j), ({mu }_{1}) and ({mu }_{0}) as well as ({sigma }_{1}) and ({sigma }_{0}) are the log-normal distribution parameters (respectively log-mean and log-sd).
The parameters ({mu }_{1}), ({sigma }_{1}) and ({rho }_{1}) capture the geometry of intrusion-animal topology for animals that exhibited a predicted behavioral response to the intrusion. Similarly, ({mu }_{0}), ({sigma }_{0}) and ({rho }_{0}) are the corresponding parameters for animals that were predicted to be undisturbed. The parameters ({mu }_{1}), (mathrm{log}({sigma }_{1})) and (mathrm{log}({rho }_{1})) were fitted to the data assuming a 3rd order polynomial relationship to ({t}_{s}): the time (in minutes) since the start of the predicted behavioral response (using the maximum F1 classification score). Since the behavioral response signature is lost over time, we truncated ({t}_{s}) to 45 min (thus ({t}_{s} >45) min was set to ({t}_{s}=45)). The parameters ({mu }_{0}), ({sigma }_{0}) and ({rho }_{0}) were estimated using the data of the controls and with randomly generated intrusion locations in the study area, in order to correct for the effects of geometry of the study area on the predicted response surfaces. The probability surface ({P}_{i}) was then calculated as:

$${P}_{i}=alpha sum_{j}{O}_{i,j}$$

where (alpha) is a normalization constant so that ({P}_{i}) integrates to 1 over the area covered by the rectangular axis-aligned bounding box around the study area.
To measure the prediction accuracy of each localization surface, we simplified each surface to a point coordinate located at the location of maximum probability, and computed the Euclidian distance to the known true position of the intrusion. We then summarized each experimental intrusion by selecting the 10 prediction surfaces with the most condense highest probability density, i.e., those in which the top 5% probability density is contained in the smallest, most condense, area. The spatial error of the localization prediction associated with these selected predictions was further summarized by taking the average Euclidian distance over the 10 selected predictions. More

1.
Andela, N. et al. A human-driven decline in global burned area. Science 356, 1356–1362 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 
2.
Westerling, A. L., Hidalgo, H. G., Cayan, D. R. & Swetnam, T. W. Warming and earlier spring increase western US forest wildfire activity. Science 313, 940–943 (2006).
CAS  PubMed  Article  Google Scholar 

3.
Turner, M. G. Disturbance and landscape dynamics in a changing world. Ecology 91, 2833–2849 (2010).
PubMed  Article  PubMed Central  Google Scholar 

4.
Higgins, S. I. & Scheiter, S. Atmospheric CO2 forces abrupt vegetation shifts locally, but not globally. Nature 488, 209–212 (2012).
CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
van der Werf, G. R. G. R. et al. Global fire emissions estimates during 1997–2016. Earth Syst. Sci. Data 9, 697–720 (2017).
Article  Google Scholar 

6.
Schoennagel, T. et al. Adapt to more wildfire in western North American forests as climate changes. Proc. Natl Acad. Sci. USA 114, 4582–4590 (2017).
CAS  PubMed  Article  PubMed Central  Google Scholar 

7.
Westerling, A. L., Turner, M. G., Smithwick, E. A. H., Romme, W. H. & Ryan, M. G. Continued warming could transform Greater Yellowstone fire regimes by mid-21st century. Proc. Natl Acad. Sci. USA 108, 13165–13170 (2011).
CAS  PubMed  Article  PubMed Central  Google Scholar 

8.
Johnstone, J. F. et al. Changing disturbance regimes, ecological memory, and forest resilience. Front. Ecol. Environ. 14, 369–378 (2016).
Article  Google Scholar 

9.
Lewis, T. Very frequent burning encourages tree growth in sub-tropical Australian eucalypt forest. Forest Ecol. Manag. 459, 117842 (2020).
Article  Google Scholar 

10.
Peterson, D. W. & Reich, P. B. Prescribed fire in oak savanna: fire frequency effects on stand structure and dynamics. Ecol. Appl. 11, 914–927 (2001).
Article  Google Scholar 

11.
Tilman, D. et al. Fire suppression and ecosystem carbon storage. Ecology 81, 2680–2685 (2000).
Article  Google Scholar 

12.
Pellegrini, A. F. A., Hedin, L. O., Staver, A. C. & Govender, N. Fire alters ecosystem carbon and nutrients but not plant nutrient stoichiometry or composition in tropical savanna. Ecology 96, 1275–1285 (2015).
PubMed  Article  PubMed Central  Google Scholar 

13.
Russell-Smith, J., Whitehead, P. J., Cook, G. D. & Hoare, J. L. Response of eucalyptus-dominated savanna to frequent fires: lessons from Munmarlary, 1973–1996. Ecol. Monogr. 73, 349–375 (2003).
Article  Google Scholar 

14.
Uhl, C. & Kauffman, J. B. Deforestation, fire susceptibility, and potential tree responses to fire in the eastern Amazon. Ecology 71, 437–449 (1990).
Article  Google Scholar 

15.
Case, M. F., Wigley‐Coetsee, C., Nzima, N., Scogings, P. F. & Staver, A. C. Severe drought limits trees in a semi‐arid savanna. Ecology 100, e02842 (2019).
PubMed  Article  PubMed Central  Google Scholar 

16.
Keeley, J. E., Pausas, J. G., Rundel, P. W., Bond, W. J. & Bradstock, R. A. Fire as an evolutionary pressure shaping plant traits. Trends Plant Sci. 16, 406–411 (2011).
CAS  PubMed  Article  PubMed Central  Google Scholar 

17.
Schoennagel, T., Turner, M. G. & Romme, W. H. The influence of fire interval and serotiny on postfire lodgepole pine density in Yellowstone National Park. Ecology 84, 2967–2978 (2003).
Article  Google Scholar 

18.
Higgins, S. I. et al. Which traits determine shifts in the abundance of tree species in a fire-prone savanna? J. Ecol. 100, 1400–1410 (2012).
Article  Google Scholar 

19.
Lehmann, C. E. R. et al. Savanna vegetation–fire–climate relationships differ among continents. Science 343, 548–552 (2014).
CAS  PubMed  Article  PubMed Central  Google Scholar 

20.
Staver, A. C., Archibald, S. & Levin, S. A. The global extent and determinants of savanna and forest as alternative biome states. Science 334, 230–232 (2011).
CAS  PubMed  Article  PubMed Central  Google Scholar 

21.
Higgins, S. I., Bond, J. I. & Trollope, W. S. Fire, resprouting and variability: a recipe for grass–tree coexistence in savanna. J. Ecol. 88, 213–229 (2000).
Article  Google Scholar 

22.
Pellegrini, A. F. A. et al. Fire frequency drives decadal changes in soil carbon and nitrogen and ecosystem productivity. Nature 553, 194–198 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

23.
Reich, P. B., Peterson, D. W., Wedin, D. A. & Wrage, K. Fire and vegetation effects on productivity and nitrogen cycling across a forest–grassland continuum. Ecology 82, 1703–1719 (2001).
Google Scholar 

24.
Phillips, R., Brzostek, E. & Midgley, M. The mycorrhizal‐associated nutrient economy: a new framework for predicting carbon–nutrient couplings in temperate forests. New Phytol. 99, 41–51 (2013).
Article  CAS  Google Scholar 

25.
Hobbie, S. E. Plant species effects on nutrient cycling: revisiting litter feedbacks. Trends Ecol. Evol. 30, 357–363 (2015).
PubMed  Article  PubMed Central  Google Scholar 

26.
Read, D. J. & Perez‐Moreno, J. Mycorrhizas and nutrient cycling in ecosystems – a journey towards relevance? New Phytol. 157, 475–492 (2003).
Article  Google Scholar 

27.
Dixon, R. K. et al. Carbon pools and flux of global forest ecosystems. Science 263, 185–190 (1994).
CAS  PubMed  Article  PubMed Central  Google Scholar 

28.
Jackson, R. B. et al. Trading water for carbon with biological carbon sequestration. Science 310, 1944–1947 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

29.
Whitman, E., Parisien, M. A., Thompson, D. K. & Flannigan, M. D. Short-interval wildfire and drought overwhelm boreal forest resilience. Sci. Rep. 9, 18796 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

30.
Hart, S. J. et al. Examining forest resilience to changing fire frequency in a fire-prone region of boreal forest. Glob. Change Biol. 25, 869–884 (2019).
Article  Google Scholar 

31.
Stephens, S. L. et al. Managing forests and fire in changing climates. Science 342, 41–42 (2013).
CAS  PubMed  Article  Google Scholar 

32.
Steel, Z. L., Safford, H. D. & Viers, J. H. The fire frequency–severity relationship and the legacy of fire suppression in California forests. Ecosphere 6, 1–23 (2015).
Article  Google Scholar 

33.
Scott, J. & Burgan, R. Standard Fire Behavior Fuel Models: A Comprehensive Set for Use with Rothermel’s Surface Fire Spread Model General Technical Report RMRS-GTR-153 (USDA, Forest Service and Rocky Mountain Research Station, 2005).

34.
Liu, Y. Y. et al. Recent reversal in loss of global terrestrial biomass. Nat. Clim. Change 5, 470–474 (2015).
Article  Google Scholar 

35.
Brandt, M. et al. Satellite passive microwaves reveal recent climate-induced carbon losses in African drylands. Nat. Ecol. Evol. 2, 827–835 (2018).
PubMed  Article  PubMed Central  Google Scholar 

36.
Butler, O. M., Elser, J. J., Lewis, T., Mackey, B. & Chen, C. The phosphorus-rich signature of fire in the soil–plant system: a global meta-analysis. Ecol. Lett. 21, 335–344 (2018).
PubMed  Article  PubMed Central  Google Scholar 

37.
Raison, R. J., Khanna, P. K. & Woods, P. V. Transfer of elements to the atmosphere during low-intensity prescribed fires in three Australian subalpine eucalypt forests. Can. J. Forest Res. 15, 657–664 (1985).
CAS  Article  Google Scholar 

38.
Averill, C., Bhatnagar, J. M., Dietze, M. C., Pearse, W. D. & Kivlin, S. N. Global imprint of mycorrhizal fungi on whole-plant nutrient economics. Proc. Natl. Acad. Sci. USA https://doi.org/10.1073/pnas.1906655116 (2019).

39.
Shah, F. et al. Ectomycorrhizal fungi decompose soil organic matter using oxidative mechanisms adapted from saprotrophic ancestors. New Phytol. 209, 1705–1719 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

40.
Woinarski, J. C. Z., Risler, J. & Kean, L. Response of vegetation and vertebrate fauna to 23 years of fire exclusion in a tropical eucalyptus open forest, Northern Territory, Australia. Austral Ecol. 29, 156–176 (2004).
Article  Google Scholar 

41.
Steidinger, B. S. et al. Climatic controls of decomposition drive the global biogeography of forest–tree symbioses. Nature 569, 404–408 (2019).
CAS  PubMed  Article  PubMed Central  Google Scholar 

42.
Pellegrini, A. F. A. et al. Repeated fire shifts carbon and nitrogen cycling by changing plant inputs and soil decomposition across ecosystems. Ecol. Monogr. 90, e01409 (2020).
Article  Google Scholar 

43.
Newland, J. A. & DeLuca, T. H. Influence of fire on native nitrogen-fixing plants and soil nitrogen status in ponderosa pine – Douglas-fir forests in western Montana. Can. J. Forest Res. 30, 274–282 (2000).
Article  Google Scholar 

44.
Johnson, D. W. & Curtis, P. S. Effects of forest management on soil C and N storage: meta analysis. Forest Ecol. Manag. 140, 227–238 (2001).
Article  Google Scholar 

45.
Pellegrini, A. F. A. Nutrient limitation in tropical savannas across multiple scales and mechanisms. Ecology 97, 313–324 (2016).
PubMed  Article  PubMed Central  Google Scholar 

46.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).
Article  Google Scholar 

47.
Harrison, X. A. et al. A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ 2018, e4794 (2018).
Article  Google Scholar 

48.
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).
Article  Google Scholar 

49.
Jackson, J. F., Adams, D. C. & Jackson, U. B. Allometry of constitutive defense: a model and a comparative test with tree bark and fire regime. Am. Nat. 153, 614–632 (1999).
PubMed  Article  PubMed Central  Google Scholar 

50.
Chave, J. et al. Towards a worldwide wood economics spectrum. Ecol. Lett. 12, 351–366 (2009).
PubMed  Article  PubMed Central  Google Scholar 

51.
Hoffmann, W. A., Marchin, R. M., Abit, P. & Lau, O. L. Hydraulic failure and tree dieback are associated with high wood density in a temperate forest under extreme drought. Glob. Change Biol. 17, 2731–2742 (2011).
Article  Google Scholar 

52.
Harmon, M. E. Decomposition of standing dead trees in the southern Appalachian Mountains. Oecologia 52, 214–215 (1982).
PubMed  Article  PubMed Central  Google Scholar 

53.
Hedges, L. V., Gurevitch, J. & Curtis, P. S. The meta-analysis of response ratios in experimental ecology. Ecology 80, 1150–1156 (1999).
Article  Google Scholar 

54.
Gurevitch, J., Morrow, L. L., Wallace, A. & Walsh, J. S. A meta-analysis of competition in field experiments. Am. Nat. 140, 539–572 (1992).
Article  Google Scholar 

55.
Zanne, A. E. et al. Three keys to the radiation of angiosperms into freezing environments. Nature 506, 89–92 (2014).
CAS  PubMed  Article  PubMed Central  Google Scholar 

56.
Pearse, W. D. et al. pez: phylogenetics for the environmental sciences. Bioinformatics 31, 2888–2890 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

57.
Kembel, S. W. et al. Picante: R tools for integrating phylogenies and ecology. Bioinformatics 26, 1463–1464 (2010).
CAS  PubMed  Article  PubMed Central  Google Scholar 

58.
Brockway, D. G. & Lewis, C. E. Long-term effects of dormant-season prescribed fire on plant community diversity, structure and productivity in a longleaf pine wiregrass ecosystem. Forest Ecol. Manag. 96, 167–183 (1997).
Article  Google Scholar 

59.
Lewis, T. & Debuse, V. J. Resilience of a eucalypt forest woody understorey to long-term (34–55 years) repeated burning in subtropical Australia. Int. J. Wildl. Fire 21, 980–991 (2012).
Article  Google Scholar 

60.
Scudieri, C. A., Sieg, C. H., Haase, S. M., Thode, A. E. & Sackett, S. S. Understory vegetation response after 30 years of interval prescribed burning in two ponderosa pine sites in northern Arizona, USA. Forest Ecol. Manag. 260, 2134–2142 (2010).
Article  Google Scholar 

61.
Lewis, T., Reif, M., Prendergast, E. & Tran, C. The effect of long-term repeated burning and fire exclusion on above- and below-ground blackbutt (Eucalyptus pilularis) forest vegetation assemblages. Austral Ecol. 37, 767–778 (2012).
Article  Google Scholar 

62.
Stratton, R. Effects of Long-Term Late Winter Prescribed Fire on Forest Stand Dynamics, Small Mammal Populations, and Habitat Demographics in a Tennessee Oak Barrens. MSc thesis, Univ. Tennessee (2007).

63.
Wade, D. D. Long-Term Site Responses to Season and Interval of Underburns on the Georgia Piedmont (Forest Service Research Data Archive, 2016).

64.
Pellegrini, A. F. A., Hoffmann, W. A. & Franco, A. C. Carbon accumulation and nitrogen pool recovery during transitions from savanna to forest in central Brazil. Ecology 95, 342–352 (2014).
PubMed  Article  PubMed Central  Google Scholar 

65.
Nesmith, C. B., Caprio, A. C., Pfaff, A. H., McGinnis, T. W. & Keeley, J. E. A comparison of effects from prescribed fires and wildfires managed for resource objectives in Sequoia and Kings Canyon National Parks. Forest Ecol. Manag. 261, 1275–1282 (2011).
Article  Google Scholar 

66.
Haywood, J. D., Harris, F. L., Grelen, H. E. & Pearson, H. A. Vegetative response to 37 years of seasonal burning on a Louisiana longleaf pine site. South. J. Appl. For. 25, 122–130 (2001).
Article  Google Scholar 

67.
Higgins, S. I. et al. Effects of four decades of fire manipulation on woody vegetation structure in savanna. Ecology 88, 1119–1125 (2007).
PubMed  Article  PubMed Central  Google Scholar 

68.
Gignoux, J., Lahoreau, G., Julliard, R. & Barot, S. Establishment and early persistence of tree seedlings in an annually burned savanna. J. Ecol. 97, 484–495 (2009).
Article  Google Scholar 

69.
Tizon, F. R., Pelaez, D. V. & Elia, O. R. The influence of controlled fires on a plant community in the south of the Caldenal and its relationship with a regional state and transition model. Int. J. Exp. Bot. 79, 141–146 (2010).
Google Scholar 

70.
Neill, C., Patterson, W. A. & Crary, D. W. Responses of soil carbon, nitrogen and cations to the frequency and seasonality of prescribed burning in a Cape Cod oak–pine forest. Forest Ecol. Manag. 250, 234–243 (2007).
Article  Google Scholar 

71.
Ryan, C. M., Williams, M. & Grace, J. Above‐ and belowground carbon stocks in a miombo woodland landscape of Mozambique. Biotropica 43, 423–432 (2011).
Article  Google Scholar 

72.
Scharenbroch, B. C., Nix, B., Jacobs, K. A. & Bowles, M. L. Two decades of low-severity prescribed fire increases soil nutrient availability in a midwestern, USA oak (Quercus) forest. Geoderma 183–184, 80–91 (2012).
Article  CAS  Google Scholar 

73.
Burton, J. A., Hallgren, S. W., Fuhlendorf, S. D. & Leslie, D. M. Jr. Understory response to varying fire frequencies after 20 years of prescribed burning in an upland oak forest. Plant Ecol. 212, 1513–1525 (2011).
Article  Google Scholar 

74.
Stewart, J. F., Will, R. E., Robertson, K. M. & Nelson, C. D. Frequent fire protects shortleaf pine (Pinus echinata) from introgression by loblolly pine (P. taeda). Conserv. Genet. 16, 491–495 (2015).
Article  Google Scholar 

75.
Knapp, B. O., Stephan, K. & Hubbart, J. A. Structure and composition of an oak–hickory forest after over 60 years of repeated prescribed burning in Missouri, U.S.A. Forest Ecol. Manag. 344, 95–109 (2015).
Article  Google Scholar 

76.
Olson, M. G. Tree regeneration in oak–pine stands with and without prescribed fire in the New Jersey Pine Barrens: management implications. North. J. Appl. For. 28, 47–49 (2011).
Article  Google Scholar  More

1.
Shipley, B. et al. Reinforcing loose foundation stones in trait-based plant ecology. Oecologia 180, 923–931 (2016).
PubMed  Article  PubMed Central  Google Scholar 
2.
McGill, B. J., Enquist, B. J., Weiher, E. & Westoby, M. Rebuilding community ecology from functional traits. Trends Ecol. Evol. 21, 178–185 (2006).
PubMed  Article  PubMed Central  Google Scholar 

3.
Vellend, M. Conceptual synthesis in community ecology. Q. Rev. Biol. 85, 183–206 (2010).
PubMed  Article  PubMed Central  Google Scholar 

4.
Bjorkman, A. D. et al. Plant functional trait change across a warming tundra biome. Nature 562, 57–62 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
Billings, W. D. Arctic and Alpine vegetations: similarities, differences, and susceptibility to disturbance. BioScience 23, 697–704 (1973).
Article  Google Scholar 

6.
Graae, B. J. et al. Stay or go – how topographic complexity influences alpine plant population and community responses to climate change. Perspect. Plant Ecol. Evol. Syst. 30, 41–50 (2018).
Article  Google Scholar 

7.
Bruelheide, H. et al. Global trait–environment relationships of plant communities. Nat. Ecol. Evol. 2, 1906–1917 (2018).
PubMed  Article  PubMed Central  Google Scholar 

8.
Choler, P. Consistent shifts in alpine plant traits along a mesotopographical gradient. Arct. Antarct. Alp. Res. 37, 444–453 (2005).
Article  Google Scholar 

9.
Wullschleger, S. D. et al. Plant functional types in Earth system models: past experiences and future directions for application of dynamic vegetation models in high-latitude ecosystems. Ann. Bot. 114, 1–16 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

10.
Pearson, R. G. et al. Shifts in Arctic vegetation and associated feedbacks under climate change. Nat. Clim. Change 3, 673–677 (2013).
Article  Google Scholar 

11.
Myers-Smith, I. H., Thomas, H. J. D. & Bjorkman, A. D. Plant traits inform predictions of tundra responses to global change. New Phytol. 221, 1742–1748 (2019).
PubMed  Article  PubMed Central  Google Scholar 

12.
Robinson, S. A. et al. Rapid change in East Antarctic terrestrial vegetation in response to regional drying. Nat. Clim. Change 8, 879–884 (2018).
CAS  Article  Google Scholar 

13.
Post, E. et al. Ecological dynamics across the Arctic associated with recent climate change. Science 325, 1355–1358 (2009).
CAS  PubMed  Article  PubMed Central  Google Scholar 

14.
Saros, J. E. et al. Arctic climate shifts drive rapid ecosystem responses across the West Greenland landscape. Environ. Res. Lett. 14, 074027 (2019).
Article  Google Scholar 

15.
Lavorel, S. & Garnier, E. Predicting changes in community composition and ecosystem functioning from plant traits: revisiting the Holy Grail. Funct. Ecol. 16, 545–556 (2002).
Article  Google Scholar 

16.
Chapin, F. S. III et al. Consequences of changing biodiversity. Nature 405, 234–242 (2000).
CAS  PubMed  Article  PubMed Central  Google Scholar 

17.
Díaz, S. et al. The global spectrum of plant form and function. Nature 529, 167–171 (2016).
PubMed  Article  CAS  PubMed Central  Google Scholar 

18.
Wright, I. J. et al. The worldwide leaf economics spectrum. Nature 428, 821–827 (2004).
CAS  Article  Google Scholar 

19.
Thomas, H. J. D. et al. Global plant trait relationships extend to the climatic extremes of the tundra biome. Nat. Commun. 11, 1351 (2020).
CAS  PubMed  PubMed Central  Article  Google Scholar 

20.
Billings, W. D. & Bliss, L. C. An alpine snowbank environment and its effects on vegetation, plant development, and productivity. Ecology 40, 388–397 (1959).
Article  Google Scholar 

21.
Myers-Smith, I. H. & Hik, D. S. Shrub canopies influence soil temperatures but not nutrient dynamics: an experimental test of tundra snow–shrub interactions. Ecol. Evol. 3, 3683–3700 (2013).
PubMed  PubMed Central  Article  Google Scholar 

22.
Chapin, F. S. III et al. Role of land-surface changes in Arctic summer warming. Science 310, 657–660 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

23.
Cahoon, S. M. P. et al. Interactions among shrub cover and the soil microclimate may determine future Arctic carbon budgets. Ecol. Lett. 15, 1415–1422 (2012).
PubMed  Article  PubMed Central  Google Scholar 

24.
Reich, P. B. The world-wide ‘fast–slow’ plant economics spectrum: a traits manifesto. J. Ecol. 102, 275–301 (2014).
Article  Google Scholar 

25.
Diaz, S. et al. The plant traits that drive ecosystems: evidence from three continents. J. Veg. Sci. 15, 295–304 (2004).
Article  Google Scholar 

26.
Cornelissen, J. H. C. et al. Global negative vegetation feedback to climate warming responses of leaf litter decomposition rates in cold biomes. Ecol. Lett. 10, 619–627 (2007).
PubMed  Article  PubMed Central  Google Scholar 

27.
Steinbauer, M. J. et al. Accelerated increase in plant species richness on mountain summits is linked to warming. Nature 556, 231–234 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

28.
Myers-Smith, I. H. et al. Climate sensitivity of shrub growth across the tundra biome. Nat. Clim. Change 5, 887–891 (2015).
Article  Google Scholar 

29.
Post, E. et al. The polar regions in a 2 °C warmer world. Sci. Adv. 5, eaaw9883 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

30.
IPCC Special Report on Global Warming of 1.5 °C (eds Masson-Delmotte, V. et al.) (WMO, 2018).

31.
Bintanja, R. & Andry, O. Towards a rain-dominated Arctic. Nat. Clim. Change 7, 263–267 (2017).
Article  Google Scholar 

32.
Bromwich, D. H. et al. Central West Antarctica among the most rapidly warming regions on Earth. Nat. Geosci. 6, 139–145 (2013).
CAS  Article  Google Scholar 

33.
Turner, J. et al. Absence of 21st century warming on Antarctic Peninsula consistent with natural variability. Nature 535, 411–415 (2016).
CAS  PubMed  Article  Google Scholar 

34.
Sonesson, M., Wielgolaski, F. E. & Kallio, P. in Fennoscandian Tundra Ecosystems. Ecological Studies (Analysis and Synthesis) Vol. 16 (ed. Wielgolaski, F. E.) 3–28 (Springer, 1975); https://doi.org/10.1007/978-3-642-80937-8_1

35.
Niittynen, P., Heikkinen, R. K. & Luoto, M. Snow cover is a neglected driver of Arctic biodiversity loss. Nat. Clim. Change 8, 997–1001 (2018).
Article  Google Scholar 

36.
Klikoff, L. G. Photosynthetic response to temperature and moisture stress of three timberline meadow species. Ecology 46, 516–517 (1965).
Article  Google Scholar 

37.
Oberbauer, S. F. & Billings, W. D. Drought tolerance and water use by plants along an alpine topographic gradient. Oecologia 50, 325–331 (1981).
PubMed  Article  Google Scholar 

38.
Eskelinen, A., Stark, S. & Männistö, M. Links between plant community composition, soil organic matter quality and microbial communities in contrasting tundra habitats. Oecologia 161, 113–123 (2009).
PubMed  Article  Google Scholar 

39.
Ernakovich, J. G. et al. Predicted responses of Arctic and alpine ecosystems to altered seasonality under climate change. Glob. Change Biol. 20, 3256–3269 (2014).
Article  Google Scholar 

40.
Galen, C. & Stanton, M. L. Responses of snowbed plant species to changes in growing-season length. Ecology 76, 1546–1557 (1995).
Article  Google Scholar 

41.
Starr, G., Oberbauer, S. F. & Ahlquist, L. E. The photosynthetic response of Alaskan tundra plants to increased season length and soil warming. Arct. Antarct. Alp. Res. 40, 181–191 (2008).
Article  Google Scholar 

42.
Happonen, K. et al. Snow is an important control of plant community functional composition in oroarctic tundra. Oecologia 191, 601–608 (2019).
PubMed  PubMed Central  Article  Google Scholar 

43.
Niittynen, P. & Luoto, M. The importance of snow in species distribution models of Arctic vegetation. Ecography 41, 1024–1037 (2018).
Article  Google Scholar 

44.
le Roux, P. C., Aalto, J. & Luoto, M. Soil moisture’s underestimated role in climate change impact modelling in low-energy systems. Glob. Change Biol. 19, 2965–2975 (2013).
Article  Google Scholar 

45.
Lembrechts, J. J. et al. SoilTemp: a global database of near-surface temperature. Glob. Change Biol. 26, 6616–6629 (2020).

46.
Bjorkman, A. D. et al. Tundra Trait Team: a database of plant traits spanning the tundra biome. Glob. Ecol. Biogeogr. 27, 1402–1411 (2018).
Article  Google Scholar 

47.
Maitner, B. S. et al. The bien r package: a tool to access the Botanical Information and Ecology Network (BIEN) database. Methods Ecol. Evol. 9, 373–379 (2018).
Article  Google Scholar 

48.
Kattge, J. et al. TRY – a global database of plant traits. Glob. Change Biol. 17, 2905–2935 (2011).
Article  Google Scholar 

49.
Pedersen, E. J., Miller, D. L., Simpson, G. L. & Ross, N. Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7, e6876 (2019).
PubMed  PubMed Central  Article  Google Scholar 

50.
Niittynen, P. et al. Fine-scale tundra vegetation patterns are strongly related to winter thermal conditions. Nat. Clim. Change 10, 1143–1148 (2020).

51.
Belluau, M. & Shipley, B. Predicting habitat affinities of herbaceous dicots to soil wetness based on physiological traits of drought tolerance. Ann. Bot. 119, 1073–1084 (2017).
PubMed  PubMed Central  Article  Google Scholar 

52.
Kemppinen, J., Niittynen, P., Riihimäki, H. & Luoto, M. Modelling soil moisture in a high-latitude landscape using LiDAR and soil data. Earth Surf. Proc. Land. 43, 1019–1031 (2018).
Article  Google Scholar 

53.
Kemppinen, J., Niittynen, P., Aalto, J., le Roux, P. C. & Luoto, M. Water as a resource, stress and disturbance shaping tundra vegetation. Oikos 128, 811–822 (2019).
Article  Google Scholar 

54.
Giblin, A. E., Nadelhoffer, K. J., Shaver, G. R., Laundre, J. A. & McKerrow, A. J. Biogeochemical diversity along a riverside toposequence in Arctic Alaska. Ecol. Monogr. 61, 415–435 (1991).
Article  Google Scholar 

55.
le Roux, P. C., Virtanen, R. & Luoto, M. Geomorphological disturbance is necessary for predicting fine-scale species distributions. Ecography 36, 800–808 (2013).
Article  Google Scholar 

56.
Finger Higgens, R., Hicks Pries, C. & Virginia, R. A. Trade-offs between wood and leaf production in Arctic shrubs along a temperature and moisture gradient in West Greenland. Ecosystems https://doi.org/10.1007/s10021-020-00541-4 (2020).

57.
Porporato, A. & Rodriguez-Iturbe, I. Ecohydrology-a challenging multidisciplinary research perspective / Ecohydrologie: une perspective stimulante de recherche multidisciplinaire. Hydrol. Sci. J. 47, 811–821 (2002).
Article  Google Scholar 

58.
Legates, D. R. et al. Soil moisture: a central and unifying theme in physical geography. Prog. Phys. Geogr. 35, 65–86 (2011).
Article  Google Scholar 

59.
McLaughlin, B. C. et al. Hydrologic refugia, plants, and climate change. Glob. Change Biol. 23, 2941–2961 (2017).
Article  Google Scholar 

60.
Choler, P. Winter soil temperature dependence of alpine plant distribution: implications for anticipating vegetation changes under a warming climate. Perspect. Plant Ecol. Evol. Syst. 30, 6–15 (2018).
Article  Google Scholar 

61.
Happonen, K. et al. Snow is an important control of plant community functional composition in oroarctic tundra. Oecologia 191, 601–608 (2019).

62.
Doran, P. T. et al. Antarctic climate cooling and terrestrial ecosystem response. Nature 415, 517–520 (2002).
CAS  PubMed  Article  PubMed Central  Google Scholar 

63.
French, D. D. & Smith, V. R. A comparison between Northern and Southern Hemisphere tundras and related ecosystems. Polar Biol. 5, 5–21 (1985).
Article  Google Scholar 

64.
le Roux, P. C. in The Prince Edward Islands: Land–Sea Interactions in a Changing Ecosystem (eds Chown, S. L. & Froneman, P. W.) 39–64 (African Sun Media, 2008).

65.
Devau, N., Le Cadre, E., Jaillarda, B. & Gérarda, F. Soil pH controls the environmental availability of phosphorus: experimental and mechanistic modelling approaches. Appl. Geochem. 24, 2163–2174 (2009).

66.
Stevens, R. J., Laughlin, R. J. & Malone, J. P. Soil pH affects the processes reducing nitrate to nitrous oxide and di-nitrogen. Soil Biol. Biochem. 30, 1119–1126 (1998).
CAS  Article  Google Scholar 

67.
Freschet, G. T., Cornelissen, J. H. C., Van Logtestijn, R. S. P. & Aerts, R. Evidence of the ‘plant economics spectrum’ in a subarctic flora. J. Ecol. 98, 362–373 (2010).
Article  Google Scholar 

68.
Bergholz, K. et al. Fertilization affects the establishment ability of species differing in seed mass via direct nutrient addition and indirect competition effects. Oikos 124, 1547–1554 (2015).
Article  Google Scholar 

69.
Curtin, D., Campbell, C. A. & Jalil, A. Effects of acidity on mineralization: pH-dependence of organic matter mineralization in weakly acidic soils. Soil Biol. Biochem. 30, 57–64 (1998).
CAS  Article  Google Scholar 

70.
Blondeel, H. et al. Light and warming drive forest understorey community development in different environments. Glob. Change Biol. 26, 1681–1696 (2020).
Article  Google Scholar 

71.
Dahlgren, J. P., Eriksson, O., Bolmgren, K., Strindell, M. & Ehrlén, J. Specific leaf area as a superior predictor of changes in field layer abundance during forest succession. J. Veg. Sci. 17, 577–582 (2006).
Article  Google Scholar 

72.
Lembrechts, J. J. et al. Comparing temperature data sources for use in species distribution models: from in‐situ logging to remote sensing. Glob. Ecol. Biogeogr. 28, 1578–1596 (2019).
Article  Google Scholar 

73.
Körner, C. & Hiltbrunner, E. The 90 ways to describe plant temperature. Perspect. Plant Ecol. Evol. Syst. 30, 16–21 (2018).
Article  Google Scholar 

74.
Maclean, I. M. D. Predicting future climate at high spatial and temporal resolution. Glob. Change Biol. 26, 1003–1011 (2019).

75.
Aalto, J., Scherrer, D., Lenoir, J., Guisan, A. & Luoto, M. Biogeophysical controls on soil–atmosphere thermal differences: implications on warming Arctic ecosystems. Environ. Res. Lett. 13, 074003 (2018).
Article  Google Scholar 

76.
Aalto, J., le Roux, P. C. & Luoto, M. Vegetation mediates soil temperature and moisture in Arctic-alpine environments. Arct. Antarct. Alp. Res. 45, 429–439 (2013).
Article  Google Scholar 

77.
Moles, A. T. et al. Which is a better predictor of plant traits: temperature or precipitation? J. Veg. Sci. 25, 1167–1180 (2014).
Article  Google Scholar 

78.
Taylor, R. V. & Seastedt, T. R. Short- and long-term patterns of soil moisture in alpine tundra. Arct. Alp. Res. 26, 14–20 (1994).
Article  Google Scholar 

79.
Lembrechts, J. J. & Lenoir, J. Microclimatic conditions anywhere at any time! Glob. Change Biol. https://doi.org/10.1111/gcb.14942 (2019).

80.
Zellweger, F., De Frenne, P., Lenoir, J., Rocchini, D. & Coomes, D. Advances in microclimate ecology arising from remote sensing. Trends Ecol. Evol. 34, 327–341 (2019).
PubMed  Article  PubMed Central  Google Scholar 

81.
Bramer, I. et al. Advances in monitoring and modelling climate at ecologically relevant scales. Adv. Ecol. Res. 58, 101–161 (2018).

82.
Halbritter, A. H. et al. The handbook for standardized field and laboratory measurements in terrestrial climate change experiments and observational studies (ClimEx). Methods Ecol. Evol. 2, 16147 (2019).
Google Scholar 

83.
Wild, J. et al. Climate at ecologically relevant scales: a new temperature and soil moisture logger for long-term microclimate measurement. Agr. Forest Meteorol. 268, 40–47 (2019).
Article  Google Scholar 

84.
Aalto, J., Riihimäki, H., Meineri, E., Hylander, K. & Luoto, M. Revealing topoclimatic heterogeneity using meteorological station data. Int. J. Climatol. 37, 544–556 (2017).
Article  Google Scholar 

85.
Kearney, M. R., Gillingham, P. K., Bramer, I., Duffy, J. P. & Maclean, I. M. D. A method for computing hourly, historical, terrain‐corrected microclimate anywhere on Earth. Methods Ecol. Evol. https://doi.org/10.1111/2041-210x.13330 (2019).

86.
Bjorkman, A. D. et al. Status and trends in Arctic vegetation: evidence from experimental warming and long-term monitoring. Ambio 49, 678–692 (2020).
PubMed  Article  PubMed Central  Google Scholar 

87.
Vandvik, V., Halbritter, A. H. & Telford, R. J. Greening up the mountain. Proc. Natl Acad. Sci. USA 115, 833–835 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

88.
Bonfils, C. J. W. et al. On the influence of shrub height and expansion on northern high latitude climate. Environ. Res. Lett. 7, 015503 (2012).
Article  Google Scholar 

89.
Zwieback, S., Chang, Q., Marsh, P. & Berg, A. Shrub tundra ecohydrology: rainfall interception is a major component of the water balance. Environ. Res. Lett. 14, 055005 (2019).
Article  Google Scholar 

90.
Robinson, D. A. et al. Global environmental changes impact soil hydraulic functions through biophysical feedbacks. Glob. Change Biol. 25, 1895–1904 (2019).
Article  Google Scholar 

91.
Loranty, M. M. et al. Reviews and syntheses: changing ecosystem influences on soil thermal regimes in northern high-latitude permafrost regions. Biogeosciences 15, 5287–5313 (2018).
CAS  Article  Google Scholar 

92.
Parker, T. C., Subke, J.-A. & Wookey, P. A. Rapid carbon turnover beneath shrub and tree vegetation is associated with low soil carbon stocks at a subarctic treeline. Glob. Change Biol. 21, 2070–2081 (2015).
Article  Google Scholar 

93.
DeMarco, J., Mack, M. C. & Bret-Harte, M. S. Effects of Arctic shrub expansion on biophysical vs. biogeochemical drivers of litter decomposition. Ecology 95, 1861–1875 (2014).
PubMed  Article  PubMed Central  Google Scholar 

94.
Qian, H., Joseph, R. & Zeng, N. Enhanced terrestrial carbon uptake in the northern high latitudes in the 21st century from the Coupled Carbon Cycle Climate Model Intercomparison Project model projections. Glob. Change Biol. 16, 641–656 (2010).
Article  Google Scholar 

95.
Sistla, S. A. et al. Long-term warming restructures Arctic tundra without changing net soil carbon storage. Nature 497, 615–618 (2013).
CAS  PubMed  Article  PubMed Central  Google Scholar 

96.
Climate in Svalbard 2100 – A Knowledge Base for Climate Adaptation (Norwegian Centre for Climate Services, 2019); https://go.nature.com/3tFTKAr

97.
Weather Observations from Greenland 1958–2018 – Observation Data with Description DMI Report 19-08 (Danish Meteorological Institute, 2019); https://go.nature.com/36RkdBk

98.
Enontekiö Kilpisjärvi Saana. Daily Climate Observations (Finnish Meteorological Institute, 2019); https://en.ilmatieteenlaitos.fi/download-observations

99.
Enontekiö Kilpisjärvi Kyläkeskus. Daily Climate Observations (Finnish Meteorological Institute, 2019); https://en.ilmatieteenlaitos.fi/download-observations

100.
Smith, V. R. & Steenkamp, M. Classification of the terrestrial habitats on Marion Island based on vegetation and soil chemistry. J. Veg. Sci. 12, 181–198 (2001).
Article  Google Scholar 

101.
Beck, H. E. et al. Present and future Köppen–Geiger climate classification maps at 1-km resolution. Sci. Data 5, 180214 (2018).
PubMed  PubMed Central  Article  Google Scholar 

102.
Canadell, J. et al. Maximum rooting depth of vegetation types at the global scale. Oecologia 108, 583–595 (1996).
CAS  PubMed  Article  PubMed Central  Google Scholar 

103.
Iversen, C. M. et al. The unseen iceberg: plant roots in Arctic tundra. New Phytol. 205, 34–58 (2015).
PubMed  Article  PubMed Central  Google Scholar 

104.
Kern, R. et al. Comparative vegetation survey with focus on cryptogamic covers in the high Arctic along two differing catenas. Polar Biol. 42, 2131–2145 (2019).
Article  Google Scholar 

105.
Miller, R. O. & Kissel, D. E. Comparison of soil pH methods on soils of North America. Soil Sci. Soc. Am. J. 74, 310–316 (2010).
Article  CAS  Google Scholar 

106.
McCune, B. & Keon, D. Equations for potential annual direct incident radiation and heat load. J. Veg. Sci. 13, 603 (2002).
Article  Google Scholar 

107.
McCune, B. Improved estimates of incident radiation and heat load using non- parametric regression against topographic variables. J. Veg. Sci. 18, 751 (2007).
Article  Google Scholar 

108.
Karger, D. N. et al. Climatologies at high resolution for the earth’s land surface areas. Sci. Data 4, 170122 (2017).
PubMed  PubMed Central  Article  Google Scholar 

109.
NASA/METI/AIST/Japan Spacesystems, and U.S./Japan ASTER Science Team ASTER Global Digital Elevation Model (GDEM) V003 (NASA EOSDIS Land Processes DAAC, 2018); https://doi.org/10.5067/ASTER/ASTGTM.003

110.
Hengl, T. et al. SoilGrids250m: global gridded soil information based on machine learning. PLoS ONE 12, e0169748 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

111.
Wood, S. N. Generalized Additive Models: An Introduction with R 2nd edn (Chapman and Hall/CRC, 2017).

112.
R Core Team R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2018).

113.
Husson, F., Le, S. & Pagès, J. Exploratory Multivariate Analysis by Example Using R (CRC, 2017).

114.
Lê, S., Josse, J. & Husson, F. FactoMineR: an R package for multivariate analysis. J. Stat. Softw. 25, 31844 (2008).
Article  Google Scholar 

115.
Kemppinen, J. et al. Data from: Consistent trait–environment relationships within and across tundra plant communities. Zenodo https://doi.org/10.5281/zenodo.4362216 (2020). More

Figure 2

Schematic figure of the labeling process. Participants usually upload several images in a single report. The best photo is picked by the validator who first marks the harassing or non-appropriate photos as hidden. All the non-best photos are marked as not classified. In some rare events, two or three images are annotated from the same report. The mosquito images are classified into four different categories (Aedes albopictus, Aedes aegypti, other species or can not tell) and also the confidence of the label is marked as probable or confirmed. In this paper we excluded the not classified, the hidden and the can not tell images.

Full size image

Between 2014 and 2019, 7686 citizen-made mosquito photos were labeled through Mosquito Alert by entomology experts, with labels indicating whether Ae. albopictus appear in the photos. The photos were included in reports that Mosquito Alert participants uploaded, and each report could contain several photos, see Fig. 2. The entomology experts usually labeled the best photo of the report, but sometimes they labeled two (420 times) or three (49 times) for a single report, meaning that the dataset consisted of 7168 reports. For 6699 reports, only one image was labeled by the experts; for 420 reports two were labeled; for 49 reports three were labeled. Although these reports usually contain several photos, only the ones with expert labels were used in the analysis, as cannot be assumed that all of the photos in a report would have been given the same label.
The main goals of Mosquito Alert during this 6 year period were to monitor Ae. albopictus spreading and provide early detection of Ae. aegypti in Spain. Although people participate in Mosquito Alert all over the world, the majority of the participants and the majority of the photos are in Spain (see Fig. 1). As Ae. aegypti has not been reported in Spain in recent times, most Mosquito Alert participants lived in areas where Ae. aegypti is not present, so most of the photos are of Ae. albopictus. For the detailed yearly distribution of the photos, see Table 1.
Table 1 The collected and expert validated dataset for the period 2014–2019.
Full size table

A popular deep learning model, ResNet5026 was trained and evaluated on the collected dataset with yearly cross-validation. ResNet50 was used because of its wide popularity and its proven classification power in various datasets. As presenting infinitesimal increments of the classification power is not a goal of this paper, we do not report various ImageNet state-of-the-art model performances. Yearly cross-validation was used to rule out any possibility of information leakage (possibility of a user submitting multiple reports for the same mosquito).
The trained model is not only capable of generating highly accurate predictions, but it can also ease the human annotator workload by auto-marking the images where the neural network is confident and more accurate, leaving more uncertain cases for the entomology experts. Moreover, while visualizing the erroneous predictions a few re-occurring patterns were identified, which can serve as a proposal for how to make images that can be best processed by the model.
Several aspects of the dataset were explored as follows.
Classification
Since Mosquito Alert was centered around Ae. albopictus during the relevant time period (2014–2019), the collected dataset is biased towards this species (Table 1). We explored training classifiers on the Mosquito Alert dataset alone and also tied training on a balanced dataset, where 3896 negative samples were added from the IP10227 dataset of various non-mosquito insects as negative samples. From the IP102 dataset, images similar to mosquitoes, and images of striped insects were selected. Although the presented mosquito alert dataset is filtered to contain only mosquito images, in later use, non-mosquito images might be uploaded by the citizens. Training the CNN on a combination of mosquito and non-mosquito images can improve the model to make correct predictions, classifying non-tiger mosquitoes for those cases too. For testing, in each fold, only the Mosquito Alert dataset was used.
The trained classifiers achieved an extremely high area under the receiver operating characteristic curve (ROC AUC) score of 0.96 (see Fig. 3). The fact that the ROC AUC score for each fold was always over 0.95 proves the consistency of our classifier. Inspecting the confusion matrix shows us that the model tends to make more false positive predictions (assuming tiger mosquito is defined as the positive outcome) than false negatives, resulting in high sensitivity. The augmentation of the Mosquito Alert dataset with various insects from IP102 images to make it more balanced resulted in a slight performance boost and narrowed the gap between the number of false positive and false negative samples as expected, see Table 2.
Figure 3

Left: ROC curve calculated on the prediction of the 7686 images in the Mosquito Alert dataset with yearly cross-validation. The blue line shows the case when only the Mosquito Alert dataset was used for training, the orange when the training dataset was balanced out with the addition of non-tiger mosquito insect images from the IP102 dataset. Also a zoom into the part of the ROC curve, where the two methods differ the most is highlighted. Right: the confusion matrix was calculated on the same predictions when only the Mosquito Alert dataset was used for training. For both, a positive label means tiger mosquito is present.

Full size image

Table 2 Yearly cross-validation results with using the Mosquito Alert dataset alone and its IP102 augmented version.
Full size table

How to take a good picture?
Inspection of the weaknesses of a machine learning model is a fruitful way to gain a deeper understanding of the underlying problems and mechanisms. In our case, a careful review of the mispredicted images led us to useful insights into what makes a photo hard to classify for the deep learning model. On Fig. 8, a few selected examples are presented. Unlike humans, deep learning models rely more on textures than on shapes28. As a consequence, grid-like background patterns or striped objects may easily confuse the machine classifier. A larger rich training set can help to avoid these pitfalls, but we also have the option to advise the participants. If participants avoid confusing setups when taking photos, this can improve the accuracy of the automated classification. These guidelines can be added to the Mosquito Alert application to help participants make good images of mosquitoes.

Do not use striped structure (e.g. mosquito net or fly-flap) as a background.

Avoid complex backgrounds when possible. A few examples: patterned carpet, different nets, reflecting/shiny background, bumpy wallpaper.

Use clear, white background (e.g. a sheet of plain paper is perfect if possible) or hold the mosquito with finger pads.

Make sure that as much as possible the mosquito is in focus and covers a large area of the photo.

In general, it is desirable to have a clean white background with the mosquito centered, and with the image containing as little background as possible.
Dataset size impact on model performance
Modern deep CNNs tend to generate better predictions when trained on larger datasets. In this experiment, we trained a ResNet50 model on 10–20–(cdots )–90–100% of 6686 images and evaluated the model on the remaining 1000 images. The 1000 images were selected from the same year (2019) and all of them came from reports with only one photo. There were 709 tiger mosquitoes out of the 1000 test images. ROC AUC and accuracy were calculated with a 500 round bootstrapping of the 1000 test images.
Figure 4

Training a ResNet50 model on a subsampled training dataset. The model was tested against the same 1000 test images for all the steps and statistics of the test metric was calculated with a 500 round bootstrapping. The curve proves the diversity of the Mosquito Alert dataset and also suggests that in the future when the dataset will be even larger, the classification performance will increase.

Full size image

The mean and the standard deviation of the 500 rounds are shown in Fig. 4 for each training data size. From the figure, we can conclude that the predictive power of the model increases as more data are used. The shape of the curves also suggests that the dataset did not reach its plateau. In the upcoming years, as the dataset size increases, ROC AUC and accuracy enhancement is expected.
On measuring image quality
Through the examined period, Mosquito Alert outreach was promoting a mosquito-targeted data collection strategy. Participants were expected to report two mosquito species (Ae. aegypti and Ae. albopictus). By defining these species as positive samples and all the other potential species of mosquito as negative, the submission decision by participants becomes a binary classification problem. In the majority of cases, when participants submit an image we should expect them to think of having a positive sample. Later, based on entomological expert validation, the true label for the image was obtained.
The main goal of such a surveillance system is to keep the sensitivity of the users as high as possible while keeping their specificity at an acceptable level. Therefore, measuring the sensitivity and specificity of the users would be a plausible quality measure. Unfortunately, there is no available information regarding the non-submitted mosquitoes (the true negative and false negative ones), meaning it is impossible to measure sensitivity. The specificity can be measured only in a special case, when there are no false positive images submitted by the user, resulting in a specificity of 1. Based on the latter argument, focusing on metrics derived from the ratio of the submitted tiger mosquito images vs. all submitted images is not meaningful. Instead, the quality can be measured by the usefulness of the photos from the viewpoint of the expert validator or a CNN, as presented in the next chapter.
Quality evolution of the images through time and space
The Mosquito Alert dataset is a unique collection of mosquito images, because, among other things, it is built from 5 consecutive years (not counting 2014, where less than 100 reports were submitted) and it also provides geolocation tags. This uniqueness of the dataset provides potential identification of time and spatial evolution and dependence of the citizen-based mosquito image quality. To explore such an evolution, we performed two different experiments. Geolocation tags were converted to country, region, and city-level information via the geopy Python package. It was found, that the vast majority (95% of all) of the reports were coming from Spain so we performed the analysis only for the Spanish data.
Figure 5

Number of submitted reports and the fraction of their ratio where the entomology expert annotator could tell if tiger mosquito was presented on the photo or not. The charts are shown for the four cities, where Mosquito Alert was the most popular.

Full size image

First, we explored the fraction of the photos, where the entomology expert marked “can not tell”, because the photo was not descriptive enough to decide which species were presented. Figure 5 shows the ratio of the useful mosquito reports, when mosquito decision was possible, compared to all the mosquito reports. The chart shows the above-mentioned ratio for four Spanish cities, which have the most reports submitted (the same information is showed on Supplementary Fig. S1 as a heatmap over Spain). The Mann–Kendall test on the fraction of useful reports shows p-values of 0.09, 0.09, 0.81, 0.22 for Barcelona, Valencia, Málaga, and Girona, which does not justify the presence of a significant trend in image quality, although any conclusions drawn from five data points must be handled with a pinch of salt. It does not mean anything about the individual participants’ quality progression, because Mosquito Alert is highly open and dynamic, and active participants can constantly change. Of note, through these years, the tiger mosquitoes have widely spread from the east coast to the southern and western regions of Spain29. New (and naive) citizen scientists living in the newly colonized regions have been systematically called to action and participation, thus, limiting the overall learning rate of the Mosquito Alert participants’ population. Our results suggest, that either a dynamic balance exists between naive and experienced participants over the period of data recollection, or mosquito photographing skills are independent of the user experience level. The expectation would be that as the population in Spain became more aware of the presence of tiger mosquitoes and their associated public health risks, the system should experience an increase in the useful report ratio, at least for tiger mosquitoes, and most tiger mosquito photos maybe classified automatically.
Figure 6

1000 random samples were selected for each years data. Separated ResNet50 models were trained on each of the years and each model was tested on the rest of the years data. Metrics were calculated with a 500 round random sampling with replacement from the test data. Left: mean of the 500 round bootstrapped accuracy calculations. Right: mean of the 500 round bootstrapped ROC AUC calculations.

Full size image

Second, we subsampled randomly 1000 images from all years between 2015 and 2019. Then we trained a different ResNet50 on data from the different years and generated predictions for the rest of the data, for each year separately. This way we can explore if data from any year is a “better training material” than the others. The results see Fig. 6, shows that 2015 is the worst training material, providing 0.83–0.84 ROC AUC score for the test period, while the rest (period 2016–2019) is similar, ROC AUC varies between 0.90 and 0.93. The reason why the 2015 data found to be the least favourable for training is its class imbalance, meaning that data from 2015 is extremely biased towards tiger mosquitoes (94%), so when training on 2015 data, the model does not see enough non-tiger mosquito samples, while for the other years lower class imbalance was found (70–80%), see Table 1. In general, machine learning models for classification require a substantial amount of examples for each possible class, in our case tiger and non-tiger mosquitoes, therefore worse performance is expected when training on the 2015 data.
Other than the varying class imbalance, we can conclude that the Mosquito Alert dataset quality is consistent, we did not find any concerning difference between training and testing our model for any of the 2016–2017–2018–2019 data pairs.
Pre-filtering the images before expert validation
Generating human annotations for an image classification task is a labour-intensive and expensive part of any project especially if the annotation requires expert knowledge. Therefore, having a model that generates accurate predictions for a well-defined subset of the data saves a lot of time and cost. We assume that the trained classifier is more accurate when the prediction probability is whether high or low and more inaccurate when it is close to 0.5. With this assumption in mind one can tune the (p_{low}) and (p_{high}) probabilities, in a way that images with a prediction probability (p_{low}< p < p_{high}) are discarded and sent to human validation. Figure 7 Randomly selecting 100,000 (p_{low}) and (p_{high}) thresholds on the predictions which were created via yearly cross-validation. Each time only samples were kept where the predicted probability were out of the ([p_{low};p_{high}]) interval. Each point shows the kept data fraction and the prediction accuracy. Varying the lower and upper predicted probability almost 98% of the images are correctly predicted while keeping 80% of all the images. Full size image Varying (p_{low}) and (p_{high}) provides a trade-off between prediction accuracy and the portion of images sent to human validation. Based on Fig. 7 sending 20% of the images to human validation while having an almost 98% accurate prediction for 80% of the dataset is a fruitful way to combine human labour-power and machine learning together. More

Variation in wood density
The values of Chinese fir’s wood physical properties varied considerably among different geographic sources and Tukey-HSD testing showed that some of these differences were statistically significant (Fig. 1). The maximum value (HNYX-T) of wood all-dry density (WDD) was 62.70% higher than the minimum (FJYK-P). The WDD of each source was consistent with the classification and performance indexes of conifer trees in the timber strength grade for structural use, a standard in China’s forestry industry39: FJYK-P was at level S10 ( HNZJJ-P  > FJYK-P, for which the maximum 58.0% higher than the minimum value. According to the wood grading standards in the grain compression index, HNZJJ-P and FJYK-P were at level II (29.1–44.0 MPa) and the rest of geographic sources were at level III (44.1–59.0 MPa) (Table 3).
The compression strength perpendicular to the grain of total tensile (CPG.TT) among geographic sources was ranked as follows: HNYX-T  > JXCS-R  > HNYX-P  > HNZJJ-P  > FJYK-P (Table 4). Its maximum value (HNYX-T) was 29.3% higher than the minimum (FJYK-P). The ranking for compression strength perpendicular to the grain of total radial (CPG.TR) was slightly different: HNYX-T  > JXCS-R  > HNYX-P  > HNZJJ-P  > FJYK-P, for which the maximum was 42.1% higher than the minimum value. Compression strength perpendicular to the grain of part radial (CPG.PR) had the same rank order as CPG.TT, with a maximum value (HNYX-T) 35.0% higher than the minimum (FJYK-P). Finally, compression strength perpendicular to the grain of part tensile (CPG.PT) was ranked as HNYX-T  > JXCS-R  > HNZJJ-P  > HNYX-P  > FJYK-P for the five geographic sources of Chinese fir.
Table 4 The statistical analysis of wood mechanical properties of Chinese fir.
Full size table

Factors influencing wood physical properties
Climate factors effect on wood physical properties
The influence of precipitation on the three kinds of density was consistent. Pre in January, October, November, and December was positively related to wood density, while it was negatively correlated with density in others months, especially in May (r = − 0.39), June (r = − 0.59), and August (r = − 0.64). On a seasonal scale, Pre in summer was negatively correlated with density (r = − 0.77), but it was positively correlated with autumn (r = 0.22). MaxT was positively correlated with density during the whole year, except in May (r = − 0.34), and likewise with wood density but most strongly in summer (r = 0.75). MinT was positively correlated with density, especially in Jan (r  > 0.7), though it was not significantly so in February and October (r  0.45). Pre showed no significant correlation with TSR.LD, RSR.LD, DDS.LD, and DDS.RD, whose correlation coefficients were 0.1–0.3. But Pre was negatively correlated with VSR.LD most of the year (except July, October). AveT was negatively correlation with TSR.RD, RSR.RD, and VSR.RD in January, February, March, and winter; however, AveT showed no significant correlation with DDS.RD. AveT was negatively correlated with TSR.LD, RSR.LD, DDS.LD, and VSR.LD during the whole year. In general, MinT had a significant positive relationship to TSR.RD (r = 0.47), RSR.RD (r = 0.48), and VSR.RD (r = 0.52), except in October, and it was negatively correlated with DDS.RD. MinT was positively related to RSR.LD, VSR.LD, yet negative related to DDS.LD. MaxT was negatively correlated with TSR.RD, RSR.RD, VSR.RD in January, February, May, and December, and winter. MaxT showed no significant correlation with DDS.RD, RSR.LD, DDS.LD or VSR.LD (Fig. 2c).
Pre had significant negative correlations with all of the mechanical properties in May, June, August, and summer, as evince by Fig. 2b, which also showed positive correlations in October. As we can seen, the effects of Pre on wood density and mechanical properties have the same tendency. Pre in all other months was not significantly correlated with mechanical properties (r  0.75), while it was showed no significant correlation in Feb and Oct (r  1000. Through stepwise regression modeling, 14 variables without multicollinearity were retained (i.e., MOE, MOR, TSG, CSG, CPG.TT, CPG.TR, CPG.PT, CPG.PR, DDS.RD, WDD, DDS.LD, TSR.RD, RSR.RD, VSR.LD).
PCA was applied to the above 14 selected physical variables. These results showed that the physical properties of wood loaded strongly on the first axis of the PCA, explaining 51.8% of variation in the 14 tested properties, while the second axis explained 11.0% of it. MOE, MOR, TSR.RD, RSR.RD, and VSR.LD loaded on the positive axis of PC1 and PC2. Both DDS.LD and DDS.RD loaded on the negative axis of PC1 and PC2, while TSG, CSG, CPG.TT, CPG.TR, CPG.PT, CPG.PR, and WDD loaded on the positive axis of PC1 and the negative axis of PC2 (Fig. 3). For a comprehensive evaluation of Chinese fir’s wood physical properties, we calculated the comprehensive scores of five geographic sources via the PCA. In this respect, significant differences were detected among the five geographic sources. Among them, the comprehensive score of HNYX-T was the highest whereas that of FJYK-P was the lowest (Fig. 4).
Figure 3

Sequence diagram plot of PCA analysis showing the relationship among physical properties of wood.

Full size image

Figure 4

Mean comprehensive score of PCA plot with 95% CI. Different letters (a, b, c, d, e) mean significant difference at 0.05 level.

Full size image More

1.
Brown, J. H. Why are there so many species in the tropics? J. Biogeogr. 41, 8–22 (2014).
PubMed  Article  PubMed Central  Google Scholar 
2.
Classen, A. et al. Temperature versus resource constraints: Which factors determine bee diversity on Mount Kilimanjaro, Tanzania? Glob. Ecol. Biogeogr. 24, 642–652 (2015).
Article  Google Scholar 

3.
Rahbek, C. et al. Humboldt’s enigma: What causes global patterns of mountain biodiversity? Science 365, 1108–1113 (2019).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

4.
Toranza, C. & Arim, M. Cross-taxon congruence and environmental conditions. BMC Ecol. 10, 18 (2010).
PubMed  PubMed Central  Article  Google Scholar 

5.
Gioria, M., Bacaro, G. & Feehan, J. Evaluating and interpreting cross-taxon congruence: Potential pitfalls and solutions. Acta Oecol. 37, 187–194 (2011).
ADS  Article  Google Scholar 

6.
Graham, C. H. et al. The origin and maintenance of montane diversity: Integrating evolutionary and ecological processes. Ecography (Cop.) 37, 711–719 (2014).
Article  Google Scholar 

7.
Westgate, M. J., Tulloch, A. I. T., Barton, P. S., Pierson, J. C. & Lindenmayer, D. B. Optimal taxonomic groups for biodiversity assessment: A meta-analytic approach. Ecography (Cop.) 40, 539–548 (2017).
Article  Google Scholar 

8.
Lomolino, M. V. Elevation gradients of species-density: Historical and prospective views. Glob. Ecol. Biogeogr. 8, 1–2 (2001).
Google Scholar 

9.
McCain, C. M. Global analysis of bird elevational diversity. Glob. Ecol. Biogeogr. 18, 346–360 (2009).
Article  Google Scholar 

10.
Peters, M. K. et al. Predictors of elevational biodiversity gradients change from single taxa to the multi-taxa community level. Nat. Commun. 7, 13736 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

11.
Sundqvist, M. K., Sanders, N. J. & Wardle, D. A. Community and ecosystem responses to elevational gradients: Processes, mechanisms, and insights for global change. Annu. Rev. Ecol. Evol. Syst. 44, 261–280 (2013).
Article  Google Scholar 

12.
Ruggiero, A. & Hawkins, B. A. Why do mountains support so many species of birds? Ecography (Cop.) 31, 306–315 (2008).
Article  Google Scholar 

13.
Mccain, C. M. & Colwell, R. K. Assessing the threat to montane biodiversity from discordant shifts in temperature and precipitation in a changing climate. Ecol. Lett. 14, 1236–1245 (2011).
PubMed  Article  PubMed Central  Google Scholar 

14.
Hawkins, B. A. et al. Energy, water, and broad-scale geographic patterns of species richness. Ecology 84, 3105–3117 (2003).
Article  Google Scholar 

15.
Currie, D. J. Energy and large-scale patterns of animal and plant species richness. Am. Nat. 137, 27–49 (1991).
Article  Google Scholar 

16.
Stein, A., Gerstner, K. & Kreft, H. Environmental heterogeneity as a universal driver of species richness across taxa, biomes and spatial scales. Ecol. Lett. 17, 866–880 (2014).
PubMed  Article  PubMed Central  Google Scholar 

17.
Costanza, J. K., Moody, A. & Peet, R. K. Multi-scale environmental heterogeneity as a predictor of plant species richness. Landsc. Ecol. 26, 851–864 (2011).
Article  Google Scholar 

18.
Vetaas, O. R., Paudel, K. P. & Christensen, M. Principal factors controlling biodiversity along an elevation gradient: Water, energy and their interaction. J. Biogeogr. https://doi.org/10.1111/jbi.13564 (2019).
Article  Google Scholar 

19.
Lande, R. Risks of population extinction from demographic and environmental stochasticity and random catastrophes. Am. Nat. 142, 911–927 (1993).
PubMed  Article  PubMed Central  Google Scholar 

20.
Kaspari, M., Alonso, L. & O’Donnell, S. Three energy variables predict ant abundance at a geographical scale. Proc. R. Soc. B Biol. Sci. 267, 485–489 (2000).
CAS  Article  Google Scholar 

21.
Pianka, E. R. Latitudinal gradients in species diversity: A review of concepts. Am. Nat. 100, 33–46 (1966).
Article  Google Scholar 

22.
Werenkraut, V. & Ruggiero, A. The richness and abundance of epigaeic mountain beetles in north-western Patagonia, Argentina: Assessment of patterns and environmental correlates. J. Biogeogr. 41, 561–573 (2014).
Article  Google Scholar 

23.
R Core Team. R version 3.6.2 ‘Dark and Stormy Night’ (2019). (Accessed 12 December 2019). https://www.r-project.org. 

24.
Blanchet, F. G., Cazelles, K. & Gravel, D. Co-occurrence is not evidence of ecological interactions. Ecol. Lett. 23, 1050–1063 (2020).
PubMed  Article  Google Scholar 

25.
Hodkinson, I. D. Terrestrial insects along elevation gradients: Species and community responses to altitude. Biol. Rev. 80, 489–513 (2005).
PubMed  Article  Google Scholar 

26.
Kampmann, D. et al. Mountain grassland biodiversity: Impact of site conditions versus management type. J. Nat. Conserv. 16, 12–25 (2008).
Article  Google Scholar 

27.
Janzen, D. H. et al. Changes in the arthropod community along an elevational transect in the Venezuelan Andes. Biotropica 8, 193–203 (1976).
Article  Google Scholar 

28.
Sirin, D., Eren, O. & Ciplak, B. Grasshopper diversity and abundance in relation to elevation and vegetation from a snapshot in Mediterranean Anatolia: Role of latitudinal position in altitudinal differences. J. Nat. Hist. 44, 1343–1363 (2010).
Article  Google Scholar 

29.
Alexander, G. & Hilliard, J. R. Altitudinal and seasonal distribution of Orthoptera in the Rocky Mountains of northern Colorado. Ecol. Monogr. 39, 385–432 (1969).
Article  Google Scholar 

30.
Mojica, A. S. & Fagua, G. Estructura de las comunidades de orthoptera (insecta) en un gradiente altitudinal de un bosque andino. Rev. Colomb. Entomol. 32, 200–213 (2006).
Google Scholar 

31.
Grytnes, J. A. Species-richness patterns of vascular plants along seven altitudinal transects in Norway. Ecography (Cop.) 26, 291–300 (2003).
Article  Google Scholar 

32.
McCain, C. M. & Grytnes, J.-A. Elevational gradients in species richness. Encyl. Life Sci. https://doi.org/10.1002/9780470015902.a0022548 (2010).
Article  Google Scholar 

33.
Xu, X. et al. Altitudinal patterns of plant species richness in the Honghe region of China. Pak. J. Bot. 49, 1039–1048 (2017).
Google Scholar 

34.
Kerr, J. T. & Packer, L. Habitat heterogeneity as a determinant of mammal species richness in high-energy regions. Nature 385, 252–254 (1997).
ADS  CAS  Article  Google Scholar 

35.
Röder, J. et al. Heterogeneous patterns of abundance of epigeic arthropod taxa along a major elevation gradient. Biotropica 49, 217–228 (2017).
Article  Google Scholar 

36.
Evans, K. L., Warren, P. H. & Gaston, K. J. Species-energy relationships at the macroecological scale: A review of the mechanisms. Biol. Rev. Camb. Philos. Soc. 80, 1–25 (2005).
PubMed  Article  Google Scholar 

37.
Kissling, W. D., Rahbek, C. & Böhning-Gaese, K. Food plant diversity as broad-scale determinant of avian frugivore richness. Proc. R. Soc. Biol. Sci. 274, 799–808 (2007).
Article  Google Scholar 

38.
Kissling, W. D., Field, R. & Böhning-Gaese, K. Spatial patterns of woody plant and bird diversity: Functional relationships or environmental effects? Glob. Ecol. Biogeogr. 17, 327–339 (2008).
Article  Google Scholar 

39.
Chown, S. L. & Gaston, K. J. Exploring links between physiology and ecology at macro scales: The role of respiratory metabolism in insects. Biol. Rev. 74, 87–120 (1999).
Article  Google Scholar 

40.
de Araújo, W. S. Different relationships between galling and non-galling herbivore richness and plant species richness: A meta-analysis. Arthropod. Plant. Interact. 7, 373–377 (2013).
Article  Google Scholar 

41.
Qian, H. & Kissling, W. D. Spatial scale and cross-taxon congruence of terrestrial vertebrate and vascular plant species richness in China. Ecology 91, 1172–1183 (2010).
PubMed  Article  PubMed Central  Google Scholar 

42.
Burrascano, S. et al. Congruence across taxa and spatial scales: Are we asking too much of species data? Glob. Ecol. Biogeogr. 27, 980–990 (2018).
Article  Google Scholar 

43.
Field, R. et al. Spatial species-richness gradients across scales: A meta-analysis. J. Biogeogr. 36, 132–147 (2009).
Article  Google Scholar 

44.
Giorgis, M. A. et al. Composición florística del Bosque Chaqueño Serrano de la provincia de Córdoba, Argentina. Kurtziana 36, 9–43 (2011).
Google Scholar 

45.
Cabido, M., Funes, G., Pucheta, E., Vendramani, F. & Díaz, S. A chorological analysis of the mountains from Central Argentina. Is all what we call Sierra Chaco really Chaco? Contribution to the study of the flora and vegetation of the Chaco: 12. Candollea 53, 321–331 (1998).
Google Scholar 

46.
Giorgis, M. A. et al. Changes in floristic composition and physiognomy are decoupled along elevation gradients in central Argentina. Appl. Veg. Sci. 20, 553–571 (2017).
Article  Google Scholar 

47.
Cabrera, A. L. Fitogeografia de la República Argentina. In Enciclopedia Argentina de Agricultura y Jardinería Vol. 14 (ed. Kugler, W. F.) 1–42 (ACME, New York, 1976).
Google Scholar 

48.
Giorgis, M. A. et al. Diferencias en la estructura de la vegetación del sotobosque entre una plantación de Pinus taedaL. (Pinaceae) y un matorral serrano (Cuesta Blanca, Córdoba). Kurtziana 31, 39–49 (2005).
Google Scholar 

49.
Martínez, G. A., Arana, M. D., Oggero, A. J. & Natale, E. S. Biogeographical relationships and new regionalisation of high-altitude grasslands and woodlands of the central Pampean Ranges (Argentina), based on vascular plants and vertebrates. Aust. Syst. Bot. 29, 473–488 (2016).
Article  Google Scholar 

50.
QGIS Development Team. QGIS Geographic Information System (Accessed 19 April 2019). (2019).

51.
Kent, M. The description of vegetation in the field. In Vegetation Description and Data Analysis: A Practical Approach (ed. Kent, M.) 65–116 (Wiley-Blackwell, Hoboken, 2012).
Google Scholar 

52.
Catálogo de las plantas vasculares del Cono Sur : (Argentina, Sur de Brasil, Chile, Paraguay y Uruguay). (Missouri Botanical Garden Press, 2008).

53.
Haddad, N., Tilman, D., Haarstad, J., Ritchie, M. & Knops, J. M. N. Contrasting effects of plant richness and composition on insect communities: a field experiment. Am. Nat. 158, 17–35 (2001).
CAS  PubMed  Article  PubMed Central  Google Scholar 

54.
Braun, H. & Zubarán, G. Tettigoniidae (Orthoptera) Species from Argentina and Uruguay (2019).

55.
Carbonell, C. S., Cigliano, M. M. & Lange, C. E. Acridomorph (Orthoptera) Species of Argentina and Uruguay. Version II [2019]. https://biodar.unlp.edu.ar/acridomorph/.

56.
Cigliano, M. M., Braun, H., Eades, D. C. & Otte, D. Orthoptera Species File. Version 5.0/5.0 (2018). http://orthoptera.speciesfile.org.

57.
Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).
Article  Google Scholar 

58.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).
Article  Google Scholar 

59.
Carrara, R., Silvestro, V. A., Cheli, G. H., Campón, F. F. & Flores, G. E. Disentangling the effect of climate and human influence on distribution patterns of the darkling beetle Scotobius pilularius Germar, 1823 (Coleoptera: Tenebrionidae). Ann. Zool. 66, 693–701 (2016).
Article  Google Scholar 

60.
Aisen, S., Werenkraut, V., Márquez, M. E. G., Ramírez, M. J. & Ruggiero, A. Environmental heterogeneity, not distance, structures montane epigaeic spider assemblages in north-western Patagonia (Argentina). J. Insect Conserv. 21, 1–12 (2017).
Article  Google Scholar 

61.
Bilskie, J. Soil Water Status: Content and Potential (Campbell Scientific Inc., Logan, 2001).
Google Scholar 

62.
Tucker, C. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ. 8, 127–150 (1979).
ADS  Article  Google Scholar 

63.
Wang, J., Rich, P. M., Price, K. P. & Dean-Kettle, W. Relations between NDVI, grassland production, and crop yield in the central great plains. Geocarto Int. 20, 5–11 (2005).
Article  Google Scholar 

64.
Oindo, B. O., de By, R. A. & Skidmore, A. K. Interannual variability of NDVI and bird species diversity in Kenya. Int. J. Appl. Earth Obs. Geoinf. 2, 172–180 (2000).
ADS  Article  Google Scholar 

65.
IGN. Modelo Digital de Elevaciones de la República Argentina. (Instituto Geográfico Nacional—Dirección General de Servicios Geográficos—Dirección de Geodesia, 2016).

66.
Riley, S. J., DeGloria, S. D. & Elliot, R. A terrain ruggedness index that quantifies topographic heterogeneity. Intermt. J. Sci. 5, 23–27 (1999).
Google Scholar 

67.
Stein, A. & Kreft, H. Terminology and quantification of environmental heterogeneity in species-richness research. Biol. Rev. 90, 815–836 (2015).
PubMed  Article  Google Scholar 

68.
Tilman, D. & Pacala, S. W. The maintenance of species richness in plant communities. In Species Diversity in Ecological Communities (eds Ricklefs, R. E. & Schulter, D.) 13–25 (University of Chicago Press, Chicago, 1993).
Google Scholar 

69.
Cleveland, W. S., Grosse, E. & Shyu, W. M. Local regresion models. In Statistical Models in S (eds Chambers, J. M. & Hastie, T. J.) 227 (Chapman and Hall, London, 1993).
Google Scholar 

70.
Szewczyk, T. & Mccain, C. M. A systematic review of global drivers of ant elevational diversity. PLoS ONE 11, e0155404 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

71.
Beck, J. et al. Elevational species richness gradients in a hyperdiverse insect taxon: A global meta-study on geometrid moths. Glob. Ecol. Biogeogr. 26, 412–424 (2017).
Article  Google Scholar 

72.
Bolker, B. M. Ecological Statistics: Contemporary Theory and Application (Oxford University Press Inc., Oxford, 2015).
Google Scholar 

73.
Grace, J. B. Structural Equation Modeling and Natural Systems (Cambridge University Press, Cambridge, 2006).
Google Scholar 

74.
Lefcheck, J. S. piecewiseSEM: Piecewise structural equation modelling in r for ecology, evolution, and systematics. Methods Ecol. Evol. 7, 573–579 (2016).
Article  Google Scholar 

75.
Jiménez-Alfaro, B., Chytrý, M., Mucina, L., Grace, J. B. & Rejmánek, M. Disentangling vegetation diversity from climate-energy and habitat heterogeneity for explaining animal geographic patterns. Ecol. Evol. 6, 1515–1526 (2016).
PubMed  PubMed Central  Article  Google Scholar  More

We have introduced and investigated a spatially explicit evolutionary food web model that allows us to explore the distribution of species in space and time, as well as the waxing and waning of species ranges with time. This is the first model that makes it possible to explore these features in the context of trophic networks, showing how they are influenced by competition as well as by predators and prey. Indeed, our model produces empirically well-known patterns in space and time, such as lifetime distributions, species-area relationships, distance decay of similarity, and temporal change of geographic range. On top, we obtain a variety of additional result and gain insights into the mechanisms that generate these patterns. While one might argue that some patterns must emerge trivially in a model like ours, the multitude of patterns emerging together is remarkable.
We find that most species only appear for short times and have small ranges. For species that conquer larger portions of the web we analysed the shape of the range evolution and found that basal species show the empirically observed “hat” pattern more often than species in higher trophic levels. This indicates that the trophic position of a species plays a major role for its range expansion success as well as the shape of its range expansion trajectory over time. To our knowledge, this has not been discussed in the literature so far.
Recently Zliobaite et al.11 analyzed which factors are more correlated to the rise and fall of the range expansion trajectory in fossil data sets of basal mammals. The range expansion follows the so called “hat pattern” that consists of five phases in a species lifetime: origination, expansion, peak, decline and extinction. They found that the temporal location of the peak of the hat pattern is more impacted by competition while the brims are more influenced by abiotic environmental factors. They also compared the range curves of different random walk models with the shape of empirical data and found that a random walk model with competition and environmental factor provides the most realistic looking curves.
Our results for basal species provide an illustration of the mechanism that might lead to these findings. There is one major difference in assumptions: the “environment” in our case is the trophic environment (network structure and abundance distribution). We do not model an abiotic surrounding, yet basal range curves look strikingly hat shaped. A species needs to fit into this trophic environment (network) to first establish a viable population (origination). To successfully increase its range it needs to disperse to neighbouring habitats and be a viable competitor there as well. This leads to the extinction of another species as the dispersing competitor takes its place. As neighbouring basal communities are similar in our systems, the chances are high that the species can spread on a large portion on the grid replacing other species (expansion). This continues until the species has reached the maximum range (peak). It is only a matter of time then until this process is repeated with the species having become the inferior competitor. The species is then successively replaced by a better adapted species (decline). The species thus ages as the network structure changes. At some point the species has vanished on all habitats (extinction). This means that we observe the same dynamics as suggested by Zliobaite et al., but with the trophic and not the abiotic environment as the main driver of the initial increase and later decrease of the range. The truth is probably that both the trophic and abiotic environment are important, as both play an important role in real ecosystems.
Regarding higher trophic species in our system the range expansion curves look more diverse and often do not resemble the hat shape. These species depend on the composition on the layer below. As this layer changes the fate of the higher species changes as well. The emergence of a new prey species that spreads over the network can save a consumer species from extinction. This cannot happen for basal species as these are more or less completely controlled by competition. The studies that we know often deal with basal species, so we are not convinced that the hat pattern is ubiquitous for all species. Future empirical work could focus on predator range expansion and try to find a case where a predator species could regain its range after the emergence of a new prey.
A qualitative difference between the basic trophic layer and higher layers occurred in our data also with respect to the similarity of networks in nearby habitats. The similarity index decays particularly slowly for the basal layer. Theory on distance decay suggests that spatial heterogeneity is a main driver in community turnover in two ways: (1) competitive species sorting along environmental gradients and (2) topological influences that let species with different dispersal abilities experience different landscapes6. As we use a homogeneous landscape we expect the first driver to be non-existent for the basal layer, as all habitats hold the same type and amount of resource. As species do not fundamentally differ in their dispersal abilities and we do not have a heterogeneous spatial topology, the second point is only weakly relevant for the basal species. They have slightly different chances of being chosen for dispersal as we choose the next disperser depending on the biomass density. As we have seen, biomass densities are quite similar for basal species. What remains is a temporal aspect: species that are older can reside on more habitats and have thus a higher chance of being chosen to disperse. To put the cart before the horse, this confirms the theory on distance decay: we expect a much faster decay in a heterogeneous environment, and this is exactly what we observe for higher trophic layers, which experience heterogeneity due to the spatial turnover in basal species composition. A trend to faster decay rates in higher trophic levels was also found in a meta-study7.
A model is always a simplification of reality. Some of the assumptions underlying our model are worth discussing. Species in nature are not restricted to the one-dimensional niche space that we assume. In fact, the original Webworld model characterised species by a large vector of traits31. We, in contrast, characterized species by three traits, all of which are based on body mass. We think that this is the reason why we do not observe super abundant species, but species densities are all of a similar order of magnitude. With a higher-dimensional trait space there must exist more diverse species types and probably also super abundant species, which have a globally optimal trait vector. Nevertheless, our simplification leads to an overall shape of the rank abundance curves that is realistic, as one would expect from a niche apportionment model40. Harpole and Tilman showed that diversity and evenness of grassland communities decreased when niche dimensionality was reduced by adding limiting nutrients to plot experiments41. In turn this indicates that rank abundance curves for less dimensional communities will be flatter. This is in line with the shape of our rank abundance curves.
The probably most interesting trait to add to each species would be its dispersal rate, and to let the dispersal rate evolve. We would expect that this would lead to higher-level species dispersing faster than lower-level species, thus making the differences in range between the different trophic levels smaller. In fact, a previous predator-prey model showed that the predator’s dispersal ability evolved in accordance to spatio-temporal fluctuations of the prey; with higher dispersal rates evolving for larger fluctuations in prey42. However, in the context of food webs the evolution of dispersal is poorly understood and a model like ours would be a good starting point.
Our choice of parameters is guided by the aim to make the model feasible. To be able to perform computer simulations on a large number of habitats we chose a relatively small value for the amount of resource R, so that the number of species of a local food web remained below 25. As we wanted to simulate food webs and not only basal communities we needed to choose a value of the efficiency (lambda ) that allows for the emergence of several trophic layers. In the original Webworld Model, (lambda ) was identified with the proportion of biomass that is passed from one trophic layer to the next. The value of 0.65 that we use here is much larger than the empirically established value of 0.143. However, in the original Webworld model no distinction was made between resident biomass and biomass fluxes. Therefore the variable B was identified with biomass, while its occurrence in Eq. (4) in fact suggests that it represents biomass flux. A further reason for the difference between the model value and the empirical value of (lambda ) is that the model does not take into account energy input due to the below-ground ecological processes. Due to all these simplifications of the model, we do not consider it important to provide an empirical justification of the precise value of the parameter (lambda ), but base its value on the condition that the model yields food webs with several trophic levels.
In contrast to other models, the model used here does not rely on an extrinsic extinction rate that randomly extirpates species that might be well adapted to the network. All extinction events are driven by the trophic dynamics, yet we observe an ongoing species turn over. We thus study the pure food web dynamics without a heterogeneous or fluctuating environment and still observe ecological reasonable species distributions. This indicates that incorporating abiotic environments and their fluctuations is not necessarily needed to study food web dynamics.
Rogge et al. analysed lifetime distributions and SAR curves in a model that is simpler than ours as it does not include population sizes22. The lifetime distributions that we find are considerably steeper (slope (-2.4)) compared to their value around (-1.7). It is also larger than the values reported for empirical findings, which lie between 1 and 2; Newman and Palmer pin them down to (1.7pm 0.3)10. Data for contemporary lifetime distributions show a power-law like shape23,44 with exponents that are in agreement with the exponent of paleological data10. It is noteworthy that there is no consensus whether lifetime distributions follow a power-law or an exponential law, as data often allow for both types of fit, due to (large) uncertainties in fossil data10,45. Exponentials, of course, have a changing slope in a double-logarithmic plot and can thus also be compatible with the exponent observed by us.
Curiously, our model shows steep lifetime distributions even though there is no external random extinction implemented as in other models22. One implication of our value (alpha > 2) is that our distribution has a well-defined mean. This features is shared with exponential distributions.
McPeek argues that lifetime distributions depend on the number and survival time of “transient” species, i.e. species that are on their way to extinction25. He reasons that the time to extinction is elongated for species that are similar, because the inferior competitor holds out longer when it competes with more similar species. If this applies to our type of model, this indicates that species in our system are, despite the one-dimensional niche axis, not as similar as species in the model of Rogge et al.22 that uses the same niche axis, as we observe shorter lifetimes. The difference is that interaction links in22 are binary (presence versus absence), whilst we use Gaussian feeding kernels. The fact that this difference affects the lifetime distributions emphasises the importance of considering details of the trophic interactions. The SAR curves on the contrary are flatter in our model than in the model by Rogge et al.  and in better agreement with empirical data.
O’Sullivan et al.18 found in a competitive metacommunity assembly models a similar collection of macroecological patterns (SAD, range size distribution RSD and SAR) as we did, when regional diversity was near equilibrium. They refer to the work of McGill16 who analysed the assumptions underlying models of macroecological patterns and found that three key ingredients seem to be sufficient for such patterns to emerge. Those are a left skewed SAD, clumping of populations in space, and species distributions in space that are uncorrelated from other species spatial distributions. O’Sullivan et al.18 report that all three ingredients occur in their model and are shaped by regional diversity equilibrium. The closer the system to regional equilibrium the stronger are the observed key patterns (SAD, RSD, Spatial non-correlation). They relate their finding with the theory of ecological structural stability, which revolves around the dynamics on a regional scale. Our communities, in contrast, are trophic communities, operate always near local and regional species equilibrium, i.e. in the regime where O’Sullivan and coauthors18 find the most prominent form of the basic patterns. Comparing the patterns we observe, we also see SADs that are left skewed, and a local clumping of species. We did not analyse the spatial correlation between species. As we have trophic layers of species there will be some correlation between predators and their prey as they can only persist in a habitat if prey is present. In addition to the results obtained by O’Sullivan et al.18, we also derive liefetime distributions, i.e., a paleoecological pattern that also seem to be connected to the metacommunity dynamics. This might indicate that spatial non-correlation is not the most important factor in the mechanisms producing macroecological patterns.
To conclude, our evolutionary food web model produces empirically well studied ecological and paleological patterns. We thus are armed with a valuable tool to broaden our understanding of the mechanisms behind those patterns. Our findings that trophic position influences geographic range and lifetime of a species might motivate further work regarding the interplay of abiotic and trophic factors on range expansion on evolutionary time scales.
More generally, evolutionary models can assist us in forming a deeper knowledge of the processes that lead to what is remnant in fossils. As recently pointed out by Marshall46 in his fifth law of paleobiology, extinction erases information. It is a strength of evolutionary food web models that they allow us to study processes whose extent eludes direct observations. More