Philippa Kaur

1.Clarke, A. R., Powell, K. S., Weldon, C. W. & Taylor, P. W. The ecology of Bactrocera tryoni (Diptera: Tephritidae): What do we know to assist pest management?. Ann. Appl. Biol. 158, 26–54. https://doi.org/10.1111/j.1744-7348.2010.00448.x (2011).Article 

Google Scholar 
2.Jessup, A. J. et al. in Area-Wide Control of Insect Pests: From Research to Field Implementation (eds M. J. B. Vreysen, A. S. Robinson, & J. Hendrichs) 685–697 (Springer Netherlands, 2007).3.HIA. Australian Horticulture Statistics Handbook 2019/20. Horticulture Innovation Australia Limited (2020).4.Fanson, B. G., Sundaralingam, S., Jiang, L., Dominiak, B. C. & D’Arcy, G. A review of 16 years of quality control parameters at a mass-rearing facility producing Queensland fruit fly, Bactrocera tryoni. Entomol. Exp. Appl. 151, 152–159. https://doi.org/10.1111/eea.12180 (2014)5.Knipling, E. F. Possibilities of Insect Control or Eradication Through the Use of Sexually Sterile Males. J. Econ. Entomol. 48, 459–462. https://doi.org/10.1093/jee/48.4.459 (1955).Article 

Google Scholar 
6.Shelly, T. & McInnis, D. Sterile Insect Technique and Control of Tephritid Fruit Flies: Do Species With Complex Courtship Require Higher Overflooding Ratios?. Ann. Entomol. Soc. Am. 109, 1–11. https://doi.org/10.1093/aesa/sav101 (2015).CAS 
Article 

Google Scholar 
7.Hendrichs, J. & Robinson, A. in Encyclopedia of Insects (Second Edition) (eds Vincent H. Resh & Ring T. Cardé) 953–957 (Academic Press, 2009).8.Dominiak, B., Clifford, C. S. & Nielsen, S. G. Queensland fruit fly (Bactrocera tryoni Froggatt) attraction to and chemical analysis of male annihilation blocks using three concentratios of cuelure at Dubbo, NSW, Australia. Plant Prot. Q. 24, 157–160 (2009).9.Dominiak, B. C., Ekman, J. & Broughton, S. Mass trapping and other management option for mediterranean fruit fly and Queensland fruit fly in Australia. Gen. Appl. Entomol. 44, 1–8 (2016).
Google Scholar 
10.Khan, M. A. M. et al. Semiochemical mediated enhancement of males to complement sterile insect technique in management of the tephritid pest Bactrocera tryoni (Froggatt). Sci. Rep. 7, 13366. https://doi.org/10.1038/s41598-017-13843-w (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
11.Vargas, R. I., Shelly, T. E., Leblanc, L. & Piñero, J. C. in Vitamins & Hormones Vol. 83 (ed Gerald Litwack) Ch. 23, 575–595 (Academic Press, 2010).12.Benelli, G. et al. Sexual communication and related behaviours in Tephritidae: current knowledge and potential applications for Integrated Pest Management. J. Pest. Sci. 87, 385–405. https://doi.org/10.1007/s10340-014-0577-3 (2014).Article 

Google Scholar 
13.Shelly, T. E. Consumption of methyl eugenol by male Bactrocera dorsalis (Diptera: Tephritidae): low incidence of repeat feeding. Florida Entomologist 77, 201–208 (1994).CAS 
Article 

Google Scholar 
14.Shelly, T. E. Trapping male oriental fruit flies (Diptera: Tephritidae): does feeding on a natural source of methyl eugenol reduce capture probability?. Florida Entomologist 83, 109–111 (2000).Article 

Google Scholar 
15.Tan, K. H., & Toong, Y. C. Floral synomone of a wild orchid, Bulbophyllum cheiri, lures Bactrocera fruit flies for pollination. J. Chem. Ecol. 28, 1161–1172 (2002).16.Shelly, T. E. Evaluation of a Genetic Sexing Strain of the Oriental Fruit Fly as a Candidate for Simultaneous Application of Male Annihilation and Sterile Insect Techniques (Diptera: Tephritidae). J. Econ. Entomol. 113, 1913–1921. https://doi.org/10.1093/jee/toaa099 (2020).CAS 
Article 
PubMed 

Google Scholar 
17.Joseph, R. M. & Carlson, J. R. Drosophila Chemoreceptors: A Molecular Interface Between the Chemical World and the Brain. Trends in genetics : TIG 31, 683–695. https://doi.org/10.1016/j.tig.2015.09.005 (2015).CAS 
Article 
PubMed 

Google Scholar 
18.Venthur, H., & Zhou, J.-J. Odorant receptors and odorant-binding proteins as insect pest control targets: A comparative analysis. Front. Physiol. 9, 1. https://doi.org/10.3389/fphys.2018.01163 (2018).19.Leal, W. S. Odorant reception in insects: roles of receptors, binding proteins, and degrading enzymes. Annu Rev. Entomol. 58, 373–391. https://doi.org/10.1146/annurev-ento-120811-153635 (2013).CAS 
Article 
PubMed 

Google Scholar 
20.Zheng, W. et al. Identification and Expression Profile Analysis of Odorant Binding Proteins in the Oriental Fruit Fly Bactrocera dorsalis. Int. J. Mol. Sci. 14, 14936 (2013).Article 

Google Scholar 
21.Siciliano, P. et al. Identification of pheromone components and their binding affinity to the odorant binding protein CcapOBP83a-2 of the Mediterranean fruit fly, Ceratitis capitata. Insect Biochem. Mol. Biol. 48, 51–62. https://doi.org/10.1016/j.ibmb.2014.02.005 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
22.Siciliano, P. et al. Sniffing Out Chemosensory Genes from the Mediterranean Fruit Fly, Ceratitis capitata. Plos One 9, e85523https://doi.org/10.1371/journal.pone.0085523 (2014)23.Campanini, E. B. & de Brito, R. A. Molecular evolution of Odorant-binding proteins gene family in two closely related Anastrepha fruit flies. BMC Evol. Biol. 16, 198–198. https://doi.org/10.1186/s12862-016-0775-0 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
24.Park, K., et al. Expression patterns of two putative odorant-binding proteins in the olfactory organs of Drosophila have different implications for their functions. Vol. 300 (2000).25.Shanbhag, S. R. et al. Expression mosaic of odorant-binding proteins in Drosophila olfactory organs. Microsc. Res. Tech. 55, 297–306. https://doi.org/10.1002/jemt.1179 (2001).CAS 
Article 
PubMed 

Google Scholar 
26.Wu, Z. et al. Discovery of chemosensory genes in the oriental fruit fly, Bactrocera dorsalis. PloS One 10, e0129794-e0129794.https://doi.org/10.1371/journal.pone.0129794 (2015)27.Liu, Z., Smagghe, G., Lei, Z. & Wang, J.-J. Identification of male- and female-specific olfaction genes in antennae of the Oriental Fruit Fly (Bactrocera dorsalis). PLoS ONE 11, e0147783–e0147783. https://doi.org/10.1371/journal.pone.0147783 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
28.Zhang, J. et al. Identification and expression profiles of novel odorant binding proteins and functional analysis of OBP99a in Bactrocera dorsalis. Arch. Insect Biochem. Physiol. 98, e21452. https://doi.org/10.1002/arch.21452 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
29.Cheng, J. et al. Identification and expression analysis of chemosensory genes in the citrus fruit fly Bactrocera (Tetradacus) minax. PeerJ Preprints 6, e27297v1. https://doi.org/10.7287/peerj.preprints.27297v1 (2018).30.Kumaran, N., Prentis, P. J., Mangalam, K. P., Schutze, M. K. & Clarke, A. R. Sexual selection in true fruit flies (Diptera: Tephritidae): transcriptome and experimental evidences for phytochemicals increasing male competitive ability. Mol. Ecol. 23, 4645–4657. https://doi.org/10.1111/mec.12880 (2014).CAS 
Article 
PubMed 

Google Scholar 
31.Liu, H. et al. BdorOBP2 plays an indispensable role in the perception of methyl eugenol by mature males of Bactrocera dorsalis (Hendel). Sci. Rep. 7, 15894. https://doi.org/10.1038/s41598-017-15893-6 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
32.Kumaran, N. et al. Plant-Mediated Female Transcriptomic Changes Post-Mating in a Tephritid Fruit Fly, Bactrocera tryoni. Genome biology and evolution 10, 94–107.https://doi.org/10.1093/gbe/evx257 (2017)33.Idrees, A. et al. Protein baits, volatile compounds and irradiation influence the expression profiles of odorant-binding protein genes in Bactrocera dorsalis (Diptera: Tephritidae). Appl. Ecol. Environ. Res. 15, 1883–1899 (2017).Article 

Google Scholar 
34.Pavlidi, N. et al. Transcriptomic responses of the olive fruit fly Bactrocera oleae and its symbiont Candidatus Erwinia dacicola to olive feeding. Sci. Rep. 7, 42633. https://doi.org/10.1038/srep42633 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
35.Arya, G. H. et al. The genetic basis for variation in olfactory behavior in Drosophila melanogaster. Chem Senses 40, 233–243. https://doi.org/10.1093/chemse/bjv001 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
36.Akter, H., Adnan, S., Morelli, R., Rempoulakis, P. & Taylor, P. W. Suppression of cuelure attraction in male Queensland fruit flies provided raspberry ketone supplements as immature adults. PLoS ONE 12, e0184086. https://doi.org/10.1371/journal.pone.0184086 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
37.Gilchrist, A. S. et al. The draft genome of the pest tephritid fruit fly Bactrocera tryoni: resources for the genomic analysis of hybridising species. BMC Genomics 15, 1153. https://doi.org/10.1186/1471-2164-15-1153 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
38.Fleischer, J. & Krieger, J. Insect pheromone receptors—Key elements in sensing intraspecific chemical signals. Front. Cell. Neurosci. 12. https://doi.org/10.3389/fncel.2018.00425 (2018).39.Liu, Z. et al. An antennae-specific odorant-binding protein is involved in Bactrocera dorsalis Olfaction. Front. Ecol. Evol. 8.https://doi.org/10.3389/fevo.2020.00063 (2020).40.Tan, K. H. Recaptures of feral Bactrocera dorsalis and B. umbrosa (Diptera: Tephritidae) males after feeding on methyl eugenol. Bull. Entomol. Res. 110, 15–21, https://doi.org/10.1017/S0007485319000208 (2019).41.Liu, L. et al. Contribution of Drosophila DEG/ENaC genes to salt taste. Neuron 39, 133–146. https://doi.org/10.1016/s0896-6273(03)00394-5 (2003).CAS 
Article 
PubMed 

