Ecology

In the current study we utilize an extensive network and data from citizen science in order to test for among taxa variation in biases and value of information (VoI) in image recognition training data. We use data from the Norwegian Species Observation Service as an example dataset due to the generic nature of this citizen science platform, where all multicellular taxa from any Norwegian region can be reported both with and without images. The platform is open to anyone willing to report under their full real name, and does not record users’ expertise or profession. The platform had 6,205 active contributors in 2021 out of its 17,655 registered users, and currently publishes almost 27 million observations through GBIF, of which 1.08 million with one or more images. Observations have been bulk-verified by experts appointed by biological societies receiving funding for this task, with particular focus on red listed species, invasive alien species, and observations out of range or season. Observations containing pictures receive additional scrutiny, as other users can alert reporters and validators to possible mistaken identifications. An advantage of this particular platform is that no image recognition model has been integrated. This ensures that the models trained in this experiment are not trained on the output resulting from the use of any model, but with identifications and taxonomic biases springing from the knowledge and interest of human observers. Moreover, the platform’s compliance with the authoritative Norwegian taxonomy allows for analyses on taxonomic coverage.In an exploration procedure we determined the taxonomic level of orders to be suitable examples of taxa with a sufficiently wide taxonomic diversity, and enough data in the dataset to be evaluated for models in this experiment. Data collection was done by acquiring taxon statistics and observation data from the Global Biodiversity Information Facility (GBIF), the largest aggregator of biodiversity observations in the world37 for the selected orders, as well as the classes used by Troudet et al.5. The authoritative taxonomy for Norway was downloaded from the Norwegian Biodiversity Information Centre38. In the experimental procedure, models were trained for 12 distinct orders (listed in Fig. 4), artificially restricting these models to different amounts of data. In the data analysis stage, model performances relative to the amount of training data were fitted for each order, allowing the estimation of a VoI. Using the number of observations per species on GBIF, and the number of species known to be present in Norway from the Norwegian Species Nomenclature Database, we calculated relative taxonomic biases.ExplorationInitial pilot runs were done on 8 taxa (see Supplementary Information), using different subset sizes of observations for each species, and training using both an Inception-ResNet-v239 as well as an EfficientNetB340 architecture for each of these subsets. These initial results indicated that the Inception-ResNet-v2 performance (F(_1)) varied less between replicate runs and was generally higher, so subsequent experiments were done using this architecture. The number of observations which still improved the accuracy of the model was found to be between 150 and 200 in the most extreme cases, so the availability of at least 220 observations with images per species was chosen as an inclusion criteria for the further experiment. This enabled us to set aside at least 20 observations per species as a test dataset for independent model analysis.From a Darwin Core Archive file of Norwegian citizen science observations from the Species Observation Service with at least one image33, a tally of the number of such observations per species was generated. We then calculated how many species, with a minimum of 220 such observations, would, at a minimum, be available per taxon if a grouping was made based on each taxon rank level with the constraint of resulting in at least 12 distinct taxa. For each taxonomic level, we calculated how many species having at least 220 such observations were available per taxon when dividing species based on that taxon level. When deciding on the appropriate taxon level to use, we limited the options to taxon levels resulting in at least 12 different taxa.A division by order was found to provide the highest minimum number of species (17) per order within these constraints, covering 12 of the 96 eligible orders. The next best alternative was the family level, which would contain 15 species per family, covering 12 of the 267 eligible families.Data collectionWe retrieved the number of species represented in the Norwegian data through the GBIF API, for all observations, all citizen science observations, and all citizen science observations with images for the 12 selected orders and the classes used by Troudet et al.5. We also downloaded the Norwegian Species Nomenclature Database38 for all kingdoms containing taxa included in these datasets. Observations with images were collected from the Darwin Core Archive file used in the exploration phase, filtering on the selected orders. For these orders, all images were downloaded and stored locally. The average number of images per observation in this dataset was 1.44, with a maximum of 17 and a median of 1.Experimental procedureFor each selected order, a list of all species with at least 220 observations with images was generated from the Darwin Core Archive file33. Then, runs were generated according to the following protocol (Fig. 5):Figure 5Data selection and subdivision. Each run is generated by selecting 17 taxonomically adjacent species per order, and randomly assigning all available images of each selected species to that run’s test-, train- or validation set. Training data are used as input during training, using the validation data to evaluate performance after each training round in order to adjust training parameters during training. The test set is used to measure model performance independently after the model is finalized28. For each subsequent model in that run, training and validation data are reduced by 25% (or slightly less than 25% if not divisible by 4). The test set is not reduced, and used for all models within a run.Full size image

1.

From a list sorted alphabetically by the full taxonomy of the species, a subset of 17 consecutive species starting from a random index was selected. If the end of the list was reached with fewer than 17 species selected, selection continued from the start of the list. The taxonomic sorting ensures that closely related species (belonging to the same family or genus), bearing more similarity, are more likely to be part of the same experimental set. This ensures that the classification task is not simplified for taxa with many eligible species.

2.

Each of the 220+ observations for each species were tagged as being either test, training or validation data. A random subset of all but 200 were assigned to the test set. The remaining 200 observations were, in a 9:1 ratio, randomly designated as training or validation data, respectively. In all cases, images from the same observation were assigned to the same subset, to keep the information in each subset independent from the others. The resulting lists of images are stored as the test set and 200-observation task.

3.

The 200 observations in the training and validation sets were then repeatedly reduced by discarding a random subset of 25% of both, maintaining a validation data proportion of (le)10%. The resulting set was saved as the next task, and this step was repeated as long as the resulting task contained a minimum of 10 observations per species. The test set remained unaltered throughout.

Following this protocol results in a single run of related training tasks with 200, 150, 113, 85, 64, 48, 36, 27, 21, 16 and 12 observations for training and validation per species. The seeds for the randomization for both the selection of the species and for the subsetting of training- and validation datasets were stored for reproducibility. The generation of runs was repeated 5 times per order to generate runs containing tasks with different species subsets and different observation subsetting.Then, a Convolutional Neural Network based on Inception-ResNet-v239 (see the Supplementary Information for model configuration) was trained using each predesignated training/validation split. When the learning rate had reached its minimum and accuracy no longer improved on the validation data, training was stopped and the best performing model was saved. Following this protocol, each of the 12 orders were trained in 5 separate runs containing 11 training tasks each, thus producing a total of 660 recognition models. After training, each model was tested on all available test images for the relevant run.Data analysisThe relative representation of species within different taxa were generated using the number of species present in the GBIF data for Norway within each taxon and the number of accepted species within that taxon present in the Norwegian Species Nomenclature Database38, in line with Troudet et al.5: (R_x = n_x – (n frac{s_x}{s})) where (R_x) is the relative representation for taxon (x), (n_x) is the number of observations for taxon (x), (n) is the total number of observations for all taxa, (s_x) is the number of species within taxon (x), and (s) is the total number of species within all taxa.As a measure of model performance, we use the F(_1) score, the harmonic mean of the model’s precision and recall, given by$$begin{aligned} F_1 = frac{tp}{tp + frac{1}{2}(fp + fn)} end{aligned}$$where (tp), (fp) and (fn) stand for true positives, false positives and false negatives, respectively. The F(_1) score is a commonly used metric for model evaluation, as it is less susceptible to data imbalance than model accuracy28.The value of information (VoI) can be generically defined as “the increase in expected value that arises from making the best choice with the benefit of a piece of information compared to the best choice without the benefit of that same information”32. In the current context, we define the VoI as the expected increase in model performance (F(_1) score) when adding one observation with at least one image. To estimate this, for every order included in the experiment, the increase in average F(_1) score over increasing training task sizes were fitted using the Von Bertalanffy Growth Function, given by$$begin{aligned} L = L_infty (1 – e^{-k(t-t_0)}) end{aligned}$$where (L) is the average F(_1) score, (L_infty) is the asymptotic maximum F(_1) score, (k) is the growth rate, (t) is the number of observations per species, and (t_0) is a hypothetical number of observations at which the F(_1) score is 0. The Von Bertalanffy curve was chosen as it contains a limited number of parameters which are intuitive to interpret, and fits the growth of model performance well.The estimated increase in performance at any given point is then given by the slope of this function, i.e. the result of the differentiation of the Von Bertalanffy Growth Curve, given41 by$$begin{aligned} frac{dL}{dt} = bke^{-kt} end{aligned}$$where$$begin{aligned} b = L_infty e^{kt_0} end{aligned}$$Using this derivative function, we can estimate the expected performance increase stemming from one additional observation with images for each of the species within the order. Filling in the average number of citizen science observations with images per Norwegian species in that order for t, and dividing the result by the total number of Norwegian species within the order, provides the VoI of one additional observation with images for that order, expressed as an average expected F(_1) increase. More

