Philippa Kaur

Classification of snow algae in the ice core based on ITS2 sequencesWe used high-throughput sequencing to obtain DNA sequences of algae from 19 layers of an ice core drilled on a glacier in central Asia, dated from present time to 8000 years ago (Fig. 1 and Table S1). In total, 17,016 unique sequences (phylotypes) for the fast-evolving algal nuclear rDNA internal transcribed spacer 2 (ITS2) region were determined in the ice core, from which 290 OTUs were defined with ≥98% nt sequence identity among all OTUs.The ITS2 sequences were classified at the species level according to the genetic species concept based on secondary structural differences in the ITS2 region, which correlate with the boundaries of most biological species [38]. The ITS2 sequences from ice core samples were classified into 24 subgroups consisting of 17 chlorophycean, 5 trebouxiophycean, and 2 ulvophycean groups based on their secondary structures and BLASTn results (Fig. S1 and Tables S3–S4). The 17 subgroups of Chlorophyceae were subdivided into 10 subgroups of the Chloromonadinia clade, 1 subgroup of the Monadinia clade (recently treated as the genus Microglena [54]), 3 subgroups of the Reinhardtinia clade, 2 subgroups of the Stephanosphaerinia clade, and 1 subgroup corresponding to an unnamed group (which is related to Ploeotila sp. CCCryo 086-99) (for the clade names, see [55]). Although the Chloromonadinia clade contains several snow species belonging to Chloromonas or Chlainomonas, the 10 subgroups of the Chloromonadinia clade were considered to be Chloromonas. The 5 trebouxiophycean subgroups were composed of 2 subgroups of the Chlorella group, 1 subgroup of the Raphidonema group, 1 subgroup of the Trebouxia group, and 1 subgroup of the Neocystis group. The 2 subgroups of Ulvophyceae were closely related to the genus Chamaetrichon and Planophila, respectively. It is noted that Sanguina (‘Chlamydomonas’-snow group B [6]), Ancylonema, and Mesotaenium, which are snow algal genera found throughout the world [56, 57], were not detected in the ice core samples (Tables S3–S4).Global distribution of the Raphidonema groupTo understand the process by which snow algae form geographically specific population structures and how they migrate globally across the glaciers and snow fields, it is necessary to focus on the microbial species that inhabit the global cryosphere. Previous work elucidated that the Raphidonema group and ‘Chlamydomonas’-snow group B (Sanguina) are the cosmopolitans at both poles [6], but the latter was not detected in ice core samples examined in this study. Therefore, to elucidate the evolutionary history of the Raphidonema group, we further analyzed the ITS2 sequences obtained from the ice core sample as well as the glacier-surface samples from both poles [6] and from the mid-latitudes (samples from 10 sites, obtained in this study) (hereafter, surface samples; Table S2). Members of the Raphidonema group were detected in the older (deep core) layers of the ice core and at the glacier surface of central Asia (Fig. S1 and Tables S3–S4), as well as in the red snow samples from both poles [6]. In central Asia, the Raphidonema group was found in the Russian, Chinese, and Kyrgyz samples but was not detected in the Japanese and Tajik samples (Fig. S1 and Tables S3–S4). Combining these sequences yielded 893,649 reads and 22,389 unique sequences for subsequent detailed analysis (Tables S5–S6). The taxonomic composition of the Raphidonema communities differed among the mid-latitude, ice core, Arctic, and Antarctic samples as determined by PERMANOVA (Table S7). Most of the unique sequences in the Raphidonema group were consistent with an endemic distribution (Tables S8–S10). An average of 77% of the unique sequences of the Raphidonema group were endemic to a specific region (mid-latitude, 96%; Antarctic, 66%; Arctic, 79%), accounting for 40% of the total sequencing reads (mid-latitude, 77%; Antarctic, 74%; Arctic, 22%) (Fig. 2a, b and Tables S9–S10). This result suggested that most of the unique sequences are endemic, indicating that their dispersal has been limited to their respective regions [58,59,60,61].Fig. 2: Distribution types of the Raphidonema group obtained from each region and the ice core based on ITS2 unique sequences.Proportions of unique sequence and number of sequencing reads are shown. a Unique sequences from surface snow and ice-core samples. b Number of sequencing reads from surface snow and ice core samples. c Unique sequences from the indicated locations within the ice core. d Number of sequencing reads of the unique sequences from the indicated locations within the ice core.Full size imageNext, we analyzed the global distribution of the cosmopolitan phylotypes of the Raphidonema group, because a previous study analyzed their distribution only at the poles [6]. Only a limited number of unique sequences were distributed in all regions (mid-latitude, 1.4%; Antarctic, 5.6%; Arctic, 3.1%), accounting for a large proportion of the sequencing reads in polar regions but for only a small proportion in the mid-latitudes (mid-latitude, 2.8%; Antarctic, 20%; Arctic, 55%) (Figs. 2a, b, S2–S3, and Tables S9–S10). The distribution types of the Raphidonema group obtained from each region and the ice core were similar between the USEARCH and DADA2 analyses (Figs. 2, S4). In addition, we note that in ancient samples, post-mortem nt substitutions, such as cytosine to thymine, accumulate over many years of deposition [62], and these are not included in the DADA2 error model, which leads to the elimination of minor sequences in the DADA2 analysis. Therefore, we based our analysis on the results of the USEARCH unique sequences. These results suggested that only a few snow algae in the Raphidonema group were detected in samples from the mid-latitude regions.Snow algae of the Raphidonema group were detected in different ice core layers, corresponding to different time periods. The ice core records revealed that the distribution types of the Raphidonema group have not changed significantly for the last 8000 years, with p = 0.1924 based on a PERMANOVA between the newer (1800–2001 AD) and the older (6000–8000 years before present) layers (Fig. 2c, d). In ice core samples, 77% of the unique sequences of the Raphidonema group were detected only in the ice core samples, accounting for 23% of the total sequencing reads (Fig. S5). Although some of these unique sequences may be artifacts of the post-mortem nt substitution or sequencing errors, because we conducted sequence quality filtering and removed the majority of artifact sequences by removing the singleton clusters, most of the unique sequences in the ice core are not likely to be artifacts, but they could represent endemic phylotypes (Figs. 2a, b, S5).The cosmopolitan phylotypes were detected over a broad period as represented by ice core samples. They were present in approximately similar ratios in the newer and older layers (Fig. 2c, d). The cosmopolitan phylotypes were relatively abundant in the ice core samples (average, 4.0%; range, 0.2–13%), accounting for 13% (0.9–81% in the samples) of the total sequencing reads (Figs. 2c, d and S5).Microevolution of cosmopolitan and endemic phylotypesWe analyzed the evolutionary relationship between cosmopolitan and endemic phylotypes of the Raphidonema group among all snow surface and ice core samples. In total, 22,389 unique sequences of the Raphidonema group were clustered into 170 OTUs that were defined with ≥98% nt sequence identity among sequences within OTUs. The OTU sequences were subdivided into five subgroups (Groups A–E) based on phylogenetic analysis (Figs. S6–S11 and Tables S11–S12). Based on a previous study [63], Groups A–C and Group E were assigned to R. sempervirens and R. nivale, respectively, but Group D was not consistent with any species examined in that study (Fig. S6).The phylotypes were categorized into three subsets: the cosmopolitan phylotypes found at both poles and the mid-latitude regions; the multi-region phylotypes found in any two of the Antarctic, Arctic, and mid-latitude regions; and the endemic phylotypes found in only one of the three regions. Cosmopolitan phylotypes were found in Groups A, B, and C and accounted for 64.6% of the unique sequences. We then analyzed the dispersal of the three groups in detail.MJ networks [47] for the ITS2 sequences in each subgroup revealed that the cosmopolitan phylotypes were located at the center of the networks in Groups A and C that contained any types (endemics, multi-regions, and cosmopolitans) of the phylotypes, whereas the endemic phylotypes were considered to be derived from the cosmopolitan phylotypes (Figs. 3 and S12–S13). Moreover, the outgroup phylotypes were directly connected to the cosmopolitan phylotypes. These findings clearly showed that the cosmopolitan phylotypes were ancestral, whereas the endemic phylotypes were derived. In contrast, there were remarkable differences in the shape of the networks between Group B and the others (Groups A and C). In Group B, the Antarctic endemic phylotypes formed a distinct clade, and multi-region phylotypes seemed to be recently derived from this clade. In addition, the Arctic endemic phylotypes formed another distinct clade. These two Group-B clades split directly from a cosmopolitan phylotype (5.3% of the total sequencing reads). For Groups A and C, however, major portions of the total sequencing reads belonged to cosmopolitan phylotypes in Groups A (48.2%) and C (62.4%), and the endemic and multi-region phylotypes were directly connected to these major cosmopolitan phylotypes in a radial manner—the so-called “star-like” pattern [64]. These contrasting network shapes seem to have been formed as a consequence of the unique evolution of each of these groups. We also found that sequences from ice cores did not represent a basal position (Figs. 3 and S12–S13). This is because the haplotypes found in the modern samples have existed from times earlier than the ice core ages, due to the very small mutation numbers expected to have occurred since the ice core ages. Therefore, detected ice core ages were not included in the molecular evolution calculations of our demographic model. However, the phylogenetic networks themselves do not provide information on the evolutionary time scale. Hence, the ice core samples provide further direct evidence that Raphidonema, especially cosmopolitans belonging to this genus, persistently grew on snow and ice at least during the Holocene, and their ITS2 sequences have not changed over the last 8000 years.Fig. 3: Phylogenetic relationships among phylotypes of the Raphidonema groups.Phylotype networks for ITS2 sequences within Groups A (a), B (b), and C (c) of the Raphidonema group that include the cosmopolitan phylotypes in this study. The median-joining method was used. Circles indicate phylotypes; the size of each circle is proportional to the number of unique sequences. Each notch on the edges represents a mutation. Phylotypes are colored according to geographic region. The arrow represents the phylotype in the outgroup (see Fig. S6).Full size imageReferring to “ancestral” phylotypes as those having a longer history than other, more recently derived phylotypes, it is possible that individuals not closely related can share the same ancestral phylotype. In such cases, if genetically far-related individuals from various geographical regions share the same ancestral phylotype, they appear to be “cosmopolitan” (Fig. S14a). In order to distinguish between these “apparent cosmopolitans”, and “true cosmopolitans” that migrate globally, it is necessary to show that the cosmopolitan and endemic phylotypes have distinct demographic histories rather than being part of a continuous population sharing certain demographic dynamics (Fig. S14). Because phylotype networks are not useful for quantifying the rate(s) of microevolution, we used the coalescent model to quantify phylotype demographics [65]. As numerous phylotypes must be analyzed with this approach, we concentrated on statistical inference based on pairwise comparisons of phylotypes, for which the likelihood can be determined in a practical manner (see Materials and Methods). Histograms for the number of mismatched sites between two phylotypes chosen from a set of phylotypes, which will be called the pairwise mismatch distribution, are shown in Figs. 4 and S15. For Groups A and C, the distribution among cosmopolitans, multi-regions, and endemics was unimodal, in which the modes align from left to right with the order cosmopolitans, multi-regions, and endemics. Rogers and Harpending [48] noted that this “wave” propagation results from the expansion in size of a population, which leads to large mismatches, and the mode shifts to the right (see Fig. 2 of [48]). As time passes, the mode shifts to the left and eventually returns to the origin, i.e., representing a population that has not undergone an expansion event. Rogers and Harpending obtained an approximate solution for the wave and fitted the solution to human mitochondrial sequence data. We improved upon their method based on the coalescent model (see Materials and Methods) and applied it to the ITS2 sequence data for snow algae.Fig. 4: Mismatch distribution based on the number of pairwise differences in each distribution type in Raphidonema groups.The lines represent the observed number of pairwise differences in each distribution type (cosmopolitan, multi-region, endemic) within the Raphidonema Groups A (a), B (b) and C (c). Calculations were performed for all distribution types of Raphidonema Groups A and C, for which various cosmopolitan phylotypes were detected. On the other hand, calculations for only multi-region and endemic phylotypes were performed for Raphidonema group B, because no variation was found in cosmopolitan phylotypes.Full size imageFor Group A, when we fit the single demographic model to all phylotypes, the log-likelihood was –414,487. In contrast, when we fit the demographic model to each subset, that is, cosmopolitans, multi-regions, and endemics, separately, the log-likelihood was –341,964. Because the latter is larger than the former, we fit the model to each subset of phylotypes separately. For Group C, when we fit the demographic model to the cosmopolitans, multi-regions, and endemics separately, the log-likelihood was –142,106, which is larger than the log-likelihood, –218,080, when we fit the single demographic model to all phylotypes. In contrast to Groups A and C then, we fit the single demographic model to all phylotypes of Group B because the log-likelihood, –196,070, was larger than the log-likelihood, –220,145, when we fit the demographic model to the cosmopolitans, multi-regions, and endemics separately. These results suggested that cosmopolitans, multi-regions, and endemics experienced different demographic histories in Groups A and C, whereas they had the same demographic history in Group B (Table S13). These results indicate the cosmopolitans in Group A and C are true cosmopolitans, whereas the those in Group B can be regarded as an apparent cosmopolitan.The ML estimates of (tau = 2ut_0), (theta _0 = 2N_0u), and (theta _1 = 2N_1u) are shown in Table S13 with standard deviation values. The population expanded t years ago, with the size before and after the expansion being represented by N0 and N1, respectively. The mutation rate (u) was assumed to be 7.9 × 10–8/ sequence/generation, and the generation interval was assumed to be 24 days (Materials and Methods). In Group A, for the cosmopolitans, the estimates of t, N0, and N1 were (33.8/(2 times 7.9) times 10^8 times {textstyle{{24} over {365}}} = 1.4 times 10^7) years, ((0.108 – 0.010)/(2 times 7.9) times 10^8 = (6.8 – 0.63) times 10^5), and ((0.217)/(2 times 7.9) times 10^8 = 1.4 times 10^6), respectively. In the same way, we computed estimates of t, N0, and N1 of other phylotypes and other groups (Table S14). For the endemics, the respective values were 9.2 × 106 years, 80, and 2.1 × 107, and the values were 4.6 × 106 years, 139, and 1.5 × 107 for the multi-regions. Taking into account the minimum and maximum ranges of the mutation rates per generation as well as the generation intervals, t for cosmopolitans was 3.6 × 106–4.0 × 107 years ago, and t for endemics was 2.3 × 106–2.6 × 107 years ago (Table S14). These results suggested that the cosmopolitans existed at least 1.4 × 107 years ago, and the endemics were derived from the cosmopolitans 9.2 × 106 years ago. The size of the endemics expanded 2.6 × 105-fold, which may have resulted from extensive dispersal. The multi-regions tended to mimic the endemics. Note that our demographic model was simplified to avoid overparameterization. In reality, considering the branching patterns of the MJ network, it is plausible that the endemic phylotypes have been repetitively and continuously derived from the cosmopolitans in multiple lineages—from 9.2 × 106 years ago to the present. In the same way, as for Group C, our results suggested that the cosmopolitan population expanded 3.9-fold ~3.2 × 106 years ago, and the endemics were derived from the cosmopolitans 1.9 × 105 years ago. The size of the endemics expanded 59-fold. In contrast to the phylotypes of Groups A and C, those of Group B experienced no significant expansion (Supplementary Results). In Groups A and C, the derived endemics (and multi-regions) expanded greatly as compared with the ancestral cosmopolitans (Table S14). These extraordinary expansions constitute evidence for local adaptation by the endemic/multi-region populations. In contrast, there was no evidence of local adaptation in Group B. The mismatch distribution of the entire Group B (multi-regions + endemics) showed a multimodal pattern (Fig. 4), which is present in the populations with stable sizes for a long period. When the populations finally reach equilibrium, the mismatch distributions show the exponential distribution [48]. Based on our ML estimates (Table S14), the historical population of Group B has been stable. More