Google Scholar 
42.Swarup, S., Huang, W., Mackay, T. F. C. & Anholt, R. R. H. Analysis of natural variation reveals neurogenetic networks for Drosophila olfactory behavior. Proc Natl Acad Sci USA 110, 1017–1022. https://doi.org/10.1073/pnas.1220168110 (2013).ADS 
Article 
PubMed 

Google Scholar 
43.Vijayan, V., Thistle, R., Liu, T., Starostina, E. & Pikielny, C. W. Drosophila pheromone-sensing neurons expressing the ppk25 Ion channel subunit stimulate male courtship and female receptivity. PLoS Genet. 10, e1004238. https://doi.org/10.1371/journal.pgen.1004238 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
44.Thistle, R., Cameron, P., Ghorayshi, A., Dennison, L. & Scott, K. Contact chemoreceptors mediate male-male repulsion and male-female attraction during Drosophila courtship. Cell 149, 1140–1151. https://doi.org/10.1016/j.cell.2012.03.045 (2012).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
45.Khan, M. A. M. et al. Raspberry ketone accelerates sexual maturation and improves mating performance of sterile male Queensland fruit fly, Bactrocera tryoni (Froggatt). Pest Manag. Sci. 75, 1942–1950. https://doi.org/10.1002/ps.5307 (2019).CAS 
Article 
PubMed 

Google Scholar 
46.Weldon, C. W., Perez-Staples, D. & Taylor, P. W. Feeding on yeast hydrolysate enhances attraction to cue-lure in Queensland fruit flies, Bactrocera tryoni. Entomol. Experim. et Applicata 129, 200–209. https://doi.org/10.1111/j.1570-7458.2008.00768.x (2008)47.Najar-Rodriguez, A. J., Galizia, C. G., Stierle, J. & Dorn, S. Behavioral and neurophysiological responses of an insect to changing ratios of constituents in host plant-derived volatile mixtures. J. Exp. Biol. 213, 3388 (2010).CAS 
Article 

Google Scholar 
48.Bertschy, C., Turlings, T. C., Bellotti, A. C. & Dorn, S. Chemically-mediated attraction of three parasitoid species to mealybug-infested cassava leaves. Florida Entomologist 80, 383–395 (1997).Article 

Google Scholar 
49.Butler, D. G., Cullis, B. R., Gilmour, A. R. & Gogel, B. J. (Queensland Department of Primary Industries and Fisheries, Australia, 2009).50.Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. https://doi.org/10.1093/bioinformatics/btu170 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
51.Zhao, Q.-Y. et al. Optimizing de novo transcriptome assembly from short-read RNA-Seq data: a comparative study. BMC Bioinf. 12, S2. https://doi.org/10.1186/1471-2105-12-S14-S2 (2011).CAS 
Article 

Google Scholar 
52.Fu, L., Niu, B., Zhu, Z., Wu, S. & Li, W. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28, 3150–3152. https://doi.org/10.1093/bioinformatics/bts565 (2012).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
53.Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210–3212. https://doi.org/10.1093/bioinformatics/btv351 (2015).CAS 
Article 
PubMed 

Google Scholar 
54.Conesa, A. et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. https://doi.org/10.1093/bioinformatics/bti610 (2005).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
55.Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760. https://doi.org/10.1093/bioinformatics/btp324 (2009).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
56.Patro, R., Duggal, G., Love, M. I., Irizarry, R. A. & Kingsford, C. Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 14, 417–419. https://doi.org/10.1038/nmeth.4197 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
57.Risso, D., Ngai, J., Speed, T. P. & Dudoit, S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat. Biotechnol. 32, 896–902. https://doi.org/10.1038/nbt.2931 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
58.Al-Shahrour, F., Díaz-Uriarte, R. & Dopazo, J. FatiGO: a web tool for finding significant associations of Gene Ontology terms with groups of genes. Bioinformatics 20, 578–580. https://doi.org/10.1093/bioinformatics/btg455 (2004).CAS 
Article 
PubMed 

Google Scholar 
59.Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucl. Acids Res. 43, e47–e47. https://doi.org/10.1093/nar/gkv007 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar  More