Reference data collectionIn February 2019, we involved forest experts from different regions around the world and organized a workshop to (1) discuss the variety of forest management practices that take place in various parts of the world; (2) explore what types of forest management information could be collected by visual interpretation of very high-resolution images from Google Maps and Microsoft Bing Maps, in combination with Sentinel time series and Normalized Difference Vegetation Index (NDVI) profiles derived from Google Earth Engine (GEE); (3) generalize and harmonize the definitions at global scale; (4) finalize the Geo-Wiki interface for the crowdsourcing campaigns; and (5) build a data set of control points (or the expert data set), which we used later to monitor the quality of the crowdsourced contributions by the participants. Based on the results of this analysis, we launched the crowdsourcing campaigns by involving a broader group of participants, which included people recruited from remote sensing, geography and forest research institutes and universities. After the crowdsourcing campaigns, we collected additional data with the help of experts. Hence, the final reference data consists of two parts: (1) a randomly stratified sample collected by crowdsourcing (49,982 locations); (2) a targeted sample collected by experts (176,340 locations, at those locations where the information collected from the crowdsourcing campaign was not large enough to ensure a robust classification).DefinitionsTable 1 contains the initial classification used for visual interpretation of the reference samples and the aggregated classes presented in the final reference data set. For the Geo-Wiki campaigns, we attempted to collect information (1) related to forest management practices and (2) recognizable from very high-resolution satellite imagery or time series of vegetation indices. The final reference data set and the final map contain an aggregation of classes, i.e., only those that were reliably distinguishable from visual interpretation of satellite imagery.Table 1 Forest management classes and definitions.Full size tableSampling design for the crowdsourcing campaignsInitially, we generated a random stratified sample of 110,000 sites globally. The total number of sample sites was chosen based on experiences from past Geo-Wiki campaigns12, a practical estimation of the potential number of volunteer participants that we could engage in the campaign, and the expected spatial variation in forest management. We used two spatial data sets for the stratification of the sample: World Wildlife Fund (WWF) Terrestrial Ecoregions13 and Global Forest Change14. The samples were stratified into three biomes, based on WWF Terrestrial Ecoregions (Fig. 2): boreal (25 000 sample sites), temperate (35,000 sample sites) and tropical (50,000 sample sites). Within each biome, we used Hansen’s14 Global Forest Change maps to derive areas with “forest remaining forest” 2000–2015, “forest loss or gain”, and “permanent non-forest” areas.Fig. 2Biomes for sampling stratification (1 – boreal, 2 – temperate, 3 – sub-tropical and tropical).Full size imageThe sample size was determined from previous experiences, taking into account the expected spatial variation in forest management within each biome. Tropical forests had the largest sample size because of increasing commodity-driven deforestation15, the wide spatial extent of plantations, and slash and burn agriculture. Temperate forests had a larger sample compared to boreal forests due to their higher fragmentation. Each sample site was classified by at least three different participants, thus accounting for human error and varying expertise16,17,18. At a later stage, following a preliminary analysis of the data collected, we increased the number of sample sites to meet certain accuracy thresholds for every mapped class (aiming to exceed 75% accuracy).The Geo‐Wiki applicationGeo‐Wiki.org is an online application for crowdsourcing and expert visual interpretation of satellite imagery, e.g., to classify land cover and land use. This application has been used in several data collection campaigns over the last decade16,19,20,21,22,23. Here, we implemented a new custom branch of Geo‐Wiki (‘Human impact on Forest’), which is devoted to the collection of forest management data (Fig. 3). Various map overlays (including satellite images from Google Maps, Microsoft Bing Maps and Sentinel 2), campaign statistics and tools to aid interpretation, such as time series profiles of NDVI, were provided as part of this Geo‐Wiki branch, giving users a range of options and choices to facilitate image classification and general data collection. Google Maps and Microsoft Bing Maps include mosaics of very high-resolution satellite and aerial imagery from different time periods and multiple image providers, including the Landsat satellites operated by NASA and USGS as base imagery to commercial image providers such as Digital Globe. More information on the spatial and temporal distribution of very high-resolution satellite imagery can be found in Lesiv et al.24. This collection of images was supplied as guidance for visual interpretation16,20. Participants could analyze time series profiles of NDVI from Landsat, Sentinel 2 and MODIS images, which were derived from Google Earth Engine (GEE). More information on tools can be found in Supplementary file 1.Fig. 3Screenshot of the Geo‐Wiki interface showing a very high-resolution image from Google Maps and a sample site as a 100 mx100 m blue square, which the participants classified based on the forest management classes on the right.Full size imageThe blue box in Fig. 3 corresponds to 100 m × 100 m pixels aligned with the Sentinel grid in UTM projection. It is the same geometry required for the classification workflow that is used to produce the Copernicus Land Cover product for 201511.Before starting the campaign, the participants were shown a series of slides designed to help them gain familiarity with the interface and to train them in how to visually determine and select the most appropriate type of land use and forest management classes at each given location, thereby increasing both consistency and accuracy of the labelling tasks among experts. Once completed, the participants were shown random locations (from the random stratified sample) on the Geo‐Wiki interface and were then asked to select one of the forest management classes outlined in the Definition section (see Table 1 above).Alternatively, if there was either insufficient quality in the available imagery, or if a participant was unable to determine the forest management type, they could skip such a site (Fig. 3). If a participant skipped a sample site because it was too difficult, other participants would then receive this sample site for classification, whereas in the case of the absence of high-resolution satellite imagery, i.e., Google Maps and Microsoft Bing Maps, this sample site was then removed from the pool of available sample sites. The skipped locations were less than 1% of the total amount of locations assigned for labeling. Table 2 shows the distribution of the skipped locations by countries, based on the subset of the crowdsourced data where all the participants agreed.Table 2 Distribution of the skipped locations by countries.Full size tableQuality assurance and data aggregation of the crowdsourced dataBased on the experience gained from previous crowdsourcing campaigns12,19, we invested in the training of the participants (130 persons in total) and overall quality assurance. Specifically, we provided initial guidelines for the participants in the form of a video and a presentation that were shown before the participants could start classifying in the forest management branch (Supplementary file 1). Additionally, the participants were asked to classify 20 training samples before contributing to the campaign. For each of these training samples, they received text‐based feedback regarding how each location should be classified. Summary information about the participants who filled in the survey at the end of the campaign (i.e., gender, age, level of education, and their country of residence) is provided in the Supplementary file 2. We would like to note that 130 participants is a high number, especially taking the complexity of the task into consideration.Furthermore, during the campaign, sample sites that were part of the “control” data set were randomly shown to the participants. The participants received text-based feedback regarding whether the classification had been made correctly or not, with additional information and guidance. By providing immediate feedback, our intention was that participants would learn from their mistakes, increasing the quality and classification accuracy over time. If the text‐based feedback was not sufficient to provide an understanding of the correct classification, the participants were able to submit a request (“Ask the expert”) for a more detailed explanation by email.The control set was independent of the main sample, and it was created using the same random stratified sampling procedure within each biome and the stratification by Global Forest Change maps14 (see “Sample design” section). To determine the size of the control sample, we considered two aspects: (a) the maximum number of sample sites that one person could classify during the entire campaign; (b) the frequency at which control sites would appear among the task sites (defined at 15%, which is a compromise between the classification of as many unknown locations as possible and a sufficient level of quality control, based on previous experience). Our control sample consisted of 5,000 sites. Each control sample site was classified twice by two different experts. When the two experts agreed, these sample sites were added to the final control sample. Where disagreement occurred (in 25% of cases), these sample sites were checked again by the experts and revised accordingly. During the campaign, participants had the option to disagree with the classification of the control site and submit a request with their opinion and arguments. They received an additional quality score in the situation when they were correct, but the experts were not. This procedure also ensured an increase in the quality of the control data set.To incentivize participation and high-quality classifications, we offered prizes as part of the campaign design. The ranking system for the prize competition considered both the quality of the classifications and the number of classifications provided by a participant. The quality measure was based on the control sample discussed above. The participants randomly received a control point, which was classified in advance by the experts. For every control point, a participant could receive a maximum of +30 points (fully correct classification) to a minimum of −30 points (incorrect classification). In the case where the answer was partly correct (e.g., the participant correctly classified that the forest is managed, but misclassified the regeneration type), they received points ranging from 5 to 25.The relative quality score for each participant was then calculated as the total sum of gained points divided by the maximum sum of points that this participant could have earned. For any subsequent data analysis, we excluded classifications from those participants whose relative quality score was less than 70%. This threshold corresponds to an average score of 10 points at each location (out of a maximum of 30 points), i.e., where participants were good at defining the aggregated forest management type but may have been less good at providing the more detailed classification.Unfortunately, we observed some imbalance in the proportion of participants coming from different countries, e.g. there were not so many participants from the tropics. This could have resulted in interpretation errors, even when all the participants agreed on a classification. To address this, we did an additional quality check. We selected only those sample sites where all the participants agreed and then randomly checked 100 sample sites from each class. Table 3 summarizes the results of this check and explains the selection of the final classes presented in Table 1.Table 3 Qualitative analysis of the reference sample sites with full agreement.Full size tableAs a result of the actions outlined in Table 3, we compiled the final reference data set, which consisted of 49,982 consistent sample sites.Additional expert data collectionWe used the reference data set to produce a test map of forest management (the classification algorithm used is described in the next section). By checking visually and comparing against the control data set, we found that the map was of insufficient quality for many locations, especially in the case of heterogeneous landscapes. While several reasons for such an unsatisfactory result are possible, the experts agreed that a larger sample size would likely increase the accuracy of the final map, especially in areas of high heterogeneity and for forest management classes that only cover a small spatial extent. To increase the amount of high-quality training data and hence to improve the map, we collected additional data using a targeted approach. In practice, the map was uploaded to Geo-Wiki, and using the embedded drawing tools, the experts randomly checked locations on the map, focusing on their region of expertise and added classified polygons in locations where the forest management was misclassified. To limit model overfitting and oversampling of certain classes, the experts also added points for correctly mapped classes to keep the density of the points the same. This process involved a few iterations of collecting additional points and training the classification algorithm until the map accuracy reached 75%. In total, we collected an additional 176,340 training points. With the 49,982 consistent training points from the Geo-Wiki campaigns, this resulted in 226,322 (Fig. 4). This two-pronged approach would not have been possible without the exhaustive knowledge obtained from running the initial Geo-Wiki campaigns, including numerous questions raised by the campaign participants. Figure 4 also highlights in yellow the areas of very high sampling density, I.e., those collected by the experts. The sampling intensity of these areas is much higher in comparison with the randomly distributed crowdsourced locations, and these are mainly areas with very mixed forest classes or small patches, in most cases, including plantations.Fig. 4Distribution of reference locations.Full size imageClassification algorithmTo produce the forest management map for the year 2015, we applied a workflow that was developed as part of the production of the Copernicus Global Land Services land cover at 100 m resolution (CGLS-LC100) collection 2 product11. A brief description of the workflow (Fig. 5), focusing on the implemented changes, is given below. A more thorough explanation, including detailed technical descriptions of the algorithms, the ancillary data used, and the intermediate products generated, can be found in the Algorithm Theoretical Basis Document (ATBD) of the CGLS-LC100 collection 2 product25.Fig. 5Workflow overview for the generation of the Copernicus Global Land Cover Layers. Adapted from the Algorithm Theoretical Basis Document25.Full size imageThe CGLS-LC100 collection 2 processing workflow can be applied to any satellite data, as it is unspecific to different sensors or resolutions. While the CGLS-LC100 Collection 2 product is based on PROBA-V sensor data, the workflow has already been tested with Sentinel 2 and Landsat data, thereby using it for regional/continental land cover (LC) mapping applications11,26. For generating the forest management layer, the main Earth Observation (EO) input was the PROBA-V UTM Analysis Ready Data (ARD) archive based on the complete PROBA-V L1C archive from 2014 to 2016. The ARD pre-processing included geometric transformation into a UTM coordinate system, which reduced distortions in high northern latitudes, as well as improved atmospheric correction, which converted the Top-of-Atmosphere reflectance to surface reflectance (Top-of-Canopy). In a further processing step, gaps in the 5-daily PROBA-V UTM multi-spectral image data with a Ground Sampling Distance (GSD) of ~0.001 degrees (~100 m) were filled using the PROBA-V UTM daily multi-spectral image data with a GSD of ~0.003 degrees (~300 m). This data fusion is based on a Kalman filtering approach, as in Sedano et al.27, but was further adapted to heterogonous surfaces25. Outputs from the EO pre-processing were temporally cleaned by using the internal quality flags of the PROBA-V UTM L3 data, a temporal cloud and outlier filter built on a Fourier transformation. This was done to produce consistent and dense 5-daily image stacks for all global land masses at 100 m resolution and a quality indicator, called the Data Density Indicator (DDI), used in the supervised learning process of the algorithm.Since the total time series stack for the epoch 2015 (a three-year period including the reference year 2015 +/− 1 year) would be composed of too many proxies for supervised learning, the time and spectral dimension of the data stack had to be condensed. The spectral domain was condensed by using Vegetation Indices (VIs) instead of the original reflectance values. Overall, ten VIs based on the four PROBA-V reflectance bands were generated, which included: Normalized Difference Vegetation Index (NDVI); Enhanced Vegetation Index (EVI); Structure Intensive Pigment Index (SIPI); Normalized Difference Moisture Index (NDMI); Near-Infrared reflectance of vegetation (NIRv); Angle at NIR; HUE and VALUE of the Hue Saturation Value (HSV) color system transformation. The temporal domain of the time series VI stacks was then condensed by extracting metrics, which are used as general descriptors to enable distinguishing between the different LC classes. Overall, we extracted 266 temporal, descriptive, and textual metrics from the VI times series stacks. The temporal descriptors were derived through a harmonic model, fitted through the time series of each of the VIs based on a Fourier transformation28,29. In addition to the seven parameters of the harmonic model that describe the overall level and seasonality of the VI time series, 11 descriptive statistics (mean, standard deviation, minimum, maximum, sum, median, 10th percentile, 90th percentile, 10th – 90th percentile range, time step of the first minimum appearance, and time step of the first maximum appearance) and one textural metric (median variation of the center pixel to median of the neighbours) were generated for each VI. Additionally, the elevation, slope, aspect, and purity derived at 100 m from a Digital Elevation Model (DEM) were added. Overall, 270 metrics were extracted from the PROBA-V UTM 2015 epoch.The main difference to the original CGLS-LC100 collection 2 algorithms is the use of forest management training data instead of the global LC reference data set, as well as only using the discrete classification branch of the algorithm. The dedicated regressor branch of the CGLS-LC100 collection 2 algorithm, i.e., outputting cover fraction maps for all LC classes, was not needed for generating the forest management layer.In order to adapt the classification algorithm to sub-continental and continental patterns, the classification of the data was carried out per biome cluster, with the 73 biome clusters defined by the combination of several global ecological layers, which include the ecoregions 2017 dataset30, the Geiger-Koeppen dataset31, the global FAO eco-regions dataset32, a global tree-line layer33, the Sentinel-2 tiling grid and the PROBA-V imaging extent;30,31 this, effectively, resulted in the creation of 73 classification models, each with its non-overlapping geographic extent and its own training dataset. Next, in preparation for the classification procedure, the metrics of all training points were analyzed for outliers, as well as screened via an all-relevant feature selection approach for the best metric combinations (i.e., best band selection) for each biome cluster in order to reduce redundancy between parameters used in the classification. The best metrics are defined as those that have the highest separability compared to other metrics. For each metric, the separability is calculated by comparing the metric values of one class to the metric values of another class; more details can be found in the ATBD25. The optimized training data set, together with the quality indicator of the input data (DDI data set) as a weight factor, were used in the training of the Random Forest classifier. Moreover, a 5-fold cross-validation was used to optimize the classifier parameters for each generated model (one per biome).Finally, the Random Forest classification was used to produce a hard classification, showing the discrete class for each pixel, as well as the predicted class probability. In the last step, the discrete classification results (now called the forest management map) are modified by the CGLS-LC100 collection 2 tree cover fraction layer29. Therefore, the tree cover fraction layer, showing the relative distribution of trees within one pixel, was used to remove areas with less than 10% tree cover fraction in the forest management layer, following the FAO definition of forest. Figure 6 shows the class probability layer that illustrates the model behavior, highlighting the areas of class confusion. This layer shows that there is high confusion between forest management classes in heterogeneous landscapes, e.g., in Europe and the Tropics while homogenous landscapes, such as Boreal forests, are mapped with high confidence. It is important to note that a low probability does not mean that the classification is wrong.Fig. 6The predicted class probability by the Random Forest classification.Full size image More