Plaisance, L., Caley, M. J., Brainard, R. E. & Knowlton, N. The diversity of coral reefs: what are we missing? PLoS One. 6, e25026 (2011).Article 
ADS 
CAS 

Google Scholar 
Hoegh-Guldberg, O., Poloczanska, E. S., Skirving, W. & Dove, S. Coral reef ecosystems under climate change and ocean acidification. Front. Mar. Sci. 4, 158 (2017).Article 

Google Scholar 
Mora, C. et al. Global human footprint on the linkage between biodiversity and ecosystem functioning in reef fishes. PLoS Biol. 9, e1000606 (2011).Article 
CAS 

Google Scholar 
Burke, L., Reytar, K., Spalding, M. & Perry, A. Reefs at Risk Revisited. 130 pp. (World Resources Institute, Washington, D.C., 2011).Hicks, C. C., Graham, N. A. J., Maire, E. & Robinson, J. P. W. Secure local aquatic food systems in the face of declining coral reefs. One Earth. 4, 1214–1216 (2021).Article 
ADS 

Google Scholar 
Cinner, J. E. et al. Gravity of human impacts mediates coral reef conservation gains. PNAS 115, E6116–E6125 (2018).Article 
CAS 

Google Scholar 
Eddy, T. D. et al. Global decline in capacity of coral reefs to provide ecosystem services. One Earth. 4, 1278–1285 (2021).Article 
ADS 

Google Scholar 
Graham, N. A. J. et al. Human disruption of coral reef trophic structure. Curr. Biol. 27, 231–236 (2017).Article 
CAS 

Google Scholar 
Sherman, C. S., Heupel, M. R., Moore, S. K., Chin, A. & Simpfendorfer, C. A. When sharks are away rays will play: effects of top predator removal in coral reef ecosystems. Mar. Ecol. Prog. Ser. 641, 145–157 (2020).Article 
ADS 

Google Scholar 
Ruppert, J. L. W., Travers, M. J., Smith, L. L., Fortin, M.-J. & Meekan, M. G. Caught in the middle: combined Impacts of shark removal and coral loss on the fish communities of coral reefs. PLoS One. 8, e74648 (2013).Article 
ADS 
CAS 

Google Scholar 
Last, P. R. et al. Rays of the World. (CSIRO Publishing, 2016).Ebert, D. A., Dando, M. & Fowler, S. Sharks of the World. 2nd edn, 608 (Princeton University Press, 2021).Heupel, M. R., Lédée, E. J. I. & Simpfendorer, C. A. Telemetry reveals spatial separation of co-occurring reef sharks. Mar. Ecol. Prog. Ser. 589, 179–192 (2018).Article 
ADS 

Google Scholar 
Heupel, M. R., Papastamatiou, Y. P., Espinoza, M., Green, M. E. & Simpfendorfer, C. A. Reef shark science – key questions and future directions. Front. Mar. Sci. 6, 12 (2019).Article 

Google Scholar 
Roff, G., Brown, C. J., Priest, M. A. & Mumby, P. J. Decline of coastal apex shark populations over the past half century. Commun. Biol. 1, 223 (2018).Article 

Google Scholar 
Williams, J. J., Papastamatiou, Y. P., Caselle, J. E., Bradley, D. & Jacoby, D. M. P. Mobile marine predators: an understudied source of nutrients to coral reefs in an unfished atoll. Proc. R. Soc. B. 285, 20172456 (2018).Article 

Google Scholar 
Heithaus, M. R., Wirsing, A. J. & Dill, L. M. The ecological importance of intact top-predator populations: a synthesis of 15 years of research in a seagrass ecosystem. Mar. Freshw. Res. 63, 1039–1050 (2012).Article 

Google Scholar 
Peel, L. R. et al. Stable isotope analyses reveal unique trophic role of reef manta rays (Mobula alfredi) at a remote coral reef. R. Soc. Open Sci. 6, 190599 (2019).Article 
ADS 
CAS 

Google Scholar 
O’Shea, O. R., Thums, M., van Keulen, M. & Meekan, M. Bioturbation by stingrays at Ningaloo Reef, Western Australia. Mar. Freshw. Res. 63, 189–197 (2012).Article 

Google Scholar 
Takeuchi, S. & Tamaki, A. Assessment of benthic disturbance associated with stingray foraging for ghost shrimp by aerial survey over an intertidal sandflat. Continental Shelf Res. 84, 139–157 (2014).Article 
ADS 

Google Scholar 
Burkholder, D. A., Heithaus, M. R., Fourqurean, J. W., Wirsing, A. & Dill, L. M. Patterns of top-down control in a seagrass ecosystem: could a roving apex predator induce a behaviour-mediated trophic cascade? J. Anim. Ecol. 82, 1192–1202 (2013).Article 

Google Scholar 
Creel, S. & Christianson, D. Relationships between direct predation and risk effects. TRENDS Ecol. Evolution. 23, 194–201 (2008).Article 

Google Scholar 
Ward-Paige, C. A. et al. Large-scale absence of sharks on reefs in the greater-Caribbean: a footprint of human presence. PLoS One. 5, e11968 (2010).Article 
ADS 

Google Scholar 
Espinoza, M., Cappo, M., Heupel, M. R., Tobin, A. J. & Simpfendorfer, C. A. Quantifying shark distribution patterns and species-habitat associations: implications of marine park zoning. PLoS One. 9, e106885 (2014).Article 
ADS 

Google Scholar 
Graham, N. A., Spalding, M. D. & Sheppard, C. R. Reef shark declines in remote atolls highlight the need for multi-faceted conservation action. Aquat. Conserv.: Mar. Freshw. Ecosyst. 20, 543–548 (2010).Article 

Google Scholar 
Nadon, M. O. et al. Re-creating missing population baselines for Pacific reef sharks. Conserv. Biol. 26, 493–503 (2012).Article 

Google Scholar 
MacNeil, M. A. et al. Global status and conservation potential of reef sharks. Nature 583, 801–806 (2020).Article 
ADS 
CAS 

Google Scholar 
Dulvy, N. K. et al. Overfishing drives over one-third of all sharks and rays toward a global extinction crisis. Curr. Biol. 31, 1–15 (2021).Article 

Google Scholar 
Walls, R. H. L. & Dulvy, N. K. Eliminating the dark matter of data deficiency by predicting the conservation status of Northeast Atlantic and Mediterranean Sea sharks and rays. Biol. Conserv. 246, 108459 (2020).Article 

Google Scholar 
Yan, H. F. et al. Overfishing and habitat loss drives range contraction of iconic marine fishes to near extinction. Science Adv. 7, eabb6026, (2021).Butchart, S. H. M. et al. Using Red List Indices to measure progress towards the 2010 target and beyond. Philos. Trans. R. Soc. B 360, 255–268 (2005).Article 
CAS 

Google Scholar 
Sherman, C. S. et al. Taeniura lymma. The IUCN Red List of Threatened Species, eT116850766A116851089 (2021). 10.2305/IUCN.UK.2021-1.RLTS.T116850766A116851089.enPacoureau, N. et al. Half a century of global decline in oceanic sharks and rays. Nature 589, 567–571 (2021).Article 
ADS 
CAS 

Google Scholar 
Cardeñosa, D. et al. Small fins, large trade: a snapshot of the species composition of low-value shark fins in the Hong Kong markets. Anim. Conserv. 23, 203–211 (2019).Article 

Google Scholar 
Haque, A. B. & Spaet, J. L. Y. Trade in threatened elasmobranchs in the Bay of Bengal, Bangladesh. Fish. Res. 243, 106059 (2021).Article 

Google Scholar 
Alcala, A. C. & Russ, G. R. A direct test of the effects of protective management on abundance and yield of tropical marine resources. ICES J. Mar. Sci. 47, 40–47 (1990).Article 

Google Scholar 
Serrano, A. et al. Effects of anti-trawling artificial reefs on ecological indicators of inner shelf fish and invertebrate communities in the Cantabrian Sea (southern Bay of Biscay). J. Mar. Biol. Assoc. U. Kingd. 91, 623–633 (2011).Article 

Google Scholar 
Cortés, E. Perspectives on the intrinsic rate of population growth. Methods Ecol. Evolution. 7, 1136–1145 (2016).Article 

Google Scholar 
McClenachan, L., Cooper, A. B. & Dulvy, N. K. Rethinking trade-driven extinction risk in marine and terrestrial megafauna. Curr. Biol. 26, 1–7 (2016).Article 

Google Scholar 
Tamburello, N., Cote, I. M. & Dulvy, N. K. Energy and the scaling of animal space use. Am. Naturalist 186, 196–211 (2015).Article 

Google Scholar 
Dulvy, N. K. et al. Challenges and priorities in shark and ray conservation. Curr. Biol. 27, R565–R572 (2017).Article 
CAS 

Google Scholar 
Davidson, L. N. K. & Dulvy, N. K. Global marine protected areas to prevent extinctions. Ecol. Evolution. 1, 1–6 (2017).
Google Scholar 
Pauly, D., Zeller, D. & Palomares, M. L. D. Sea Around Us Concepts, Design and Data, (2021).Simpfendorfer, C. A. & Dulvy, N. K. Bright spots of sustainable shark fishing. Curr. Biol. 27, R83–R102 (2017).Article 

Google Scholar 
Booth, H., Squires, D. & Milner-Gulland, E. J. The mitigation hierarchy for sharks: a risk-based framework for reconciling trade-offs between shark conservation and fisheries objectives. Fish. Fish. 21, 269–289 (2019).Article 

Google Scholar 
Grorud-Colvert, K. et al. The MPA Guide: A framework to achieve global goals for the ocean. Science 373, eabf0861 (2021).Article 
CAS 