Data collectionData on the birds’ reproductive success and the number of humans present at nest sites were collected as part of a long-term, ongoing monitoring project in Hungary, in which we investigate the impacts of urbanization on populations of great tits. The great tit is an insectivorous passerine bird that is widespread across the Western Palearctic, occupies both urban and forest habitats, readily accepts nestboxes, and shares many important ecological traits with other tit or chickadee species also occurring in urban habitats27. These traits make this species an ideal model organism for studying the effects of the anthropause on wildlife in different environments.Study sitesWe monitored breeding great tit populations and also collected human presence data in two urban areas and at one forest study site. In one of the urban sites, Veszprém (47°05′17.29″N, 17°54′29.66″E; human population: c. 56,000; the monitoring scheme started in 2013), the nestboxes were placed in public green spaces (public parks, university campuses, a bus station, and a cemetery) that are surrounded by built-up areas and roads, and experience frequent anthropogenic disturbance. At the other urban size, Budapest (47°30′27.4″N, 19°01′03.4″E; the capital city of Hungary, human population: c. 1.75 million; the monitoring scheme started in 2019), the nestboxes were placed in two public urban parks, located c. 400 m from each other in the city core area and separated by high-traffic roads. The parks are freely accessible to residents and are heavily embedded within the urban matrix. At both urban sites, most of the nestboxes are distributed along paths or walking trails. Even though the two cities greatly differ in their size and human population, our urban study plots in both cities have similar general characteristics: these are surrounded by built-up areas, are at a similar distance (c. 3–4 km) from the nearest forested areas (for Veszprém, this is the forest at Vilma-puszta: 47°05′06.7″N, 17°51′51.4″E; for Budapest, this is the forest at Normafa: 47°30′27.7″N 18°57′51.1″E), and nests also experienced a similar level of human disturbance in the pre-COVID reference period (Fig. 1b). The forest site, Szentgál (47°06′39.75″N, 17°41′17.94″E; the monitoring scheme started in 2013), is a mature woodland, dominated by beech (Fagus sylvatica) and hornbeam (Carpinus betulus), located 3 km away from the nearest human habitation (Szentgál, human population: c. 2.800), c. 20 km and 110 km away from Veszprém and Budapest, respectively. There are no paved roads in the forest, and the area is relatively free from human disturbance although it experiences occasional hunting and logging activity.Human presence around nestsTo quantify human presence at our study sites for 2020 and the reference years we counted the number of humans (motorized vehicles excluded) during each nest check, for 30 s, in the proximity of the nestboxes (for similar approach see Corsini et al. 2019). The number of humans was recorded within a 50-m radius of the nestboxes between 2013 and 2018 (Veszprém, Szentgál), and within a 15-m radius distance in 2019–2020 (all sites). We changed the counting distance in 2019 due to methodological reasons following28. However, to be able to compare the human presence data of 2020 in Veszprém and Szentgál to that recorded in earlier years, in 2020 we performed the counts with both the 15-m and the 50-m radius distances at these two sites. Thus, for 2020 in Veszprém, we have human presence data both for the 50-m and the 15-m radius areas that were used in the forest-city and the between-cities comparisons, respectively (see below). For each year and study site, we used human presence data only from seasonally first broods (defined below), and only from nests where there were either already eggs or nestlings in the nest, resulting in 9.4 ± 3.6 (mean ± SD) observations per brood which is a reliable indicator of human presence28.Birds’ reproductive successWe monitored nestboxes each year at least twice a week from mid-March to early June to record laying and hatching dates, clutch size, hatching success, and the number of nestlings in active great tit nests. We ringed nestlings at day 14–16 post-hatch (i.e. a few days before fledging; hatching day of the first chick = day 1) with a numbered metal ring and also recorded their body mass (to the nearest 0.1 g), tarsus length (to the nearest 0.1 mm and following Svensson’s ‘alternative’ method29) and wing length (from the bend of the wing to the longest primary; to nearest 1 mm). Shortly after the expected date of fledging we carefully examined the nest material to identify and count the number of chicks that died after ringing (due to e.g. starvation, predation) that we included in the calculation of nestling survival (detailed below). The aim of this is to get a more accurate estimate for the number of offspring that could indeed fledge from the nest. The number of broods (nestlings) that suffered partial or complete mortality between ringing and fledging were: n = 6 (13) in Budapest (2019–2020), n = 70 (152) in Veszprém (2013–2020), and n = 25 (83) in the Szentgál forest.From these data we determined clutch size (the maximum number of eggs observed in a brood), hatching success (the proportion of chicks hatched / eggs laid), the number of fledglings, and nestling survival (the proportion of fledged young / hatched chicks). The number of fledglings (i.e. the number of young fledged successfully) was calculated as the number of chicks ringed minus the number of chicks found dead in the nest after the ringing. We involved only seasonally first breeding attempts (as this period overlapped with the lockdown period; detailed at the Statistical analyses), and defined first broods as follows. In our study system breeding great tits are captured on their nests and receive a unique combination of colour rings. Active nests are also routinely equipped with a small, concealed video camera enabling us to reliably identify over 80% of breeding individuals each year30. Thus, relying on this setup, we considered a clutch as a first breeding attempt of a pair if it was initiated before the date of the first egg laid in the earliest second clutch at that site by an individually identifiable (i.e. colour-ringed) female that successfully raised her first clutch (i.e. fledged at least one young) in that year.Air pollution and meteorological conditionsTo describe the levels of traffic-related air pollution (nitrogen dioxide [NO2], nitrogen oxides [NOX] and ozone [O3]) and the meteorological conditions (temperature and precipitation) at the two urban study sites (Veszprém and Budapest), we used data provided by the Hungarian Air Quality Monitoring Network and the Hungarian Meteorological Service, respectively. To better understand which aspect of the anthropause might have affected great tits’ breeding success we thus assessed if the lockdown affected air pollution levels differently at the two urban study sites (compared to 2019), or if weather conditions showed different fluctuations between 2019 and 2020 at the two cities. For more details on the statistical analyses and results, see ESM: Sect. 1.Statistical analysesThe duration of the official restrictions on human mobility (lockdown) spanned between 28 March–4 May in Veszprém (calendar date: 88–125; 01 January = 1) and 28 March–18 May (88–139) in Budapest. During this period people were allowed to leave their homes e.g. to run essential errands including individual sport and recreational activities in public green spaces, although with keeping at least 1.5 m from each other (social distancing). Very importantly, from the point of view of our study, the period of movement restrictions had almost completely overlapped with the seasonally first breeding attempts (from egg-laying to fledging) of great tits at both urban sites. The date of laying the 1st egg (calendar date, mean ± SD) in Veszprém was 94.2 ± 6.4, while in Budapest 97 ± 7.8; the date of chick ringing and measuring in Veszprém was 128 ± 5.3, while in Budapest: 133 ± 9.1. Thus, we decided not to exclude any first broods based on the date in order to maximize our sample size. Similarly, the period from which we involved human presence data was also strongly overlapped with the duration of the movement restrictions in both cities. Therefore, in Veszprém, the calendar dates of the first and the last human count at each nest were 87–108 (median: 100) and 121–142 (median: 132), respectively, while in Budapest 87–128 (median: 98) and 118–155 (median: 128).Human presence around nestsIn accordance with our first objective (forest-city comparisons), we explored if the lockdown in 2020 caused any changes in human disturbance around the great tit nests. To do so, we compared the number of humans (50-m radius of the nests) between 2020 and the 2013–2018 reference period, separately for the forest (Szentgál) and urban (Veszprém) study sites. Note that in 2019, we did not collect data on human presence within a 50-m radius at Veszprém and Szentgál (see above: Data collection), therefore 2019 was not included in the reference period of this analysis. We, however, also compared human presence in Veszprém between 2019 and 2020 using the 15-m radius data which indicates a change that is consistent with the differences found using the 50-m radius data (detailed below).First, we built generalized linear mixed-effects (GLM, lme4 R package) models with Poisson error distribution with the number of humans as the response variable, including year as a fixed factor and nestbox ID as random factor to control for non-independence of the data. Next, we extracted the mean values (least-squares means; package emmeans31) and associated standard errors for each year as estimated by the model. We computed the mean of these yearly mean estimates for the 2013–2018 reference period (i.e. calculated a single overall mean describing the whole reference period) and compared this long-term mean to the mean estimate of 2020 by calculating the linear contrast between them (with the ‘contrast’ function of the emmeans package), and expressed linear contrasts as 2020 minus the reference period.For our second objective (between-cities comparisons), we compared the changes in human disturbance around the nestboxes at the two urban study sites, Veszprém and Budapest, using the number of humans recorded within the 15-m radius of the active nests in 2019 and 2020. We analysed the data from Budapest and Veszprém separately and built generalized linear mixed-effects models with Poisson error distribution with the number of humans (15-m radius of the nests) as the response variable, including year as a fixed factor and nestbox ID as random factor to control for non-independence.Birds’ reproductive successWe used data from 2019 (reference; for justification see below in this section) and 2020 (lockdown). First, we constructed separate linear models to analyse each component of reproductive success (response variables), and for the forest-city and the between-cities comparisons. We used linear models (LM) for clutch size and the number of fledglings, generalized linear models (GLM, with quasi-binomial error distribution) for hatching success and nestling survival, and linear-mixed effects models (LME) for nestling body size traits (body mass, tarsus length, and wing length). Models on nestling body size traits contained nestlings’ age at ringing as a confounding variable (three-level factor: 14, 15, or 16 d of age) and brood ID as a random factor to control for the non-independence of chicks raised in the same brood. Finally, these models always contained a habitat (Veszprém or Szentgál) × year (2019 or 2020) interaction term for forest-city comparisons and a city (Budapest or Veszprém) × year (2019 or 2020) interaction term for between-cities comparisons. We checked assumptions of residuals’ normality and homogeneity of variance by inspecting the residuals plots which were respected for all models.Next, to test the prediction for our first objective (forest-city comparisons), we extracted the mean values (least-squares means) and associated standard errors of each response variable for each habitat × year combination as estimated by the linear model’s interaction. Then, from these estimates, we calculated habitat contrasts, i.e. the mean forest-city difference (forest minus urban) for each year (i.e. for 2019 and 2020), and compared the mean habitat contrast for the 2019 reference year to the mean habitat contrast of 2020; for similar approach see14,32,33.For our second objective (between-cities comparisons), we followed the same procedure as for the forest-city comparisons (detailed above) except that here we compared the differences between cities (Budapest minus Veszprém) in 2020 and 2019. These full models (i.e. for the forest-city and between-cities comparisons) are presented in Table S1–S2 (ESM: Sect. 2).In our study, we chose 2019 as a reference year for multiple reasons. First, because this was temporally the closest year without a lockdown. Second, because for Budapest we have monitoring data only from 2019 to 2020, using 2019 and 2020 in all analyses makes the results more comparable. Finally, although we have monitoring data from a total of eight years (2013–2020) for Szentgál (forest site) and Veszprém (urban site), for the forest-city comparisons we did not include years before 2019 in the reference because we noticed a negative trend in birds’ reproductive success throughout the study years (Fig. S3). This trend was especially apparent in the forest population, and may have reduced the forest-city difference by the end of the study period. Indeed, 2019 and 2020 were amongst the poorest years and resulted in a very similar reproductive success between both years within both habitats (Fig. S3). Because such temporal trend may have confounded the comparisons of 2020 with earlier years, to take account for its effect, and to further justify our approach of using 2019 as the reference year, we conducted additional analyses on the birds’ reproductive success by comparing both 2019 and 2020 (separately) to the 2013–2018 long-term reference period. We predict that if 2019 and 2020 are similarly affected by the decreasing trend in reproductive success than then the differences between the long-term reference period and 2019 and 2020, respectively, should be similar. For the details of these long-term forest-city comparisons see ESM: Sect. 3 and Table S3).Finally, we did not conduct the forest-city comparisons (first objective) between the forest site (Szentgál) and the other urban site (Budapest) for two reasons. First, because unlike to the Szentgál vs. Veszprém setup, we did not have an appropriate forest (control) location which is close to Budapest. Second, because conducting comparisons between the long-term data and 2019 and 2020, respectively (see: ESM Sect. 3) was not possible for Budapest because we do not have similar long-term data for the latter site.Clutches that failed before reaching the incubation stage (due to predation or desertion; i.e. final clutch size was uncertain), suffered complete mortality due to weather (e.g. nestbox fall from the tree due to strong wind), and cases when complete or partial clutch or brood loss may have occurred due to the monitoring process (e.g. when a nestbox was dropped or when complete brood failure occurred soon after capturing a parent on the nest) were excluded from all analyses. In the analyses investigating the number of fledglings, fledging success, and nestling body size traits we involved nests only where at least one nestling hatched, and excluded broods that were involved in a food-supplementation experiment (as treatment group) during the nestling rearing period in 201714. We used the R 4.0.5 software environment for statistical analysis and creating figures34.Ethical statementAll procedures were in accordance with Hungarian laws, and adhered to the ASAB/ABS guidelines for the use of animals in behavioural research and teaching. Permit to the use of animals in this study was provided by the National Scientific Ethical Committee on Animal Experimentation (permit number: PE-06/KTF/997–8/2018, FPH061/1329–5/2018, PE-06/KTF/06,543–7/2020 and FPH061/3036–4/2020). Permits to study protected species and access to protected areas were provided by the Middle Transdanubian Inspectorate for Environmental Protection, Natural Protection and Water Management (permit numbers: 31559/2011, 24,861/2014 and VE-09Z/03,454–8/2018, for working in Veszprém and Szentgál) and the Environment Protection and Nature Conservation Department of the Pest County Bureau of the Hungarian Government and the Mayor’s Office of Budapest (permit numbers: PE-06/KTF/997–8/2018, FPH061/1329–5/2018, PE-06/KTF/06,543–7/2020 and FPH061/3036–4/2020, for working in Budapest). More

ScopeWe compared disposable face masks that were used once with face masks that were sterilized and used five more times (six times in total). Sterilisation and PFE test data of the Aura 1862+ (3M, Saint Paul, Minnesota, USA) face mask indicate that this type of face mask shows good performance after multiple sterilisation cycles10,11,12. In a previous pilot study, the company CSA Services (Utrecht, the Netherlands), a sterilization facility for cleaning, disinfection and sterilization of medical instruments, was rebuild to process FFP2 face masks. In total, 18,166 single use FFP2 masks were sterilised after use in a medical autoclave. As the majority (n = 7993) were Aura 1862+ (3M, Saint Paul, Minnesota, USA), this particular type of face mask was chosen for the LCA.The total weight of the face masks and packaging together during end-of-life consists of incineration for the face masks (97%) and landfill for the carton box packaging of new face masks (3%). There is no recycling potential used in our model since the materials coming from the operating room and its packaging is commonly disposed as medical waste. In the Netherlands, no energy recovery takes place at the incineration of regulated medical waste. Therefore, no co-function was applicable for the end-of-life scenario.Recycling is often a multi-functional process that produces two or more goods. To deal with the multi-functionality in the background processes, the cut-off approach was applied to exclude the allocation of the greenhouse gas emissions to additional goods. This means that potential rest materials such as energy gained during incineration are cut-off and that the greenhouse gas emissions are fully allocated to the waste treatment processes itself.In the LCA, the ‘functional unit’ defines the primary function that is fulfilled by the investigational products and indicates how much of this function is considered18. In this study, we pragmatically chose as a definition for the protection of 100 health care workers against airborne viruses, using one FFP2 certified face mask, each during one working shift of an average of 2 h in a hospital in the Netherlands.Table 1 shows the differences between the two scenarios:

1.

100 masks including packaging, transported from production to the hospital, used and disposed.

2.

100 times use of reprocessed masks. We calculated that 27.1 masks are being produced and transported from production to the hospital. The 27.1 are being reprocessed five times, taking into account that 20% of the batch cannot be reprocessed. Therefore 80% of the batch could be used for reprocessing after each step resulting in: 27.1 (new) + 21.7 (repro 1) + 17.3 (repro 2) + 13.9 (repro 3) + 11.1 (repro 4) + 8.9 (repro 5) = 100 times of use. For each time of reprocessing the batch is transported from the hospital to the (hospital) Central Sterilization Services Department (CSSD) and disposed after five times of reprocessing.