Data collectionThe acoustic recordings were collected during 2017 off the Changhua coast (24° 4.283 N/120° 19.102 E) (Fig. 5) by deploying a passive acoustic monitoring (PAM) device from Wildlife Acoustics, which was an SM3M recorder moored at a depth of 18–20 m. The hydrophone recorded continuously with a sampling frequency of 48 kHz and a sensitivity of −164.2 dB re:1 v/µPa. The acoustic files were recorded in the.WAV format with a duration of 60 minutes. The hydrophone setup prior to deployment is shown in Fig. 6. Table 2 contains the details for the monitoring period with the corresponding season and the number of hours of recordings each season used in this study. Studies have shown that the presence of seasonal chorusing at this monitoring site in the frequency range of 500–2500 Hz is caused by two types of chorusing15,38, with chorusing starting in early spring, reaching a peak in summer, and starting to decline late autumn, and silencing in winter6. Previous studies6,15,38 at this monitoring site have derived the details of two types of fish calls responsible for chorusing (Type 1 and Type 2); Supplementary Fig. 1 shows the spectrogram, waveform, and power spectrum density of the individual calls. Supplementary Table 1 tabulated are the acoustic features of the two call types. The monitoring region, Changhua, lies in the Eastern Taiwan Strait (ETS). The ETS is ~350 km in length and ~180 km wide64. The ETS experiences diverse oceanographic and climatic variations influenced by monsoons in summer and winter65 and extreme events caused by tropical storms, typhoons in summer, and wind/cold bursts occurring in winter66,67,68.Fig. 5: Study area located off the Taiwan Strait.Map of the Changhua coast located in Taiwan Strait, Taiwan depicting the deployed passive acoustic monitoring recorder at site A1. The map was produced in Matlab 9.11 (The Mathworks, Natick, MA; http://www.mathworks.com/) using mapping toolbox function geobasemap(). Full global basemap composed of high-resolution satellite imagery hosted by Esri (https://www.esri.com/).Full size imageFig. 6: Setup of the SM3M submersible recorder.SM3M recorder fastened to the steel frame (length and breadth = 1.22 m, height = 0.52 m) with plastic cable zip ties prior to deployment.Full size imageTable 2 Passive acoustic monitoring device specifications and monitoring duration during different seasons.Full size tableAcoustic data analysisThe acoustic data were analyzed using the PAMGuide toolbox in Matlab60. The seasonal spectrograms were computed with an FFT size of 1024 points and a 1 s time segment averaged to a 60 s resolution. The sound pressure levels (SPL) were computed in the frequency band of 500–3500 Hz and programmed to provide a single value every hour, thus resulting in 984, 1344, and 1440 data points in spring, summer, and winter, respectively (Supplementary Data 1).Determining the regularity and complexity with the complexity-entropy planeThe complexity-entropy plane was utilized in this study to quantify the structural statistical complexity and the regularity in the hourly acoustical and seasonal SPL time series data. The C-H plane is a 2D plane representation of the permutation entropy on the horizontal axis that quantifies the regularity, and the vertical axis is represented by the statistical complexity quantifying the correlation structure in the temporal series.For a given time series ({{x(t)}}_{t=1}^{N}), the N’ ≡ N − (m − 1) the values of the vectors for the length m  > 1 are ranked as$${X}_{s}=left({x}_{s-(m-1)},{x}_{s-(m-2)},ldots ,{x}_{s}right),s=1,ldots ,,{N}^{{prime} }$$
(1)
Within each vector, the values are reordered in the ascending order of their amplitude, yielding the set of ordering symbols (patterns) ({r}_{0},{r}_{1},ldots ,{r}_{m-1}) such that$${x}_{s-{r}_{0}}le {x}_{s-{r}_{1}}le ..,..le {x}_{s-{r}_{(m-1)}}$$
(2)
This symbolization scheme was introduced by Bandt and Pompe69. The scheme performs the local ordering of a time series to construct a probability mass function (PMF) of the ordinal patterns of the vector. The corresponding vectors (pi ={r}_{0},{r}_{1},ldots ,{r}_{(m-1)}) may presume any of the m! possible permutations of the set ({{{{{mathrm{0,1}}}}},ldots ,m-1}) and symbolically represent the original vector. For instance, for a given time series {9, 4, 5, 6, 1,…} with length m = 3, provides 3! different order patterns with six mutually exclusive permutation symbols are considered. The first three-dimensional vector is (9, 4, 5), following the Eq. (1), this vector corresponds to ((,{x}_{s-2},{x}_{s-1},{x}_{s})). According to Eq. (2), it yields ({x}_{s-1}le {x}_{s}le {x}_{s-2}). Then, the ordinal pattern satisfying the Eq. (2) will be (1, 0, 2). The second 3-dimensional vector is (4, 5, 6), and (2, 1, 0) will be its associated permutation, and so on.The permutation entropy (H) of order m ≥ 2 is defined as the Shannon entropy of the Brandt-Pompe probability distribution p(π)69$$Hleft(mright)=,-{mathop{sum}limits _{{pi }}}pleft(pi right){{{log }}}_{2}p(pi )$$
(3)
where ({pi }) represents the summation over all possible m! permutations of order m, (p(pi )) is the relative frequency of each permutation (pi), and the binary logarithm (base of 2) is evaluated to quantify the entropy in bits. H(m) attains the maximum ({{log }}m!) for (p(pi )=1/m!). Then the normalized Shannon entropy is given by$$0le H(m)/{{{{{rm{ln}}}}}},m!le 1$$
(4)
where the lower bound H = 0 corresponds to more predictable signals with fewer fluctuations, an either strictly increasing or decreasing series (with a single permutation), and the upper bound H = 1 corresponds to an unpredictable random series for which all the m! possible permutations are equiprobable. Thus, H quantifies the degree of disorder inherent in the time series. The choice of the pattern length m is essential for calculating the appropriate probability distribution, particularly for m, which determines the number of accessible states given by m!70,71. As a rule of thumb, the length T of the time series must satisfy the condition T (gg) m! in order to obtain reliable statistics, and for practical purposes, Bandt and Pompe suggested choosing m = 3,…,7 69.The statistical complexity measure is defined with the product form as a function of the Bandt and Pompe probability distribution P associated with the time series. (Cleft[Pright]) is represented as33$$Cleft[Pright]=frac{J[P,U]}{{J}_{{max }}}{H}_{s}[P]$$
(5)
where ({H}_{s}left[Pright]=Hleft[Pright]/{{log }}m!) is the normalized permutation entropy. (J[P,U]) is the Jensen divergence$$Jleft[P,Uright]=left{Hleft[frac{P+U}{2}right]-frac{H[P]}{2}-frac{H[U]}{2}right}$$
(6)
which quantifies the difference between the uniform distributions U and P, and ({J}_{{max }})is the maximum possible value of (Jleft[P,Uright]) that is obtained from one of the components of P = 1, with all the other components being zero:$$Jleft[P,Uright]=-frac{1}{2}left[frac{m!+1}{m!}{{log }}left(m!+1right)-2{{log }}left(2m!right)+{{log }}(m!)right]$$
(7)
For each value of the normalized permutation entropy (0le {H}_{s}[P]le 1) there is a corresponding range of possible statistical complexity (Cleft[Pright]) values. Thus, the upper (({C}_{{max }})) and lower ((C_{{min }})) complexity bounds in the C-H plane are formed. The periodic sequences such as sine and series with increasing and decreasing (with ({H}_{s}[P]=0)) and completely random series such as white noise (for which (Jleft[P,Uright]=0) and ({H}_{s}[P]=1)) will have zero complexity. Furthermore, for each given value of the (0le {H}_{s}[P]le 1), there exists a range of possible values of the statistical complexity, ({C}_{{min }}le C[P]le {C}_{{max }}). The procedure for evaluating the complexity bounds ({C}_{{min }}) and ({C}_{{max }}) is given in Martin et al.72. We utilized the R package ‘statcomp’73 to evaluate the statistical complexity (C) and the permutation entropy (H) using the command global-complexity() for the order m = 6, and the command limit_curves(m, fun = ‘min/max’) was utilized to evaluate the complexity boundaries ({C}_{{min }}) and ({C}_{{max }}). In this study, we constructed two C-H planes: (1) C and H was computed for each hourly acoustic file during different seasons. The resulting lengths of C and H during spring, summer, and autumn-winter are similar to the number of hours in the particular season (Table 2). (2) C and H was computed every 4–5 days for the seasonal SPL. The resulting length of C and H obtained during spring was 9 points (each value of C and H for every 109 h), and in summer and autumn-winter was 12 points (each value of C and H for every 112 and 120 h).Determining predictability and dynamics (linear/nonlinear) using EDMIn this study, we utilized EDM to quantify the predictability (forecasting) and distinguish between the linear stochastic and nonlinear dynamics in the seasonal soundscape SPL. EDM involves phase-space reconstruction via delay coordinate embeddings to make forecasts and to determine the ‘predictability portrait’ of time series data40. The first step in EDM is to determine the optimal embedding dimension (E), and this is obtained using a method based on simplex projection41. The simplex projection is carried out by dividing the dataset into two equal parts, of which the first part is called the library and the other part is called the target. The library set is used to build a series of non-parametric models (known as predictors) for the one step ahead predictions for the E varying between 1 and 10. Then the model’s accuracies are determined when the model is applied to the target dataset and the prediction skill (⍴) for the actual and predicted datasets is measured. The embedding dimension with the highest predictive skill is the optimal E.For the appropriate optimal E chosen, the predictability profile of the time series data is obtained by evaluating ⍴ at Tp = 1, 2, 3, … steps ahead. The flat prediction profile of the ⍴–Tp curve indicates that the time series is purely random (low ⍴) or regularly oscillating (high ⍴). In contrast, a decreasing ⍴ as Tp increases may reject the possibility of an underlying uncorrelated stochastic process and indicate the presence of low-dimensional deterministic dynamics. However, the concern with the predictability profile is that it may exhibit predictability even if time series are purely stochastic (such as autocorrelated red noise). Hence, a nonlinear test can be performed by using S-maps (sequential locally weighted global linear maps) to distinguish between linear stochastic and nonlinear dynamics in the time series dataset by fitting a local linear map. S-maps similar to simplex projects provide the forecasts in phase-space by quantifying the degree to which points are weighted when fitting the local linear map, which is given by the nonlinear localization parameter θ. When θ = 0, the entire library set will exhibit equal weights regardless of the target prediction, which mathematically resembles the model of a linear autoregressive process. In contrast, if θ  > 0, the forecasts of the library provided by the S-map depend on the local state of the target prediction, thus producing large weights, and the unique local fittings can vary in phase-space to incorporate nonlinear behavior. Consequently, if the (⍴–θ) dynamics profile shows the highest ⍴ at θ = 0 that is reduced as θ increases, it represents linear stochastic dynamics. If the ⍴ achieves the highest value at θ  > 0, then the dynamics are represented by nonlinear dynamics.In this study, the R package “rEDM”74 was used to evaluate the optimal E, prediction profile (⍴–Tp), and dynamics profile (⍴–θ) for the seasonal SPL dataset. While evaluating these entities, the data points are equally into two parts, and sequentially the first half is chosen as the library set and the other as the target set. The length of the library and the target set for spring, summer, and autumn-winter are 492, 672, and 720. The command EmbedDimension() was used to determine the forecast skill for the E ranging from 1 to 10 and the optimal E with the highest forecast skill (Supplementary Fig. 2) was chosen. In this study, we found that for all seasons, the optimal E was 2. The (⍴–Tp) was evaluated for Tp varying between 1 and 100 using the command PredictInterval() and the (⍴–θ) was evaluated using the command PredictNonlinear() for θ = 0, 0.0001, 0.0003, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5,0.75, 1.0, 1.5, 2, and 3 to 20.StatisticsThe nonparametric Kruskal–Wallis test, followed by post hoc Bonferroni’s multiple comparisons, was used to test differences in the seasonal H and C that were obtained directly from the hourly acoustic data during chorusing hours, as well as the H and C obtained for the seasonal SPL and the seasonal forecast skill. The statistical calculations were performed using the R package “agricolae” 75. More

Our model focused on the trophic interactions among CoTS and two groups of coral within a feedback loop with natural and anthropogenic forcing. Our model draws on accepted features of the published dynamics described by Morello et al.37, Condie et al.28 and Condie et al.17, but is a substantial advance in terms of adding spatial structure and coupling with climate variables. Here we have resolved a fine spatiotemporal model structure, developed a novel recruitment formulation for CoTS, integrated tactical management control dynamics and incorporated the impact of broad-scale drivers upon the population dynamics of corals and CoTS at the local scale. Our model is formally fitted to a subset of the CoTS control program data described by Westcott et al.12. We operationalised our model as a tactical and strategic tool to inform how CoTS management strategies interact with alternative disturbance and ecological realisations at the sub-reef scale, the scale at which management operates.DataWe fitted our model to a subset of four reefs from the dataset described by Westcott et al.12, which were consistently and intensively managed (for a map with reef locations see Fig. 2). We restricted our focus to a subset to avoid parametrisation of reef and management site dynamics. Thus, ~39% of site visits were concentrated over the 13 management sites we considered, with a mean of 20.73 ± 5.5 (mean ± standard deviation) visits across the time series relative to a mean visitation rate of 12.23 ± 4.7 (mean ± standard deviation) for the rest of the sites. Each reef in the subset contained two or more management sites where each site was visited at least 18 times. The subset was used because it contained sufficient data for estimating the 11 model parameters for each management site. Across included sites were a range of CoTS densities, coral abundances and disturbance histories12,72,73. Given the intensity with which these sites were managed, they therefore provided us with a valuable opportunity to formally fit the interactions between management intervention, coral abundance and CoTS dynamics in the presence of regional sequential bleaching events.Model spatial structure and ecological componentsSpatially, we considered a circular 300 km region of the Great Barrier Reef centred between Cairns and Cape Tribulation, and resolved at a daily timescale and a sub-reef spatial scale, matching the scale at which observed data were resolved12,19. Reefs were randomly generated as points to capture possible spatial correlation in disturbance impacts between nearby reefs, as well as to allow variability in reef locations. Coral, CoTS and disturbance dynamics within the management sites of each reef were resolved relative to a 1 ha focal region. That is, each management site was captured as a 1 ha area representative of the whole site. In the Pacific, Acanthaster spp. disproportionally target faster-growing corals, predominantly Acropora, Pocillopora and Montipora22. Coral taxa characterised by slow growth rates and massive morphologies, such as Porites, are generally consumed less than expected based on their abundance22 and are thus non-preferred prey. The two modelled coral groups were the fast-growing favoured prey items of CoTS, and the slower-growing non-preferred prey. Processes resolved in the model included reproduction, density dependence, the effect of bleaching and cyclonic disturbances on corals and the impact of manual control (culling) upon CoTS and coral dynamics.CoTS population structureWe used an age-structured approach to model CoTS population dynamics. We defined our age classes to encapsulate plausible size-at-age variation due to plastic growth. This was achieved through linking catch size classes of the management control program19 to age classes through size-age relationships developed from observations spanning multiple environmental realisations, manipulated scenarios and methodologies55,70,74,75. Delayed growth in juvenile CoTS due to deferral of their switch to coral prey or composition of their pre-coral diet, may induce variability in the size-at-age of juveniles52,53. However, the population-level consequences of prolonged juvenile phases are not easily observed nor understood. For example, juveniles are subject to high mortality rates in situ, delayed growth may reduce lifetime fitness and there have been no observations of juveniles during spawning periods that would indicate protracted juvenile phases55,56,57. Consequently, suggests size-at-age is—due to an early life history mortality bottleneck or otherwise—predominantly concordant with growth curves of the literature55,70,74,75 and the size classes we have used here. Age classes comprised annual 0, 1, 2 and 3+ groups, with 3+ being an absorbing class – once there, they stay there. Age-0 ( ; 32.5)). This induced a slope change in the relationship between maximum wind velocity and its radius at a wind velocity of 32.5 m.s−1 (≥ category 3 intensities). However, whilst maximum wind velocity was modelled to determine ({d}_{{{{{{rm{m}}}}}}}), the overall size of the cyclone was uncorrelated with its intensity. The overall size was uniformly sampled from 130 to 460 km diameter which allowed for the potential of complete focal area coverage and for a range of intensity-size relationships to be captured. Given a cyclone footprint of radius ({d}_{0}) (km), wind velocity, (V) (m.s−1), at a distance, (d) (km), was interpolated104 through:$$Vleft(dright)=left{begin{array}{c}{V}_{0}+left({V}_{{{{{{rm{m}}}}}}}-{V}_{0}right){left(frac{sqrt{{d}_{0}}-sqrt{d}}{sqrt{{d}_{0}}-sqrt{{d}_{{{{{{rm{m}}}}}}}}}right)}^{alpha },,dge {d}_{{{{{{rm{m}}}}}}}\ {V}_{{{{{{rm{m}}}}}}}, , d ; < ; {d}_{{{{{{rm{m}}}}}}}end{array}right.$$ (32) The distance from the cyclone centre to the reef perimeter, (D) (km), is calculated through:$$D=sqrt{{left({x}_{{{{{{rm{rf}}}}}}}-{r}_{1}-{x}_{{{{{{rm{cyc}}}}}}}right)}^{2}+{left({x}_{{{rm{rf}}}}-{r}_{1}-{y}_{{{{{{rm{cyc}}}}}}}right)}^{2}}$$ (33) Thus, given a reef strike occurs ((sqrt{{d}_{0}}-sqrt{d}ge 0) required from non-integer (alpha)), the wind velocity experienced at said reef due to the tropical cyclone was calculated as (Vleft(Dright)). Wind velocity was subsequently categorised and damage to reef zone corals calculated as per Supplementary Table 4.We resolved stochasticity in cyclone dynamics in projection scenarios. In projected scenarios cyclone arrivals, locations and intensities were probabilistically sampled and their inflicted damage upon coral communities sampled from damage ranges. Cyclone locations, their footprints, intensity ranges and corresponding damage ranges were sampled from uniform distributions. Cyclone arrivals were sampled from a Poisson distribution and considered in scenarios from 2018 to 2029. Projections were averaged over 80 simulations to capture mean dynamics and bound trajectory uncertainty due to said stochasticity.Our cyclone model was calibrated to parameters sourced from the literature (Supplementary Tables 4-5). This was necessary since our data time series did not encompass a cyclone event and/or impacts upon a reef and cyclone-induced mortality is typically a key coral mortality source30. Consequently, we were unable to validate the impacts of cyclones through formal estimation in our model. However, our endeavours to source parameters from empirical and modelling studies in conjunction with our formulation allowed us to plausibly capture the cumulative outcomes of a cyclone event at discrete locations. Our cyclone model offers a limited complexity approach that is empirically grounded to simply resolve cyclone impacts in local-scale models without the need to be coupled to a regional-scale model.Cyclones, induced thermal stress and tactical managementThe occurrence of cyclone events was modelled to directly interact with both management interventions and thermal stress events. Cyclones were assumed to realistically preclude co-occurring co-located management interventions. This was such that a management site control visit was abandoned if a cyclone preceded or was forecast within five days of a control voyage. The later interaction of cyclones with thermal stress events operated through an induced thermal cooling of sea surface temperatures (SST) at impacted locations.In the case of the overlapping cyclone and thermally induced bleaching events, we first accounted for cyclone impacts. This was because, in addition to physical damage to corals, cyclones have the potential for regional-scale cooling of SST which can reduce coral bleaching43,107. To capture this interaction, we resolved the duration108,109 and amplitude107 of tropical cyclone-induced cooling. We captured this interaction through Degree Heating Weeks (DHW) which is a useful metric for the accumulated thermal stress experienced by corals94.The duration of tropical cyclone-induced cooling was modelled through a temporal-SST response curve consistent with the work of Lloyd and Vecchi108 and Vincent et al.109. Cooling rapidly occurs once a tropical cyclone arrives at a location and decays in an asymptotic manner over a period of ~40–60 days108,109. Temperatures however do not return to pre-cyclone levels and plateau at ~1/4 of the cooling signal amplitude below pre-cyclone levels108,109. We expressed this cooling response curve as it related to bleaching-induced coral mortality through DHWs.We based the average expected DHW cooling signal on the work of Carrigan and Puotinen107. This was achieved through scaling the difference in amplitude of overlapping thermal stress-tropical cyclone events and thermal stress only events—a cooling signal amplitude of ({{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}} sim 1.5) DHW. Consistent with the model of Carrigan and Puotinen107, we then resolved cooling within the radius of gale-force winds (category 1, 17 m.s−1) to model tropical cyclone-induced cooling. Depending on the size of the tropical cyclone, this meant that an individual cyclone would not necessarily cool all reefs within the model region. However, the culmination of multiple cyclones may have limited bleaching exposure for corals across the region107.We did not treat the cooling consequences of multiple cyclones additively nor the complex interplay of oceanic feedbacks upon cyclone intensity and cooling. Such processes were beyond the scope of our study and model. If multiple cyclones occurred within our model, then the cooling signal timeline was re-initialised at impacted reefs for the last tropical cyclone at said location. Non-impacted reefs maintained the timeline for the decay of the cooling signal originating from their previous tropical cyclone interaction.Once a tropical cyclone impacted a reef, the duration of the induced cooling signal was modelled. Price et al.110 found that cooling decays exponentially which is reflective of the recovery of SST following tropical cyclones as demonstrated by Lloyd and Vecchi108 and Vincent et al.109. We operationalised the exponential functional form in conjunction with the decay timelines of Lloyd and Vecchi108 and Vincent et al.109 and the DHW amplitude of Carrigan and Puotinen107. We modelled the level of cooling ({{{{{{rm{DHW}}}}}}}_{{{{{{rm{cool}}}}}}}) after ({d}_{{{{{{rm{postTC}}}}}}}) days post-cyclone event by:$${{{{{{rm{DHW}}}}}}}_{{{{{{rm{cool}}}}}}}left({d}_{{{{{{rm{postTC}}}}}}}right)=frac{1}{4}{{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}}+frac{frac{3}{4}{{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}}}{{e}^{{d}_{{{{{{rm{postTC}}}}}}}/10}}$$ (34) This ensured that once a reef experienced a tropical cyclone event, the cooling signal initialised at ({{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}}) and decayed to (sim frac{1}{4}{{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}}) after 40–60 days108,109. The rate of decay was given by the e-folding time (days required for the cooling signal to be reduced by a factor of (e)) which we took to be 10. This is consistent with the results of Price et al.110, Lloyd and Vecchi108 and Vincent et al.109 who found e-folding times ranging from 5 through to 20 days. Thermally induced bleaching mortality of corals was computed after cyclone physical damage and cooling had been accounted for.Formal model fittingWe formally fitted our coral-CoTS model simultaneously to coral cover data, catch-per-unit-effort data and catch numbers obtained from the management control program with dive effort (minutes) treated as an input (visits summarised in Supplementary Table 7)12. Simultaneously fitting CoTS and coral dynamics at concurrent locations was useful here as it allowed for coral cover trajectories to help inform local CoTS abundance (sensu CoTS feeding vs. coral trajectories63,79 and local site fidelity24). Our model also used Long Term Monitoring Program (LTMP) data (based on manta tows and provided by the Australian Institute of Marine Science) which provides an independent index of relative abundance of CoTS. This was such that our model here was developed and parametrised based on an earlier version37,111 which did not use CPUE information but was fitted to the LTMP data on CoTS relative abundance, as well as the corresponding coral cover, to estimate a number of CoTS-coral interaction parameters used in the present model (Supplementary Table 3).Fitting and estimation of our model were achieved through Maximum Likelihood Estimation (MLE). Our objective function was the outcome of combining the negative log-likelihood contributions arising from fitting the model to multiple sets of location-specific data, across a range of environmental and ecological realisations, in conjunction with penalty terms. Specifically, we fitted coral cover (data series ({x}^{{{{{{rm{Coral}}}}}}})) and CoTS CPUEs (data series ({x}^{{{{{{rm{CoTS}}}}}}})) at each management site which contained ({n}_{{{{{{rm{Coral}}}}}}}) and ({n}_{{{{{{rm{CoTS}}}}}}}) data points respectively. This involved fitting parameters that were specific to management sites (e.g. thermal stress - DHW), reefs (e.g. recruitment variability) as well as those that were common amongst reefs (e.g. CoTS consumption rates). A parametrisation that optimised one contribution was unlikely to optimise all contributions and hence we obtained a parametrisation across all reefs and sub-regions. For a modelled catch of (N) (sum of catches across age classes), a catchability coefficient (a constant of proportionality) of ({q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}), and data standard deviation of ({sigma }_{{{{{{rm{LL}}}}}}}) our likelihood contribution arising from a management site CPUEs was given by:$$-{{log }}{{{{{rm{L}}}}}}left({q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}N,{{sigma }_{{{{{{rm{LL}}}}}}}}^{2}{{{{{rm{|}}}}}}{x}_{i}^{{{{{{rm{CoTS}}}}}}}right) = {n}_{{{{{{rm{CoTS}}}}}}},{{{{{rm{ln}}}}}}left({sigma }_{{{{{{rm{LL}}}}}}}right)+{sum }_{i=1}^{{n}_{{{{{{rm{CoTS}}}}}}}}frac{{left({{{{{rm{ln}}}}}}left({x}_{i}^{{{{{{rm{CoTS}}}}}}}right)-{{{{{rm{ln}}}}}}left({q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}{N}_{i}right)right)}^{2}}{2{{sigma }_{{{{{{rm{LL}}}}}}}}^{2}}$$ (35) From which the data series variance and catchability coefficient were computed for the maximum likelihood estimate. The derived variance and the catchability were respectively computed as per:$${sigma }_{{{{{{rm{LL}}}}}}}=sqrt{frac{1}{{n}_{{{{{{rm{CoTS}}}}}}}}{sum }_{i=1}^{{n}_{{{{{{rm{CoTS}}}}}}}}{left({{{{{rm{ln}}}}}}left({x}_{i}^{{{{{{rm{CoTS}}}}}}}right)-{{{{{rm{ln}}}}}}left({q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}right)right)}^{2}}$$ (36) and$${q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}=frac{1}{{n}_{{{{{{rm{CoTS}}}}}}}}{sum }_{i=1}^{{n}_{{{{{{rm{CoTS}}}}}}}}left({{{{{rm{ln}}}}}}left({x}_{i}^{{{{{{rm{CoTS}}}}}}}right)-{{{{{rm{ln}}}}}}left({N}_{i}right)right)$$ (37) Similarly, the likelihood contribution arising from fitting to a management site coral cover with standard deviation ({sigma }_{{Coral}}) was described by:$$-{{log }}{{{{{rm{L}}}}}}left(frac{{C}_{y,d}^{{{{{{rm{f}}}}}}}+{C}_{y,d}^{{{{{{rm{s}}}}}}}}{{K}^{{{{{{rm{coral}}}}}}}},{{sigma }_{{{{{{rm{Coral}}}}}}}}^{2}{{{{{rm{|}}}}}}{x}_{i}^{{{{{{rm{Coral}}}}}}}right) = {n}_{{{{{{rm{Coral}}}}}}},{{{{{rm{ln}}}}}}left({sigma }_{{{{{{rm{Coral}}}}}}}right)+{sum }_{i=1}^{{n}_{{{{{{rm{Coral}}}}}}}}frac{{left({ln}left({x}_{i}^{{{{{{rm{Coral}}}}}}}right)-left(frac{{C}_{y,d}^{{{{{{rm{f}}}}}},i}+{C}_{y,d}^{{{{{{rm{s}}}}}},i}}{{K}^{{{{{{rm{coral}}}}}}}}right)right)}^{2}}{2{{sigma }_{{{{{{rm{Coral}}}}}}}}^{2}}$$ (38) Where the standard deviation was given by:$${sigma }_{{{{{{rm{Coral}}}}}}}=sqrt{frac{1}{{n}_{{{{{{rm{Coral}}}}}}}}{sum }_{i=1}^{{n}_{{{{{{rm{Coral}}}}}}}}{left({{{{{rm{ln}}}}}}left({x}_{i}^{{{{{{rm{Coral}}}}}}}right)-{{{{{rm{ln}}}}}}left(frac{{C}_{y,d}^{{{{{{rm{f}}}}}},i}+{C}_{y,d}^{{{{{{rm{s}}}}}},i}}{{K}^{{{{{{rm{coral}}}}}}}}right)right)}^{2}}$$ (39) We computed the negative log-likelihood objective function by summing the contributions from all management sites across considered reefs.Fitting was conducted through the modelling language Automatic Differentiation Model Builder (ADMB) which implements a Quasi-Newton optimisation algorithm for estimation of parameters and provides Hessian based estimation of standard errors112. Penalty terms were added to our likelihood function to integrate a prior understanding of system dynamics and to reduce model variability. Penalty terms encompassed recruitment variability and the magnitude of catches observed in the data.Recruitment was expressed through recruitment deviations, ({r}_{y}), given a standard deviation of ({sigma }_{{{{{{rm{R}}}}}}}) about underlying modelled recruitment (sum of self-recruitment and immigration sources described previously). The recruitment variability negative log-likelihood penalty contribution was given by:$$-{{log }}{{{{{rm{L}}}}}}left(0,{sigma }_{{{{{{rm{R}}}}}}}^{2}{{{{{rm{|}}}}}}{r}^{{{{{{rm{rec}}}}}}}right)={sum }_{y=1}^{{{{{{rm{#Years}}}}}}}{sum }_{{{{{{rm{reef}}}}}}=1}^{{{{{{rm{#Reefs}}}}}}}{r}_{y,{{{{{rm{reef}}}}}}}^{{rec}}/2{sigma }_{{{{{{rm{R}}}}}}}^{2}$$ (40) An additional penalty term for model deviations from the magnitude of observed catches was encompassed. This was such that a constant of proportionality relating modelled catches to observed catches tended to one. For an allowed standard deviation of ({sigma }_{{{{{{rm{CM}}}}}}}), the likelihood function was penalised for deviations from unity proportionality, ({r}^{{{{{{rm{CM}}}}}}}), through:$$-{{log }}{{{{{rm{L}}}}}}left(0,{sigma }_{{{{{{rm{CM}}}}}}}^{2}{{{{{rm{|}}}}}}{r}^{{{{{{rm{CM}}}}}}}right)={sum }_{{{{{{rm{zone}}}}}}=1}^{{{{{{rm{#Zones}}}}}}}{r}_{{{{{{rm{zone}}}}}}}^{{{{{{rm{CM}}}}}}}/2{sigma }_{{{{{{rm{CM}}}}}}}^{2}$$ (41) Model simulations were conducted in ADMB with output analysis and visualisation conducted in MATLAB.Sensitivity to CoTS controlTo test whether our projected scenarios were consistent with the period over which data were collected, we conducted a model-based before and after comparison to the impact of control. Specifically, we used the fitted trajectory for sites, including both the coral data and CoTS control data (voyages and time spent), and compared this to the model-suggested coral trajectories if CoTS control had not taken place. These were modelled over the fitted period (2013–2018) and, unlike the projected scenarios (2019–2029), were variable in terms of the timing of control (amount of time between visits was variable), the amount of time spent at sites (not a consistent number of dive minutes per visit), CoTS dynamics (recruitment was fitted and hence different annually and between reefs), and in the level of thermal stress they experienced (different sites experienced different effective levels and some sites experience back-to-back events).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Scheffers, B. R. et al. The broad footprint of climate change from genes to biomes to people. Science 354, aaf7671 (2016).PubMed 
Article 
CAS 

Google Scholar 
Pecl, G. T. et al. Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being. Science 355, eaai214 (2017).Article 
CAS 

Google Scholar 
Nogués-Bravo, D. et al. Cracking the code of biodiversity responses to past climate change. Trends Ecol. Evol. 33, 765–776 (2018).PubMed 
Article 

Google Scholar 
Forsman, A., Betzholtz, P.-E. & Franzén, M. Faster poleward range shifts in moths with more variable colour patterns. Sci. Rep. 6, 36265 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Voelkl, B. et al. Reproducibility of animal research in light of biological variation. Nat. Rev. Neurosci. 21, 384–393 (2020).CAS 
PubMed 
Article 

Google Scholar 
Rödder, D., Schmitt, T., Gros, P., Ulrich, W. & Habel, J. C. Climate change drives mountain butterflies towards the summits. Sci. Rep. 11, 14382 (2021).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Habel, J. C., Teucher, M., Gros, P., Schmitt, T. & Ulrich, W. Land use and climate change affects butterfly diversity across northern Austria. Landscape Ecol. 36, 1741–1754 (2021).Article 

Google Scholar 
Hill, J. K. et al. Responses of butterflies to twentieth century climate warming: implications for future ranges. Proc. Biol. Sci. 269, 2163–2171 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Chen, I. C. et al. Elevation increases in moth assemblages over 42 years on a tropical mountain. Proc. Natl. Acad. Sci. 106, 1479–1483 (2009).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Cohen, J. M., Lajeunesse, M. J. & Rohr, J. R. A global synthesis of animal phenological responses to climate change. Nat. Clim. Change 8, 224–228 (2018).ADS 
Article 

Google Scholar 
Bell, J. R. et al. Spatial and habitat variation in aphid, butterfly, moth and bird phenologies over the last half century. Glob. Change Biol. 25, 1982–1994 (2019).ADS 
Article 

Google Scholar 
Hällfors, M. H. et al. Shifts in timing and duration of breeding for 73 boreal bird species over four decades. Proc. Natl. Acad. Sci. 117, 18557–18565 (2020).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Pruett, J. E. & Warner, D. A. Spatial and temporal variation in phenotypes and fitness in response to developmental thermal environments. Funct. Ecol. 35, 2635–2646 (2021).Article 

Google Scholar 
Hall, M., Nordahl, O., Larsson, P., Forsman, A. & Tibblin, P. Intra-population variation in reproductive timing covaries with thermal plasticity of offspring performance in perch Perca fluviatilis. J. Animal Ecol 90, 2236–2347 (2021).Article 

Google Scholar 
Ehrlich, P. R. & Hanski, I. On the Wings of Checkerspots: A Model System for Population Biology (Oxford University Press, 2004).
Google Scholar 
Warren, M. S. et al. The decline of butterflies in Europe: Problems, significance, and possible solutions. Proc. Natl. Acad. Sci. 118, e2002551117 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kristensen, N. P. Lepidoptera: Moths and Butterflies. 1. Evolution, Systematics, and Biogeography. Handbook of Zoology Vol. IV, Part 35 (De Gruyter, 1999).
Google Scholar 
Forsman, A. & Wennersten, L. Inter-individual variation promotes ecological success of populations and species: Evidence from experimental and comparative studies. Ecography 39, 630–648 (2016).Article 

Google Scholar 
Zografou, K. et al. Species traits affect phenological responses to climate change in a butterfly community. Sci. Rep. 11, 3283 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Stevens, C. J. et al. Nitrogen deposition threatens species richness of grasslands across Europe. Environ. Pollut. 158, 2940–2945 (2010).CAS 
PubMed 
Article 

Google Scholar 
Heinrich, B. The Thermal Warriors (Harvard University Press, 2013).
Google Scholar 
Bladon, A. J. et al. How butterflies keep their cool: Physical and ecological traits influence thermoregulatory ability and population trends. J. Anim. Ecol. 89, 2440–2450 (2020).PubMed 
Article 

Google Scholar 
Tsai, C.-C. et al. Physical and behavioral adaptations to prevent overheating of the living wings of butterflies. Nat. Commun. 11, 551 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ahnesjö, J. & Forsman, A. Differential habitat selection by pygmy grasshopper color morphs; interactive effects of temperature and predator avoidance. Evol. Ecol. 20, 235–257 (2006).Article 

Google Scholar 
Ma, C.-S., Ma, G. & Pincebourde, S. Survive a warming climate: Insect responses to extreme high temperatures. Annu. Rev. Entomol. 66, 163–184 (2021).CAS 
PubMed 
Article 

Google Scholar 
Forsman, A. Rethinking phenotypic plasticity and its consequences for individuals, populations and species. Heredity 115, 276–284 (2015).CAS 
PubMed 
Article 

Google Scholar 
Hill, G. M., Kawahara, A. Y., Daniels, J. C., Bateman, C. C. & Scheffers, B. R. Climate change effects on animal ecology: Butterflies and moths as a case study. Biol. Rev. Camb. Philos. Soc. 96, 2113–2126 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Gilbert, A. L. & Miles, D. B. Natural selection on thermal preference, critical thermal maxima and locomotor performance. Proc. R. Soc. B Biol. Sci. 284, 20170536 (2017).Article 

Google Scholar 
Eliasson, C. U., Ryrholm, N., Holmér, M., Gilg, K. & Gärdenfors, U. Nationalnyckeln till Sveriges flora och fauna. Fjärilar: Dagfjärilar. Hesperidae – Nymphalidae. (ArtDatabanken, SLU, 2005).Thomas, J. A. & Wardlaw, J. C. The capacity of a Myrmica ant nest to support a predacious species of Maculinea butterfly. Oecologia 91, 101–109 (1992).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Vilbas, M. et al. Habitat use of the endangered parasitic butterfly Phengaris arion close to its northern distribution limit. Insect Conserv. Divers. 8, 252–260 (2015).Article 

Google Scholar 
Johansson, V., Kindvall, O., Askling, J. & Franzén, M. Extreme weather affects colonization–extinction dynamics and the persistence of a threatened butterfly. J. Appl. Ecol. 57, 1068–1077 (2020).Article 

Google Scholar 
Johansson, V., Kindvall, O., Askling, J. & Franzén, M. Intense grazing of calcareous grasslands has negative consequences for the threatened marsh fritillary butterfly. Biol. Cons. 239, 108280 (2019).Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical. R version 4.1.1. (2021).Eubank, R. L. & Speckman, P. Curve fitting by polynomial-trigonometric regression. Biometrika 77, 1–9 (1990).MathSciNet 
MATH 
Article 

Google Scholar 
Allen, J. C. A modified sine wave method for calculating degree days. Environ. Entomol. 5, 388–396 (1976).Article 

Google Scholar 
Wickham, H. & Wickham, M. H. The ggplot package. Google Scholar. http://ftp.uni-bayreuth.de/math/statlib/R/CRAN/doc/packages/ggplot.pdf, (2007).Lüdecke, D. ggeffects: Tidy data frames of marginal effects from regression models. J. Open Source Softw. 3, 772 (2018).ADS 
Article 

Google Scholar 
Forsman, A. Some like it hot: Intra-population variation in behavioral thermoregulation in color-polymorphic pygmy grasshoppers. Evol. Ecol. 14, 25–38 (2000).Article 

Google Scholar 
Forsman, A., Ringblom, K., Civantos, E. & Ahnesjo, J. Coevolution of color pattern and thermoregulatory behavior in polymorphic pygmy grasshoppers Tetrix undulata. Evolution 56, 349–360 (2002).PubMed 
Article 

Google Scholar 
Ahnesjö, J. & Forsman, A. Correlated evolution of colour pattern and body size in polymorphic pygmy grasshoppers, Tetrix undulata. J. Evol. Biol. 16, 1308–1318 (2003).PubMed 
Article 

Google Scholar 
Zeuss, D., Brandl, R., Brändle, M., Rahbek, C. & Brunzel, S. Global warming favours light-coloured insects in Europe. Nat. Commun. 5, 3874 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Heidrich, L. et al. The dark side of Lepidoptera: colour lightness of geometrid moths decreases with increasing latitude. Glob. Ecol. Biogeogr. 27, 407–416 (2018).MathSciNet 
Article 

Google Scholar 
Porter, K. Basking behaviour in larvae of the butterfly Euphydryas aurinia. Oikos 38, 308–312 (1982).Article 

Google Scholar 
Rolff, J., Johnston, P. R. & Reynolds, S. Complete metamorphosis of insects. Philos. Trans. R. Soc. B 374, 20190063 (2019).Article 

Google Scholar 
Thomas, J. A., Simcox, D. J. & Clarke, R. T. Successful conservation of a threatened Maculinea butterfly. Science 325, 80–83 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Nilsson, M. & Forsman, A. Evolution of conspicuous colouration, body size and gregariousness: A comparative analysis of lepidopteran larvae. Evol. Ecol. 17, 51–66 (2003).Article 

Google Scholar 
Mappes, J., Kokko, H., Ojala, K. & Lindström, L. Seasonal changes in predator community switch the direction of selection for prey defences. Nat. Commun. 5, 5016 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Bale, J. S. et al. Herbivory in global climate change research: direct effects of rising temperature on insect herbivores. Glob. Change Biol. 8, 1–16 (2002).ADS 
Article 

Google Scholar 
Otaki, J. M., Hiyama, A., Iwata, M. & Kudo, T. Phenotypic plasticity in the range-margin population of the lycaenid butterfly Zizeeria maha. BMC Evol. Biol. 10, 1–13 (2010).Article 

Google Scholar 
Galarza, J. A. et al. Evaluating responses to temperature during pre-metamorphosis and carry-over effects at post-metamorphosis in the wood tiger moth (Arctia plantaginis). Philos. Trans. R. Soc. B 374, 20190295 (2019).CAS 
Article 

Google Scholar 
Kingsolver, J. G. The well-temperatured biologist: (American Society of Naturalists Presidential Address). Am. Nat. 174, 755–768 (2009).PubMed 
Article 

Google Scholar 
Lafuente, E. & Beldade, P. Genomics of developmental plasticity in animals. Front. Genet. 10, 720 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Angilletta, M. J. Jr., Niewiarowski, P. H. & Navas, C. A. The evolution of thermal physiology in ectotherms. J. Therm. Biol 27, 249–268 (2002).Article 

Google Scholar 
Posledovich, D., Toftegaard, T., Wiklund, C., Ehrlén, J. & Gotthard, K. Phenological synchrony between a butterfly and its host plants: Experimental test of effects of spring temperature. J. Anim. Ecol. 87, 150–161 (2018).PubMed 
Article 

Google Scholar 
Adams, A. Succisa pratensis Moench. J. Ecol. 43, 709–718 (1955).Article 

Google Scholar 
Lawton, J. H. & Strong, D. J. Community patterns and competition in folivorous insects. Am. Nat. 118, 317–338 (1981).Article 

Google Scholar 
Forsman, A. Effects of genotypic and phenotypic variation on establishment are important for conservation, invasion, and infection biology. Proc. Natl. Acad. Sci. 111, 302–307 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Forsman, A., Betzholtz, P.-E. & Franzén, M. Variable coloration is associated with dampened population fluctuations in noctuid moths. Proc. R. Soc. B 282, 1–9 (2015).Article 

Google Scholar 
Betzholtz, P. E., Franzén, M. & Forsman, A. Colour pattern variation can inform about extinction risk in moths. Anim. Conserv. 20, 72–79 (2017).Article 

Google Scholar 
Klemme, I. & Hanski, I. Heritability of and strong single gene (Pgi) effects on life-history traits in the Glanville fritillary butterfly. J. Evol. Biol. 22, 1944–1953 (2009).CAS 
PubMed 
Article 

Google Scholar 
Mattila, A. L. Thermal biology of flight in a butterfly: genotype, flight metabolism, and environmental conditions. Ecol. Evol. 5, 5539–5551 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Russell, B. D. et al. Predicting ecosystem shifts requires new approaches that integrate the effects of climate change across entire systems. Biol. Let. 8, 164–166 (2012).Article 

Google Scholar 
van Bergen, E. et al. The effect of summer drought on the predictability of local extinctions in a butterfly metapopulation. Conserv. Biol. 34, 1503–1511 (2020).PubMed 
Article 

Google Scholar 
Thomas, J. A., Clarke, R. T., Elmes, G. W. & Hochberg, M. E. in Insect Populations in theory and in practice: 19th Symposium of the Royal Entomological Society 10–11 September 1997 at the University of Newcastle (eds J. P. Dempster & I. F. G. McLean) 261–290 (Springer Netherlands, 1998).Nakonieczny, M., Kedziorski, A. & Michalczyk, K. Apollo butterfly (Parnassius apollo L.) in Europe—Its history, decline and perspectives of conservation. Funct. Ecosyst. Communities 1, 56–79 (2007).
Google Scholar 
Schweiger, O., Harpke, A., Wiemers, M. & Settele, J. CLIMBER: Climatic niche characteristics of the butterflies in Europe. ZooKeys 367, 65–84 (2014).Article 

Google Scholar 
Ashton, S., Gutierrez, D. & Wilson, R. J. Effects of temperature and elevation on habitat use by a rare mountain butterfly: Implications for species responses to climate change. Ecological Entomology 34, 437–446 (2009).Article 

Google Scholar 
Klockmann, M. & Fischer, K. Effects of temperature and drought on early life stages in three species of butterflies: Mortality of early life stages as a key determinant of vulnerability to climate change?. Ecol. Evol. 7, 10871–10879 (2017).PubMed 
PubMed Central 
Article 

Google Scholar  More

Biological control in IPM systems in Africa. (CABI, 2002). https://doi.org/10.1079/9780851996394.0000Kvakkestad, V., Sundbye, A., Gwynn, R. & Klingen, I. Authorization of microbial plant protection products in the Scandinavian countries: A comparative analysis. Environ. Sci. Policy 106, 115–124 (2020).Article 

Google Scholar 
Barzman, M. et al. Eight principles of integrated pest management. Agron. Sustain. Dev. 35, 1199–1215 (2015).Article 

Google Scholar 
Popp, J., Pető, K. & Nagy, J. Pesticide productivity and food security. A review. Agron. Sustain. Dev. 33, 243–255 (2013).Article 

Google Scholar 
Bale, J., van Lenteren, J. & Bigler, F. Biological control and sustainable food production. Philos. Trans. R. Soc. B Biol. Sci. 363, 761–776 (2008).CAS 
Article 

Google Scholar 
Vacante, V. & Bonsignore, C. P. Natural enemies and pest control. In Handbook of Pest Management in Organic Farming 60–77 (CABI, 2018). https://doi.org/10.1079/9781780644998.0060Eilenberg, J., Hajek, A. & Lomer, C. Suggestions for unifying the terminology in biological control. Biocontrol 46, 387–400 (2001).Article 

Google Scholar 
Lacey, L. A. et al. Insect pathogens as biological control agents: Back to the future. J. Invertebr. Pathol. 132, 1–41 (2015).CAS 
PubMed 
Article 

Google Scholar 
Hatting, J. L., Moore, S. D. & Malan, A. P. Microbial control of phytophagous invertebrate pests in South Africa: Current status and future prospects. J. Invertebr. Pathol. 165, 54–66 (2019).PubMed 
Article 

Google Scholar 
Karimi, S., Askari Seyahooei, M., Izadi, H., Bagheri, A. & Khodaygan, P. Effect of arsenophonus endosymbiont elimination on fitness of the date palm hopper, ommatissus lybicus (Hemiptera: Tropiduchidae). Environ. Entomol. 48, 614–622 (2019).CAS 
PubMed 
Article 

Google Scholar 
Kumar, K. K. et al. Microbial biopesticides for insect pest management in India: Current status and future prospects. J. Invertebr. Pathol. 165, 74–81 (2019).CAS 
PubMed 
Article 

Google Scholar 
Mascarin, G. M. et al. Current status and perspectives of fungal entomopathogens used for microbial control of arthropod pests in Brazil. J. Invertebr. Pathol. 165, 46–53 (2019).PubMed 
Article 

Google Scholar 
Shapiro-Ilan, D. I., Bruck, D. J. & Lacey, L. A. Principles of epizootiology and microbial control. Insect Pathol. https://doi.org/10.1016/B978-0-12-384984-7.00003-8 (2012).Article 

Google Scholar 
Hawkins, B. A. & Cornell, H. V. Theoretical Approaches to Biological Control. https://doi.org/10.1017/CBO9780511542077 (Cambridge University Press, 2009).Tonnang, H. E. Z., Nedorezov, L. V., Ochanda, H., Owino, J. & Löhr, B. Assessing the impact of biological control of Plutella xylostella through the application of Lotka—Volterra model. Ecol. Model. 220, 60–70 (2009).Article 

Google Scholar 
Hesketh, H., Roy, H. E., Eilenberg, J., Pell, J. K. & Hails, R. S. Challenges in modelling complexity of fungal entomopathogens in semi-natural populations of insects. Biocontrol 55, 55–73 (2010).Article 

Google Scholar 
Fuxa, J. R. & Tanada, Y. Epizootiology of Insect Diseases (Wiley, 1987).
Google Scholar 
Lacey, L. A. Manual of Techniques in Insect Pathology. Manual of Techniques in Insect Pathology (Academic, 1997). https://doi.org/10.1016/b978-0-12-432555-5.x5000-3.Book 

Google Scholar 
Lomer, C. J., Bateman, R. P., Johnson, D. L., Langewald, J. & Thomas, M. Biological control of locusts and grasshoppers. Annu. Rev. Entomol. 46, 667–702 (2001).CAS 
PubMed 
Article 

Google Scholar 
Arthurs, S. & Thomas, M. B. Effects of a mycoinsecticide on feeding and fecundity of the brown locust Locustana pardalina. Biocontrol Sci. Technol. 10, 321–329 (2000).Article 

Google Scholar 
Jiang, W. et al. Effects of the entomopathogenic fungus Metarhizium anisopliae on the mortality and immune response of Locusta migratoria. Insects 11, 36 (2020).Article 

Google Scholar 
Thomas, M. B. & Blanford, S. Thermal biology in insect-parasite interactions. Trends Ecol. Evol. 18, 344–350 (2003).Article 

Google Scholar 
Douthwaite, M. B. Development and Commercialization of the Green Muscle Biopesticide 21 (2001).Douthwaite, B., Langewald, J., & Harris, J. Development and commercialization of the Green Muscle biopesticide. (International Institute of Tropical Agriculture, 2002).CABI. Green Muscle providing strength against devastating locusts in the horn of Africa—CABI.org. CABI.org https://www.cabi.org/news-article/green-muscle-providing-strength-against-devastating-locusts-in-the-horn-of-africa/ (2020).Geoff, G. & Steve, W. Biological Control (Springer, 1996). https://doi.org/10.1007/978-1-4613-1157-7.Book 

Google Scholar 
Fargues, J., Ouedraogo, A., Goettel, M. S. & Lomer, C. J. Effects of temperature, humidity and inoculation method on susceptibility of Schistocerca gregaria to Metarhizium flavoviride. Biocontrol Sci. Technol. 7, 345–356 (1997).Article 

Google Scholar 
Aragón, P., Coca-Abia, M. M., Llorente, V. & Lobo, J. M. Estimation of climatic favourable areas for locust outbreaks in Spain: Integrating species’ presence records and spatial information on outbreaks. J. Appl. Entomol. 137, 610–623 (2013).Article 

Google Scholar 
Arthurs, S. & Thomas, M. B. Effect of dose, pre-mortem host incubation temperature and thermal behaviour on host mortality, mycosis and sporulation of Metarhizium anisopliae var. acridum in Schistocerca gregaria. Biocontrol Sci. Technol. 11, 411–420 (2001).Article 

Google Scholar 
van der Valk, H. Review of the efficacy of Metarhizium anisopliae var. acridum. FAO—U.N. Publ. (2007).Klass, J. I., Blanford, S. & Thomas, M. B. Development of a model for evaluating the effects of environmental temperature and thermal behaviour on biological control of locusts and grasshoppers using pathogens. Agric. For. Entomol. 9, 189–199 (2007).Article 

Google Scholar 
Devi, K. U., Sridevi, V., Mohan, C. M. & Padmavathi, J. Effect of high temperature and water stress on in vitro germination and growth in isolates of the entomopathogenic fungus Beauveria bassiana (Bals.) Vuillemin. J. Invertebr. Pathol. 88, 181–189 (2005).PubMed 
Article 

Google Scholar 
Dimbi, S., Maniania, N. K., Lux, S. A. & Mueke, J. M. Effect of constant temperatures on germination, radial growth and virulence of Metarhizium anisopliae to three species of African tephritid fruit flies. Biocontrol 49, 83–94 (2004).Article 

Google Scholar 
Ekesi, S., Maniania, N. K. & Ampong-Nyarko, K. Effect of temperature on germination, radial growth and virulence of Metarhizium anisopliae and Beauveria bassiana on Megalurothrips sjostedti. Biocontrol Sci. Technol. 9, 177–185 (1999).Article 

Google Scholar 
Thomas, M. B. & Jenkins, N. E. Effects of temperature on growth of Metarhizium flavoviride and virulence to the variegated grasshopper Zonocerus variegatus. Mycol. Res. 101, 1469–1474 (1997).Article 

Google Scholar 
Klass, J. I., Blanford, S. & Thomas, M. B. Use of a geographic information system to explore spatial variation in pathogen virulence and the implications for biological control of locusts and grasshoppers. Agric. For. Entomol. 9, 201–208 (2007).Article 

Google Scholar 
Castro, T., Moral, R., Demétrio, C., Delalibera, I. & Klingen, I. Prediction of sporulation and germination by the spider mite pathogenic fungus Neozygites floridana (Neozygitomycetes: Neozygitales: Neozygitaceae) based on temperature, humidity and time. Insects 9, 69 (2018).PubMed Central 
Article 

Google Scholar 
Hajek, A. E., Larkin, T. S., Carruthers, R. I. & Soper, R. S. Modelling the dynamics of Entomophaga maimaga (Zygomycetes: Entomophtorales) epizootics in gypsy moth (Lepidoptera: Lymantridae) populations. Environ. Entomol. 22, 1172–1187 (1993).Article 

Google Scholar 
Gul, H. T., Saeed, S. & Khan, F. A. Z. Entomopathogenic fungi as effective insect pest management tactic: A review. Appl. Sci. Bus. Econ. 1, 10–18 (2014).
Google Scholar 
Davidson, G. et al. Study of temperature—Growth interactions of entomopathogenic fungi with potential for control of Varroa destructor (Acari: Mesostigmata) using a nonlinear model of poikilotherm development. J. Appl. Microbiol. 94, 816–825 (2003).CAS 
PubMed 
Article 

Google Scholar 
Hallsworth, J. E. & Magan, N. Water and temperature relations of growth of the entomogenous fungi Beauveria bassiana, Metarhizium anisopliae, and Paecilomyces farinosus. J. Invertebr. Pathol. 74, 261–266 (1999).CAS 
PubMed 
Article 

Google Scholar 
Fargues, J. et al. Climatic factors on entomopathogenic hyphomycetes infection of Trialeurodes vaporariorum (Homoptera: Aleyrodidae) in Mediterranean glasshouse tomato. Biol. Control 28, 320–331 (2003).Article 

Google Scholar 
Boulard, T. et al. Effect of greenhouse ventilation on humidity of inside air and in leaf boundary-layer. Agric. For. Meteorol. 125, 225–239 (2004).ADS 
Article 

Google Scholar 
Mishra, S., Kumar, P. & Malik, A. Effect of temperature and humidity on pathogenicity of native Beauveria bassiana isolate against Musca domestica L. J. Parasit. Dis. 39, 697–704 (2015).PubMed 
Article 

Google Scholar 
Klingen, I., Westrum, K. & Meyling, N. V. Effect of Norwegian entomopathogenic fungal isolates against Otiorhynchus sulcatus larvae at low temperatures and persistence in strawberry rhizospheres. Biol. Control 81, 1–7 (2015).Article 

Google Scholar 
Thaochan, N., Benarlee, R., Shekhar Prabhakar, C. & Hu, Q. Impact of temperature and relative humidity on effectiveness of Metarhizium guizhouense PSUM02 against longkong bark eating caterpillar Cossus chloratus Swinhoe under laboratory and field conditions. J. Asia. Pac. Entomol. 23, 285–290 (2020).Article 

Google Scholar 
Kryukov, V. et al. Ecological preferences of Metarhizium spp. from Russia and neighboring territories and their activity against Colorado potato beetle larvae. J. Invertebr. Pathol. 149, 1–7 (2017).PubMed 
Article 

Google Scholar 
Saldarriaga Ausique, J. J., D’Alessandro, C. P., Conceschi, M. R., Mascarin, G. M. & Delalibera Júnior, I. Efficacy of entomopathogenic fungi against adult Diaphorina citri from laboratory to field applications. J. Pest Sci. 2017 903 90, 947–960 (2017).
Google Scholar 
Dwyer, G. Density dependence and spatial structure in the dynamics of insect pathogens. Am. Nat. 143, 533–562 (1994).ADS 
Article 

Google Scholar 
Dwyer, G., Elkinton, J. & Hajek, A. Spatial scale and the spread of a fungal pathogen of gypsy moth. Am. Nat. 152, 485–494 (1998).CAS 
PubMed 
Article 

Google Scholar 
Knudsen, G. R. & Schotzko, D. J. Spatial simulation of epizootics caused by Beauveria bassiana in Russian wheat aphid populations. Biol. Control 16, 318–326 (1999).Article 

Google Scholar 
Weseloh, R. M. Effect of conidial dispersal of the fungal pathogen Entomophaga maimaiga (Zygomycetes: Entomophthorales) on survival of its gypsy moth (Lepidoptera: Lymantriidae) host. Biol. Control 29, 138–144 (2004).Article 

Google Scholar 
Meynard, C. N. et al. Climate-driven geographic distribution of the desert locust during recession periods: Subspecies’ niche differentiation and relative risks under scenarios of climate change. Glob. Chang. Biol. 23, 4739–4749 (2017).ADS 
PubMed 
Article 

Google Scholar 
Anderson, R. M. & May, R. M. Infectious diseases of humans: Dynamics and control. Aust. J. Public Health 16, 208–212 (1991).
Google Scholar 
Cáceres, C. E. et al. Complex Daphnia interactions with parasites and competitors. Math. Biosci. 258, 148–161 (2014).MathSciNet 
PubMed 
MATH 
Article 

Google Scholar 
Briggs, C. J. & Godfray, H. C. J. The dynamics of insect-pathogen interactions stage-structured populations c. J. Am. Nat. 145, 855–887 (1995).Article 

Google Scholar 
Rapti, Z. & Cáceres, C. E. Effects of intrinsic and extrinsic host mortality on disease spread. Bull. Math. Biol. 78, 235–253 (2016).MathSciNet 
CAS 
PubMed 
MATH 
Article 

Google Scholar 
Hartemink, N. A., Randolph, S. E., Davis, S. A. & Heesterbeek, J. A. P. The basic reproduction number for complex disease systems: Defining R0 for tick-borne infections. Am. Nat. 171, 743–754 (2014).Article 

Google Scholar 
Arthur, F. H. Toxicity of diatomaceous earth to red flour beetles and confused flour beetles (Coleoptera: Tenebrionidae): Effects of temperature and relative humidity. J. Econ. Entomol. 93, 526–532 (2000).CAS 
PubMed 
Article 

Google Scholar 
Arthurs, S. & Thomas, M. B. Effects of temperature and relative humidity on sporulation of Metarhizium anisopliae var. acridum in mycosed cadavers of Schistocerca gregaria. J. Invertebr. Pathol. 78, 59–65 (2001).CAS 
PubMed 
Article 

Google Scholar 
Whipps, J. M. & Davies, K. G. Success in Biological Control of Plant Pathogens and Nematodes by Microorganisms. In Biological Control: Measures of Success 1st edn, (eds Gurr, G. & Wratten, S.) 429. https://doi.org/10.1007/978-94-011-4014-0_8 (Springer, Dordrecht, 2000).Gilchrist, M. A., Sulsky, D. L. & Pringle, A. Identifying fitness and optimal life-history strategies for an asexual filamentous fungus. Evolution 60, 970–979 (2006).PubMed 
Article 

Google Scholar 
Frank, S. A. Spatial processes in host-parasite genetics. In Metapopulation Biology, 1st edn, (eds Hanski, I. A. & Gilpin, M. E.) 325–352. https://doi.org/10.1016/B978-012323445-2/50018-3 (Elsevier, 1997).Yan, Y., Wang, Y.-C., Feng, C.-C., Wan, P.-H.M. & Chang, K.T.-T. Potential distributional changes of invasive crop pest species associated with global climate change. Appl. Geogr. 82, 83–92 (2017).Article 

Google Scholar 
Inglis, G. D., Johnson, D. L. & Goettel, M. S. Effects of temperature and thermoregulation on mycosis by Beauveria bassianain grasshoppers. Biol. Control 7, 131–139 (1996).Article 

Google Scholar 
Lactin, D. J. & Johnson, D. L. Temperature-dependent feeding rates of Melanoplus sanguinipes nymphs (Orthoptera: Acrididae) laboratory trials. Environ. Entomol. 24, 1291–1296 (1995).Article 

Google Scholar 
FAO. Biopesticides for locust control | FAO Stories | Food and Agriculture Organization of the United Nations. Food and Agriculture Organisation of the UN http://www.fao.org/fao-stories/article/en/c/1267098/ (2021).Kimathi, E. et al. Prediction of breeding regions for the desert locust Schistocerca gregaria in East Africa. Sci. Rep. 10, 11937 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Cordovez, J. M., Rendon, L. M., Gonzalez, C. & Guhl, F. Using the basic reproduction number to assess the effects of climate change in the risk of Chagas disease transmission in Colombia. Acta Trop. 129, 74–82 (2014).PubMed 
Article 

Google Scholar 
Hartemink, N. A. et al. Mapping the basic reproduction number ( R 0) for vector-borne diseases: A case study on bluetongue virus. EPIDEM 1, 153–161 (2009).CAS 
Article 

Google Scholar 
Jamison, A., Tuttle, E., Jensen, R., Bierly, G. & Gonser, R. Spatial ecology, landscapes, and the geography of vector-borne disease: A multi-disciplinary review. Appl. Geogr. 63, 418–426 (2015).Article 

Google Scholar 
Moukam Kakmeni, F. M. et al. Spatial panorama of malaria prevalence in Africa under climate change and interventions scenarios. Int. J. Health Geogr. 17, 2 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Ngarakana-Gwasira, E. T., Bhunu, C. P., Masocha, M. & Mashonjowa, E. Transmission dynamics of schistosomiasis in Zimbabwe: A mathematical and GIS approach. Commun. Nonlinear Sci. Numer. Simul. 35, 137–147 (2016).ADS 
MathSciNet 
MATH 
Article 

Google Scholar 
Ogden, N. H. & Radojevic, M. Estimated effects of projected climate change on the basic reproductive number of the Lyme disease vector ixodes scapularis. Environ. Health Perspect. 122, 631–639 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Parham, P. E. & Michael, E. Modeling the effects of weather and climate change on malaria transmission. Environ. Health Perspect. 118, 620–626 (2010).PubMed 
Article 

Google Scholar 
Phillips, J. Climate change and surface mining: A review of environment-human interactions & their spatial dynamics. Appl. Geogr. 74, 95–108 (2016).Article 

Google Scholar 
Rogers, D. J. & Randolphz, S. E. The global spread of malaria in a future. Warmer World Sci. 2, 1763–1766 (2000).
Google Scholar 
Wu, X. et al. Developing a temperature-driven map of the basic reproductive number of the emerging tick vector of Lyme disease Ixodes scapularis in Canada. J. Theor. Biol. 319, 50–61 (2013).ADS 
MathSciNet 
PubMed 
MATH 
Article 

Google Scholar 
CABI. Green Muscle providing strength against devastating locusts in the horn of Africa. https://www.cabi.org/news-article/green-muscle-providing-strength-against-devastating-locusts-in-the-horn-of-africa/ (2020).Piou, C. et al. Mapping the spatiotemporal distributions of the Desert Locust in Mauritania and Morocco to improve preventive management. Basic Appl. Ecol. 25, 37–47 (2017).Article 

Google Scholar 
FAO. FAO Locust Hub. https://locust-hub-hqfao.hub.arcgis.com/ (2021).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 
DeJesus, E. X. & Kaufman, C. Routh-Hurwitz criterion in the examination of eigenvalues of a system of nonlinear ordinary differential equations. Phys. Rev. A 35, 5288–5290 (1987).ADS 
MathSciNet 
CAS 
Article 

Google Scholar 
QGIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.osgeo.org. Qgisorg (2014).RCoreTeam. R: A language and environment for statistical computing. The R Foundation for Statistical Computing. (2020).Marino, S., Hogue, I. B., Ray, C. J. & Kirschner, D. E. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J. Theor. Biol. 254, 178–196 (2008).ADS 
MathSciNet 
PubMed 
PubMed Central 
MATH 
Article 

Google Scholar  More

Soulsbury, C. D. & White, P. C. L. Human–wildlife interactions in urban areas: A review of conflicts, benefits and opportunities. Wildl. Res. 42, 541 (2015).Article 

Google Scholar 
Tucker, M. A. et al. Moving in the Anthropocene: Global reductions in terrestrial mammalian movements. Science 359, 466–469 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Wilson, M. W. et al. Ecological impacts of human-induced animal behaviour change. Ecol. Lett. 23, 1522–1536 (2020).PubMed 
Article 

Google Scholar 
Silva-Rodríguez, E. A., Gálvez, N., Swan, G. J. F., Cusack, J. J. & Moreira-Arce, D. Urban wildlife in times of COVID-19: What can we infer from novel carnivore records in urban areas?. Sci. Total Environ. https://doi.org/10.1016/j.scitotenv.2020.142713 (2020).Article 
PubMed 

Google Scholar 
Joshi, Y. V. & Musalem, A. Lockdowns lose one third of their impact on mobility in a month. Sci Rep 11, 22658 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Chung, P.-C. & Chan, T.-C. Impact of physical distancing policy on reducing transmission of SARS-CoV-2 globally: Perspective from government’s response and residents’ compliance. PLoS ONE 16, e0255873 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Corlett, R. T. et al. Impacts of the coronavirus pandemic on biodiversity conservation. Biol. Conserv. 246, 108571 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Connellan, I. The ‘anthropause’ during COVID-19. Cosmos Magazine https://cosmosmagazine.com/nature/animals/the-anthropause-during-covid-19/ (2020).Rutz, C. et al. COVID-19 lockdown allows researchers to quantify the effects of human activity on wildlife. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-020-1237-z (2020).Article 
PubMed 

Google Scholar 
Derryberry, E. P., Phillips, J. N., Derryberry, G. E., Blum, M. J. & Luther, D. Singing in a silent spring: Birds respond to a half-century soundscape reversion during the COVID-19 shutdown. Science 370, 575–579 (2020).CAS 
PubMed 
Article 

Google Scholar 
Gordo, O., Brotons, L., Herrando, S. & Gargallo, G. Rapid behavioural response of urban birds to COVID-19 lockdown. Proc. R. Soc. B Biol. Sci. 288, 20202513 (2021).CAS 
Article 

Google Scholar 
Gaynor, K. M. et al. Anticipating the impacts of the COVID-19 pandemic on wildlife. Front. Ecol. Environ. 18, 542–543 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Humphrey, C. Under cover of COVID-19, loggers plunder Cambodian wildlife sanctuary. Mongabay Environmental News https://news.mongabay.com/2020/08/under-cover-of-covid-19-loggers-plunder-cambodian-wildlife-sanctuary/ (2020).Bates, A. E., Primack, R. B., Moraga, P. & Duarte, C. M. COVID-19 pandemic and associated lockdown as a “Global Human Confinement Experiment” to investigate biodiversity conservation. Biol. Cons. 248, 108665 (2020).Article 

Google Scholar 
Nickel, B. A., Suraci, J. P., Allen, M. L. & Wilmers, C. C. Human presence and human footprint have non-equivalent effects on wildlife spatiotemporal habitat use. Biol. Cons. 241, 108383 (2020).Article 

Google Scholar 
Zellmer, A. J. et al. What can we learn from wildlife sightings during the COVID-19 global shutdown?. Ecosphere 11, e03215 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Jägerbrand, A. K., Antonson, H. & Ahlström, C. Speed reduction effects over distance of animal-vehicle collision countermeasures – a driving simulator study. Eur. Transp. Res. Rev. 10, 40 (2018).Article 

Google Scholar 
Abra, F. D. et al. Pay or prevent? Human safety, costs to society and legal perspectives on animal-vehicle collisions in São Paulo state. Brazil. PLoS One 14, e0215152 (2019).CAS 
PubMed 
Article 

Google Scholar 
Canal, D., Martín, B., de Lucas, M. & Ferrer, M. Dogs are the main species involved in animal-vehicle collisions in southern Spain: Daily, seasonal and spatial analyses of collisions. PLoS One 13, e0203693 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Visintin, C., van der Ree, R. & McCarthy, M. A. Consistent patterns of vehicle collision risk for six mammal species. J. Environ. Manage. 201, 397–406 (2017).PubMed 
Article 

Google Scholar 
Kreling, S. E. S., Gaynor, K. M. & Coon, C. A. C. Roadkill distribution at the wildland-urban interface. J. Wildl. Manag. 83, 1427–1436 (2019).Article 

Google Scholar 
Bíl, M. et al. COVID-19 related travel restrictions prevented numerous wildlife deaths on roads: A comparative analysis of results from 11 countries. Biol. Cons. 256, 109076 (2021).Article 

Google Scholar 
Langbein, J., Putman, R. & Pokorny, B. Traffic collisions involving deer and other ungulates in Europe and available measures for mitigation. Ungulate management in Europe: problems and practices 215–259 (2010).Filonchyk, M., Hurynovich, V. & Yan, H. Impact of Covid-19 lockdown on air quality in the Poland, Eastern Europe. Environ. Res. 198, 110454 (2021).CAS 
PubMed 
Article 

Google Scholar 
Porębska, A. et al. Lockdown in a disneyfied city: Kraków Old Town and the first wave of the Covid-19 pandemic. Urban Des Int 26, 315–331 (2021).Article 

Google Scholar 
Tarkowski, M., Puzdrakiewicz, K., Jaczewska, J. & Połom, M. COVID-19 lockdown in Poland – changes in regional and local mobility patterns based on Google Maps data. Prace Komisji Geografii Komunikacji PTG 2020, 46–55 (2020).Article 

Google Scholar 
Dean, W. R. J., Seymour, C. L., Joseph, G. S. & Foord, S. H. A review of the impacts of roads on wildlife in semi-arid regions. Diversity 11, 81 (2019).Article 

Google Scholar 
Saint-Andrieux, C., Calenge, C. & Bonenfant, C. Comparison of environmental, biological and anthropogenic causes of wildlife–vehicle collisions among three large herbivore species. Popul. Ecol. 62, 64–79 (2020).Article 

Google Scholar 
Grosman, P. D., Jaeger, J. A. G., Biron, P. M., Dussault, C. & Ouellet, J.-P. Trade-off between road avoidance and attraction by roadside salt pools in moose: An agent-based model to assess measures for reducing moose-vehicle collisions. Ecol. Model. 222, 1423–1435 (2011).Article 

Google Scholar 
Barbosa, P., Schumaker, N. H., Brandon, K. R., Bager, A. & Grilo, C. Simulating the consequences of roads for wildlife population dynamics. Landsc. Urban Plan. 193, 103672 (2020).PubMed 
Article 

Google Scholar 
Silva, C., Simões, M. P., Mira, A. & Santos, S. M. Factors influencing predator roadkills: The availability of prey in road verges. J Environ Manage 247, 644–650 (2019).PubMed 
Article 

Google Scholar 
Sullivan, J. M. Trends and characteristics of animal-vehicle collisions in the United States. J. Safety Res. 42, 9–16 (2011).PubMed 
Article 

Google Scholar 
Morelle, К, Lehaire, F. & Lejeune, P. Spatio-temporal patterns of wildlife-vehicle collisions in a region with a high-density road network. Nature Conservation 5, 53–73 (2013).Article 

Google Scholar 
Bartonička, T., Andrášik, R., Duľa, M., Sedoník, J. & Bíl, M. Identification of local factors causing clustering of animal-vehicle collisions. J. Wildl. Manag. 82, 940–947 (2018).Article 

Google Scholar 
Saxena, A., Chatterjee, N., Rajvanshi, A. & Habib, B. Integrating large mammal behaviour and traffic flow to determine traversability of roads with heterogeneous traffic on a Central Indian Highway. Sci Rep 10, 18888 (2020).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Basak, S. M. et al. Human-wildlife conflicts in Krakow City, Southern Poland. Animals 10, 1014 (2020).PubMed Central 
Article 

Google Scholar 
Gil-Fernández, M., Harcourt, R., Newsome, T., Towerton, A. & Carthey, A. Adaptations of the red fox (Vulpes vulpes) to urban environments in Sydney, Australia. J. Urban Ecol. https://doi.org/10.1093/jue/juaa009 (2020).Article 

Google Scholar 
Podgórski, T. et al. Spatiotemporal behavioral plasticity of wild boar (Sus scrofa) under contrasting conditions of human pressure: primeval forest and metropolitan area. J Mammal 94, 109–119 (2013).Article 

Google Scholar 
Steiner, W., Schöll, E. M., Leisch, F. & Hackländer, K. Temporal patterns of roe deer traffic accidents: Effects of season, daytime and lunar phase. PLoS ONE 16, e0249082 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Cagnacci, F. et al. Partial migration in roe deer: migratory and resident tactics are end points of a behavioural gradient determined by ecological factors. Oikos 120, 1790–1802 (2011).Article 

Google Scholar 
Kämmerle, J.-L. et al. Temporal patterns in road crossing behaviour in roe deer (Capreolus capreolus) at sites with wildlife warning reflectors. PLoS One 12, e0184761 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Romanowski, J. Vistula river valley as the ecological corridor for mammals. Pol. J. Ecol. 55, 805–819 (2007).
Google Scholar 
Abraham, J. O. & Mumma, M. A. Elevated wildlife-vehicle collision rates during the COVID-19 pandemic. Sci Rep 11, 20391 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Gunson, K. E., Mountrakis, G. & Quackenbush, L. J. Spatial wildlife-vehicle collision models: A review of current work and its application to transportation mitigation projects. J. Environ. Manage. 92, 1074–1082 (2011).PubMed 
Article 

Google Scholar 
Leblond, M., Dussault, C. & Ouellet, J.-P. Avoidance of roads by large herbivores and its relation to disturbance intensity. J. Zool. 289, 32–40 (2013).Article 

Google Scholar 
Bissonette, J. A. & Kassar, C. A. Locations of deer–vehicle collisions are unrelated to traffic volume or posted speed limit. Human-Wildlife Conflicts 2, 122–130 (2008).
Google Scholar 
Steiner, W., Leisch, F. & Hackländer, K. A review on the temporal pattern of deer–vehicle accidents: Impact of seasonal, diurnal and lunar effects in cervids. Accid. Anal. Prev. 66, 168–181 (2014).PubMed 
Article 

Google Scholar 
Kušta, T., Keken, Z., Ježek, M., Holá, M. & Šmíd, P. The effect of traffic intensity and animal activity on probability of ungulate-vehicle collisions in the Czech Republic. Saf. Sci. 91, 105–113 (2017).Article 

Google Scholar 
Shilling, F. et al. A Reprieve from US wildlife mortality on roads during the COVID-19 pandemic. Biol. Cons. 256, 109013 (2021).Article 

Google Scholar 
Yasin, Y. J., Grivna, M. & Abu-Zidan, F. M. Global impact of COVID-19 pandemic on road traffic collisions. World J Emerg Surg 16, 51 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Seiler, A. & Helldin, J. O. Mortality in wildlife due to transportation. In The Ecology of Transportation: Managing Mobility for the Environment (eds Davenport, J. & Davenport, J. L.) (Springer, 2006).
Google Scholar 
Smits, R., Bohatkiewicz, J., Bohatkiewicz, J. & Hałucha, M. A Geospatial Multi-scale Level Analysis of the Distribution of Animal-Vehicle Collisions on Polish Highways and National Roads. In Vision Zero for Sustainable Road Safety in Baltic Sea Region (eds Varhelyi, A. et al.) (Springer International Publishing, 2020).
Google Scholar 
Sozański, B. et al. Psychological responses and associated factors during the initial stage of the coronavirus disease (COVID-19) epidemic among the adult population in Poland – a cross-sectional study. BMC Public Health 21, 1929 (2021).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Sidor, A. & Rzymski, P. Dietary choices and habits during COVID-19 lockdown: Experience from Poland. Nutrients 12, E1657 (2020).PubMed 
Article 
CAS 

Google Scholar 
Vingilis, E. et al. Coronavirus disease 2019: What could be the effects on Road safety?. Accid. Anal. Prev. 144, 105687 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Kioko, J. et al. Driver knowledge and attitudes on animal vehicle collisions in Northern Tanzania. Trop. Conserv. Sci. 8, 352–366 (2015).Article 

Google Scholar 
Stokstad, E. Pandemic lockdown stirs up ecological research. Science 369, 893–893 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Dandy, N. Behaviour, lockdown and the natural world. Environ. Values 29, 253–259 (2020).Article 

Google Scholar 
Baścik, M. & Degórska, B. Środowisko przyrodnicze Krakowa. Zasoby – Ochrona – Kształtowanie. vol. 2 (2015).Borcard, D., Gillet, F. & Legendre, P. Numerical Ecology with R (Springer, 2011).MATH 
Book 

Google Scholar 
R Core Team. R: a language and environment for statistical computing. Vienna (Austria): R Foundation for Statistical Computing. https://www.r-project.org/ (2020).Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).MATH 
Book 

Google Scholar 
Oksanen, J. et al. vegan: Community Ecology Package (2019).Hervé, M. RVAideMemoire: Testing and Plotting Procedures for Biostatistics (2020).Hancock, J. M. Jaccard Distance (Jaccard Index, Jaccard Similarity Coefficient). in Dictionary of Bioinformatics and Computational Biology (American Cancer Society, 2014). https://doi.org/10.1002/9780471650126.dob0956 More

Calow, P. Towards a definition of functional ecology. Funct. Ecol. 1, 57–61 (1987).Article 

Google Scholar 
Violle, C. et al. Let the concept of trait be functional! Oikos 116, 882–892 (2007).Article 

Google Scholar 
McGill, B. J., Enquist, B. J., Weiher, E. & Westoby, M. Rebuilding community ecology from functional traits. Trends Ecol. Evol. 21, 178–185 (2006).Article 
PubMed 

Google Scholar 
Cadotte, M. W. Functional traits explain ecosystem function through opposing mechanisms. Ecol. Lett. 20, 989–996 (2017).Article 
PubMed 

Google Scholar 
Dehling, D. M. & Stouffer, D. B. Bringing the Eltonian niche into functional diversity. Oikos 127, 1711–1723 (2018).Article 

Google Scholar 
Schleuter, D., Daufresne, M., Massol, F. & Argillier, C. A user’s guide to functional diversity indices. Ecol. Monogr. 80, 469–484 (2010).Article 

Google Scholar 
Leinster, T. & Cobbold, C. A. Measuring diversity: the importance of species similarity. Ecology 93, 477–489 (2012).Article 
PubMed 

Google Scholar 
Carmona, C. P., de Bello, F., Mason, N. W. H. & Lepš, J. Traits without borders: integrating functional diversity across scales. Trends Ecol. Evol. 31, 382–394 (2016).Article 
PubMed 

Google Scholar 
Chao, A. et al. An attribute-diversity approach to functional diversity, functional beta diversity, and related (dis)similarity measures. Ecol. Monogr. 89, e01343 (2019).ADS 
Article 

Google Scholar 
Morris, E. K. et al. Choosing and using diversity indices: insights for ecological applications from the German biodiversity exploratories. Ecol. Evol. 4, 3514–3524 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Kattge, J., Bönisch, G. & D’iaz, S. et al. TRY plant trait database – enhanced coverage and open access. Glob. Change Biol. 26, 119–188 (2020).ADS 
Article 

Google Scholar 
Fajardo, A. & Siefert, A. Intraspecific trait variation and the leaf economics spectrum across resource gradients and levels of organization. Ecology 99, 1024–1030 (2018).Article 
PubMed 

Google Scholar 
Mayfield, M. M. et al. What does species richness tell us about functional trait diversity? Predictions and evidence for responses of species and functional trait diversity to land-use change. Glob. Ecol. Biogeogr. 19, 423–431 (2010).
Google Scholar 
Wieczynski, D. J. et al. Climate shapes and shifts functional biodiversity in forests worldwide. Proc. Natl. Acad. Sci. 116, 587–592 (2019).CAS 
Article 
PubMed 

Google Scholar 
Tilman, D., Isbell, F. & Cowles, J. M. Biodiversity and ecosystem functioning. Annu. Rev. Ecol. Evol. Syst. 45, 471–493 (2014).Article 

Google Scholar 
Petchey, O. L. & Gaston, K. J. Functional diversity: back to basics and looking forward. Ecol. Lett. 9, 741–758 (2006).Article 
PubMed 

Google Scholar 
Bolnick, D. I. et al. Why intraspecific trait variation matters in community ecology. Trends Ecol. Evol. 26, 183–192 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Violle, C. et al. The return of the variance: intraspecific variability in community ecology. Trends Ecol. Evol. 27, 244–252 (2012).Article 
PubMed 

Google Scholar 
Díaz, S. & Cabido, M. Vive la différence: plant functional diversity matters to ecosystem processes. Trends Ecol. Evol. 16, 646–655 (2001).Article 

Google Scholar 
Stuart-Smith, R. D. et al. Integrating abundance and functional traits reveals new global hotspots of fish diversity. Nature 501, 539–542 (2013).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Hillebrand, H., Bennett, D. M. & Cadotte, M. W. Consequences of dominance: a review of evenness effects on local and regional ecosystem processes. Ecology 89, 1510–1520 (2008).Article 
PubMed 

Google Scholar 
Loreau, M. The Challenges of Biodiversity Science. Excellence in Ecology Series (International Ecology Institute, 21385 Oldendorf/Luhe, Germany, 2010).Hulshof, C. M. et al. Intra-specific and inter-specific variation in specific leaf area reveal the importance of abiotic and biotic drivers of species diversity across elevation and latitude. J. Veg. Sci. 24, 921–931 (2013).Article 

Google Scholar 
Siefert, A. et al. A global meta-analysis of the relative extent of intraspecific trait variation in plant communities. Ecol. Lett. 18, 1406–1419 (2015).Article 
PubMed 

Google Scholar 
Dall, S. R. X., Bell, A. M., Bolnick, D. I. & Ratnieks, F. L. W. An evolutionary ecology of individual differences. Ecol. Lett. 15, 1189–1198 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
Bolnick, D. I. & Ballare, K. M. Resource diversity promotes among individual diet variation, but not genomic diversity, in lake stickleback. Ecol. Lett. 23, 495–505 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Lande, R. & Arnold, S. J. The measurement of selection on correlated characters. Evolution 37, 1210–1226 (1983).Article 
PubMed 

Google Scholar 
Mullon, C. & Lehmann, L. An evolutionary quantitative genetics model for phenotypic (co)variance under limited dispersal, with an application to socially synergistic traits. Evolution 73, 1695–1728 (2019).Article 
PubMed 

Google Scholar 
Taper, M. L. & Case, T. J. Quantitative genetic models for the coevolution of character displacement. Ecology 66, 355–371 (1985).Article 

Google Scholar 
Engen, S., Grotan, V., Saether, B.-E. & Coste, C. F. D. An evolutionary and ecological community model for distribution of phenotypes and abundances among competing species. Am. Natur. 198, 1 (2021). https://doi.org/10.1086/714529.Kohyama, T. & Takada, T. The stratification theory for plant coexistence promoted by one-sided competition. J. Ecol. 97, 463–471 (2009).Article 

Google Scholar 
Kinzig, A. P., Levin, S. A., Dushoff, J. & Pacala, S. W. Limiting similarity, species packing, and system stability for hierarchical competition-colonization models. Am. Nat. 153, 371–383 (1999).CAS 
Article 
PubMed 

Google Scholar 
Adler, F. R. & Mosquera, J. Is space necessary? Interference competition and limits to biodiversity. Ecology 81, 3226–3232 (2000).Article 

Google Scholar 
Suding, K. N. et al. Scaling environmental change through the community-level: a trait-based response-and-effect framework for plants. Glob. Change Biol. 14, 1125–1140 (2008).ADS 
Article 

Google Scholar 
Loreau, M. & Hector, A. Partitioning selection and complementarity in biodiversity experiments. Nature 412, 72–76 (2001).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Parent, C. E. & Crespi, B. J. Ecological opportunity in adaptive radiation of Galápagos endemic land snails. Am. Nat. 174, 898–905 (2009).Article 
PubMed 

Google Scholar 
Geist, D. J., Snell, H., Snell, H., Goddard, C. & Kurz, M. D. A. Paleogeographic Model of the Galápagos Islands and Biogeographical and Evolutionary Implications. In Geophysical Monograph Series, (eds Harpp, K. S., Mittelstaedt, E., d’Ozouville, N. & Graham, D. W.), chap. 8, 145–166 (2014).Parent, C. E. & Crespi, B. J. Sequential colonization and diversification of Galápagos endemic land snail genus Bulimulus (Gastropoda, Stylommatophora). Evolution 60, 2311–2328 (2006).CAS 
PubMed 

Google Scholar 
Parent, C. E. Diversification on islands: bulimulid land snails of Galápagos. Ph.D. thesis, Simon Fraser University, Burnaby, Canada (2008).Kraemer, A. C., Roell, Y. E., Shoobs, N. F. & Parent, C. E. Does island ontogeny dictate both the accumulation of species richness and functional diversity? Glob. Ecol. Biogeogr. 31, 123–137 (2021).Kraemer, A. C., Philip, C. W., Rankin, A. M. & Parent, C. E. Trade-offs direct the evolution of coloration in Galápagos land snails. Proc. R. Soc. B 286, 20182278 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Barabás, G. & D’Andrea, R. The effect of intraspecific variation and heritability on community pattern and robustness. Ecol. Lett. 19, 977–986 (2016).Article 
PubMed 

Google Scholar 
Barton, N. H., Etheridge, A. M. & Véber, A. The infinitesimal model: definition, derivation, and implications. Theor. Popul. Biol. 118, 50–73 (2017).CAS 
MATH 
Article 
PubMed 

Google Scholar 
Govaert, L. et al. Eco-evolutionary feedbacks—theoretical models and perspectives. Funct. Ecol. 33, 13–30 (2019).Article 

Google Scholar 
Keddy, P. A. & Shipley, B. Competitive hierarchies in herbaceous plant communities. Oikos 54, 234–241 (1989).Article 

Google Scholar 
Allesina, S. et al. Predicting the stability of large structured food webs. Nat. Commun. 6, 7842 (2015).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Kandlikar, G. S., Johnson, C. S., Yan, X., Kraft, N. J. B. & Levine, J. M. Winning and losing with microbes: how microbially mediated fitness differences influence plant diversity. Ecol. Lett. 22, 1178–1191 (2019).PubMed 

Google Scholar 
Reich, P. B. et al. Impacts of biodiversity loss escalate through time as redundancy fades. Science 336, 589–592 (2012).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Spaak, J. W. & De Laender, F. Effects of pigment richness and size variation on coexistence, richness and function in light limited phytoplankton. J. Ecol. 109, 2385–2394 (2021).Article 

Google Scholar 
Parain, E. C., Rohr, R. P., Gray, S. M. & Bersier, L.-F. Increased temperature disrupts the biodiversity–ecosystem functioning relationship. Am. Nat. 193, 227–239 (2019).Article 
PubMed 

Google Scholar 
Gonzalez, A. et al. Scaling-up biodiversity-ecosystem functioning research. Ecol. Lett. 23, 757–776 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Leibold, M. A., Urban, M. C., De Meester, L., Klausmeier, C. A. & Vanoverbeke, J. Regional neutrality evolves through local adaptive niche evolution. Proc. Natl Acad. Sci. USA 116, 2612–2617 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Dieckmann, U. & Doebeli, M. On the origin of species by sympatric speciation. Nature 400, 354–357 (1999).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Edwards, K. F. et al. Evolutionarily stable communities: a framework for understanding the role of trait evolution in the maintenance of diversity. Ecol. Lett. 21, 1853–1868 (2018).Article 
PubMed 

Google Scholar 
Bolnick, D. I., Svanbäck, R., Araujo, M. S. & Persson, L. Comparative support for the niche variation hypothesis that more generalized populations also are more heterogeneous. Proc. Natl Acad. Sci. USA 104, 10075–10079 (2007).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Van Valen, L. Morphological variation and width of ecological niche. Am. Nat. 99, 377–390 (1965).Article 

Google Scholar 
Goodfriend, G. A. Variation in land-snail shell form and size and its causes: a review. Syst. Biol. 35, 204–223 (1986).Article 

Google Scholar 
Machin, J. Structural adaptation for reducing water-loss in three species of terrestrial snail. J. Zool. 152, 55–65 (1967).Article 

Google Scholar 
McMahon, R. F. Thermal tolerance, evaporative water loss, air-water oxygen consumption and zonation of intertidal prosobranchs: a new synthesis. In Progress in Littorinid and Muricid Biology, 241–260 (Springer, Dordrecht, The Netherlands, 1990).Rees, B. B. & Hand, S. C. Heat dissipation, gas exchange and acid-base status in the land snail oreohelix during short-term estivation. J. Exp. Biol. 152, 77–92 (1990).Article 

Google Scholar 
Newkirk, G. F. & Doyle, R. W. Genetic analysis of shell-shape variation in Littorina saxatilis on an environmental cline. Mar. Biol. 30, 227–237 (1975).Article 

Google Scholar 
Seeley, R. H. Intense natural selection caused a rapid morphological transition in a living marine snail. Proc. Natl Acad. Sci. USA 83, 6897–6901 (1986).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Schmidt-Nielsen, K., Taylor, C. R. & Shkolnik, A. Desert snails: problems of heat, water and food. J. Exp. Biol. 55, 385–398 (1971).CAS 
Article 
PubMed 

Google Scholar 
Xavier Jordani, M. et al. Intraspecific and interspecific trait variability in tadpole metacommunitiees from the Brazilian Atlantic rainforest. Ecol. Evol. 9, 4025–4037 (2019).Article 
PubMed 
PubMed Central 

Google Scholar  More