Google Scholar 
Enright, S. R., Meneses-Orellana, R. & Keith, I. The Eastern Tropical Pacific Marine Corridor (CMAR): The emergence of a voluntary regional cooperation mechanism for the conservation and sustainable use of marine biodiversity within a fragmented regional ocean governance landscape. Front. Mar. Sci. 8, 674825 (2021).Article 

Google Scholar 
Chin, A., Kyne, P. M., Walker, T. I. & McAuley, R. B. An integrated risk assessment for climate change: analysing the vulnerability of sharks and rays on Australia’s Great Barrier Reef. Glob. Change Biol. 16, 1936–1953 (2010).Article 
ADS 

Google Scholar 
Dwyer, R. G. et al. Individual and population benefits of marine reserves for reef sharks. Curr. Biol. 30, 480–489 (2020).Article 
CAS 

Google Scholar 
Speed, C. W., Cappo, M. & Meekan, M. G. Evidence for rapid recovery of shark populations within a coral reef marine protected area. Biol. Conserv. 220, 308–319 (2018).Article 

Google Scholar 
Mizrahi, M. I., Diedrich, A., Weeks, R. & Pressey, R. L. A systematic review of the socioeconomic factors that influence how marine protected areas impact on ecosystems and livelihoods. Soc. Nat. Resour. 32, 4–20 (2019).Article 

Google Scholar 
IPBES. Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. 1148 (Bonn, Germany, 2019).Butchart, S. H. M. et al. Global biodiversity: indicators of recent declines. Science 328, 1164–1168 (2010).Article 
ADS 
CAS 

Google Scholar 
Hanh, T. T. H. & Boonstra, W. J. What prevents small-scale fishing and aquaculture households from engaging in alternative livelihoods? A case study in the Tam Giang lagoon, Viet Nam. Ocean Coast. Manag. 182, 104943 (2019).Article 

Google Scholar 
Ahmed, N., Troell, M., Allison, E. H. & Muir, J. F. Prawn postlarvae fishing in coastal Bangladesh: challenges for sustainable livelihoods. Mar. Policy. 34, 218–227 (2010).Article 

Google Scholar 
Prasetyo, A. P. et al. Shark and ray trade in and out of Indonesia: addressing knowledge gaps on the path to sustainability. Mar. Policy. 133, 104714 (2021).Article 

Google Scholar 
McClanahan, T., Polunin, N. & Done, T. Ecological states and the resilience of coral reefs. Conserv. Ecol. 6, 18 (2002).
Google Scholar 
Bellwood, D. R., Hughes, T. P. & Hoey, A. S. Sleeping functional group drives coral-reef recovery. Curr. Biol. 16, 2434–2439 (2006).Article 
CAS 

Google Scholar 
Cinner, J. E. et al. Vulnerability of coastal communities to key impacts of climate change on coral reef fisheries. Glob. Environ. Change. 22, 12–20 (2012).Article 
ADS 

Google Scholar 
Víe, J.-C., Hilton-Taylor, C. & Stuart, S. N. Wildlife in a Changing World – An analysis of the 2008 IUCN Red List of Threatened Species. 180 (Gland, Switzerland, 2009).Mace, G. M. et al. Quantification of extinction risk: IUCN’s system for classifying threatened species. Conserv. Biol. 22, 1424–1442 (2008).Article 

Google Scholar 
Sherley, R. B. et al. Estimating IUCN Red List population reduction: JARA – A decision-support tool applied to pelagic sharks. Conserv. Lett. 13, e12688 (2019).
Google Scholar 
IUCN Red List. Threats Classification Scheme (Version 3.2), (2021).Salafsky, N. et al. A standard lexicon for biodiversity conservation: unified classifications of threats and actions. Conserv. Biol. 22, 897–911 (2008).Article 

Google Scholar 
Moore, A. Chiloscyllium arabicum. The IUCN Red List of Threatened Species 2017, e.T161426A109902537 (2017). 10.2305/IUCN.UK.2017-2.RLTS.T161426A109902537.enSadovy de Mitcheson, Y. J. et al. Valuable but vulnerable: Over-fishing and under-management continue to threaten groupers so what now? Mar. Policy. 116, 103909 (2020).Article 

Google Scholar 
R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2021).Regression Models for Ordinal Data v. 2019.12.10 (CRAN, 2019).Econometric Tools for Performance and Risk Analysis v. 2.0.4 (2020).Dormann, C. F. et al. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46 (2013).Article 

Google Scholar 
Akinwande, M. O., Dikko, H. G. & Samson, A. Variance inflation factor: As a condition for the inclusion of suppressor variable(s) in regression analysis. Open J. Stat. 5, 754–767 (2015).Article 

Google Scholar 
Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods Res. 33, 261–304 (2004).Article 
MathSciNet 

Google Scholar 
Plots Coefficients from Fitted Models v. 1.2.8 (2022).Fisheries and Aquaculture Software. FishStatJ – Software for Fishery and Aquaculture Statistical Time Series., (2020).Reynolds, R. W., Rayner, N. A., Smith, T. M., Stokes, D. C. & Wang, W. An improved in situ and satellite SST analysis for climate. J. Clim. 15, 1609–1625 (2002).Article 
ADS 

Google Scholar 
NASA Ocean Biology (OB.DAAC). Mean annual sea surface chlorophyll-a concentration for the period 2009-2013 (composite dataset created by UNEP-WCMC). Data obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua Ocean Colour website (NASA OB.DAAC, Greenbelt, MD, USA), (2014).General Bathymetric Chart of the Oceans. GEBCO_2014 Grid. version 20150318. www.gebco.net (2015).XGBoost: A Scalable Tree Boosting System v. 1.4.1.1 (In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (pp. 785–794). New York, NY, USA: ACM, 2016).Elith, J., Leathwick, J. R. & Hastie, T. A working guide to boosted regression trees. J. Anim. Ecol. 77, 802–813 (2008).Article 
CAS 

Google Scholar 
ArcGIS Pro 2.7.0 (Environmental Systems Research Institute) (2020).Ferreira, L. C. & Simpfendorer, C. Galeocerdo cuvier. The IUCN Red List of Threatened Species 2019, e.T39378A2913541 (2019).Beta Regression v. 3.1-4 (2021).Butchart, S. H. et al. Improvements to the Red List Index. PLoS ONE. 2, e140 (2007).Article 
ADS 

Google Scholar 
Sherman, C. S. et al. Half a century of rising extinction risk of coral reef sharks and rays, sammsherman27/CoralReefSharkRayIUCN: Data and Code Used in Sherman et al. Half a century of rising extinction risk of coral reef sharks and rays v1.0.0. https://doi.org/10.5281/zenodo.7267904 (2022). More

With human movement and globalization, invasive container breeding vectors responsible for dengue, Zika, chikungunya and now malaria, with An. stephensi, are being introduced and establishing populations in new locations. They are bringing with them the threat of increasing or novel cases of vector-borne diseases to new locations where health systems may not be prepared.Anopheles stephensi was first detected on the African continent in Djibouti in 2012 and has since been confirmed in Ethiopia, Somalia, and Sudan. Unlike most malaria vectors, An. stephensi is often found in artificial containers and in urban settings. This unique ecology combined with its initial detection in seaports in Djibouti, Somalia, and Sudan has led scientists to believe that the movement of this vector is likely facilitated through maritime trade.By modeling inter- and intra-continental maritime connectivity in Africa we identified countries with higher likelihood of An. stephensi introduction if facilitated through maritime movement and ranked them based on this data. Anopheles stephensi was not detected in Africa (Djibouti) until 2012. To determine whether historical maritime data would have identified the first sites of introduction, 2011 maritime data were analyzed to determine whether the sites with confirmed An. stephensi would rank highly in connectivity to An. stephensi endemic countries. Using 2011 data on maritime connectivity alone, Djibouti and Sudan were identified as the top two countries at risk of An. stephensi introduction if it is facilitated by marine cargo shipments. In 2021, these are two of the three African coastal nations where An. stephensi is confirmed to be established.When 2011 maritime data were combined with the HSI for An. stephensi establishment, the top five countries remain the same as with maritime data alone: Sudan, Djibouti, Egypt, Kenya and Tanzania, in that order. The maritime data show likelihood of introduction and HSI shows likelihood of establishment. When combined, the analyses show a likelihood of being able to establish and survive once introduced. Interestingly, the results of the combined analyses align with the detection data being reported in the Horn of Africa. The 2011 maritime data reinforces the validity of the model as it points to Sudan and Djibouti, where An. stephensi established in the following years. Similarly, the HSI data for Ethiopia has aligned closely with detections of the species to date15. Interestingly, around this time of initial detection in Djibouti, Djibouti City port underwent development and organizational change. The government of Djibouti took back administrative control of the port as early as 201230.Following this method, maritime trade data from 2020 could point to countries at risk of An. stephensi introduction from endemic countries as well as from the coastal African countries with newly introduced populations. Here we provide a prioritization list and heat map of countries for the early detection, rapid response, and targeted surveillance of An. stephensi in Africa based on this data and the HSI (Fig. 4). Further invasion of An. stephensi on the African continent has the potential to reverse progress made on malaria control in the last century. Anopheles stephensi thrives in urban settings and in containers, in contrast to the rural settings and natural habitats where most Anopheles spp. are found20. The situation in Djibouti may be a harbinger for what is to come if immediate surveillance and control strategies are not initiated18.Figure 4Prioritization Heat Map of African Countries. These 2020 heat maps rank African countries using (A) the Likelihood of An. stephensi through Maritime Trade Index (LASIMTI) data alone and (B) LASIMTI and HSI combined, based on maritime connectivity to countries where An. stephensi is endemic. Higher ranking countries which are at greater risk of An. stephensi introduction are darker in red color than those that are lower ranking (lighter red). Countries which are shaded grey are inland countries that do not have a coast and therefore no data on maritime movement into ports. Countries which are grey and checkered have established or endemic An. stephensi populations and are considered source locations for potential An. stephensi introduction in this analysis. Map was generated using MapChart (mapchart.net).Full size imageMaritime data from 2020, with Djibouti and Sudan considered as potential source populations for intracontinental introduction of An. stephensi, indicate the top five countries at risk for maritime introduction are Egypt, Kenya, Mauritius, Tanzania, and Morocco, suggesting that targeted larval surveillance in these countries near seaports may provide a better understanding of whether there are maritime introductions. When the data from 2020 data is combined with HSI for An. stephensi, the top five countries are instead Egypt, Kenya, Tanzania, Morocco, and Libya. Interestingly, historical reports of An. stephensi in Egypt exist; however, following further identification these specimens were determined to be An. ainshamsi31. With several suitable habitats both along the coast and inland of Egypt, revisiting surveillance efforts there would provide insight into how countries that are highly connected to An. stephensi locations through maritime traffic may experience introductions.Further field validation of this prioritization list is necessary, because it is possible that An. stephensi is being introduced through other transportation routes, such as dry ports or airports32, or may even be dispersed through wind facilitation33. However, countries highlighted here with high levels of connectivity to known An. stephensi locations should be considered seriously at risk and surveillance urgently established to determine whether An. stephensi introduction has already occurred or to enable early detection. Primary vector surveillance for both Ae. aegypti and An. stephensi are through larval surveys, and the two mosquitoes are commonly detected in the same breeding habitats. It could therefore be beneficial to coordinate with existing Aedes surveillance efforts to be able to simultaneously gather data on medically relevant Aedes vectors while seeking to determine whether An. stephensi is present. Similarly, in locations with known An. stephensi and not well established Aedes programs, coordinating surveillance efforts provides an opportunity to conduct malaria and arboviral surveillance by container breeding mosquitoes simultaneously.Efforts to map pinch points or key points of introduction based on the movement of goods and populations could provide high specificity for targeted surveillance and control efforts. For example, participatory mapping or population mobility data collection methods, such as those used to determine routes of human movement for malaria elimination, may simultaneously provide information on where targeted An. stephensi surveillance efforts should focus. Several methods have been proposed in the literature for modeling human movement and one in particular, PopCAB, which is often used for communicable diseases, combined quantitative and qualitative data with geospatial information to identify points of control34.Data on invasive mosquito species has shown that introduction events are rarely a one-time occurrence. Population genetics data on Aedes species indicate that reintroductions are very common and can facilitate the movement of genes between geographically distinct populations, raising the potential for introduction of insecticide resistance, thermotolerance, and other phenotypic and even behavioral traits which may be facilitated by gene flow and introgression35. Djibouti, Sudan, Somalia, and Ethiopia, countries with established invasive populations of An. stephensi, should continue to monitor invasive populations and points of introduction to control and limit further expansion and adaptation of An. stephensi. Work by Carter et al. has shown that An. stephensi populations in Ethiopia in the north and central regions can be differentiated genetically, potentially indicating that these populations are a result of more than one introduction into Ethiopia from South Asia, further emphasizing the potential role of anthropogenic movement on the introduction of the species17.One major limitation of this work is that Somalia is the third coastal nation where An. stephensi has been confirmed; however, marine traffic data were not available for Somalia so it could not be included in this analysis. The potential impact of Somalia on maritime trade is unknown and it should not be excluded as a potential source population. Additionally, this model does not account for the possibility of other countries with An. stephensi populations that have not been detected yet. As new data on An. stephensi expansion becomes available, more countries will be at higher risk. Other countries with An. stephensi populations, such as Iran, Myanmar, and Iraq, constitute lower relative percentages of trade with these countries so were not included in the analysis. However, genetic similarities were noted from An. stephensi in Pakistan, so this nation was included10.Due to the nature of maritime traffic, inland countries were also not included in this prioritization ranking. Countries which are inland but share borders with high-risk countries according to the LASTIMI index should also be considered with high priority. For example, the ranking from 2011 highlights Sudan and Djibouti, both which border Ethiopia, and efforts to examine key land transportation routes between bordering nations where humans and goods travel may provide additional insight into the expansion routes of this invasive species.In Ethiopia, An. stephensi was detected in 2016. It has largely been detected along major transportation routes although further data is needed to understand the association between movement and An. stephensi introductions and expansion since most sampling sites have also been located along transport routes. Importantly, Ethiopia relies heavily on the ports of Djibouti and Somalia for maritime imports and exports. Surveillance efforts have revealed that the species is also frequently associated with livestock shelters and An. stephensi are frequently found with livestock bloodmeals15. Interestingly, the original detection of An. stephensi was found in a livestock quarantine station in the port of Djibouti. Additionally, livestock constitutes one of the largest exports of maritime trade from this region. For countries with high maritime connectivity to An. stephensi locations, surveillance efforts near seaports, in particular those with livestock trade, may be targeted locations for countries without confirmed An. stephensi to begin larval surveillance.As Ae. aegypti and Culex coronator were detected in tires or Ae. albopictus through tire and bamboo (Dracaena sanderiana) trade, An. stephensi could be carried through maritime trade of a specific good36,37,38. Future examination of the movement of specific goods would be beneficial in interpreting potential An. stephensi invasion pathways. Additionally, the various types of vessels used to transport certain cargo such as container, bulk, and livestock ships could affect An. stephensi survivability during transit. Sugar and grain are often shipped in bulk or break bulk vessels which store cargo in large unpackaged containers. Container ships transport products stored in containers sized for land transportation via trucks and carry goods such as tires. Livestock vessels are often multilevel, open-air ships which require more hands working on deck and water management39.Using LSBCI index data from 2020, we developed a network to highlight how coastal African nations are connected through maritime trade (Fig. 4). The role of this network analysis is two-fold, (1) it demonstrates an understanding of intracontinental maritime connectivity; and (2) it highlights the top three countries connected via maritime trade through an interactive html model (Supplemental File). For example, if An. stephensi is detected and established in a specific coastal African nation such as Djibouti, selecting the Djibouti node reveals the top three locations at risk of introduction from that source country (Djibouti links to Sudan, Egypt and Kenya). This can be used as an actionable prioritization list for surveillance if An. stephensi is detected in any given country and highlights major maritime hubs in Africa which could be targeted for surveillance and control. For example, since the development of this model, An. stephensi has been detected in Nigeria. Through the use of this interactive model, Ghana, Cote d’Ivoire, and Benin have been identified as countries most connected to Nigeria through maritime trade and therefore surveillance prioritization activities could consider these locations.The network analysis reveals the significance of South African trade to the rest of the continent. Due to the distance, South Africa did not appear to be high in risk of An. stephensi introduction. However, this analysis does reveal that if An. stephensi were to enter nearby countries, it could very easily be introduced because of its high centrality. Western African countries such as Ghana, Togo, and Morocco are also heavily connected to other parts of Africa. Interestingly, Mauritius appears to be highly significant to this network of African maritime trade. Based on 2020 maritime data, Mauritius is ranked as the country with the third greatest likelihood of introduction of An. stephensi and has the second highest centrality rank value of 0.159. Considering these factors, Mauritius could serve as an important port of call connecting larger ports throughout Africa or other continents. With long standing regular larval surveillance efforts across the island for Aedes spp., this island nation is well suited to look for Anopheles larvae as part of Aedes surveillance efforts for early detection and rapid response to prevent the establishment of An. stephensi. If An. stephensi were to become established in countries with high centrality ranks, further expansion on the continent could be accelerated drastically. These ports could serve as important watchpoints and indicators of An. stephensi’s incursion into Africa. Anopheles stephensi is often found in shared habitats with Aedes spp. and a great opportunity exists to leverage Aedes arboviral surveillance efforts to initiate the search for An. stephensi, especially in countries that have high potential of introduction through maritime trade. More