Table 1 Comparison between reference flow 1 and 2.Full size tableCombining the functional unit with the two alternative scenarios results in the reference flows for the protection of 100 health care workers against airborne viruses, either using a face mask one single time (100 virgin masks produced for the 1st scenario), or reusing a face mask for five additional times (27.1 virgin masks produced for the 2nd scenario). For both reference flows, only FFP2 certified face masks are considered. For the calculations each mask is used for a single two hours working shift in an average hospital in the Netherlands.Life cycle inventory (LCI) analysisThe inventory data includes all phases from production (including material production and part production), transport, sterilisation to end-of-life of the life cycle of the single use and reprocessed face masks. We disassembled one face mask to obtain the weight of each individual component on a precision scale (Fit Evolve, Bangosa Digital, Groningen, the Netherlands) with a calibrated inaccuracy of 1.5%. Component information and materials were obtained from the data fact sheet provided by the manufacturer. We conducted a separate validation experiment to establish the material composition in the filtering fabric (Supplement file).This LCA with the Aura 3M masks was based on steam sterilization by means of a hospital autoclave and therefore part of this study. Therefore, face masks were placed in a sterilization bag that contained up to five masks. A total of 1000 masks were placed into an autoclave (Getinge, GSS6713H-E, Sweden) per cycle. After sterilization, the masks were transported to the hospital. Masks were reprocessed for a maximum of five times before final disposal10,11.The assessment of climate change impact is done following as closely as possible the internationally accepted Life Cycle Assessment (LCA) method following the ISO 14040 and 14044 standards19,20. The LCA examines all the phases of the product’s life cycle from raw material extraction to production, packaging, transport, use and reprocessing until final disposal19. The LCA was modelled using SimaPro 9.1.0.7 (PRé Sustainability, Amersfoort, The Netherlands). The background life cycle inventory data were retrieved from the ecoinvent database (Ecoinvent version 3.6, Zürich, Switzerland)21.To make a valid comparison between the disposable and reprocessing face masks, the system boundaries should be equal in both scenarios. The system boundaries in this study consisted of the production, the use and the disposal and waste treatment of the masks. For the reprocessed face masks, the lifecycle is extended due to the sterilisation process (Fig. 1). Therefore, the additional PPE’s and materials needed to safely process the masks (e.q. masks, gloves and protective sheets) are included in the production phase. The production of machinery for the manufacturing of the face masks and the autoclave were not included in this study.Figure 1System boundary overview of new and reprocessed face masks including waste treatment by incineration.Full size imageThe production facility for the face masks is located in Shanghai, China22,23. Further distribution took place from Bracknell, UK to Neuss, Germany and the final destination was set in Rotterdam, the Netherlands.The packaging materials were disposed in the hospital where the face masks are used primarily. After first use, face masks were transported to the sterilisation department. All masks were manually checked before reprocessing by personnel wearing PPE. Of all used Aura 1862+ facemasks that entered the CSA, approximately 10% was discarded. To remain conservative, the LCA was conducted based on a 20% rejection rate as a result of face masks which could not be reused anymore due to deformities, lipstick, and broken elastic bands.A full overview of the life cycle inventory table for the two scenarios and details on model assumptions are added in the Supplemental file (Supplemental file, Part B).Life cycle impact assessmentThe carbon footprint (kg CO2 eq) was chosen as the primary unit in the impact category. ReCiPe was applied at midpoint level and used to translate greenhouse gas emissions into climate change impact16.Uncertainty analysisThe final LCA model contains several uncertainties based on assumptions and measurement inaccuracies24. The included uncertainties were based on weighted components of the masks as well as the packaging which were measured with 1.5% inaccuracy of the precision scale apparatus. A Monte Carlo sampling25 was conducted for both alternatives (disposable and reprocessing) where input parameters for the LCA were sampled randomly from their respective statistical distributions in for 10,000 ‘runs’. Because input parameters between scenarios were partly overlapping, we compared these two scenarios directly using a discernibility analysis. This technique, establishes which scenario is beneficial for each of 10,000 Monte Carlo runs. We report the percentage of instances where the reprocessing scenario has a lower carbon footprint than the disposable scenario.Sensitivity analysisA sensitivity analysis was conducted to check the sensitivity of the outcome measures to variation in the input parameters. To determine which parameters are interesting to investigate, three aspects were considered: the variations in number of face masks per sterilization cycle (autoclave capacity), rejection rate (number of losses per cycle) and transport distance to the CSSD. Finally, we included the relative contribution of these variations. The following three parameter variations were chosen for the sensitivity analysis:

1.

Rejection percentage. The rejection rate was defined based on experiences from the participating sterilisation department and studies that show that sterilisation of the face masks up to 5 times is possible. Masks were re-used for 5 times, approximately 10% was discarded during the total life cycle. Out of this experience and to remain conservative, the total rejection rate was set on 20%. Therefore it is interesting to investigate whether variation in PFE testing outcomes or differences in user protocols influence the outcomes. This should indicate if masks from higher or lower quality can also be suitable candidates for reprocessing.

2.

Autoclave capacity, which largely depends on the loading of the autoclave. To mimic different loads of the autoclave, it is interesting to know the influence of sterilizing fewer masks per run on the model.

3.

Transport. As it is likely that many hospitals have a Central Sterilisation Services Department (CSSD) it is interesting to know the effect of having zero transportation. Moreover, in case hospitals are not willing to change the routing in their CSSD it is interesting to observe how outcomes are influenced if transportation is set on the maximal realistic value of 200 km.

The parameters have been varied with 250 and 500 face masks per sterilisation batch. A rate varying with 10% and 30% of the face masks being rejected due to quality reasons and variation in transport kilometres of 0–200 km.There is a small difference between the baselines of the sensitivity, LCIA and contribution analyses because all these are performed using separate Monte-Carlo simulations. The output of the different simulations may show minor differences due to statistical distribution.Cost price comparisonA cost analysis was made to give insight in costing from a procurement perspective. The cost analysis is conducted with five face masks that were steam sterilized per batch in a permeable laminate bag, Halyard type CLFP150X300WI-S20 and includes the expenses of energy, depreciation, water consumption, cost of personnel, overhead and compared to the prices for a new disposable 3M Aura face mask during the first and second Corona waves. Five pieces per bag were chosen in order to have enough space between the masks to sterilise each mask properly. The cost analysis is based on actual sterilization as well as associated costs compared to the prices of new disposable face masks. The costs were then related to the functional unit of protecting 100 health care workers by calculating the difference in the amount of Euros per 100 face masks. More

Isolation of marine sediment-derived actinomycetes from west-central PhilippinesThe geographical sites identified in this study were evaluated to explore the actinomycete distribution in west-central Philippines (Fig. 1a). A total of 16 sediment cores were recovered from the 11 sampling sites and were processed in the laboratory using cultured-dependent actinomycete isolation (Supplementary Table S2). The seawater’s physicochemical conditions in all sampling sites were uniform with pH 7.0 and salinity ranging from 3.1 to 3.2. The characteristics of marine sediments and the distance of the actual collection sites identified using the given criteria varied per sampling location. The sediment characteristics vary from coarse to fine sand with mixture of broken corals and pebbles (Supplementary Table S3). Culture-dependent isolation revealed that actinomycete-like colonies and spores were observed in the minimal marine media after 30 to 60 days of incubation. Actinomycete isolates were repeatedly sub-cultured in enriched marine medium 1 (MM1) to obtain pure cultures as shown in Fig. 1b. Actinomycete growth observed in enriched media was white to gray aerial spores with brown to yellow mycelia or without diffusible pigmentations. Notably, there were strains that produced colonies with no diffusible pigmentations.Figure 1Distribution, abundance, and antibacterial activity of marine sediment-derived actinomycetes in the Philippines. (a) The overall map showing the 11 sampling sites situated within the west-central region in the Philippines. The enlarged map showed the details of the sampling sites and their corresponding actinomycetes abundance. Sampling sites are number-coded as shown inside the circle. The color gradient indicates strain abundance. (b) Actinomycete colonies were sub-cultured repeatedly to obtain pure culture of isolates. (c) A total of 92 out of 2212 actinomycetes strains have confirmed antibacterial activities as verified using microbroth susceptibility assay. The map with Streptomyces abundance plot was generated using ggplot2 package in Rstudio ver. 1.2.5042 (https://www.rstudio.com/).Full size imageIn this work, a total of 2212 pure actinomycete strains were isolated from marine sediments collected in 11 geographically distant sampling sites across the west-central Philippines (Fig. 1a). Actinomycete strains were highly abundant in Negros Occidental with 580 isolates (26%), followed by Southern Antique with 348 isolates (16%) and Southeastern Iloilo with 228 (10%). We recovered least actinomycete strains in Occidental Mindoro and South Central Visayas with only 94 and 81 strains, respectively.Antibacterial activity profile of actinomycete strainsWe assessed the antibacterial activities of actinomycete strains against a multidrug-resistant Gram-positive bacterium (S. aureus ATCC BAA-44) and three Gram-negative bacteria (E. coli ATCC 25922, P. aeruginosa ATCC 27853, and E. aerogenes ATCC 13048) using resazurin agar overlay assay and microbroth susceptibility assay as initial and confirmatory screenings, respectively. A total of 218 (9.9%) out of the 2212 actinomycete isolates have antibacterial activities in the initial screening as indicated by positive results or retained blue resazurin color in wells containing actinomycete overlaid with the test pathogens (Supplementary Fig. S2). The 218 active isolates were fermented to produce biomass for secondary metabolite extraction and to confirm their antibacterial activities by microbroth susceptibility testing. Figure 1c showed the confirmed 92 (4.1%) antibiotic-producing actinomycete strains. The majority of the strains (71 isolates) exhibited activity against Gram-positive S. aureus ATCC BAA-44. Thirty-nine (39) strains (42%) were active against E. coli ATCC 25922. Six strains were active against P. aeruginosa ATCC 27853, while all strains tested were inactive against E. aerogenes ATCC 13048 as indicated with no or less than 50% growth inhibition. Twenty-three (23) active strains targeted 2–3 test pathogens, while 69 active strains were only active against one test pathogen (Supplementary Table S4).Phylogenetic diversity of multiple antibiotic-producing strainsThe 92 active actinomycete strains were further identified and confirmed as Streptomyces species based on genomic analysis of their 16S rRNA and rpoB gene sequences. Comparison of 16S rRNA gene sequences (ranging from 1150 to 1500 nucleotides) and rpoB (700–995 nucleotides) gene sequences with their similar matches in the GenBank verified that the 92 active strains were closely related (97 to 100%) with 19 species under the genus Streptomyces. The nearly complete 16S rRNA and rpoB gene sequences were analyzed in a phylogenetic tree using maximum likelihood algorithms. The 16S rRNA and rpoB gene sequences of active strains reported in the present study were deposited in the GenBank nucleotide database (Supplementary Table S5).Phylogenetic analysis revealed multiple strains with identical 16S rRNA gene sequences which clustered together into 13 major clusters (shown by the colored nodes in the tree) with high bootstrap values ( > 90%) in the phylogenetic tree (Supplementary Fig. S3). Thirty-three strains (36%) were highly similar to S. parvulus presented in red circle. Followed by 12 strains with high similarity to S. enissocaesilis (light blue circle), 11 S. rochei strains (dark pink), six S. mutabilis strains (dark blue), five S. diastaticus strains (light green), four S. kunmingensis strains (green) and three S. geysiriensis strains (light orange). A phylogenetic analysis of rpoB gene sequence was conducted to provide a better resolution of the evolutionary relationship among strains within and between species supporting the taxonomic identity of the phylogenetically identical strains.Similarly, the rpoB gene sequences phylogenetic tree (Fig. 2) showed 13 major clusters that were highly supported with bootstrap replicates  > 90%, except for monophyletic clusters III (Streptomyces sp. strain DSD176) and X (Streptomyces sp. strain DSD1006) with low bootstrap replicates ( 1% abundance) can be recovered in both methods. We also noted that some species were exclusively recovered using a specific method. Eight species were exclusively recovered in DSM compared to 6 species in HSM. All active strains of S. kunmingensis, S. mutabilis, S. sedi and S. olivaceus were only recovered by HSM. Contrary, the active strains of S. carpaticus and S. harbinensis were only isolated using DSM.Carbon source composition of marine minimal media, along with effective isolation techniques, was crucial for the isolation of antibiotic-producing Streptomyces. Among the five minimal marine media used, three carbon sources yielded high isolation rate: glucose, mannitol, and trehalose yielded nine species with 32 strains (35%), nine species with 18 strains (20%), and ten species with 18 strains (20%) respectively. However, only eight species (12 strains, 20%) and two species (12 strains, 20%) were isolated in raffinose, and starch-based media, respectively. In the contrary, high diversity indices were observed in trehalose (2.197), mannitol (2.0), and raffinose (1.979) (Fig. 5a). As expected, the starch-containing media had the lowest isolation rate and diversity as only two species (S. enissocaesilis and S. parvulus) were able to utilize a more complex carbon source.Figure 5Diversity of antibiotic-producing Streptomyces using five different carbon sources. (a) From the five carbon sources in the minimal marine media utilized by Streptomyces strains in this study, mannitol yielded the highest number of active strains, while high diversity was recorded in active strains that utilized glucose (n = 92). (b) Venn diagram of five carbon sources showed that two Streptomyces species can be isolated using all five carbon sources.Full size imageInterestingly, co-isolation of species in different carbon sources was shown in the Venn diagram (Fig. 5b). Eight species can be recovered from at least two different carbon sources, whereas 11 species were exclusively isolated from a specific carbon source. Bioactive S. enissocaesilis and S. parvulus strains were recovered from all of the carbon sources utilized in this study. Active S. rochei were isolated in four media but not in starch-based media. In contrast, more exclusive species were isolated in trehalose with four species (S. harbinensis strains, Streptomyces sp. strain DSD3025, S. pseudogriseolus, and S. xiamenensis). Followed by glucose with three species (Streptomyces sp. strain DSD742, S. carpaticus and S. sedi), and two species each on mannitol- and raffinose-based media (Fig. 5b). The results indicated that diverse Streptomyces species preferred simple sugars-containing one or two sugar molecule as nutrient source compared to complex sugars.
Streptomyces abundance and diversity in geographical sampling locationsBioactive Streptomyces species were widely distributed across the different sampling locations in west-central Philippines (Fig. 6a). Although Southern Antique, Negros Occidental and Negros Oriental have highest number of active strains isolated, we found that Southern Antique, Southern Iloilo, and Western Antique were the most diverse sampling sites (Fig. 6b). We have isolated the greatest number of antibiotic-producing Streptomyces species which were evenly distributed in Southern Antique. This indicates that Southern Antique is stable with many potential niches that can support highly diverse Streptomyces species.Figure 6(a) The abundance profile of 19 antibiotic-producing Streptomyces species in different geographical location across the Philippines revealed that S. parvulus was the most abundant species. (b) Diversity, richness and evenness of antibiotic-producing Streptomyces species in different sites were calculated using vegan package in RStudio ver1.2.5042 (https://www.rstudio.com/). (c) Recovery profile and abundance of antibiotic-producing strain per depth layer, indicated by colored circles and its size, revealed that strains in the bottom sediments were the most diverse as compared to other depth layer and S. parvulus was the most abundant species recovered.Full size imageStreptomyces parvulus emerged as the dominant antibiotic-producing species in this study. Out of the 33 active S. parvulus strains, the majority were isolated in Negros Occidental (11 strains) and Negros Oriental (ten strains). Although Negros Occidental and Negros Oriental have high abundance, its microbial community is highly dominated by one species, S. parvulus, supported by the low species richness, evenness, and diversity (Fig. 6b). This finding implies that these sites may have few potential niches that only a few species dominate. Bioactive S. enissocaesilis strains were recovered in four sampling locations only; specifically, Southern Antique (five strains), Western Antique (three strains), Southeastern Iloilo (three strains), and Northwestern Antique (one strain). Active S. rochei were isolated and evenly distributed in seven sampling locations, but were not present in Western Antique, Northwestern Antique, South Central Visayas, and Tubbataha Reefs. Notably, we observed that no bioactive S. parvulus, S. enissocaesilis, and S. rochei were isolated in Tubbataha Reefs, but antibiotic-producing S. cacaoi, S. psuedogriseolus and S. mutabilis strains were isolated only in Tubbataha Reefs marine sediments. Meanwhile, site-specific species such as S. sedi were recovered only in Occidental Mindoro. The isolation of site-specific species within genus Streptomyces can offer insight on the adaptive capacity of strains to inhibit locally coexisting resource competitors within and among these distinct locations.Distribution of bioactive Streptomyces species at different sediment depthsWe further investigated the distribution of antibiotic-producing Streptomyces strains along with the 110-cm sediment depth in different sampling sites. From the sediments that were partitioned according to depth with five categories at 25-cm increments, heterogeneous distributions of bioactive species were observed at deeper sediment with different dominant species in each depth (Fig. 6c). Although, S. parvulus, S. rochei, and S. enissocaesilis strains were ubiquitous in all depths, several species thrive abundantly in specific depths compared to other species. Streptomyces parvulus was the most dominant species in surface sediment. Meanwhile, Streptomyces rochei and S. enissocaesilis strains were more adapted in surface and sub-bottom sediments, respectively. Depth-specific Streptomyces strains were also identified as follows: Streptomyces sp. strain DSD3025 was isolated in subsurface sediments; S. mutabilis strains were abundant in the middle sediment layer; Streptomyces sp. strain DSD1006 and S. pseudogriseolus strain were recovered from sub-bottom sediments; and Streptomyces strain sp. DSD742, S. albus, S. sedi and S. xiamenensis strains were obtained from bottom sediments. High species diversity was positively correlated with increasing sediment depth, where surface sediments are known to be more prone to dispersal and wash-offs by environmental factors such as deep ocean currents42. Furthermore, the depth-specific species identified largely influenced the species richness in varying sediment depth. More

Study regions and farmsTen European regions from boreal to Mediterranean were selected (Supplementary Table 1). They represented major agricultural land uses such as arable crops including horticulture, mixed farming, grassland and perennial crops (vineyards and olives). Within each region, a pool of ~20–40 farms was selected from which 12–20 farms were randomly selected (169 in total) that belonged to the same farm type, produced under homogeneous climatic and environmental circumstances and fulfilled specific criteria regarding their main production branch. In case the selected farms were not willing to participate, we asked other farms from the pool till the sufficient number has been reached. The selected organic farms had all been certified for at least five years. Farmers were asked if they were willing to participate in the study. If they refused, additional random sampling was conducted. In the region NL, 11 organic farms agreed to participate but only three non-organic farms, whereas seven organic farms and 11 non-organic farms were available in the region HU. During the study, one non-organic farmer in the region CH ceased participation.Habitat maps and farm interviewsThe complete area of all selected farms was mapped, using the BioHab method36. Excluded from the farm area were woody and aquatic habitats larger than 800 m2 and summer pastures. Within the farm area, areal and linear habitats were recorded. For an areal habitat, the minimal mapping unit was 400 m2 with a width of at least 5 m. More narrow habitats, between 0.5 and 5 m wide and at least 30 m long, were mapped as linear habitats. Habitats were distinguished in habitat types according to Raunkiær life forms, environmental conditions and management evidence28. Further, a farmland class was assigned to each habitat that described whether the habitat was managed for agricultural production or other objectives such as e.g. nature conservation. In face-to-face interviews following a standardized questionnaire, farmers provided detailed information on field management and yield.Categorization as production fields and semi-natural habitatsBased on the habitat maps and available information about management intensity, we categorized all habitats as either semi-natural habitats or production fields. In agricultural landscapes, these two categories are often not clearly distinguishable. There is a gradient from more intensively managed production fields to less intensively used semi-natural habitats. In addition, a categorization at the local scale can be different from an approach at a European scale (29 and see p. 45 of37). Here, we applied the same criteria for all ten study regions.In all cases, we categorized as production fields: arable crops, intensively managed grasslands (following main plant species observed, management evidence and objectives, with fertilization and/or two or more cuts a year), horticultural crops, and vineyards.We categorized as semi-natural habitats: linear habitats, habitats that were managed for nature conservation objectives, habitats where mainly geophytes, helophytes or hydrophytes were growing, grasslands with woody vegetation (shrubs and/or trees), and extensively managed grasslands (no fertilization, no or one cut a year).Species samplingVascular plant, earthworm, spider and bee species were sampled in all different habitat types of a farm. One plot per habitat type was randomly selected per farm for species sampling. This resulted in 1402 selected habitat plots on 169 farms (Supplementary Table 2). In the selected habitats, species were sampled during one growing season, using standardized protocols19,38. Plant species were identified in squares of 10 × 10 m in areal habitats and in rectangular strips of 1 × 10 m in linear habitats. Earthworms were collected at three random locations of 30 × 30 cm per habitat. First, a solution of allyl isothiocyanate (AITC) was poured out to extract earthworms from the soil. Afterwards, a 20-cm-deep soil core from the same location was hand sorted to find additional specimens. Identification took place in the lab. Spiders were sampled on three dates at five random locations per habitat within a circle of 0.1 m2. Using a modified vacuum shredder, spiders were taken from the soil surface, transferred to a cool box, frozen, or put in ethanol, sorted and identified in the lab. Bees (wild bees and bumble bees) were sampled on three dates, during dry, sunny and warm weather conditions. They were captured with an entomological aerial net along a 100 m long and 2 m wide transect, transferred to a killing jar and identified in the lab.Grouping of species dataSpecies data were pooled per taxa, habitat and region, and three sub-communities were formed with species (1) exclusively found in semi-natural habitats, i.e. unique to semi-natural habitats, (2) exclusively found in production fields, i.e. unique to production fields, and (3) found in both habitat categories i.e. shared by production fields and semi-natural habitats. For calculations of effects over all four taxa, species richness was the sum of the individual taxa species richnesses.Estimating species richnessSpecies richness was estimated using coverage- and sample-size-based rarefaction and extrapolation curves31,39,40. Rarefaction and extrapolation, including confidence intervals (bootstrap method) and sampling coverage, were calculated in R 3.4.041 using package iNEXT42. Detailed information is provided below for each topic.Estimating richness of unique species to compare semi-natural habitats and production fieldsTo legitimately compare the richness of species unique to semi-natural habitats and to production fields, we used the coverage-based method, i.e. we standardized the samples by their completeness30. The point of comparison was determined by the so-called ‘base coverage’ identified by the following procedure31: (1) select the maximum sample coverage at reference sample size (number of sampling units) of the sub-communities under comparison, (2) select the minimum sample coverage at twice the reference sample size of the sub-communities under comparison, (3) identify the maximum of the results from step (1) and step (2) as ‘base coverage’. The species richness estimates were then read off from the species sample-size-based rarefaction and extrapolation curves at the ‘base coverage’ for each sub-community being compared. If zero or exactly one species was unique to a sub-community at the reference sample size, no sample coverage could be calculated. In this case, we set the species richness at 0 or 1, respectively. The species richness estimate of the other sub-community under comparison was then read off at twice the reference sample size on the curve.The ‘base coverage’ was individually defined for each region and each taxonomic group since the mixed effects models used to analyze the data took into account the variation among regions and taxonomic groups.Differences in species richness unique to semi-natural habitats and production fieldsThe difference between the species richness unique to semi-natural habitats and unique to production fields was tested with mixed effects models using package lme4 (Version 1.1-12) in R43. The data were (Sij | β, b, x) ~ Poisson(µij) from i = 1, …, 10 regions. The model is:$${{{rm{ln}}}}left({mu }_{{ij}}right)={beta }_{0}+{beta }_{1}{x}_{1i}+{b}_{1i}$$
(1)
$${b}_{1} sim N(0,sigma 2)$$where ({beta }_{0}) is a fixed intercept, ({beta }_{1}) a fixed effect sub-community ({x}_{1{ij}}) (species unique to semi-natural habitats versus species unique to production fields), b1i are random intercepts for region i. Random effects are normally distributed with mean 0 and variance σ2. The significance of term ({beta }_{1}) was calculated by log-likelihood ratio tests with one degree of freedom. For the models over all four taxa, an additional random intercept was included, i.e. b2j with mean 0 and variance σ2 for j = 1, …, 4 taxa (Fig. 1b).Differences in species richness between organic and non-organic systemsThe comparison between organic and non-organic systems of species unique to semi-natural habitats and to production fields, and of species shared by the two habitat categories, relied on coverage-based extrapolation as described above. Differences between management systems were tested for significance using mixed-effects models with management system ({beta }_{1}) ({x}_{1{ij}}) as fixed effect in (1).Estimating species loss due to conversion of semi-natural habitats to production fieldsTo predict the species loss due to conversion of semi-natural habitats to production fields, we relied on sample-size-based extrapolations31 with species incidence frequencies. We estimated the richness of the species pool for the total number of mapped habitats including the extrapolated species richness unique to semi-natural habitats and unique to production fields, and the observed richness of shared species for each of the four taxa. This species pool provided the basis for the calculation of the species loss or gain (Table 1 and Supplementary Table 7). To model the species richness decrease for any amount of semi-natural habitats converted to production fields, we calculated and drew backward the curve composed of the accumulation curve for species unique to semi-natural habitats, to which the estimated total species richness unique to production fields (constant) and the corresponding gain of species unique to production fields (increases with increasing area of production fields as semi-natural habitats are converted), and the richness of observed shared species (constant) were added. This is the species decrease curve (Supplementary Fig. 2). If started at the observed species richness, this curve corresponds exactly to a species richness curve calculated by a cumulative random removal of semi-natural habitats one by one from the pool of all habitats. The four taxa decrease curves were added for the curve in Fig. 2. Confidence intervals (CI, 95%) shown in Figs. 2 and 3 are calculated by bootstrapping within the calculation of the species accumulation curves (iNEXT42), upper and lower bounds of the 95% CI of the four taxa being added. From the species decrease curve, we read off the predicted species richness for a conversion of 50% and 90% of the semi-natural habitats, and a conversion required to increase production by 10%.As species were sampled in 20% of all mapped habitats on average per region (min. 8%, max. 35%), extrapolated species accumulation curves used to build the species decrease curve were calculated for more than two to three times the reference sample size, which is the suggested range for reliable extrapolation of the species richness estimator31,44. Obviously, the confidence intervals (CI) of the species richness extrapolations here became wide (Supplementary Fig. 4). As we still wanted to show the impact of a conversion of the whole semi-natural area into production fields on the production gain in the ten regions, we used the uncertainty (upper and lower bounds of the 95% CI of the four taxa added) to define two situations in addition to the average case to predict species richness for a 50% and a 90% semi-natural habitat conversion, and a conversion required to increase production by 10%: (1) a worst case situation with the upper bound of the CI of the expected species richness unique to semi-natural habitats, the lower bound of the CI of the expected species richness unique to production fields, and shared species assumed not to be able to survive without semi-natural habitats and considered like species unique to semi-natural habitats (i.e. upper bound); and (2) a best case situation with the lower bound of the CI of the expected species richness unique to semi-natural habitats, the upper bound of the CI of the expected species richness unique to production fields, and the lower bound of the CI of the expected shared species richness.Estimating production gainFarmer interviews delivered an average yield per crop type per farm for the years 2008–2010 (Supplementary Data45 shows details for organic and non-organic systems separately). Farmers indicated yield in kilograms or tons per hectare. This was transformed into energy units, i.e. mega joules per hectare (MJ ha−1) using standard values46. From this, for each region, the average yield (MJ ha−1) was calculated by first multiplying individual crop type yields by the corresponding crop type areas to obtain the production per crop type, then summing up the production of all crop types, and finally dividing this sum by the total area of the crop types. For livestock farms, the fodder production of grasslands was estimated based on the average requirements per livestock unit, accounting for the amount of feed grain, legumes, silage maize and of imported feedstuff. All yields relate to plant biomass production and do not comprise livestock products. The average yield takes into account the relative cover of the different crop types in the regions. Therefore, the conversion of the semi-natural area to production fields was region-specific. The production of certain semi-natural habitats as e.g. olive groves in Spain was not part of the production calculation. The reason is that data on production for semi-natural habitats were mainly not available and/or negligible, e.g. extensively used grassland in CH or in HU, and we decided to apply the same treatment to all the regions. Consequently, in case of olive groves in Spain the effective increase in production is overestimated. To calculate the production gain per region, the production field area added by the conversion of semi-natural habitat area was multiplied by the average yield. In practice, in many regions it may be impossible to convert semi-natural habitat to productive land due to geomorphological constraints and poor soils, and even if land were converted, yields would be much lower than these averages. The results presented here, especially the 90% scenario, are therefore over-optimistic. On the other hand, our calculations are based on the area of semi-natural habitat available for conversion on existing farms, but in some regions other sources of semi-natural land may be available for conversion, e.g. former agricultural land that has been abandoned.Species loss and production gain for three scenariosWe calculated the change of species richness and the production gain under current day production efficiency for two scenarios: (1) a conversion of 90% of the semi-natural area into production fields. The 10% of semi-natural area remaining is considered unsuitable for agricultural use or even impossible to cultivate; (2) a conversion of 50% of the semi-natural area into production fields, and (3) a necessary conversion of the semi-natural area into production fields to achieve a 10% production increase per region.Standardization for organic and non-organic systemsAlthough the overall mapped area, the number of semi-natural habitats, the number of production fields and the average habitat size did not significantly differ between the two management systems (Supplementary Table 5), we standardized the number and size of habitats to the average across both systems per region to compare the species loss and production gain at current day production efficiency in the organic and non-organic systems. The total production in organic and non-organic systems per region was calculated based on the respective yield and the average mapped area of the production fields across both systems as described in section “Estimation of production gain”. The impact on biodiversity was analyzed for the scenario that organic systems should achieve the same level of production as non-organic systems by converting semi-natural habitats to production fields. We calculated the amount of the required area to be converted into production fields and the corresponding species change.Differences between management systems were again tested for significance using mixed-effects models with management system ({{{{rm{beta }}}}}_{1}) ({{{{rm{x}}}}}_{1{{{rm{ij}}}}}) as fixed effect in (1).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Africa’s population is rapidly growing, with its share of the global population projected to increase from 17% in 2020 to 39% by 2100 (ref. 8). The continent is already grappling with low agricultural productivity and food security challenges. Tremendous efforts are needed to increase food production; however, arable land continues to undergo widespread degradation due to issues such as nutrient mining, erosion, overgrazing and pollution. Climate change and more frequent weather extremes, such as floods and droughts, further degrade land and reduce agricultural productivity.Some efforts to counteract low productivity, however, can increase greenhouse gas emissions and derail efforts to meet global climate targets. Poor water management, fertilizer application and residue burning in rice production are, for example, major sources of potent greenhouse gases such as methane and nitrous oxide9,10. To ensure that the United Nations sustainable development goals and the African Union’s Agenda 2063 for food and water security are realized at minimal environmental cost, science-based land management practices are needed to decouple agricultural productivity from greenhouse gas emissions.
Credit: majimazuri21/PixabayThe Agriculture, Forestry and Other Land Uses (AFOLU) sector contributes the largest share of greenhouse gas emissions in Africa11. Thus, developing large-scale agronomic, livestock and forest management practices that increase productivity and reduce emissions is key to achieving enhanced production and environmental sustainability. However, it is impossible to effectively manage greenhouse gas emissions if there is limited capacity to quantify them in Africa.Improved data infrastructure and research are needed to quantify emissions associated with specific land management practices under different land uses. Similarly, land use mitigation strategies should be informed by existing and potential future land use changes and their impact on greenhouse gas emissions under different climate scenarios. However, past studies that examined land use changes at various temporal scales mainly used coarse resolution satellite imagery and suffered from limited availability or poor-quality of data, partly due to cost. Such challenges have resulted in limited knowledge of land management practices that reduce greenhouse gas emissions while increasing agricultural productivity.Improved greenhouse gas observation networks and in situ measurements12 will enable the development of country-specific emission factors (IPCC tier 2/3)13 and quantification and management of land use specific greenhouse emissions. It will reduce uncertainties in emissions inventory data on Agriculture, Forestry and Other Land Uses14, which are currently estimated using emission factors extracted from default value databases (tier 1 methodologies).Free earth observation data, such as those from the European Space Agency and United States Geological Surveys, are becoming increasingly available. Together with improvements in cloud-based computing infrastructure, this presents an opportunity to advance research into current and future land use and vegetation dynamics. Coupled with accurately quantified greenhouse gas emissions, this can support current and future land management practices that contribute to mitigation and adaptation objectives of countries. More

1.von Humboldt, A. Ansichten der Natur mit Wissenschaftlichen Erlauterungen (J.G. Cotta, 1808).2.Perrigo, A., Hoorn, C. & Antonelli, A. Why mountains matter for biodiversity. J. Biogeogr. 47, 315–325 (2020).Article 

Google Scholar 
3.Badgley, C. et al. Biodiversity and topographic complexity: modern and geohistorical perspectives. Trends Ecol. Evol. 32, 211–226 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
4.Rahbek, C. et al. Building mountain biodiversity: geological and evolutionary processes. Science 365, 1114–1119 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Steinbauer, M. J. et al. Topography-driven isolation, speciation and a global increase of endemism with elevation. Glob. Ecol. Biogeogr. 25, 1097–1107 (2016).Article 

Google Scholar 
6.Fjeldså, J., Bowie, R. C. K. & Rahbek, C. The role of mountain ranges in the diversification of birds. Annu. Rev. Ecol. Evol. Syst. 43, 249–265 (2012).Article 

Google Scholar 
7.Hughes, C. & Eastwood, R. Island radiation on a continental scale: exceptional rates of plant diversification after uplift of the Andes. Proc. Natl Acad. Sci. USA 103, 10334–10339 (2006).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Antonelli, A. et al. Geological and climatic influences on mountain biodiversity. Nat. Geosci. 11, 718–725 (2018).CAS 
Article 

Google Scholar 
9.Quintero, I. & Jetz, W. Global elevational diversity and diversification of birds. Nature 555, 246–250 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Gillooly, J. F., Allen, A. P., West, G. B. & Brown, J. H. The rate of DNA evolution: effects of body size and temperature on the molecular clock. Proc. Natl Acad. Sci. USA 102, 140–145 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Martin, A. P. & Palumbi, S. R. Body size, metabolic rate, generation time, and the molecular clock. Proc. Natl Acad. Sci. USA 90, 4087–4091 (1993).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
12.Rohde, K. Latitudinal gradients in species diversity: the search for the primary cause. Oikos 65, 514–527 (1992).Article 

Google Scholar 
13.Allen, A. P., Gillooly, J. F., Savage, V. M. & Brown, J. H. Kinetic effects of temperature on rates of genetic divergence and speciation. Proc. Natl Acad. Sci. USA 103, 9130–9135 (2006).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
14.Rabosky, D. L. et al. An inverse latitudinal gradient in speciation rate for marine fishes. Nature 559, 392–395 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
15.Igea, J. & Tanentzap, A. J. Angiosperm speciation cools down in the tropics. Ecol. Lett. 23, 692–700 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
16.Schluter, D. Speciation, ecological opportunity, and latitude (American Society of Naturalists address). Am. Nat. 187, 1–18 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
17.Anderson, K. J. & Jetz, W. The broad-scale ecology of energy expenditure of endotherms. Ecol. Lett. 8, 310–318 (2005).Article 

Google Scholar 
18.Clarke, A. & Gaston, K. J. Climate, energy and diversity. Proc. R. Soc. B 273, 2257–2266 (2006).PubMed 
PubMed Central 
Article 

Google Scholar 
19.Dowle, E. J., Morgan-Richards, M. & Trewick, S. A. Molecular evolution and the latitudinal biodiversity gradient. Heredity 110, 501–510 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
20.Brown, J. H. Why are there so many species in the tropics? J. Biogeogr. 41, 8–22 (2014).PubMed 
Article 

Google Scholar 
21.Stevens, G. C. The latitudinal gradient in geographical range: how so many species coexist in the tropics. Am. Nat. 133, 240–256 (1989).Article 

Google Scholar 
22.Boucher-Lalonde, V. & Currie, D. J. Spatial autocorrelation can generate stronger correlations between range size and climatic niches than the biological signal — a demonstration using bird and mammal range maps. PLoS One 11, e0166243 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
23.Cutter, A. D. & Gray, J. C. Ephemeral ecological speciation and the latitudinal biodiversity gradient. Evolution 70, 2171–2185 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
24.Morales‐Barbero, J., Martinez, P. A., Ferrer‐Castán, D. & Olalla‐Tárraga, M. Á. Quaternary refugia are associated with higher speciation rates in mammalian faunas of the Western Palaearctic. Ecography 41, 607–621 (2018).Article 

Google Scholar 
25.Xing, Y. & Ree, R. H. Uplift-driven diversification in the Hengduan Mountains, a temperate biodiversity hotspot. Proc. Natl Acad. Sci. USA 114, E3444–E3451 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
26.Lagomarsino, L. P., Condamine, F. L., Antonelli, A., Mulch, A. & Davis, C. C. The abiotic and biotic drivers of rapid diversification in Andean bellflowers (Campanulaceae). New Phytol. 210, 1430–1442 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
27.Testo, W. L., Sessa, E. & Barrington, D. S. The rise of the Andes promoted rapid diversification in Neotropical Phlegmariurus (Lycopodiaceae). New Phytol. 222, 604–613 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
28.Dowsett, H. et al. The PRISM4 (mid-Piacenzian) paleoenvironmental reconstruction. Climate 12, 1519–1538 (2016).
Google Scholar 
29.Hartley, A. J. Andean uplift and climate change. J. Geol. Soc. 160, 7–10 (2003).Article 

Google Scholar 
30.Aron, P. G. & Poulsen, C. J. in Mountains, Climate and Biodiversity (eds Hoorn, C., Perrugi, A. & Antonelli, A.) Ch. 8 (2018).31.Hewitt, G. The genetic legacy of the Quaternary ice ages. Nature 405, 907–913 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
32.Wallis, G. P., Waters, J. M., Upton, P. & Craw, D. Transverse Alpine speciation driven by glaciation. Trends Ecol. Evol. 31, 916–926 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
33.Luebert, F. & Muller, L. A. H. Effects of mountain formation and uplift on biological diversity. Front. Genet. 6, 54 (2015).34.Huang, S., Meijers, M. J. M., Eyres, A., Mulch, A. & Fritz, S. A. Unravelling the history of biodiversity in mountain ranges through integrating geology and biogeography. J. Biogeogr. 46, 1777–1791 (2019).Article 

Google Scholar 
35.Whittaker, R. J., Triantis, K. A. & Ladle, R. J. A general dynamic theory of oceanic island biogeography. J. Biogeogr. 35, 977–994 (2008).Article 

Google Scholar 
36.Li, Y. et al. Climate and topography explain range sizes of terrestrial vertebrates. Nat. Clim. Change 6, 498–502 (2016).Article 

Google Scholar 
37.Kisel, Y. & Barraclough, T. G. Speciation has a spatial scale that depends on levels of gene flow. Am. Nat. 175, 316–334 (2010).PubMed 
Article 

Google Scholar 
38.Spooner, F. E. B., Pearson, R. G. & Freeman, R. Rapid warming is associated with population decline among terrestrial birds and mammals globally. Glob. Change Biol. 24, 4521–4531 (2018).Article 

Google Scholar 
39.Rowley, D. B. & Garzione, C. N. Stable isotope-based paleoaltimetry. Annu. Rev. Earth Planet. Sci. 35, 463–508 (2007).CAS 
Article 

Google Scholar 
40.Mulch, A. Stable isotope paleoaltimetry and the evolution of landscapes and life. Earth Planet. Sci. Lett. 433, 180–191 (2016).CAS 
Article 

Google Scholar 
41.Kuhn, T. S., Mooers, A. Ø. & Thomas, G. H. A simple polytomy resolver for dated phylogenies. Methods Ecol. Evol. 2, 427–436 (2011).Article 

Google Scholar 
42.Rolland, J., Condamine, F. L., Jiguet, F. & Morlon, H. Faster speciation and reduced extinction in the tropics contribute to the mammalian latitudinal diversity gradient. PLoS Biol. 12, e1001775 (2014).43.Meredith, R. W. et al. Impacts of the Cretaceous Terrestrial Revolution and KPg Extinction on mammal diversification. Science 334, 521–524 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
44.Britton, T., Anderson, C. L., Jacquet, D., Lundqvist, S. & Bremer, K. Estimating divergence times in large phylogenetic trees. Syst. Biol. 56, 741–752 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
45.Drummond, A. J. & Rambaut, A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7, 214 (2007).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
46.Schliep, K. P. phangorn: phylogenetic analysis in R. Bioinformatics 27, 592–593 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Jetz, W., Thomas, G. H., Joy, J. B., Hartmann, K. & Mooers, A. O. The global diversity of birds in space and time. Nature 491, 444–448 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
48.Redding, D. W. & Mooers, A. Ø. Incorporating evolutionary measures into conservation prioritization. Conserv. Biol. 20, 1670–1678 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
49.Rabosky, D. L. Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees. PLoS One 9, e89543 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
50.Moore, B. R., Höhna, S., May, M. R., Rannala, B. & Huelsenbeck, J. P. Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures. Proc. Natl Acad. Sci. USA 113, 9569–9574 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
51.Meyer, A. L. S., Román-Palacios, C. & Wiens, J. J. BAMM gives misleading rate estimates in simulated and empirical datasets. Evolution 72, 2257–2266 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
52.Rabosky, D. L., Mitchell, J. S. & Chang, J. Is BAMM flawed? Theoretical and practical concerns in the analysis of multi-rate diversification models. Syst. Biol. 66, 477–498 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
53.Mitchell, J. S., Etienne, R. S. & Rabosky, D. L. Inferring diversification rate variation from phylogenies with fossils. Syst. Biol. 68, 1–18 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
54.Title, P. O. & Rabosky, D. L. Tip rates, phylogenies and diversification: what are we estimating, and how good are the estimates? Methods Ecol. Evol. 10, 821–834 (2019).Article 

Google Scholar 
55.Louca, S. & Pennell, M. W. Extant timetrees are consistent with a myriad of diversification histories. Nature 580, 502–505 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
56.Amante, C. & Eakins, B. W. ETOPO1 Arc-minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24 (NOAA, 2009).57.Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).Article 

Google Scholar 
58.Brown, J. L., Hill, D. J., Dolan, A. M., Carnaval, A. C. & Haywood, A. M. PaleoClim, high spatial resolution paleoclimate surfaces for global land areas. Sci. Data 5, 180254 (2018).PubMed 
PubMed Central 
Article 

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

Google Scholar 
60.Bivand, R. & Piras, G. Comparing implementations of estimation methods for spatial econometrics. J. Stat. Softw. 63, v063i18 (2015).
Google Scholar  More

1.Foley, J. A. et al. Global consequences of land use. Science 309, 570–574 (2005).CAS 

Google Scholar 
2.Carpenter, S. R. et al. Science for managing ecosystem services: beyond the Millennium Ecosystem Assessment. Proc. Natl Acad. Sci. USA 106, 1305–1312 (2009).CAS 

Google Scholar 
3.Nelson, E. et al. Modeling multiple ecosystem services, biodiversity conservation, commodity production, and tradeoffs at landscape scales. Front. Ecol. Environ. 7, 4–11 (2009).
Google Scholar 
4.Maes, J. et al. An indicator framework for assessing ecosystem services in support of the EU Biodiversity Strategy to 2020. Ecosyst. Serv. 17, 14–23 (2016).
Google Scholar 
5.Maes, J. et al. Mapping ecosystem services for policy support and decision making in the European Union. Ecosyst. Serv. 1, 31–39 (2012).
Google Scholar 
6.Millennium Ecosystem Assessment. Ecosystems and Human Well-Being: Synthesis (Island Press, 2005).7.Summary for Policymakers. In Global Assessment Report on Biodiversity and Ecosystem Services (IPBES, 2019).8.Martínez-Harms, M. J. & Balvanera, P. Methods for mapping ecosystem service supply: a review. Int. J. Biodivers. Sci. Ecosyst. Serv. Manage. 8, 17–25 (2012).
Google Scholar 
9.Hauck, J. et al. ‘Maps have an air of authority’: potential benefits and challenges of ecosystem service maps at different levels of decision making. Ecosyst. Serv. 4, 25–32 (2013).
Google Scholar 
10.Balvanera, P. et al. Conserving biodiversity and ecosystem services. Science 291, 2047 (2001).CAS 

Google Scholar 
11.Dick, J., Maes, J., Smith, R. I., Paracchini, M. L. & Zulian, G. Cross-scale analysis of ecosystem services identified and assessed at local and European level. Ecol. Indic. 38, 20–30 (2014).
Google Scholar 
12.UK National Ecosystem Assessment. The UK National Ecosystem Assessment Technical Report. (UNEP-WCMC, 2011); http://uknea.unep-wcmc.org/13.Orsi, F., Ciolli, M., Primmer, E., Varumo, L. & Geneletti, D. Mapping hotspots and bundles of forest ecosystem services across the European Union. Land Use Policy 99, 104840 (2020).
Google Scholar 
14.Holland, R. A. et al. The influence of temporal variation on relationships between ecosystem services. Biodivers. Conserv. 20, 3285–3294 (2011).
Google Scholar 
15.Renard, D., Rhemtull, J. M. & Bennett, E. M. Historical dynamics in ecosystem service bundles. Proc. Natl Acad. Sci. USA 112, 13411–13416 (2015).CAS 

Google Scholar 
16.Rukundo, E. et al. Spatio-temporal dynamics of critical ecosystem services in response to agricultural expansion in Rwanda, East Africa. Ecol. Indic. 89, 696–705 (2018).
Google Scholar 
17.Stürck, J., Schulp, C. J. E. & Verburg, P. H. Spatio-temporal dynamics of regulating ecosystem services in Europe—the role of past and future land use change. Appl. Geogr. 63, 121–135 (2015).
Google Scholar 
18.Rau, A. L. et al. Temporal patterns in ecosystem services research: a review and three recommendations. Ambio 49, 1377–1393 (2020).
Google Scholar 
19.Sutherland, I. J., Bennett, E. M. & Gergel, S. E. Recovery trends for multiple ecosystem services reveal non-linear responses and long-term tradeoffs from temperate forest harvesting. For. Ecol. Manage. 374, 61–70 (2016).
Google Scholar 
20.Hansen, M. C., Stehman, S. V. & Potapov, P. V. Quantification of global gross forest cover loss. Proc. Natl Acad. Sci. USA 107, 8650–8655 (2010).CAS 

Google Scholar 
21.Vanhanen, H. et al. Making Boreal Forests Work for People and Nature (IUFRO, 2012).22.Pan, Y. et al. A large and persistent carbon sink in the world’s forests. Science 333, 988–993 (2011).CAS 

Google Scholar 
23.Moen, J. et al. Eye on the Taiga: removing global policy impediments to safeguard the boreal forest. Conserv. Lett. 7, 408–418 (2014).
Google Scholar 
24.Global Forest Industry (Swedish Forest Industries, 2019); https://www.forestindustries.se/forest-industry/statistics/global-forest-industry/25.Saastamoinen, O., Kangas, K. & Aho, H. The picking of wild berries in Finland in 1997 and 1998. Scand. J. For. Res. 15, 645–650 (2000).
Google Scholar 
26.Gamfeldt, L. et al. Higher levels of multiple ecosystem services are found in forests with more tree species. Nat. Commun. 4, 1340 (2013).
Google Scholar 
27.Hou, Y., Li, B., Müller, F., Fu, Q. & Chen, W. A conservation decision-making framework based on ecosystem service hotspot and interaction analyses on multiple scales. Sci. Total Environ. 643, 277–291 (2018).CAS 

Google Scholar 
28.Blumstein, M. & Thompson, J. R. Land-use impacts on the quantity and configuration of ecosystem service provisioning in Massachusetts, USA. J. Appl. Ecol. 52, 1009–1019 (2015).
Google Scholar 
29.Fernandez-Campo, M., Rodríguez-Morales, B., Dramstad, W. E., Fjellstad, W. & Diaz-Varela, E. R. Ecosystem services mapping for detection of bundles, synergies and trade-offs: examples from two Norwegian municipalities. Ecosyst. Serv. 28, 283–297 (2017).
Google Scholar 
30.Gissi, E., Fraschetti, S. & Micheli, F. Incorporating change in marine spatial planning: a review. Environ. Sci. Policy 92, 191–200 (2019).
Google Scholar 
31.Maxwell, S. M., Gjerde, K. M., Conners, M. G. & Crowder, L. B. Mobile protected areas for biodiversity on the high seas. Science 367, 252–254 (2020).CAS 

Google Scholar 
32.Willcock, S. et al. Do ecosystem service maps and models meet stakeholders’ needs? A preliminary survey across sub-Saharan Africa. Ecosyst. Serv. 18, 110–117 (2016).
Google Scholar 
33.Jonsson, M., Bengtsson, J., Gamfeldt, L., Moen, J. & Snäll, T. Levels of forest ecosystem services depend on specific mixtures of commercial tree species. Nat. Plants 5, 141–147 (2019).
Google Scholar 
34.Pohjanmies, T. et al. Impacts of forestry on boreal forests: an ecosystem services perspective. Ambio 46, 743–755 (2017).
Google Scholar 
35.Miina, J., Hotanen, J.-P. & Salo, K. Modelling the abundance and temporal variation in the production of bilberry (Vaccinium myrtillus L.) in Finnish mineral soil forests. Silva Fenn. 43, 577–593 (2009).
Google Scholar 
36.Hertel, A. G. et al. Berry production drives bottom–up effects on body mass and reproductive success in an omnivore. Oikos 127, 197–207 (2018).
Google Scholar 
37.Thiffault, E. Boreal forests and soils. Dev. Soil Sci. 36, 59–82 (2019).
Google Scholar 
38.Jonsson, M., Bengtsson, J., Moen, J., Gamfeldt, L. & Snäll, T. Stand age and climate influence forest ecosystem service delivery and multifunctionality. Environ. Res. Lett. 15, 0940a8 (2020).
Google Scholar 
39.Stokland, J. N. Volume increment and carbon dynamics in boreal forest when extending the rotation length towards biologically old stands. For. Ecol. Manage. 488, 119017 (2021).
Google Scholar 
40.Harmon, M. E., Ferrell, W. K. & Franklin, J. F. Effects on carbon storage of conversion of old-growth forests to young forests. Science 247, 699–702 (1990).CAS 

Google Scholar 
41.Mazziotta, A. et al. Applying a framework for landscape planning under climate change for the conservation of biodiversity in the Finnish boreal forest. Glob. Change Biol. 21, 637–651 (2015).
Google Scholar 
42.Triviño, M. et al. Optimizing management to enhance multifunctionality in a boreal forest landscape. J. Appl. Ecol. 54, 61–70 (2017).
Google Scholar 
43.Qiu, J. & Turner, M. G. Spatial interactions among ecosystem services in an urbanizing agricultural watershed. Proc. Natl Acad. Sci. USA 110, 12149–12154 (2013).CAS 

Google Scholar 
44.Felipe-Lucia, M. R. et al. Multiple forest attributes underpin the supply of multiple ecosystem services. Nat. Commun. 9, 4839 (2018).
Google Scholar 
45.Eggers, J., Räty, M., Öhman, K. & Snäll, T. How well do stakeholder-defined forest management scenarios balance economic and ecological forest values? Forests 11, 86 (2020).
Google Scholar 
46.Eyvindson, K., Repo, A. & Mönkkönen, M. Mitigating forest biodiversity and ecosystem service losses in the era of bio-based economy. For. Policy Econ. 92, 119–127 (2018).
Google Scholar 
47.Rusch, A., Bommarco, R., Jonsson, M., Smith, H. G. & Ekbom, B. Flow and stability of natural pest control services depend on complexity and crop rotation at the landscape scale. J. Appl. Ecol. 50, 345–354 (2013).
Google Scholar 
48.Schipanski, M. E. et al. A framework for evaluating ecosystem services provided by cover crops in agroecosystems. Agric. Syst. 125, 12–22 (2014).
Google Scholar 
49.Hufnagel, J., Reckling, M. & Ewert, F. Diverse approaches to crop diversification in agricultural research. A review. Agron. Sustain. Dev. 40, 14 (2020).
Google Scholar 
50.Guerry, A. D. et al. Modeling benefits from nature: using ecosystem services to inform coastal and marine spatial planning. Int. J. Biodivers. Sci. Ecosyst. Serv. Manage. 8, 107–121 (2012).
Google Scholar 
51.Wikström, P. et al. The Heureka Forestry Decision Support System: An Overview. Math. Comput. For Nat.-Resour. Sci. 3, 87–94 (2011).
Google Scholar 
52.Forest Statistics (Swedish University of Agricultural Sciences, 2020).53.Eriksson, A., Snäll, T. & Harrison, P. J. Analys av miljöförhållanden ‐ SKA 15. Report 11 (Swedish Forest Agency, 2015).54.Axelsson, A.-L. et al. in National Forest Inventories—Pathways for Common Reporting (eds Tomppo, E. et al.) 541–553 (Springer, 2010).55.Marklund, L. G. Biomass Functions for Pine, Spruce and Birch in Sweden (1988).56.Petersson, H. & Ståhl, G. Functions for below-ground biomass of Pinus sylvestris, Picea abies, Betula pendula and Betula pubescens in Sweden. Scand. J. For. Res. 21, 24–83 (2006).
Google Scholar 
57.Miina, J., Pukkala, T. & Kurttila, M. Optimal multi-product management of stands producing timber and wild berries. Eur. J. For. Res. 135, 781–794 (2016).
Google Scholar 
58.Schröter, M. & Remme, R. P. Spatial prioritisation for conserving ecosystem services: comparing hotspots with heuristic optimisation. Landsc. Ecol. 31, 431–450 (2016).
Google Scholar 
59.Wu, J., Feng, Z., Gao, Y. & Peng, J. Hotspot and relationship identification in multiple landscape services: a case study on an area with intensive human activities. Ecol. Indic. 29, 529–537 (2013).CAS 

Google Scholar 
60.Akaike, H. A new look at the statistical model identification. IEEE Trans. Automat. Contr. 19, 716–723 (1974).
Google Scholar 
61.R: A Language and Environment for Statistical Computing (R Development Core Team, 2014); https://www.R-project.org/ More