Data collectionWe obtained trade data from two different sources: the United States Fish and Wildlife Service (USFWS) Law Enforcement Management Information System (LEMIS)31 and the International Union for Conservation of Nature (IUCN) Red List32. We used the former to obtain data on the live wildlife trade in general and the latter for data on the pet trade specifically. We then matched trade data with our previously compiled global scale datasets of life history traits and introductions in mammals, reptiles and amphibians25,26.We obtained data on the US live wildlife trade from LEMIS by a Freedom of Information Act Request on 12/08/2019. We requested summary data on all US imports and exports of wildlife across all available years (1999-2019) and all trade purposes, including information on species identities and shipment contents (e.g. live individuals, meat, skins, etc.). For each species, we summed the total number of recorded shipments of live individuals (including individuals that died in transit, and live eggs) as a measure of trade frequency. We classified species as in trade if there was at least one shipment of live individuals recorded in the LEMIS database, and as not traded otherwise. The LEMIS dataset is geographically limited to trade by the US, and therefore may not capture the full diversity of species involved in the wildlife trade. For example, the LEMIS database may be missing some species involved in the substantial trade in live wildlife between South–East Asian countries50. However, the US represents one of the most dominant players in the global market for live wildlife16, and by summing both imports and exports we capture demand for species in countries beyond the US to some extent. Supplementary Fig. 2 illustrates the frequency of trade between the US and countries represented in the US LEMIS dataset. LEMIS data should be considered a minimum estimate of the diversity of species involved in the wildlife trade since they mostly record only legal trade (although confiscated shipments are recorded), and shipments are sometimes not identified to the species level16,51,53,53. The LEMIS database also contains some mis-spelled and incorrectly identified species due to human input errors52. To minimise the effect of misidentified shipments on our species level classifications of US trade status, we discarded all LEMIS records that were not identified to the species level (i.e. those identified using genus, common or generic names only), and manually checked the LEMIS data for synonyms and alternate spellings when we could not automatically match any records in LEMIS with species in our life history datasets. Species classified as traded on the basis of a single recorded live shipment in LEMIS are most vulnerable to species level misclassification due to misidentified shipments. The vast majority of traded species have multiple shipments recorded in LEMIS (259/312 [83%] of traded mammals, 265/285 [93%] of traded reptiles and 72/75 [96%] of traded amphibians), reducing the potential impact of shipment level misidentification over the reliability of species level trade classifications. However, to investigate the robustness of our findings to possible errors in species identification in LEMIS, we re-ran our key analyses excluding species classified as traded on the basis of a single live shipment. We found qualitatively the same effects of life history traits on the probability of trade when removing these species as in our full sample (Supplementary Tables 25–27). Despite its limitations, LEMIS is an invaluable resource for identifying broad scale trends in the wildlife trade since few other countries maintain such detailed records, and it is the only large-scale international trade dataset that includes both CITES- and non-CITES-listed species16,41. Including non-CITES listed species in our analyses is important because CITES-listed species represent only a small minority of those in trade14 and are likely to be a biased sample in terms of life history traits, since species vulnerable to extinction typically have slower life histories40.We obtained separate data on the pet trade from the IUCN Red List. The IUCN has assessed the vast majority of mammal, reptile and amphibian species (91%, 79% and 86% respectively54). Here, we classified a species as involved in the pet trade if the IUCN species account included at least one clear description of involvement in the pet trade. Otherwise, we considered a species as not involved in the pet trade. Although LEMIS records the purpose of trade, it uses broad categories (e.g. ‘Commercial’, ‘Personal’, ‘Breeding in captivity’), none of which refers specifically to nor necessarily equates to trade for pets. Therefore, we sought this additional data on the pet trade from the IUCN Red List instead of following the approach of some previous studies which have used LEMIS data as a proxy for the pet trade (e.g. Refs. 15,19). In contrast, the IUCN Red List contains clear textual descriptions of use and trade for many species, allowing us to identify which species are traded specifically for pets32. The IUCN data has further complementary strengths compared with LEMIS in that it is global in scope and includes both legal and illegal trade. We obtained data from the IUCN Red List by manually searching the binomial name of each species in our samples and consulting the ‘Threats’ and ‘Use and Trade’ sections of the species accounts. We classified species as in the pet trade if the information clearly stated this was the case (e.g. “It has been recorded in the pet trade”, “This species appears in the international pet trade”). We discounted descriptions where the information was uncertain (e.g. the species is described as “probably” or “possibly” traded for pets). We did not count as pets those species that the IUCN categorises as used for “Pets/display animals, horticulture” but which are used only for zoos or captive display, such as beluga whales (Delphinapterus leucas). All species described as pets by the IUCN are ‘exotic’, i.e. those without a long history of domestication14, since the IUCN does not list domesticated species.We matched trade data with our previously published global scale datasets on life history traits and introductions25,26. Internationally traded species may or not be released in the wild outside their native range: some may remain in the confines of captivity (e.g. in zoos or kept by private owners). We defined a species as introduced if there was at least one reliable record of its release, by humans, into the wild outside of its native range, either accidentally or intentionally25,26. We included only species with complete data for the same life history traits as used in our prior analyses (mammals: body mass, gestation period, weaning age, neonatal body mass, litter size, litters per year, age at first reproduction and reproductive lifespan; reptiles: body mass, hatchling mass, clutch size, clutches per year, age of sexual maturity, reproductive lifespan and parity; amphibians: snout-vent length, egg size, clutch size, age of sexual maturity and reproductive lifespan) to facilitate direct comparisons with previous results and to allow us to account for covariation between life history traits55. Species with complete life history data represent 7.8%, 3.5% and 1.6% of the total estimated number of species of mammals, reptiles and amphibians respectively56,57,58. These samples are not random as they over-represent orders containing many species of interest and utility to humans (e.g. ungulates, primates, crocodilians) (Supplementary Tables 28–30). However, these biases are unlikely to undermine our results since we examine life history effects on trade and introduction within these samples. Trade and introduction data do not necessarily cover the same time periods: the US dataset covers only the years 1999-present and the IUCN descriptions also typically refer to recent trade. In contrast, our introduction dataset includes both historical and recent introductions25,26. Therefore, the goal of our analyses is not to test causal hypotheses on the direct relationship between trade and introduction but rather to investigate whether the same life history traits predispose species towards both trade and introduction across diverse taxa, locations and circumstances. When combining the datasets and phylogenies59,60,61,62,63, we resolved species name mis-matches by referring to taxonomic information from the IUCN Red List32, the Global Biodiversity Information Facility (GBIF33) and the Integrated Taxonomic Information System (ITIS64). Table 1 summarises final sample sizes and Supplementary Table 1 the degree of overlap between the trade datasets. Most species in the pet trade are also in the general live wildlife trade, but many more species are traded by the US for general purposes than are involved in the pet trade specifically.Finally, we obtained data for a proxy measure of species detectability in order to control for a potential confounding effect on relationships between life history traits and introduction: larger bodied and longer-lived species may be more likely to be recorded by human observers when introduced compared with smaller and shorter-lived species. We obtained data on species occurrence records, geographic range size and population density, assuming that highly detectable species will have a disproportionately large number of recorded observations than expected based on the size of their geographic ranges and average population densities, following similar approaches by e.g. Refs. 65,66. We obtained occurrence records from the Global Biodiversity Information Facility (GBIF33) via the R package rgbif67 selecting only records resulting from human observation. We obtained range sizes (in decimal degrees squared) from the IUCN Red List32 and processed them for analysis using functions from the rgdal package68, excluding areas of uncertain presence (i.e. limiting range to presence code 1, ‘extant’). We obtained population density estimates from the TetraDENSITY database (version 134), a global database of population density estimates for terrestrial vertebrates. Most species in the TetraDENSITY dataset are represented by estimates from multiple different studies (median = 3, range 1–408). We collapsed density estimates to the species level by taking the median value across studies, including all estimates regardless of sampling method to maximise sample size, and converting all units to individuals/km2 to ensure comparability.Statistical analysesTo investigate relationships between life history traits and trade, we run models treating US or pet trade as the outcome variable and life history traits as the predictors. For all analyses, all life history variables were included in the same models to account for covariation among life history traits55. For US trade, where data on trade frequency are available, we run models both in which trade is treated as a binary variable (traded vs. not traded) and as a count variable (frequency of live shipments, including zero values), while for the pet trade, we have no data on trade frequency and so we treat pet trade as a binary variable only. To investigate the effects of life history traits on introduction, we run models in which introduction is the outcome variable and life history traits are the predictors. In introduction models, we only include traded species (running separate models for the set of species in US trade and the set of species in the pet trade). This approach allows us to disentangle effects associated with trade and introduction and thus identify at which stage(s) life history biases emerge. We also run introduction models in which frequency of US trade is included as an additional predictor alongside life history traits, anticipating that highly traded species are more likely to be introduced. Finally, to investigate possible confounding effects of species detectability on relationships between life history traits and introduction, we investigate effects of number of observations, geographic range size and, where sample sizes allowed, population density on the probability of introduction. If highly detectable species are more likely to be recorded as introduced, we expect to find a positive effect of the number of observations (while accounting for geographic range size and population density) on the probability of introduction. If this effect confounds relationships between body mass/lifespan and introduction, the effect of these life history traits on the probability of introduction should disappear when detectability measures are included in the models alongside life history traits. All analyses were conducted using the R statistical programming environment (Version 4.2.069). Plots were coloured using palettes from the viridis package70.To estimate effects of predictor variables, we fit generalized linear mixed models (GLMMs) using Markov chain Monte-Carlo (MCMC) estimation, implemented in the MCMCglmm package35,36. For analyses with binary outcome variables (traded vs. not traded, introduced vs. not introduced) we fit probit models, while for analyses with US trade frequency as the outcome variable we fit hurdle models. Hurdle models estimate two latent variables: the probability that the outcome is zero (on the logit scale), and the probability of the outcome modelled as a Poisson distribution for non-zero values71. This method therefore allows us to estimate effects of life history traits on the probability and frequency of trade in the same model. While the binary component of a hurdle model estimates the probability that outcomes are zero, when reporting results we reverse the sign of coefficients from the binary model for ease of interpretation, so that effects can be interpreted as the probability that the outcome is not zero. Therefore, here predictors with consistent effects on the probability and frequency of trade in hurdle models will have the same sign (so that if, for example, litter size has a positive effect on both the probability and frequency of trade, both coefficients for litter size from the hurdle model will be positive).Datasets comprising biological measures from multiple related species violate the fundamental statistical assumption that observations are independent of one another, since closely related species are more phenotypically similar than expected by chance due to their shared evolutionary history72. To account for the non-independence of species due to shared ancestry, we included a phylogenetic random effect in all models, represented by a variance-covariance (VCV) matrix derived from the phylogeny. The off-diagonal elements of the VCV matrix contain the amount of shared evolutionary history for each pair of species35,37,38 based on the branch lengths of the phylogeny (here proportional to time)59,61,62,63,63. This approach allows us to estimate phylogenetic signal using the heritability (H2) parameter, which measures the proportion of total variance in the latent variable attributable to the phylogeny35,37,38. Heritability is interpreted in the same way as Pagel’s λ in phylogenetic generalized least squares regression35,37,38,72. Specifically, phylogenetic signal is constrained between 0, indicating no phylogenetic effect so that species can be treated as independent, and 1, indicating that similarity between species is directly proportional to their amount of shared evolutionary history35,38,72. As hurdle models estimate two latent variables, for each hurdle model we report two heritability estimates, one for the binary and one for the Poisson component. All continuous independent variables were log-10 transformed due to positively skewed distributions. Although GLMMs do not require normally distributed predictor variables, log-transforming positively skewed life history predictors in phylogenetic comparative analyses allows us to model life history evolution on proportional rather than absolute scales. This is important as it facilitates biologically meaningful comparisons between species across large scales of life history variation73. Further, log-transforming positively skewed predictors helps to meet assumptions of the underlying Brownian motion model of evolutionary change, which assumes that phenotypic change along branches of the phylogeny is normally distributed74.We calculated variance inflation factors (VIFs) using functions from the car R package75 to check for multicollinearity between predictor variables. Where any model reported a variance inflation factor of 5 or above, indicating potentially problematic levels of collinearity76, we re-ran the model removing the variable with the highest VIF iteratively until all the remaining variables had VIFs of More

West-Eberhard, M. J. Developmental Plasticity and Evolution (Oxford University Press, 2003).Book 

Google Scholar 
Piersma, T. & Van Gils, J. A. The Flexible Phenotype: A Body-Centred Integration of Ecology, Physiology, and Behaviour (Oxford University Press, 2011).
Google Scholar 
Piersma, T. & Drent, J. Phenotypic flexibility and the evolution of organismal design. Trends Ecol. Evol. 18, 228–233 (2003).Article 

Google Scholar 
Tabari, H. Climate change impact on flood and extreme precipitation increases with water availability. Sci. Rep. 10, 1–10 (2020).
Google Scholar 
Beck, H. E. et al. Present and future Köppen–Geiger climate classification maps at 1-km resolution. Sci. Data 5, 180214 (2018).Article 

Google Scholar 
Bogert, C. M. Thermoregulation in reptiles, a factor in evolution. Evolution 3, 195–211 (1949).Article 
CAS 

Google Scholar 
Rensch, B. Das Prinzip geographischer Rassenkreise und das Problem der Artbildung (Gebrueder Borntraeger, 1929).
Google Scholar 
Clusella Trullas, S., van Wyk, J. H. & Spotila, J. R. Thermal melanism in ectotherms. J. Therm. Biol. 32, 235–245 (2007).Article 

Google Scholar 
Delhey, K. A review of Gloger’s rule, an ecogeographical rule of colour: Definitions, interpretations and evidence. Biol. Rev. 94, 1294–1316 (2019).
Google Scholar 
Stuart-Fox, D., Newton, E. & Clusella-Trullas, S. Thermal consequences of colour and near-infrared reflectance. Philos. Trans. R. Soc. B: Biol. Sci. 372, 20160345 (2017).Article 

Google Scholar 
Friedman, N. R. & Remês, V. Ecogeographical gradients in plumage coloration among Australasian songbird clades. Glob. Ecol. Biogeogr. 26, 261–274 (2017).Article 

Google Scholar 
Delhey, K. Darker where cold and wet: Australian birds follow their own version of Gloger’s rule. Ecography 41, 673–683 (2018).Article 

Google Scholar 
Galván, I., Rodríguez-Martínez, S. & Carrascal, L. M. Dark pigmentation limits thermal niche position in birds. Funct. Ecol. 32, 1531–1540 (2018).Article 

Google Scholar 
Medina, I. et al. Reflection of near-infrared light confers thermal protection in birds. Nat. Commun 9, 3610 (2018).Article 
ADS 

Google Scholar 
Aldrich, J. W. & James, F. C. Ecogeographic variation in the American Robin (Turdus migratorius). Auk 108, 230–249 (1991).
Google Scholar 
Morales, H. E. et al. Neutral and selective drivers of colour evolution in a widespread Australian passerine. J. Biogeogr. 44, 522–536 (2017).Article 

Google Scholar 
Griffith, S. C., Owens, I. P. & Burke, T. Environmental determination of a sexually selected trait. Nature 400, 358–360 (1999).Article 
ADS 
CAS 

Google Scholar 
Fargallo, J. A., Laaksonen, T., Korpimäki, E. & Wakamatsu, K. A melanin-based trait reflects environmental growth conditions of nestling male Eurasian kestrels. Evol. Ecol. 21, 157–171 (2007).Article 

Google Scholar 
Fargallo, J. A., Martínez, F., Wakamatsu, K., Serrano, D. & Blanco, G. Sex-dependent expression and fitness consequences of sunlight derived color phenotypes. Am. Nat. 191, 726–743 (2018).Article 

Google Scholar 
Beebe, W. Geographic variation in birds, with especial reference to the effects of humidity. Zoologica 1, 3–41 (1907).
Google Scholar 
Bieber, H. Fellverdunklung beim hauskaninchen nach kälteeinwirkung. Zeitschrift für Säugetierkunde 38, 33–38 (1972).
Google Scholar 
Johnston, R. F. & Selander, R. K. House sparrows: Rapid evolution of races in North America. Science 144, 548–550 (1964).Article 
ADS 
CAS 

Google Scholar 
Galván, I., Wakamatsu, K. & Alonso-Álvarez, C. Black bib size is associated with feather content of pheomelanin in male house sparrows. Pigment Cell Melanoma Res. 27, 1159–1161 (2014).Article 

Google Scholar 
Endler, J. A. On the measurement and classification of colour in studies of animal colour patterns. Biol. J. Linn. Soc. 41, 315–352 (1990).Article 

Google Scholar 
Montgomerie, R. Analyzing colors. In Bird Colouration I. Mechanisms and Measurements (eds Hill, E. G. & McGraw, K. J.) (Harvard University Press, 2006).
Google Scholar 
McGraw, K. J., Dale, J. & Mackillop, E. A. Social environment during molt and the expression of melanin-based plumage pigmentation in male house sparrows (Passer domesticus). Behav. Ecol. Sociobiol. 53, 116–122 (2003).Article 

Google Scholar 
Lessells, C. M. & Boag, P. T. Unrepeatable repeatabilities a common mistake. Auk 104, 116–121 (1987).Article 

Google Scholar 
Anderson, T. R. Biology of the Ubiquitous House Sparrow (Oxford University Press, 2006).Book 

Google Scholar 
Gelman, A. & Hill, J. Data Analysis Using Regression and Multilevel/Hierarchical Models (Cambridge University Press, 2006).Book 

Google Scholar 
Nakagawa, S., Ockendon, N., Gillespie, D. O., Hatchwell, B. J. & Burke, T. Assessing the function of house sparrows’ bib size using a flexible meta-analysis method. Behav. Ecol. 18, 831–840 (2007).Article 

Google Scholar 
Hill, G. E. & McGraw, K. J. Bird Coloration, Volume I: Mechanisms and Measurements (Harvard University Press, 2006).Book 

Google Scholar 
D’Alba, L. & Shawkey, M. D. Melanosomes: Biogenesis, properties, and evolution of an ancient organelle. Physiol. Rev. 99, 1–19 (2018).Article 

Google Scholar 
Delhey, K., Burger, C., Fiedler, W. & Peters, A. Seasonal changes in colour: A comparison of structural, melanin- and carotenoid-based plumage colours. PLoS ONE 5, e11582 (2010).Article 
ADS 

Google Scholar 
Galván, I., Mousseau, T. A. & Møller, A. P. Bird population declines due to radiation exposure at Chernobyl are stronger in species with pheomelanin-based coloration. Oecologia 165, 827–835 (2011).Article 
ADS 

Google Scholar 
Meunier, J., Pinto, S. F., Burri, R. & Roulin, A. Eumelanin-based coloration and fitness parameters in birds: A meta-analysis. Behav. Ecol. Sociobiol. 65, 559–567 (2011).Article 

Google Scholar 
Roulin, A., Almasi, B., Meichtry-Stier, K. S. & Jenni, L. Eumelanin- and pheomelanin-based colour advertise resistance to oxidative stress in opposite ways. J. Evol. Biol. 24, 2241–2247 (2011).Article 
CAS 

Google Scholar 
Gasparini, J. et al. Strength and cost of an induced immune response are associated with a heritable melanin-based colour trait in female tawny owls. J. Anim. Ecol. 78, 608–616 (2009).Article 

Google Scholar 
Fargallo, J. A. et al. Sex-specific phenotypic integration: Endocrine profiles, coloration, and behavior in fledgling boobies. Behav. Ecol. 25, 76–87 (2013).Article 

Google Scholar 
Wittkopp, P. J. & Beldade, P. Development and evolution of insect pigmentation: Genetic mechanisms and the potential consequences of pleiotropy. Semin. Cell Dev. Biol. 20, 65–71 (2009).Article 
CAS 

Google Scholar 
Hubbard, J. K., Uy, J. A. C., Hauber, M. E., Hoekstra, H. E. & Safran, R. J. Vertebrate pigmentation: From underlying genes to adaptive function. Trends Genet. 26, 231–239 (2010).Article 
CAS 

Google Scholar 
McKinnon, J. S. & Pierotti, M. E. Colour polymorphism and correlated characters: Genetic mechanisms and evolution. Mol. Ecol. 19, 5101–5125 (2010).Article 

Google Scholar 
Poston, J. P., Hasselquist, D., Stewart, I. R. & Westneat, D. F. Dietary amino acids influence plumage traits and immune responses of male house sparrows, Passer domesticus, but not as expected. Anim. Behav. 70, 1171–1181 (2005).Article 

Google Scholar 
McGraw, K. J. Dietary mineral content influences the expression of melanin-based ornamental coloration. Behav. Ecol. 18, 137–142 (2007).Article 

Google Scholar 
Fargallo, J. A., Martínez-Padilla, J., Toledano-Díaz, A., Santiago-Moreno, J. & Dávila, J. A. Sex and testosterone effects on growth, immunity and melanin coloration of nestling Eurasian kestrels. J. Anim. Ecol. 76, 201–209 (2007).Article 

Google Scholar 
Fitze, P. S. & Richner, H. Differential effects of a parasite on ornamental structures based on melanins and carotenoids. Behav. Ecol. 13, 401–407 (2002).Article 

Google Scholar 
Roulin, A., Altwegg, R., Jensen, H., Steinsland, I. & Schaub, M. Sex-dependent selection on an autosomal melanic female ornament promotes the evolution of sex ratio bias. Ecol. Lett. 13, 616–626 (2010).Article 

Google Scholar 
Sharma, A. Effect of ambient humidity on UV/visible photodegradation of melanin thin films. Photochem. Photobiol. 86, 852–855 (2010).Article 
CAS 

Google Scholar 
Burtt, E. H. The adaptiveness of animal colors. Bioscience 31, 723–729 (1981).Article 

Google Scholar 
Heppner, F. The metabolic significance of differential absorption of radiant energy by black and white birds. Condor 72, 50–59 (1970).Article 

Google Scholar 
Clusella-Trullas, S., Terblanche, J. S., Blackburn, T. M. & Chown, S. L. Testing the thermal melanism hypothesis: A macrophysiological approach. Funct. Ecol. 22, 232–238 (2008).Article 

Google Scholar 
Zink, R. M. & Remsen, J. V. Evolutionary processes and patterns of geographic variation in birds. Curr. Ornithol. 4, 1–69 (1986).
Google Scholar 
Burtt, E. H. & Ichida, J. M. Gloger’s rule, feather-degrading bacteria, and color variation among song sparrows. Condor 106, 681–686 (2004).Article 

Google Scholar 
Ruiz-De-Castaneda, R., Burtt, E. H. Jr., Gonzalez-Braojos, S. & Moreno, J. Bacterial degradability of an intrafeather unmelanized ornament: A role for feather-degrading bacteria in sexual selection?. Biol. J. Linn. Soc. 105, 409–419 (2012).Article 

Google Scholar 
Goldstein, G. et al. Bacterial degradation of black and white feathers. Auk 121, 656–659 (2004).Article 

Google Scholar 
Ducrest, A. L., Keller, L. & Roulin, A. Pleiotropy in the melanocortin system, coloration and behavioural syndromes. Trends Ecol. Evol. 23, 502–510 (2008).Article 

Google Scholar 
Kim, S. Y., Fargallo, J. A., Vergara, P. & Martínez-Padilla, J. Multivariate heredity of melanin-based coloration, body mass and immunity. Heredity 111, 139–146 (2013).Article 
CAS 

Google Scholar 
Horrocks, N. P. C. et al. Environmental proxies of antigen exposure explain variation in immune investment better than indices of pace of life. Oecologia 177, 281–290 (2015).Article 
ADS 

Google Scholar 
McLean, N., Van Der Jeugd, H. P. & van de Pol, M. High intra-specific variation in avian body condition responses to climate limits generalisation across species. PLoS ONE 13, e0192401 (2018).Article 

Google Scholar 
Gardner, J. L. et al. Spatial variation in avian bill size is associated with humidity in summer among Australian passerines. Clim. Change Responses 3, 11 (2016).Article 

Google Scholar 
Gerson, A. R. et al. Flight at low ambient humidity increases protein catabolism in migratory birds. Science 333, 1434–1436 (2011).Article 
ADS 
CAS 

Google Scholar  More

Here, we examine infection dynamics in the four regions of focus to learn more about how various control interventions performed in each region. Since the first documented cases emerged in these regions, the virus was able to spread freely across much of the first and second phases with a gradual increase in the control interventions. In the four regions of Makkah, Madinah, Eastern, and Riyadh, cases peaked on 12th May (6397 cases; 95% CI 5960–9697), 15th May (1967 cases; 95% CI 1625–2308), 23rd June (10367; 95% CI 8948–11785) and 11th June (11273 cases; 95% CI 11068–12491) respectively according to the fitted model shown in Fig. 2. As the epidemic progressed, more measures were adopted to contain the disease, and the disease’s infectiousness sharply decreased after the third period. There are a few factors responsible for the sudden declining trend: first our model is dependent on official data on cases that have been documented, and these data will only ever reflect a portion of the overall number of cases. Second, different regions developed different testing strategies, and some locations altered their approach to testing during the course of the time period that was investigated. It is possible that the beginning of the Hajj term (period 4) was a contributing factor in the decrease in the number of documented cases. Additionally, previous to this time period, the government indicated that it would be increasing the size of its local testing in order to detect new cases. It is possible that the efficacy of interventions would be reduced if increases are found to be occurring during the falling phase of an epidemic. This may result in measures being kept in place for a longer period of time than they would have been had more accurate data been provided.Figure 2With the use of the Delay Rejection Adaptive Metropolis method, the relevant parameters were estimated for each of the four areas of interest by fitting the data from 13th March until 25th September.Full size imageWe estimate the effective reproduction number (R_t) as an indicator of SARS-CoV-2 transmission before and after the interventions. Figure 3 depicts the dramatic shift in the rate of SARS-CoV-2 transmission as a result of decreased social contact and other control measures. At the beginning of the pandemic, (R_t) for SARS-CoV-2 in Saudi regions was between 4 and 6 as illustrated in Tables 6 and 7. In other words, on average each case spread to between four and six others. Considering that each new generation of SARS-CoV-2 cases occurs every five days, it is evident that this pandemic was rapidly expanding out of control. Moreover, we assumed that the transmission rate and the documented infection rate did not change during the first two periods since interventions were carried out gradually until a complete lockdown took place. As more measures were introduced, the spread of the disease began to decrease. Therefore, our data were based on the weekly reported number of documented SARS-CoV-2 cases broken down by region. As a result, it became clear that the reliability of the (R_t) value was relatively high for transmission.Figure 3Distribution of Rt estimates derived from 10000 MCMC samples for Makkah, Madinah, Eastern, and Riyadh, respectively. The black dot in the centre of each violin plot denotes the median, the thick bar in the plot denotes the interquartile range, and the thin bar in the plot denotes the lowest and maximum values. The mean and the credible interval for 95%, which is shown in parentheses, are labelled below or above, respectively.Full size imageThe effects of the events and interventions on the dynamics of SARS-CoV-2 in the regions of interest are considered. First, if the controls remained in phase four in Makkah, our model projects that the total number of documented cases would increase to 81047 (95% CI 79421–82672). In Al-Madinah, the cumulative number of documented cases would have increased to 22997 (95% CI 19578–26415). The number of cumulative documented cases may have reached 80520 (95% CI 78335–82704) if controls stayed steady in the Eastern region at the level they were at in phase four. If the pattern shown during the fourth period is taken into account, we estimate that there would have been 67150 (95% CI 63731–70568) documented infections in the Riyadh region. Figure 4 illustrates these findings.Figure 4The relevant parameters were estimated for each of the four regions of interest by first fitting the data of each region, and then predicting using the parameters from period 4. This was done with each of the four regions of interest separately.Full size imageWe now explore the impact of controls remaining in place at the same level as that implemented in phase three. In that case, the number of documented cases in Makkah would have increased to 116641 (95% CI 105015–128266). Similarly, the total number of documented cases in Al-Madinah would have increased to 53877 (95% CI 50458-57295) if the outbreak had been allowed to continue at the same level. If the controls had remained unchanged from how they were in phase three in the Eastern region, the total number of documented cases would have been 310459 (95% CI 298362–334981). Finally, in Riyadh this would have resulted in 665241 documented cases (95% CI 651822 to 678659). Figure 5 highlights these findings.Figure 5For each of the four areas of interest, the relevant parameters were estimated by first fitting the data of each region and then predicting using the parameters from period 3. This was carried out for each of the four areas of interest independently.Full size imageWe now investigate the impact of second-period controls remaining in place. In that case, the number of documented cases would increase to 1236642 (95% 1218314–1251626), 442865 (95% CI 439446–456283), 454031 (95% CI 441846–466215), and 2322624 (95% CI 1919206-3026042) in the regions of Makkah, Madinah, Eastern, and Riyadh, respectively (see Fig. 6).Figure 6For each of the four regions of interest, the relevant variables were identified by first fitting the data from each area and then making predictions using the parameters from period 2. Each of the four regions of interest was done separately.Full size imageThe efficacy of NPIs is dependent on when they are adopted, with earlier adoption resulting in greater success in lowering transmission rates of infectious diseases. In the early stages of COVID-19, Saudi regions made the decision to gradually implement measures in order to understand the severity of the disease and reduce the economic and social costs of lockdowns, as well as the political costs. In Fig. 6, if the government were to rely on the interventions of the second phase, then the number of cases of infection would considerably rise owing to the ineffectiveness of the measures. In the third period as in Fig. 5, the government made it possible to relax some of the control measures, but it is ultimately up to each area to decide whether they will maintain the same level of control or whether they will increase or decrease it. In comparison to the control measures carried out during the third and fourth periods, this led to significantly improved outcomes. The reason that these time periods were chosen is that there was no stiffening of the NPI response in most Saudi regions during the first two periods and control interventions were improved later on.Significant undetected infections resulted in the fast spread of new coronaviruses (SARS-CoV-2) which is illustrated in Fig. 7. The proportion of undocumented infections, including asymptomatic cases and undocumented symptomatic individuals who did not seek medical treatment or be tested for mild symptoms, was greater than that of Wuhan at the onset of the pandemic26, which may be a result of the following factors: first, the medical configuration was not optimal and public awareness was limited during the onset of the pandemic while the undocumented rate progressively increased; Second, contact tracing procedures employed in Saudi regions may have become overwhelmed if the number of early-stage cases in Saudi regions rises substantially. The discrepancy between the predicted proportions of asymptomatic (undocumented) cases may be attributable to the difficulty in the un-identifiability of parameters in epidemiological models. There were a substantial number of asymptomatic infected individuals with high infectivity in Saudi regions, where the epidemic situation escalated rapidly. Our research emphasises the frequency of asymptomatic SARS-CoV-2 cases and their role in transmission in order to increase people’s knowledge of asymptomatic cases and to serve as a guide for the prevention and control of SARS-CoV-2.Table 3 Estimated transmission rate in Saudi Regions.Full size tableTable 4 Estimated ascertainable infection rate.Full size tableIn this model, we fitted dynamic transmission rates because of varied preventable measures by the Saudi government at the level of the country or region. After a series of actions taken by the government, regions and cities went into lockdown, resulting in a decrease in the transmission rate as in Table 3. Before the interventions were introduced, in the first two periods of our study, we assumed the transmission rate did not change since individual and community responses had not effectively taken place. After severe interventions were implemented, the transmission rates were allowed to vary in later periods and reduced gradually due to the control measures that reduced the spread of disease27. Estimates of documented infection rates are presented in Table 4. Our model estimates show the documented infection rate has continued to decrease in the last two periods. Thus, the parameters we fit across periods are a measure of how effective the lockdown was in bringing down the documented infection rate28.Risk of resurgenceThe risk of resurgence in Saudi Arabia’s four regions has been examined in this section after the relaxation of intervention measures. There will be a rise in disease activity if control measures are relaxed without taking into account increases in the number of cases being detected, isolated, and/or traced. We predict the first week of no new cases of infection and the week when all current infections in Saudi Arabia will be eradicated.In the Makkah region, had the trend continued into the fourth period, the number of documented infections would have dropped to zero on average by the 6th September (23rd August to 27th September), and all infections would have been eradicated by the 26th of October (7th October to 14th November). On the 28th June, the number of weekly active infections (including presymptomatic, symptomatic, and asymptomatic cases) reached its highest point of 230,230 (95% CI 226811–234364), and on 8th September, that number dropped to 44023 (95% CI 40604–47441).Therefore, the number of documented infections would have reached zero in Al-Madinah region on average on 6th November (23rd October to 22nd November), and all infections would have been eliminated by 1st December (27th November to 14th December). On 23rd June, weekly active infections (including presymptomatic, symptomatic, and asymptomatic cases) peaked at 130,134 (95% CI 126715–133552) and then declined to 60023 (95% CI 58604–63441) on 25th September.If the trend had continued as it did in the fourth period in the Eastern region, the average number of documented infections would have reached zero on 2nd November (from 23rd October to 18th November), and the total eradication of infections would have happened on 1st December (26th November to 22nd December). The number of weekly active infections (including presymptomatic, symptomatic, and asymptomatic cases) peaked at 65000 (95% CI 61581–68418) during the week of July 23rd and subsequently decreased to 800 (95% CI 765–834) on 8th of September.Lastly, the model predicted that the number of weekly active infections in the Riyadh region (including presymptomatic, symptomatic, and asymptomatic infections) peaked on 28th June at 562332 (95% CI 513379–619542) and then decreased to 188215 (95% CI 174796–191633) on 18th September. On average, we expected that the number of documented infections would have decreased to zero on 18th October (7th October to 14th November) and that the total number of infections would have been eliminated on 1st December if the trend continued as it did in the fourth period (20th November to 23rd December). Figure 7 illustrates these findings. We found that if control measures were lifted 30 days following the first day of zero documented cases.Figure 7The estimated number of infected cases that were active (presymptomatic, symptomatic, and asymptomatic) during the research period in the areas of Makkah, Madinah, Eastern, and Riyadh respectively.Full size imageThe probability of resurgence, which we define as the number of active documented cases greater than 100 could be as high as 0.96 in Eastern, 0.95 in Madinah, 0.97 in Makkah, and 0.96 in Riyadh. If we adopt more stringent conditions of lifting controls after observing no confirmed cases for a continuous period of 30 days, the probability of resurgence decreases to 0.31, 0.28, 0.30, and 0.30, with probable resurgence occurring on 13th February, 7th February, 2nd January, and 8th January for Eastern, Makkah, Madinah and Riyadh, respectively (Fig. 8). Despite the use of a simplified model, these results emphasize the hazards of ignoring undetermined occurrences when modifying intervention techniques.Figure 8Figure demonstrating the effect of relaxing all control measures in all four regions 30 days following the first day without confirmed cases.Full size imageSensitivity analysisFor the purpose of testing the robustness of our research results, we conducted a series of sensitivity analyses by varying the durations of the latent and infectious periods, the ratio of transmissibility in asymptomatic (undocumented) cases to symptomatic (documented) cases, and the initial documented infection rate. We conduct eight sensitivity analyses (S1 to S8) within each model for each region of Saudi Arabia to assess the robustness of our model results. For instance, the sensitivity analysis performed for S1 was based on the changes of the latent period and pre-symptomatic infectious period, respectively, and other parameters remain the same. These modifications were carried out with the help of reference15,29, and the same approaches were used for the other parts of the sensitivity analysis, which is summarised in Table 5.Table 5 Description of essential model parameters that were not fitted in the MCMC, where (D_e) refers to the latent period, (D_p) refers to the pre-symptomatic infectious period, (D_i) refers to the symptomatic infectious period, (gamma _0) refers to the initial ascertain rate and (alpha) refers to the ratio of the transmission rate for P and A to I.Full size tableIn particular, for (S1), we raised the incubation period to 7 days (upper 95% CI based on ref15) and the pre-symptomatic infectious period to 3 days (upper 95% CI based on ref29). Therefore we set (D_e = 4) and (D_p=3), and modified (E_0) and (P_0) as needed. The transmissibility of the undocumented cases was assumed to be 0.46 (lower 95 % CI according to ref.31) of the infection cases for (S2); for (S3), the transmissibility of the asymptomatic (undocumented) cases was assumed to be 0.62 (upper 95 % CI according to ref31). We assumed that in (S4), the initial documented infection rate was (gamma _0) = 0.14 (lower 95 % CI according to ref13) and adjusted (A_0), (P_0) and E(0) accordingly. Similarly for (S5) we assumed the initial documented infection rate was (gamma _0) = 0.42 (upper 95 % CI according to ref13) and adjusted (P_0), (A_0), and (E_0) accordingly. In (S6) we set the variables (D_ e=3) and (D_p=1.1), and altered the values of (P_0) and (E_ 0) as necessary in accordance with13. In (S7) we assumed that the transmission rate of asymptomatic (undocumented) cases was half that of documented cases by setting 0.5. Finally, in (S8) we assumed that the infectious period ((D_i)) was double that of symptomatic cases by setting 6 days. Both (S7) and (S8) were based on30. The results of our sensitivity analysis are summarised in Tables 6 and 7. We note that the variation in the model predictions of (R_t) varies from setting to setting. However, these variations appear to be fairly small, proposing the robustness of the results to the specification of associated values in fairly realistic ranges13,13. Our sensitivity analysis provides information about the importance of each parameter to the model representing the transmission of SARS-CoV-2. An increase (or decrease) in parameter values, while other parameters’ values remain the same, contributes to an increase (or decrease) in effective reproduction numbers. For example, an increase in infectious period would result in a higher effective reproduction number at the beginning of the epidemic and a longer time required to clear all infections in Saudi regions32. Our sensitivity analysis indicates that almost all model parameters may have an important role in spreading this virus among susceptible people. In particular, the contact rate from person-to-person and the transition rate of asymptomatic (undetected) individuals play a significant role in disease spread. Our important findings, of a significant decrease in (R_t) after interventions and the existence of a substantial number of presymptomatic and asymptomatic cases, were found to be robust. This highlights that Saudi authorities should pay attention to intervention strategies in the event of a resurgence of cases and quarantining those who were in contact with active cases can effectively reduce the disease33. In Tables 6 and 7 we show the estimated effective reproduction number (R_t) associated with 95% CIs obtained from those eight sensitivity analyses for all four regions and all five time periods.Table 6 Sensitivity analysis of the effective reproduction number (R_t) for Eastern and Madinah.Full size tableTable 7 Sensitivity analysis of estimated effective reproduction number (R_t) for Makkah and Riyadh.Full size table More

Plant Ecology and Nature Conservation Group, Wageningen University, Wageningen, the NetherlandsFons van der Plas & Liesje MommerSystematic Botany and Functional Biodiversity, Life Science, Leipzig University, Leipzig, GermanyThomas Schröder-Georgi, Alexandra Weigelt, Kathryn Barry & Christian WirthGerman Centre for Integrative Biodiversity Research Halle-Jena-Leipzig, Leipzig, GermanyAlexandra Weigelt, Kathryn Barry, Adriana Alzate, Nico Eisenhauer, Anke Hildebrandt, Christiane Roscher & Christian WirthTerrestrial Ecology Research Group, School of Life Sciences Weihenstephan, Technical University of Munich, Munich, GermanySebastian Meyer & Wolfgang WeisserAquaculture and Fisheries Group, Wageningen University and Research Centre, Wageningen, the NetherlandsAdriana AlzateAgroécologie, AgroSup Dijon, Institut National de la Recherche Agronomique, Université de Bourgogne, Université de Bourgogne Franche-Comté, Dijon, FranceRomain L. BarnardEidgenössische Technische Hochschule Zürich, Zurich, SwitzerlandNina BuchmannDepartment of Experimental Plant Ecology, Institute for Water and Wetland Research, Radboud University Nijmegen, Nijmegen, the NetherlandsHans de KroonInstitute of Ecology and Evolution, University Jena, Jena, GermanyAnne Ebeling & Winfried VoigtInstitute of Biology, Leipzig University, Leipzig, GermanyNico EisenhauerHumboldt-Universität zu Berlin, Berlin, GermanyChristof EngelsInstitute of Plant Sciences, University of Bern, Bern, SwitzerlandMarkus FischerMax Planck Institute for Biogeochemistry, Jena, GermanyGerd Gleixner, Ernst-Detlef Schulze & Christian WirthHelmholtz Centre for Environmental Research, Leipzig, GermanyAnke HildebrandtFriedrich Schiller University Jena, Jena, GermanyAnke HildebrandtGeoecology, University of Tübingen, Tübingen, GermanyEva Koller-France & Yvonne OelmannInstitute of Geography and Geoecology, Karlsruhe Institute of Technology, Karlsruhe, GermanySophia Leimer & Wolfgang WilckeEcotron Européen de Montpellier, Centre National de la Recherche Scientifique, Montferrier-sur-Lez, FranceAlexandru MilcuCentre d’Ecologie Fonctionnelle et Evolutive, Unité Mixte de Recherche 5175 (Centre National de la Recherche Scientifique-Université de Montpellier-Université Paul-Valéry Montpellier-Ecole Pratique des Hautes Etudes), Montpellier, FranceAlexandru MilcuDepartment of Evolutionary Biology and Environmental Studies, University of Zurich, Zurich, SwitzerlandPascal A. NiklausUFZ, Helmholtz Centre for Environmental Research, Department Physiological Diversity, Leipzig, GermanyChristiane RoscherInstitute of Landscape Ecology, University of Münster, Münster, GermanyChristoph ScherberCentre for Biodiversity Monitoring, Zoological Research Museum Alexander Koenig, Bonn, GermanyChristoph ScherberGeobotany, Faculty of Biology, University of Freiburg, Freiburg, GermanyMichael Scherer-LorenzenCentre of Biodiversity and Sustainable Land Use, University of Göttingen, Göttingen, GermanyStefan ScheuJ.F. Blumenbach Institute of Zoology and Anthropology, Animal Ecology, University of Göttingen, Göttingen, GermanyStefan ScheuDepartment of Geography, University of Zurich, Zurich, SwitzerlandBernhard SchmidInstitute of Ecology, College of Urban and Environmental Sciences, Peking University, Beijing, ChinaBernhard SchmidLeuphana University Lüneburg, Institute of Ecology, Lüneburg, GermanyVicky TempertonAgroecology, Department of Crop Sciences, University of Göttingen, Göttingen, GermanyTeja TscharntkeF.v.d.P. wrote the initial draft of the manuscript. T.S.-G., A.W., K.B., S.M., A.A., R.L.B., N.B., H.d.K., A.E., N.E., C.E., M.F., G.G., A.H., E.K.-F., S.L., A.M., L.M., P.A.N., Y.O., C.R., C.S., M.S.-L., S.S., B.S., E.-D.S., V.T., T.T., W.V., W. Weisser, W. Wilcke and C.W. helped edit the manuscript. More

The Global Assessment Report on Biodiversity and Ecosystem Services: Summary for Policy-Makers (IPBES, 2019)DeFries, R. & Nagendra, H. Ecosystem management as a wicked problem. Science 356, 265–270 (2017).Article 
CAS 

Google Scholar 
Turkelboom, F. et al. When we cannot have it all: ecosystem services trade-offs in the context of spatial planning. Ecosyst. Serv. 29, 566–578 (2018).Article 

Google Scholar 
Lee, H. & Lautenbach, S. A quantitative review of relationships between ecosystem services. Ecol. Indic. 66, 340–351 (2016).Article 

Google Scholar 
Bennett, E. M., Peterson, G. D. & Gordon, L. J. Understanding relationships among multiple ecosystem services. Ecol. Lett. 12, 1394–1404 (2009).Article 

Google Scholar 
Goldstein, J. H. et al. Integrating ecosystem-service tradeoffs into land-use decisions. Proc. Natl Acad. Sci. USA 109, 7565–7570 (2012).Article 
CAS 

Google Scholar 
Vallet, A., Locatelli, B. & Pramova, E. Ecosystem Services and Social Equity: Who Controls, Who Benefits and Who Loses? (CIFOR, 2020); https://doi.org/10.17528/cifor/007849Neyret, M. et al. Assessing the impact of grassland management on landscape multifunctionality. Ecosyst. Serv. 52, 101366 (2021).Linders, T. E. W. et al. Stakeholder priorities determine the impact of an alien tree invasion on ecosystem multifunctionality. People Nat. 3, 658–672 (2021).Article 

Google Scholar 
Herzig, A., Ausseil, A.-G. & Dymond, J. in Ecosystem Services in New Zealand—Conditions and Trends (ed. Dymond, J. R.) 511–523 (Manaaki Whenua Press, 2014).Chan, K. M. A., Shaw, M. R., Cameron, D. R., Underwood, E. C. & Daily, G. C. Conservation planning for ecosystem services. PLoS Biol. 4, e379 (2006).Article 

Google Scholar 
Pennington, D. N. et al. Cost-effective land use planning: optimizing land use and land management patterns to maximize social benefits. Ecol. Econ. 139, 75–90 (2017).Article 

Google Scholar 
Hölting, L. et al. Including stakeholders’ perspectives on ecosystem services in multifunctionality assessments. Ecosyst. People 16, 354–368 (2020).Article 

Google Scholar 
Plieninger, T. et al. Exploring futures of ecosystem services in cultural landscapes through participatory scenario development in the Swabian Alb, Germany. Ecol. Soc. 18, 39 (2013).Article 

Google Scholar 
Tasser, E., Schirpke, U., Zoderer, B. M. & Tappeiner, U. Towards an integrative assessment of land-use type values from the perspective of ecosystem services. Ecosyst. Serv. 42, 101082 (2020).Article 

Google Scholar 
Sayer, J. et al. Ten principles for a landscape approach to reconciling agriculture, conservation, and other competing land uses. Proc. Natl Acad. Sci. USA 110, 8349–8356 (2013).Article 
CAS 

Google Scholar 
Vallet, A. et al. Linking equity, power, and stakeholders: roles in relation to ecosystem services. Ecol. Soc. 24, 14 (2019).Article 

Google Scholar 
Allan, E. et al. Land use intensification alters ecosystem multifunctionality via loss of biodiversity and changes to functional composition. Ecol. Lett. 18, 834–843 (2015).Article 

Google Scholar 
Hector, A. & Bagchi, R. Biodiversity and ecosystem multifunctionality. Nature 448, 188–190 (2007).Article 
CAS 

Google Scholar 
Manning, P. et al. Redefining ecosystem multifunctionality. Nat. Ecol. Evol. 2, 427–436 (2018).Article 

Google Scholar 
Raudsepp-Hearne, C., Peterson, G. D. & Bennett, E. M. Ecosystem service bundles for analyzing tradeoffs in diverse landscapes. Proc. Natl Acad. Sci. USA 107, 5242–5247 (2010).Article 
CAS 

Google Scholar 
Daniel, T. C. et al. Contributions of cultural services to the ecosystem services agenda. Proc. Natl Acad. Sci. USA 109, 8812–8819 (2012).Article 
CAS 

Google Scholar 
Gunton, R. M. et al. Beyond ecosystem services: valuing the invaluable. Trends Ecol. Evol. 32, 249–257 (2017).Article 

Google Scholar 
Peter, S., Le Provost, G., Mehring, M., Müller, T. & Manning, P. Cultural worldviews consistently explain bundles of ecosystem service prioritisation across rural Germany. People Nat. 4, 218–230 (2022).Article 

Google Scholar 
Haines-Young, R. & Potschin, M. in Ecosystem Ecology (eds Raffaelli, D. G. & Frid, C. L. J.) 110–139 (Cambridge Univ. Press, 2010).Fischer, M. et al. Implementing large-scale and long-term functional biodiversity research: the Biodiversity Exploratories. Basic Appl. Ecol. 11, 473–485 (2010).Article 

Google Scholar 
Wilson, E. O. Half-Earth: Our Planet’s Fight for Life (Norton, 2017).Bastin, J.-F. et al. The global tree restoration potential. Science 365, 76–79 (2019).Article 
CAS 

Google Scholar 
Clapp, J. & Moseley, W. G. This food crisis is different: COVID-19 and the fragility of the neoliberal food security order. J. Peasant Stud. 47, 1393–1417 (2020).Article 

Google Scholar 
Kirwan, J. & Maye, D. Food security framings within the UK and the integration of local food systems. J. Rural Stud. 29, 91–100 (2013).Article 

Google Scholar 
Ellis, E. C. To conserve nature in the Anthropocene, half Earth is not nearly enough. One Earth 1, 163–167 (2019).Article 

Google Scholar 
Boetzl, F. A. et al. A multitaxa assessment of the effectiveness of agri-environmental schemes for biodiversity management. Proc. Natl Acad. Sci. USA 118, e2016038118 (2021).Tyllianakis, E. & Martin-Ortega, J. Agri-environmental schemes for biodiversity and environmental protection: how we are not yet ‘hitting the right keys’. Land Use Policy 109, 105620 (2021).Article 

Google Scholar 
Arroyo-Rodríguez, V. et al. Designing optimal human-modified landscapes for forest biodiversity conservation. Ecol. Lett. 23, 1404–1420 (2020).Article 

Google Scholar 
Gilroy, J. J. et al. Cheap carbon and biodiversity co-benefits from forest regeneration in a hotspot of endemism. Nat. Clim. Change 4, 503–507 (2014).Article 

Google Scholar 
Lindenmayer, D. B. et al. Avoiding bio-perversity from carbon sequestration solutions: avoiding bio-perversity in carbon markets. Conserv. Lett. 5, 28–36 (2012).Article 

Google Scholar 
Stoll-Kleemann, S. & O’Riordan, T. in The Encyclopedia of the Anthropocene Vol. 3 (eds DellaSala, D. A. & Goldstein, M. I.) 347–353 (Elsevier, 2018).Schaich, H., Bieling, C. & Plieninger, T. Linking ecosystem services with cultural landscape research. GAIA 19, 269–277 (2010).Article 

Google Scholar 
O’Connor, L. M. J. et al. Balancing conservation priorities for nature and for people in Europe. Science 372, 856–860 (2021).Article 

Google Scholar 
Büscher, B. et al. Half-Earth or Whole Earth? Radical ideas for conservation, and their implications. Oryx 51, 407–410 (2017).Article 

Google Scholar 
van der Plas, F. et al. Towards the development of general rules describing landscape heterogeneity–multifunctionality relationships. J. Appl. Ecol. 56, 168–179 (2019).Article 

Google Scholar 
Almeida, I., Rösch, C. & Saha, S. Converting monospecific into mixed forests: stakeholders’ views on ecosystem services in the Black Forest Region. Ecol. Soc. 26, 28 (2021).Meyer, M. A. & Früh-Müller, A. Patterns and drivers of recent agricultural land-use change in southern Germany. Land Use Policy 99, 104959 (2020).Article 

Google Scholar 
Kastner, T. et al. Global agricultural trade and land system sustainability: implications for ecosystem carbon storage, biodiversity, and human nutrition. One Earth 4, 1425–1443 (2021).Rasmussen, L. V. et al. Social–ecological outcomes of agricultural intensification. Nat. Sustain. 1, 275–282 (2018).Article 

Google Scholar 
Lindborg, R. et al. How spatial scale shapes the generation and management of multiple ecosystem services. Ecosphere 8, e01741 (2017).Article 

Google Scholar 
Duarte, G. T., Santos, P. M., Cornelissen, T. G., Ribeiro, M. C. & Paglia, A. P. The effects of landscape patterns on ecosystem services: meta-analyses of landscape services. Landsc. Ecol. 33, 1247–1257 (2018).Article 

Google Scholar 
Le Provost, G. et al. The supply of multiple ecosystem services requires biodiversity across spatial scales. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-022-01918-5 (2022).Martin, D. A. et al. Land-use trajectories for sustainable land system transformations: identifying leverage points in a global biodiversity hotspot. Proc. Natl Acad. Sci. USA 119, e2107747119 (2022).Article 
CAS 

Google Scholar 
Seabloom, E. W., Borer, E. T. & Tilman, D. Grassland ecosystem recovery after soil disturbance depends on nutrient supply rate. Ecol. Lett. 23, 1756–1765 (2020).Article 

Google Scholar 
Messinger, J. & Winterbottom, B. African forest landscape restoration initiative (AFR100): restoring 100 million hectares of degraded and deforested land in Africa. Nat. Faune 30, 14–17 (2016).
Google Scholar 
Whittingham, M. J. The future of agri-environment schemes: biodiversity gains and ecosystem service delivery? J. Appl. Ecol. 48, 509–513 (2011).Article 

Google Scholar 
Le Clec’h, S. et al. Assessment of spatial variability of multiple ecosystem services in grasslands of different intensities. J. Environ. Manage. 251, 109372 (2019).Article 

Google Scholar 
Forschungsethische Grundsätze und Prüfverfahren in den Sozial‐ und Wirtschaftswissenschaften Output 9, Berufungsperiode 5 (German Data Forum, 2017).Strukturdaten Reutlingen—Statistisches Bundesamt (Bundeswahlleiter, 2020); https://www.bundeswahlleiter.de/europawahlen/2019/strukturdaten/bund-99/land-8/kreis-8415.htmlStrukturdaten Uckermark—Statistisches Bundesamt (Bundeswahlleiter, 2020); https://www.bundeswahlleiter.de/europawahlen/2019/strukturdaten/bund-99/land-12/kreis-12073.htmlStrukturdaten Unstrut-Hainich-Kreis—Statistisches Bundesamt (Bundeswahlleiter, 2020); https://www.bundeswahlleiter.de/europawahlen/2019/strukturdaten/bund-99/land-16/kreis-16064.htmlBlüthgen, N. et al. A quantitative index of land-use intensity in grasslands: integrating mowing, grazing and fertilization. Basic Appl. Ecol. 13, 207–220 (2012).Article 

Google Scholar 
Ostrowski, A., Lorenzen, K., Petzold, E. & Schindler, S. Land use intensity index (LUI) calculation tool of the Biodiversity Exploratories project for grassland survey data from three different regions in Germany since 2006, BEXIS 2 module. Zenodo https://doi.org/10.5281/zenodo.3865579 (2020).Schall, P. et al. The impact of even‐aged and uneven‐aged forest management on regional biodiversity of multiple taxa in European beech forests. J. Ecol. 55, 267–278 (2018).Statistisches Jahrbuch über Ernährung, Landwirtschaft und Forsten der Bundesrepublik Deutschland Vol. 63 (Bundesministerium für Ernährung und Landwirtschaft, 2019).Simons, N. K. & Weisser, W. W. Agricultural intensification without biodiversity loss is possible in grassland landscapes. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-017-0227-2 (2017).Zinke, O. Heupreise steigen: Futter für die Bauern knapp und teuer. Agrarheute https://www.agrarheute.com/markt/futtermittel/heupreise-steigen-futter-fuer-bauern-knapp-teuer-571946 (2020).Bois de Chez Nous (Lignum, 2021); https://www.lignum.ch/files/images/Downloads_francais/Shop/20010_Bois_de_chez_nous.pdfGerman Timber Company—Internationaler Holzhandel (German Timber Company, 2021); https://www.germantimber.company/. Accessed 2021-11-24Holzeinschlag nach Holzartengruppen, Holzsorten, ausgewählten Besitzarten (Statistisches Bundesamt, 2022); https://www.destatis.de/DE/Themen/Branchen-Unternehmen/Landwirtschaft-Forstwirtschaft-Fischerei/Wald-Holz/Tabellen/holzeinschlag-deutschland.htmlJahresjagdstrecke Bundesrepublik Deutschland, 2019–2020 (Deutsche Jagdverband, 2020); https://www.jagdverband.de/sites/default/files/2021-01/2021-01_Infografik_Jahresjagdstrecke_Bundesrepublik_Deutschland_2019_2020.jpgHeinze, E. et al. Habitat use of large ungulates in northeastern Germany in relation to forest management. For. Ecol. Manage. 261, 288–296 (2011).Article 

Google Scholar 
Conant, R. T., Cerri, C. E. P., Osborne, B. B. & Paustian, K. Grassland management impacts on soil carbon stocks: a new synthesis. Ecol. Appl. 27, 662–668 (2017).Article 

Google Scholar 
Hermes, J., Albert, C. & von Haaren, C. Mapping and Assessing Local Recreation as a Cultural Ecosystem Service in Germany. UVP-Report https://doi.org/10.17442/uvp-report.034.08 (2020).Hermes, J., Albert, C. & von Haaren, C. Assessing the aesthetic quality of landscapes in Germany. Ecosyst. Serv. 31, 296–307 (2018).Article 

Google Scholar 
Ehrhart, S. & Schraml, U. Perception and evaluation of natural forest dynamics. Allg. Forst Jagdztg. 185, 166–183 (2014).
Google Scholar 
Villanueva-Rivera, L. J. & Pijanowski, B. C. soundecology: Soundscape ecology. R package version 1.3.3 (2018).Meyer, S., Wesche, K., Krause, B. & Leuschner, C. Dramatic losses of specialist arable plants in central Germany since the 1950s/60s—a cross-regional analysis. Divers. Distrib. 19, 1175–1187 (2013).Article 

Google Scholar 
Sasaki, K., Hotes, S., Kadoya, T., Yoshioka, A. & Wolters, V. Landscape associations of farmland bird diversity in Germany and Japan. Glob. Ecol. Conserv. 21, e00891 (2020).Article 

Google Scholar 
Peña, L., Casado-Arzuaga, I. & Onaindia, M. Mapping recreation supply and demand using an ecological and a social evaluation approach. Ecosyst. Serv. 13, 108–118 (2015).Article 

Google Scholar 
Schägner, J. P., Brander, L., Paracchini, M.-L., Hartje, V. & Maes, J. Mapping recreational ecosystem services and its values across Europe: a combination of GIS and meta-analysis. In European Association of Environmental and Resource Economists 22nd Annual Conference (2016).R Core Team. R: A Language and Environment for Statistical Computing v.4.2.1 (R Foundation for Statistical Computing, 2022).Rust Programming Language https://www.rust-lang.org/ v 1.44Le Provost, G. et al. Contrasting responses of above- and belowground diversity to multiple components of land-use intensity. Nat. Commun. 12, 3918 (2021).Article 

Google Scholar 
Gini, C. On the measurement of concentration and variability of characters (English translation from Italian by Fulvio de Santis in 2005). Metron 63, 1–38 (1914). More