Ecology

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

Sample collection of Liolophura japonica and the genetic diversity of COI barcoding regionTo examine genetic lineage divergence within L. japonica on the northwestern Pacific coast, we newly collected a total of 342 L. japonica samples from 12 sampling localities in the intertidal coasts of the Korean Peninsula and Japanese Archipelago (Fig. 1; Table S1). From the collected L. japonica samples, we amplified the COI barcoding region using PCR, and then sequenced the 635-bp PCR products. As a result, a total of 75 COI haplotypes based on COI sequences obtained from 342 individuals of L. japonica were detected via the present study (Table S2). In addition to this, we extracted 31 COI haplotypes based on COI sequences from 127 individuals of L. japonica (also known as Acanthopleura japonica) previously reported in the NCBI GenBank database, consisting of two Japanese and 29 Chinese COI haplotypes. Finally, we gathered 106 COI haplotypes from 469 L. japonica individuals collected in 15 localities of South Korea, Japan, and southern China (Tables S1, S2). The average haplotype (h) and nucleotide diversities (π) were 0.808 and 0.04936, respectively; the highest haplotype diversity was observed in Tsushima (TS; h = 0.963), and the highest nucleotide diversity was found in Wando (WD; π = 0.04581), located in the South Sea of the Korean Peninsula. As shown in Table S3, the population distribution pattern of COI haplotypes revealed that all collection sites had site-specific haplotype(s) except for Busan (BS), Wando (WD), Sinan (SA), and Jeju Island (JJ). The most abundant haplotype was A1, which was found in 170 (39.4%) out of the COI sequences obtained from 469 L. japonica individuals.Figure 1A map showing sampling localities and photos of a habitat landscape and wild samples of Liolophura japonica inhabiting coastal areas of the Korean Peninsula (N = 249), the Japanese Archipelago (N = 57), and southern China (N = 125) in the northwestern Pacific Ocean. (a) A map showing twelve direct sampling localities for L. japonica in coastal areas of the northwestern Pacific Ocean. The sampling localities of one southern Chinese (ZJ) and two Japanese (EH and MY) previously catalogued haplotype sequencing studies retrieved from NCBI are also depicted. Table S1 and S2 contain more accurate information on the populations and individuals. The basic map is from a free map providing site (https://d-maps.com), which is modified with Adobe Illustrator v.25.2. (https://www.adobe.com). (b) Photos of a habitat landscape and wild samples of L. japonica, taken from Seogwipo-si, Jeju Island, South Korea, photographed by Mi Young Yeo, Bia Park, and Cho Rong Shin. The photos were edited using Adobe Photoshop v.22.2 (https://www.adobe.com).Full size imagePhylogenetic and population genetic analyses based on COI
We constructed a nucleotide sequence alignment set with 106 COI haplotypes of L. japonica (Data S1), and identified 95 polymorphic sites (15.0%, Table S4) and 68 parsimoniously informative sites (10.7%). To elucidate phylogenetic relationships among the populations of L. japonica, we performed molecular phylogenetic analyses, including maximum likelihood (ML), Bayesian inference (BI), and neighbor-joining (NJ) analyses, based on these 106 COI haplotypes with the outgroup Acanthopleura spinosa (Fig. 2a, Figs. S1, S2). The resultant phylogenetic trees clearly revealed the existence of three distinct genetic lineages within the monophyletic group of L. japonica (100 BP in ML, 1.00 BPP in BI, and 100 BP in NJ): Lineage N (91 BP, 1.00 BPP, and 100 BP), Lineage S1 (79 BP, 0.82 BPP, and 98 BP), and Lineage S2 (98 BP, 1.00 BPP, and 100 BP). Among these three genetic lineages, Lineages S1 and S2 were grouped with high node confidence values (94 BP, 1.00 BPP, and 95 BP). We additionally conducted a phylogenetic network analysis using a neighbor net algorithm without an outgroup (Fig. 2b), which confirmed that these sequences were distinctly divided into three genetic lineages, in agreement with the topology of the rooted phylogenetic trees (Fig. 2a, Fig. S1).Figure 2Phylogenetic, TCS network, and PCoA analyses based on 106 COI haplotypes from 469 individuals of Liolophura japonica inhabiting coastal areas of the northwestern Pacific Ocean, suggesting the existence of the three different genetic lineages: Lineage N, Lineage S1, and Lineage S2. (a) Maximum likelihood tree showing the three different genetic lineages for L. japonica: Lineage N members are most likely from the populations inhabiting a wide range of South Korea and Japan, Lineage S1 members from the populations inhabiting southern coastal areas of South Korea and Japan only, and Lineage S2 members from the southern Chinese population. As shown in Fig. S1, Acanthopleura spinosa was used as an outgroup. Numbers on branches indicate node confidence values: BP in ML, BPP in BI, and BP in NJ in order. (b) A phylogenetic network reconstructed using the neighbor net algorithm without an outgroup, showing three different genetic lineages for L. japonica inhabiting the northwestern Pacific coast: Lineages N, S1, and S2. The COI sequence alignment set used is shown in Data S1. Detailed information of the 106 COI haplotypes used in this phylogenetic analysis is summarized in Table S1 and S2. (c) An unrooted TCS network showing three distinct genetic clusters, corresponding to Lineages N, S1, and S2. Three different genetic groups correspond to the three genetic lineages shown in the phylogenetic tree (a), respectively. The haplotype frequency is displayed by the circle size. (d) A two-dimensional PCoA plot showing the three distinct genetic groups corresponding to Lineages N, S1, and S2 shown in the phylogenetic tree (a). The score on the first two axes (Axis 1 = 79.05% and Axis 2 = 15.32%) from the matrix of genetic distances estimated with the 106 COI haplotypes are indicated.Full size imageConsistently, the TCS network analysis (Fig. 2c) and principal coordinate analysis (PCoA) (Fig. 2d) showed the existence of three distinguished genetic groups among L. japonica, in accordance with the phylogenetic analyses (Fig. 2a − b). The TCS network (Fig. 2c) revealed that Lineages S1 and S2 were separated by 18 mutation steps, which is far shorter than the distance between Lineages N and S1 (37 mutation steps) or between Lineages N and S2 (60 mutation steps), indicating that Lineages S1 and S2 have a close affinity and had only recently diverged from each other. The overwhelming dominancy of the A1 haplotype implies a recent and rapid population expansion of Lineage N. In addition to A1, it was found that haplotypes B2 for Lineage S1 and C21 for Lineage S2 were dominant. In the PCoA plot (Fig. 2d), the three genetic groups of L. japonica were also observed, as in the phylogenetic (Fig. 2a,b) and TCS network (Fig. 2c) analyses. Lineage N was distantly located from Lineages S1 and S2, while Lineages S1 and S2 were spatially much closer.Sample collection of L. japonica and the genetic diversity of 16S rRNA
The 342 individuals of L. japonica from 12 localities in the intertidal coasts on the Korean Peninsula and Japanese Archipelago (Fig. 1) were subjected to PCR amplification of a partial region of 16S rRNA (506 bp) (Tables S5, S6). Of these, only 299 samples were successfully amplified and sequenced. Based on 299 individual 16S rRNA sequences, a total of 23 16S rRNA haplotypes of L. japonica were detected (Tables S5, S6). Combined with 11 haplotypes extracted from 16S rRNA sequences of 125 L. japonica individuals known previously in southern China, we totaled 34 16S rRNA haplotypes from 425 L. japonica individuals in 13 collection localities. The average haplotype (h) and nucleotide (π) diversities were 0.702 and 0.02093, respectively; the highest haplotype diversity was found in Geojedo (GJ; h = 0.833), and the highest nucleotide diversity in Wando (WD; π = 0.02244), located in the South Sea of the Korean Peninsula. Overall, the average haplotype and nucleotide diversities of 16S rRNA were lower than those of COI (Table S1). As shown in Table S7, the population distribution pattern of 16S rRNA haplotypes revealed that most of the collection sites had site-specific haplotype(s), except for BS, GJ, WD, SA, JJ, and TT. The most abundant haplotype was RA1, which was found in 186 (48.1%) out of the 16S rRNA sequences obtained from 425 L. japonica individuals.Phylogenetic and population genetic analyses based on 16S rRNA
We constructed a nucleotide sequence alignment set with 34 16S rRNA haplotypes of L. japonica (Data S2), and identified 35 polymorphic sites (6.9%; Table S8) and 24 parsimoniously informative sites (4.7%). Phylogenetic analyses, including ML, BI, and NJ analyses, were conducted with the outgroup Acanthopleura echinata (Table S6). The resultant phylogenetic trees (Fig. S2) and unrooted phylogenetic network (Fig. 3a) consistently supported the three distinct genetic lineages of L. japonica, with the phylogenetic relationship between Lineages S1 and S2 being much closer than those inferred from the results of COI (Fig. 2a,b; Fig. S1). The TCS network (Fig. 3b) revealed that Lineages S1 and S2 were closely connected with only with 4–5 mutation steps between them, while Lineages N and S1 or Lineages N and S2 were remotely distanced by 18 mutation steps. Also, the overwhelming dominance of the RA1 haplotype implied a recent and rapid population expansion of Lineage N. In addition to RA1, haplotypes RB1 for Lineage S1 and RC1 and RC2 for Lineage S2 were dominant (Fig. 3b; Table S7). Consistent with this, in the PCoA plot (Fig. 3c), the three genetic groups of L. japonica were spatially separated. Lineage N was distantly located apart from Lineages S1 and S2, while Lineages S1 and S2 were spatially much closer.Figure 3The results of phylogenetic and population genetic analyses based on 34 16S rRNA haplotypes from 425 individuals of Liolophura japonica inhabiting coastal areas of the northwestern Pacific Ocean. (a) Phylogenetic network reconstructed using the neighbor net algorithm, showing three different genetic lineages for L. japonica: Lineage N, Lineage S1, and Lineage S2. The 16S rRNA sequence alignment set used is shown in Data S2. Detailed information of 34 16S rRNA haplotypes used in these analyses is summarized in Table S5 and S6. (b) An unrooted TCS network. There are distinctly observed three different genetic groups, corresponding to the three genetic lineages shown in the phylogenetic network (a). The haplotype frequency is displayed by the circle size. (c) A two-dimensional PCoA plot showing the three distinct genetic groups, corresponding to Lineage N, Lineage S1, and Lineage S2. The score on the first two axes (Axis 1 = 87.77% and Axis 2 = 4.4%) from the matrix of genetic distances estimated with the 34 16S rRNA haplotypes are indicated.Full size imageExamination of species discrimination of L. japonica based on COI and 16S rRNA
Using the Automatic Barcode Gap Discovery (ABGD), we performed distribution of pairwise genetic divergences, ranked pairwise difference, and automatic partition analyses based on COI and 16S rRNA of L. japonica, respectively (Fig. 4a–c), which confirmed that there were distinct barcoding gaps between intraspecific and interspecific variations, strongly supporting the possibility of species discrimination of L. japonica. the COI-based analysis yielded two different barcoding gaps, while the 16S rRNA-based analysis revealed only a single barcoding gap (Fig. 4a–c). The results of automatic partition at each value of the prior intraspecific divergence (P) divided L. japonica into three groups by COI and two groups by 16S rRNA, respectively (Fig. 4a–c). We also implemented two DNA taxonomy approaches to evaluate the possibility of species discrimination based on COI: the general mixed Yule coalescent (GMYC) approach (Fig. S3) and a Bayesian implementation of a Poisson Tree Processes model (bPTP) (Fig. S4). The results consistently and robustly supported the possibility that L. japonica can be divided into three different species, as shown in the results of ABDG (Fig. 4a–c).Figure 4Distribution of pairwise genetic divergences, ranked pairwise difference, and automatic partition based on COI and 16S rRNA haplotypes of Liolophura japonica and a COI-based NJ tree showing the phylogenetic relationship with a congeneric species L. tenuispinosa. (a) Distribution patterns of pairwise genetic divergences observed in COI and 16S rRNA for L. japonica. The horizontal axis represents intervals of pairwise Kimura-2-parameter (K2P) genetic distance in percentage, and the vertical axis represents the number of individuals associated with each distance interval. (b) The results of ranked pairwise differences based on COI and 16S rRNA, ranked by ordered value, which is similar to the distribution of pairwise genetic divergence in (a). The horizontal axis indicates a ranked ordered value based on K2P genetic distance, and the vertical axis represents the K2P genetic distance in percentage. (c) The results of automatic partition analyses based on COI and 16S rRNA. The horizontal axis represents the prior maximum intraspecific divergence (P), and the vertical axis represents the number of groups inside the partitions (primary and recursive). (d) A COI-based NJ tree with L. tenuispinosa. Refer to Fig. S3 and Data S3.Full size imageThe molecular variance analyses using analysis of molecular variance (AMOVA), based on COI and 16S rRNA, were conducted to evaluate the degree of genetic differentiation among Lineages N, S1, and S2 (Tables S9, S10). According to the results, supposing that there are three genetic lineages (N, S1, and S2) or two genetic lineages (N and S1/S2), almost all variation in both cases is attributed to variation among groups (= among lineages), whereas variations within populations (within lineages) exhibit negative values in common. We confirmed that there was a high degree of genetic differentiation among Lineages N, S1, and S2, which supports the results of the COI barcoding gap analysis shown in Fig. 4a–c, although this was not statistically significant (P  > 0.05; Table S9). When we assumed only two genetic groups, Lineages N and S1/S2, the genetic differentiation between the two groups was statistically significant (P  0.05) (Tables S9 and S10). The discrepancy between the number of barcoding gaps inferred from COI and 16S rRNA may have been affected by different gene evolutionary rates of the molecular markers11; nucleotide substitution rate of 16S rRNA is known to be generally slower than that of COI (which is especially fast in the third codon positions: 105 out of 127 polymorphic sites). When an ML tree was constructed based on 22 polymorphic sites, which are found only in the first and second codon positions of COI that are much more conserved than the third codon position, the three genetic lineages were retained in the resultant tree (Fig. S5), but Lineage S2 was nested within Lineage S1, as in the trees inferred from 16S rRNA (Fig. S2). Reflecting the powerful resolution of the COI barcoding marker well known from animals12 and the high degree of variation among the three genetic lineages (Fig. 4, Figs. S3, S4), we suggested that L. japonica could be categorized into three different species: L. koreana, Yeo and Hwang, sp. nov. for Lineage N, L. japonica for Lineage S1, and L. sinensis Choi, Park, and Hwang, sp. nov. for Lineage S2. To examine whether it is reasonable to give these a species-level taxonomic status, as shown in Fig. 4d, we reconstructed a COI-based NJ tree with one congeneric species L. tenuispinosa13, which was originally described as a subspecies-level taxon of L. japonica14,15 and was then revised as an independent species closely related to L. japonica by Saito & Yoshioka16 in 1993. The resultant tree (Fig. S6 and Data S3) showed that L. tenuispinosa forms a sister group with L. japonica (Lineage S1). This likely indicates that L. koreana and L. sinensis have taxonomic status as independent species.Morphological comparison and geographical distribution of the three Liolophura speciesWe compared morphological characteristics among Liolophura koreana, sp. nov. (Lineage N), L. japonica (Lineage S1), and L. sinensis, sp. nov. (Lineage S2). Their morphological appearances are shown in Fig. 5a–c, which indicated that black spots on the tegmentum (Fig. 5d–e) and shapes of spicules on the perinotum (Figs. 5f–k, 6e–f) represent key morphological characteristics to distinguish them from each other. Although black dots in pleural areas, which are between the middle and lateral areas of the tegmentum on valves II–VII (or VIII), are commonly shared in all three lineages (Fig. 5a–c), other black spots on the valves exhibit a high degree of variation in morphology (Fig. 5a–c, Fig. S7). Herein, we described a new species of genus Liolophura, that is, L. koreana Yeo and Hwang from South Korea and Japan, with detailed descriptions of morphological characteristics observed by light microscopy (M205, Leica Camera AG, Germany) and FE-SEM (SU8220, Hitachi, Japan). In addition, we suggested the divergence of a new species, L. sinensis Choi, Park, and Hwang from southern China, with simple remarks based on distinct genetic difference (mainly COI barcoding gaps), with possible unique morphological characteristics as follows.Figure 5Morphological comparison of Liolophura koreana, sp. nov., L. japonica, and L. sinensis, sp. nov. (a–c) Photos of dorsal views of the individuals belonging to L. koreana (Lineage N), L. japonica (Lineage S1), and L. sinensis (Lineage S2) in order. (d,e) Morphological comparison of pleural and lateral black spots on valves III and IV of the tegmentum of L. koreana (d; holotype) and L. japonica (e). (f,g) Morphological comparison of spicules on the perinotum of L. koreana (f; holotype) and L. japonica (g). (h–k) Morphological comparison of the spicule of L. koreana (h,i; paratype) and L. japonica (j,k) in lateral and dorsal views. The scale bar marks 2.0 mm (d,e), 1.0 mm (f,g), and 0.5 mm (h–k). The photos were edited using Adobe Photoshop v.22.2 (https://www.adobe.com).Full size imageFigure 6Microstructural comparison of Liolophura koreana, sp. nov. and L. japonica using field emission scanning electron microscopy (FE-SEM). (a,b) Middle and lateral areas on the tegmentum of the holotype of L. koreana. (c,d) Middle and lateral areas on the tegmentum of L. japonica. Arrows indicate that morphological difference of the posterior valve margin of the valve II between two species. The scale bar marks 1.0 mm. (e,f) The occurrence frequency, and shape and structure differences of the spicules on the perinotum between the holotype of L. koreana (e) and L. japonica (f). The scale bar marks 1.0 mm and 0.2 mm, respectively. The photos were edited using Adobe Photoshop v.22.2 (https://www.adobe.com).Full size image
Liolophura koreana Yeo and Hwang, sp. nov. (Figs. 5, 6; Figs. S7, S8)(urn:lsid:zoobank.org:act:4418355E-F55C-44FA-B4CE-585589FDCD23).Type specimens examined[Holotype] SOUTH KOREA: 1 specimen, Jeju-do, Seogwipo-si, Seongsan-eup, Seopjikoji, 3.XI.2020, UW Hwang, B Park & CR Shin (LEGOM040501); [Paratypes] SOUTH KOREA: 1 specimen, Gyeongsangbuk-do, Pohang-si, Guryongpo-eup, Janggil-ri, 27.VII.2008, UW Hwang (LEGOM040502); 3 specimens, Gyeongsangbuk-do, Ulleung-gun, Seo-myeon, Namyang-ri, Ulleungdo Island, Namtong tunnel, 12.VI.2007, UW Hwang (LEGOM040503–0505); 1 specimen, Gyeongsangbuk-do, Ulleung-gun, Namyang-ri, Ulleungdo Island, Namyang tunnel, 5.X.2007, UW Hwang (LEGOM040506); 2 specimens, Gyeongsangnam-do, Geoje-si, Nambu-myeon, Dapo-ri, 28.IV.2009, MY Yeo (LEGOM040507,0508); 4 specimens, same data as the holotype (LEGOM040509–0512); 2 specimens, Jeollanam-do, Yeosu-si, Hwajeong-myeon, Sado-ri, Sado Island, 8.IV.2008, MY Yeo (LEGOM040513,0514); 4 specimens, same data as the holotype (LEGOM040515–0518); JAPAN: 6 specimens, Tottori Prefecture, Hakuto, 24.V.2009, UW Hwang (LEGOM040519–0524); 1 specimen, Tottori Prefecture, Iwato, 25.V.2009, UW Hwang (LEGOM040525).DescriptionBody small-sized and broad oval- to oval-shaped (Fig. 5a; Fig. S7); length 3.9 (1.9–12.3) mm and width 2.4 (1.2–7.1) mm. Tegmentum entirely brown (dark brown or black entirely, or each valve with black line anteriorly or white line laterally), with black dots on the pleural areas of valves II–VII (or VIII) (Fig. 5a,d, Fig. S7); articulamentum entirely black (dark brown); whitish and blackish spicules on the perinotum scattered irregularly, sometimes forming a band besides each valve (Fig. 5a, Fig. S7). Surface of the tegmentum in middle and lateral areas as in Fig. 6a,b and Fig. S8; posterior margin of the head valve nearly straight; dorsal shape of intermediate valves round-backed and side slopes slightly convex; the posterior valve margin with a distinct central apex, its shape subtriangular to triangular (rounded or linear), particularly valve II, mainly with a strong projection (Fig. 6a). Perinotum covered with large, solid, slightly curved, and obtusely pointed spicules (rarely with smooth and radial ribbed spicules apically), its density relatively lower than that of L. japonica (Figs. 5f,h,i, 6e).DistributionSouth Korea, Japan; avobe 33°24′ N (Seogwipo, JJ) in South Korea and TT and MY in Japan (Fig. 7).Figure 7Geographical distribution of Liolophura koreana, sp. nov., L. japonica, and L. sinensis, sp. nov. inhabiting coastal areas of the northwestern Pacific Ocean. A COI-based map showing geographical distribution of L. koreana, L. japonica, and L. sinensis on the northwestern Pacific coast. L. koreana are found in a wide range of South Korea and Japan above ca. 33°24′ N (JJ), L. japonica in mainly southern coastal areas of South Korea and Japan below ca. 35°53′ N (TT), and L. sinensis in ZJ of southern China around ca. 27°02′ N–28°00′ N. The sympatric distribution of L. koreana and L. japonica is shown between 33°24′ and 35°53′ N. Table S1–S3 contain the full names of localities and detailed haplotype information. The question mark indicates that collection of Liolophura samples from such coastal areas in Japan is required to clarify distribution patterns of L. koreana and L. japonica in the East Sea (= Sea of Japan). The basic map was obtained from a free map-providing site (https://d-maps.com), which was modified using Adobe Illustrator v.25.2. (https://www.adobe.com).Full size imageHabitatThis new species appears to be attached to rocks in coastal areas with strong waves, or a calm inner shore in the northwestern Pacific Ocean (Fig. S9).EtymologyThe species is named per the locality of the new species.RemarksWe found that Liolophura koreana, sp. nov. has no black spots on lateral areas of the tegmentum (Fig. 5d), and large, slightly curved, and obtusely pointed spicules on the perinotum compared to those of L. japonica (Figs. 5f,h,i, 6e). On the other hand, L. japonica from southern South Korea and southern Japan has two black spots on the lateral areas of valves II–VII (or VIII) (Fig. 5e), and small, almost straight, and cylindrical spicules compared to those of L. koreana (Figs. 5g,j,k, 6f).As shown in Fig. 7, L. koreana (Lineage N) was observed in all the South Korean and Japanese populations examined here, except for the EH population in Japan (refer to Tables S1, S2), which were found from JJ at 33°24′ N to MY at 38°32′ N. On the other hand, L. japonica (Lineage S1) was found from the southern coastal areas of South Korea and Japan, which were found only between JJ at 33°24′ N and TT at 35°53′ N. Interestingly, we found only L. koreana north from the latitude of 35°10′ N (BS) such as UL/DD (37°24′ N) and PH (36°02′ N) in South Korea. In Japan, there was found only L. koreana at MY (38°32′ N) too, but it remains to be explored to clarify its distribution range in Japan with much more sample collections covering northern Japanese coastal areas through further study. It was also confirmed that L. koreana and L. japonica show a sympatric distribution pattern between JJ at 33°24′ N and TT at 35°53′ N on the southern coastal area of South Korea and Japan.
Liolophura sinensis Choi, Park, and Hwang, sp. nov. (Fig. 5c).(urn:lsid:zoobank.org:act:72DF7E75-1853-4F23-AC12-3AB8CD054187).Type specimens examined[Holotype] CHINA: 1 specimen, Zhejiang Province, Dongtou Island, 27°49′57.44″ N, 121°10′19.13″ E, 2017; [Paratypes] CHINA: 1 specimen, Beiji Island, 27°37′08.82″ N, 121°11′47.82″ E, 2017; 1 specimen, Beilongshan Island, 27°40′08.56″ N, 121°58′51.56″ E, 2017; 1 specimen, Chaishi Island, 27°25′40.36″ N, 121°04′54.45″ E, 2017; 1 specimen, Daleishan Island, 27°29′39.48″ N, 121°05′24.50″ E, 2017; 1 specimen, Daqu Island, 27°47′29.92″ N, 121°05′23.97″ E, 2017; 1 specimen, Dazhushi Island, 27°49′12.87″ N, 121°12′48.74″ E, 2017; 1 specimen, Dongce Island, 27°45′32.04″ N, 121°09′01.38″ E, 2017; 1 specimen, Dongxingzai Island, 27°02′40.36″ N, 121°02′47.98″ E, 2017; 1 specimen, Houjishan Island, 27°28′26.96″ N, 121°07′40.71″ E, 2017; 1 specimen, Luxi Island, 27°59′33.43″ N, 121°12′50.70″ E, 2017; 1 specimen, Nanji Island, 27°27′30.57″ N, 121°03′06.28″ E, 2017; 1 specimen, Nanpanshan Island, 28°00′15.29″ N, 121°15′33.62″ E, 2017.DistributionSouthern China; Zhejiang Province around ca. 27°02′–28°00′ N (Fig. 7).HabitatThis new species appears to be attached to rocks in coastal areas of the northwestern Pacific Ocean.EtymologyThe species is named after its locality.RemarksL. sinensis, sp. nov. was examined and established mainly by molecular data such as the COI barcoding gap (Fig. 4a–c) presented in this study and photos provided from Prof. Yong-Pu Zhang (Wenzhou University, Zhejiang Province, China) without any direct real sample observation. The morphology of this new species is very similar to that of a previously known species, L. japonica. L. sinensis from southern China has black spots on the lateral areas of valves II–VII similar to L. japonica, but black bowtie-shaped spots anterior to valves II–VII (Fig. 5c) are a unique characteristic for L. sinensis. L. sinensis (Lineage S2) was found in Zhejiang Province (ZJ) in southern China, around ca. 27°02′–28°00′ N (Fig. 7; Tables S1, S5).Demographic history and divergence time estimation analysesMismatch distribution analyses (MDA) based on COI were performed for L. koreana, L. japonica, and L. sinensis, respectively. The MDA results (Fig. 8a) showed a unimodal curve for each of the three lineages. In addition, when neutrality tests were performed with COI and 16S rRNA (Table S11), all three showed statistically significant negative values in both Tajima’s D and Fu’s Fs, except for the Tajima’s D values in COI (L. sinensis) and 16S rRNA (L. japonensis and L. sinensis), implying that these had experienced population expansions. Bayesian skyline plot (BSP) analyses with COI (Fig. 8b) were performed to examine the fluctuation patterns in effective population sizes for L. koreana, L. japonica, and L. sinensis, respectively. The effective population sizes of L. japonica and L. sinensis had gradually grown between ca. 100 Ka and ca. 50 Ka, while those of L. japonica had grown between ca. 80–50 Ka, L. sinensis had grown between 100–60 Ka, and that of L. koreana had begun to rapidly expand ca. 85 Ka and ceased ca. 75 Ka. This indicated that population expansion had occurred more dramatic in L. koreana than in L. japonica and L. sinensis, following the last interglacial age, called the Eemian (129–116 Ka). As shown in Fig. 8c and Fig. S10, according to the molecular clock analysis by the BEAST program, it was estimated that L. japonica and L. koreana shared their most recent common ancestor about 3.37 Ma, around the mid-Pliocene warm period (3.30–3.00 Ma), before the extensive glaciation in the late Pliocene (ca. 3.00 mya). L. japonica and L. sinensis likely diverged around 1.84 Ma, around the beginning stage of the Early Pleistocene Transition (EPT; 1.85–1.66 mya). The augmentation of haplotype diversity in L. japonica, L. sinensis, and L. koreana might have intensified in the interglacial stages during the late-middle (0.35–0.126 Ma) and late Pleistocene (0.126–0.012 Ma), before the last glacial maximum (LGM: 0.026–0.019 Ma).Figure 8The results of mismatch distribution analyses (MDA), Bayesian skyline plots (BSPs), and molecular clock analysis performed with COI haplotypes for Liolophura koreana, sp. nov., L. japonica, and L. sinensis, sp. nov. (a) MDA plots resulting in a unimodal curve for L. koreana, L. japonica, and L. sinensis. Dotted lines indicate the observed distribution of mismatches, and solid lines represent the expected distribution under a demographic expansion model. (b) BSP results showing the demographic history of population expansions of L. koreana, L. japonica, and L. sinensis. The graph in gray depicts sea level changes during the last 330 Ka. (c) Time-calibrated Bayesian tree reconstructed using BEAST with the inference of ancestral areas under the Bayesian binary MCMC (BBM) model implemented in RASP ver 3.2. Ancestral areas were hypothesized based on the distribution range of the fossil records of Mopalia and the contemporary distribution of L. koreana, L. japonica, and L. sinensis. LGM indicates the last glacial maximum (0.026–0.019 Ma; blue vertical bar) and three interglacial periods are indicated by light green boxes during the late-middle and late Pleistocene. The pictures were edited using Adobe Illustrator v.25.2. (https://www.adobe.com).Full size image More

We consider public goods games played iteratively over a fixed connected network. The vertices of the network represent the players and the edges represent neighboring connections5,10,11,12. The dynamics evolve in discrete rounds. In each round, each player chooses a behavior that minimizes its cost, where the player’s cost is affected by its own behavior and the behaviors of its neighbors.Our main model includes three behavior types, namely, defection, hypocrisy, and cooperation, in which those who hardly contribute to the social welfare, i.e., defector and hypocritical players, face the risk of being caught and punished by their neighbors who are non-defectors. The level of risk together with the extent of punishment is captured by a notion that we call “social-pressure”. The main result is that adjusting the level of social-pressure employed against hypocritical players compared to the one employed against defectors can have a dramatic impact on the dynamics of the system. Specifically, letting the former level of social-pressure be within a certain range below the latter level, allows the system to quickly transform from being composed almost exclusively of defectors to being fully cooperative. Conversely, setting the level to be either too low or too high locks the system in a degenerate configuration.As mentioned, our main model assumes that non-defectors induce mild social-pressure on the defectors among their neighbors. This implicitly assumes that inducing the corresponding social-pressure is beneficial (e.g., allows for a social-upgrade), although other explanations have also been proposed21. To remove this implicit assumption we also consider a generalized model, called the two-order model, which includes costly punishments. Consistent with previous work on the second-order problem, e.g.,23,25,26,27,36,40, this model distinguishes between first-order cooperation, that corresponds to actions that directly contribute to the social welfare, and second-order cooperation, that corresponds to applying (costly) social-pressure, or punishments, on others. As in the main model, the level of punishment employed against first-order defectors may differ from that employed against second-order defectors. We identify a simple criteria for the emergence of cooperation: For networks with minimal degree (Delta), cooperation emerges when two conditions hold. The first condition states that the cost (alpha _2) of employing punishments against second-order defectors should be smaller than the corresponding punishment (beta _2) itself, i.e., (alpha _2 More

Technical criterion for ecosystem status evaluationOn March 3, 2015, the Ministry of Ecology and Environment of the People’s Republic of China approved the “Technical Criterion for Ecosystem Status Evaluation” as the national environmental protection standard. This standard is based on the former standards released in 2006, and 48 relevant documents from 2006 to 2012 were searched to propose new standards and factor weights based on actual utilization effects and expert guidance. The eco-environment assessment uses a comprehensive index (eco-environmental status index, EI) to reflect the overall state of the regional eco-environment. The indicator system includes the biological abundance index, vegetation coverage index, river density index, land stress index, and pollution loading index. These indexes reflect the abundance of organisms in the evaluated area, the level of vegetation coverage, the abundance of water, the intensity of land stress, and the extent of the pollution load. Each indicator was calculated according to its weight to obtain an eco-environment assessment map (Table 1). All parameters involved in the calculation are derived from this standard.Table 1 Weights of the evaluation indicators.Full size tableThe calculation of the eco-environment status is as follows:$$begin{aligned} EI & = 0.35*biological ; abundance ; index + 0.25*vegetation ; coverage ; index hfill \ & quad + 0.15*river ; density ; index + 0.15*(1 – land ; stress ; index) hfill \ & quad + 0.1*(1 – pollution ; loading ; index) hfill \ end{aligned}$$
(9)
Biological abundance indexThe biological abundance index refers to the number of certain organisms in this area. The calculation method is as follows:$$Biological , abundance , index , = , left( {BI , + , HQ} right)/2$$
(10)
In this formula, BI is the biodiversity index and HQ is the habitat quality index. When the biodiversity index does not have dynamic data updates, the change in the biological abundance index is equal to the change in the HQ.Biodiversity is a general term for the complexity of species and their genetic variation and ecosystems in space over time. Biodiversity plays an important role in maintaining soil fertility, ensuring water quality, regulating the climate, stabilizing the environment, and maintaining ecological balance.The BI method is as follows:$$BI = NPP_{mean} *F_{pre} *F_{tem} *(1 – F_{alt} )$$
(11)
NPPmean is the net primary productivity. Fpre is the annual average precipitation. Ftem is the temperature parameter. Falt is the altitude parameter.NPP refers to the amount of organic matter accumulated per unit area and unit time of green plants. NPP is the remainder of the total amount of organic matter produced by photosynthesis after deducting autotrophic respiration and is usually expressed as dry weight. In this study, the estimation of NPP was based on the absorbed photosynthetically active radiation (APAR) and actual light-use efficiency (LUE) (ε) of the CASA ecosystem model40. The CASA model is a process-based remote sensing model that couples ecosystem productivity and soil carbon and nitrogen fluxes, driven by gridded global climate, radiation, soil, and remote sensing vegetation index datasets41. The model can be expressed generally as follows:$$NPP(x,t) = APAR(x,t)*varepsilon (x,t)$$
(12)
The entire study area is divided into 11,303 pixels on a 30 * 30 m grid. x indicates the location of each pixel, and t indicates time; the data were collected once a month. APAR(x,t) represents the photosynthetically active radiation absorbed by pixel x in that month (gC * m−2* month−1). Ɛ(x, t) is LUE (gC * MJ−1) of the vegetation42.Estimation of the fraction of APAR using RS data is based on the reflection characteristics of the vegetation on the infrared and near-infrared bands. The value of APAR is determined by the effective radiation of the sun and the absorption ratio of the vegetation to the effective photosynthetic radiation. The formula is as follows:$$APAR(x,t) = SOL(x,t)*FPAR(x,t)*0.5$$
(13)

where SOL(x,t) represents the total amount of solar radiation at pixel x in month t, FPAR(x,t) represents the absorption ratio of the vegetation layer to the incident photosynthetically active radiation, and a constant of 0.5 indicates the ratio of the effective solar radiation that can be utilized by the vegetation to the total solar radiation.Since there is a linear relationship between FPAR and NDVI within a certain range, this relationship can be determined according to the maximum and minimum values of a certain vegetation type NDVI and the corresponding FPAR maximum and minimum values.$$FPAR(x,t) = frac{{(NDVI(x,t) – NDVI_{i,min } )}}{{(NDVI_{i,max } – NDVI_{i,min } )}}*(FPAR_{max } – FPAR_{min } ) + FPAR_{min }$$
(14)

where NDVImax and NDVImin correspond to the NDVI maximum and minimum values of the ith planting type, respectively. There is also a good linear relationship between FPAR and the simple ratio index (SR) of vegetation, which is represented by the following formula:$$FPAR(x,t) = frac{{(SR(x,t) – SR_{i,min } )}}{{(SR_{i,max } – SR_{i,min } )}}*(FPAR_{max } – FPAR_{min } ) + FPAR_{min }$$
(15)

where the values of FPARmin and FPARmax are independent of vegetation type and are 0.001 and 0.95, respectively; SRi,max and SRi,min correspond to the 95% and 5% percentiles, respectively, of the ith NDVI. SR(x,t) is represented by the following formula:$$SR(x,t) = frac{1 + NDVI(x,t)}{{1 – NDVI(x,t)}}$$
(16)
A comparison of the estimated results of FPAR-NDVI and FPAR-SR shows that the FPAR estimated by NDVI is higher than the measured value, while the FPAR estimated by SR is lower than the measured value, but the error is less than that estimated directly by NDVI. As a result, these two values can be combined, and their weighted average value is taken as an estimate of the estimated FPAR, while ɑ means weight:$$FPAR(x,t) = alpha FPAR_{NDVI} + (1 – alpha )FPAR_{SR}$$
(17)
Light use efficiency (LUE) refers to the ratio of chemical energy contained in organic dry matter produced per unit area over a certain period of time to the photosynthetically active radiation absorbed by plants projected onto the same area at the same time. Different vegetation types and the same types of vegetation have different light energy utilization rates in different living environments43. The differences are mainly due to the characteristics of the vegetation itself, temperature, moisture, and soil44. Vegetation has the highest utilization rate of light energy under ideal conditions, but the maximum light energy utilization rate in the real environment is mainly affected by temperature and moisture, which can be expressed as follows:$$varepsilon left( {x,t} right) = T_{varepsilon 1} left( {x,t} right) cdot T_{varepsilon 2} left( {x,t} right) cdot W_{varepsilon } left( {x,t} right) cdot varepsilon_{max }$$
(18)
where Tε1(x,t) and Tε2(x,t) represent the stress effects of low temperature and high temperature on light energy utilization, respectively, Wε(x,t) is the effect of water stress on the maximum light energy utilization under ideal conditions, and εmax is the maximum light energy utilization under ideal conditions (gC * MJ−1). The maximum solar energy utilization rate εmax varies depending on the vegetation type. In this study, the maximum light energy utilization rate of different land use types simulated by an improved Carnegie-Ames-Stanford Approach (CASA) model is used as the input parameter of light energy utilization in the CASA model (Table 2). The monthly maximum light energy utilization rate is determined in three steps: first, calculate the APAR, temperature, and water stress factors of all pixels; then, select the NPP measured data of the same time period in the study area; finally, simulate the εmax of vegetation according to the principle of minimum error45. Figure 5 shows the calculation process of NPP. The weight of each habitat type in the HQ is shown in Table 3. The weight value is derived from the official document30. To facilitate the calculation, this paper normalizes the calculation results from 0 to 1 (Fig. 6a).Table 2 Maximum LUE rates of different land use types.Full size tableFigure 5NPP calculation process.Full size imageTable 3 Weight of each habitat type in the HQ.Full size tableVegetation coverage indexThe vegetation coverage index was obtained from the NDVI, which is a simple, effective, and empirical measure of surface vegetation status. The vegetation index mainly describes the difference between the reflection of vegetation in the visible and near-infrared bands and the soil background. This index also reduces the solar elevation angle and noise caused by the atmosphere and is thus the most widely used and effective calculation method. Each vegetation index can be used to quantitatively describe the growth of vegetation under certain conditions. The expression is as follows:$$NDVI = frac{NIR – R}{{NIR + R}}$$
(19)

where NIR and R are reflectance values in the near-infrared and red bands, respectively.NDVI values are obtained by processing the RS images of the Landsat 8 satellite. This satellite is equipped with an operational land imager (OLI) that includes nine bands with a spatial resolution of 30 m, including a 15-m panchromatic band. To facilitate the calculation, this paper normalizes the calculation results from 0 to 1 (Fig. 6b).Figure 6Eco-environment assessment indexes and evaluation rating map (the first quarter was used as an example). (a) Biological abundance index; (b) vegetation coverage index; (c) river density index; (d) land stress index; (e) pollution loading index; and (f) environmental status classification. The Figure is created using ArcGIS ver.10.3 (https://www.esri.com/).Full size imageRiver density indexThe river density index refers to the total length of rivers, lakes, and water resources in the assessed area as a percentage of the assessed area, which is used to reflect the abundance of water in the assessed area and is calculated as follows:$$begin{gathered} River ; density ; index = (A_{riv} *river ; length / area + A_{lak} *water ; area/area hfill \ + A_{res} {*}amount ; of ; resources/area , )/3 hfill \ end{gathered}$$
(20)

where Ariv is the normalization coefficient of river length, with a reference value of 84.3704, Alak is the normalization coefficient of the lake area, with a reference value of 591.7909, and Ares is the normalization coefficient of water resources, with a reference value of 86.387. Finally, the calculation results were normalized from 0 to 1 (Fig. 6c).Land stress indexThe land stress index is the degree to which the land quality in the assessment area is under stress. The weight of the land stress index evaluation is shown in Table 4.Table 4 Weight of the land stress index evaluation.Full size tableThe calculation method is as follows:$$begin{gathered} Land ; stress ; index = A_{ero} *(0.4*severe ; erosion ; area + 0.2*{text{mod}}erate ; erosion ; area hfill \ + 0.2*construction ; land ; area + 0.2*other ; land ; stress ; area)/area hfill \ end{gathered}$$
(21)

where Aero is the normalization coefficient of the land stress index, with a reference value of 236.0436. According to the “Classification criteria for soil erosion”46, the influencing factors of soil erosion, vegetation, soil texture, landform, and precipitation are ranked according to importance. In the calculation of the land stress index, all the land is divided into three categories, in which the weight of severe erosion is 0.4, the weight of non-erosion is 0, and the other erosion types such as moderate erosion and construction land are 0.2. The areas with severe erosion include vegetation coverage less than 30% and areas of soil erosion greater than 3.7 mm/a due to human activities. These areas are generally developed on highly erosive-sensitive soils. Cinnamon soil and loess soil in the study area are highly erosive-sensitive soils. Therefore, the industrial and mining areas of cinnamon and loess soil types are regarded as severe erosion areas. Areas with vegetation coverage greater than 50% are non-eroded areas, so water bodies and woodlands are divided into non-erodible areas. All areas except these two types have a weight of 0.2. Finally, the calculation results were normalized from 0 to 1 (Fig. 6d).Pollution loading indexThe pollution loading index refers to the load of pollutants in a certain area or an environmental element. In this study, the AQI was used to calculate the pollution loading index, and the results were normalized from 0 to 1 (Fig. 6e).The eco-environmental evaluation score was calculated based on the national environmental protection standard according to the weight of each indicator (Fig. 6f).Improved evaluation system and intelligent evaluation modelImproved evaluation systemConsidering that the evaluation factors in the national environmental protection standards are applicable to ordinary areas, areas affected by mines should have more evaluation factors than those in the standards. Thus, an improved evaluation system was proposed. The improved evaluation system has added factors that affect the environment of the mine based on the factors of the original system. The impact of the mine on the environment is reflected in the pollution of the atmosphere, such as dust from open pits and industrial waste from concentrators; the occupation of land by solid waste, such as ore piles and coal piles; soil pollution, such as the diffusion of heavy metals from coal piles, coal mine concentrator plants, and mines; and the increased likelihood of geological disasters, such as collapse caused by underground mining, spontaneous coal combustion and landslides caused by open-pit mining surfaces. Therefore, the improved evaluation system adds an air pollution range, a solid waste area, a geologic hazard range, and a metallic and non-metallic mine soil pollution buffer to the national environmental protection standards.The area of air pollution in mining areas is generally near open pits and concentrator plants. Therefore, the air pollution range was selected within 50 m around the open pits and the concentrator plants. Due to the dilution and dispersion of the air itself, an estimate of the pollution is the reciprocal of the wind speed (Fig. 7a).Figure 7New factors in the improved evaluation system. (a) Air pollution range; (b) solid waste area; (c) geological hazard range; (d1) non-metallic mine soil pollution buffer; and (d2) metal mine soil pollution buffer. The Figure is created using ArcGIS ver.10.3 (https://www.esri.com/).Full size imageMine solid waste pollution includes a large amount of waste rock from open-pit mining and pit mining, coal gangue produced by coal mining, tailings from beneficiation and slag from smelting. These solid wastes are generally piled up near the mining area. They not only occupy large areas of land and induce geological disasters such as landslides and mudslides but also cause chemical pollution, spontaneous combustion, and radiation from radioactive materials due to long-term stacking. This may affect the health and safety of humans and other biological organisms. The scope of the solid waste area is determined by the ore piles, coal piles, and dumping sites (Fig. 7b).Mine geological disasters are caused by a large number of mining wells and rock and soil deformation, as well as serious changes in the geological, hydrogeological, and natural environments of the mining area, endangering human life and property and destroying mining engineering equipment and mining resources. In this study, the geologic hazard range consists of areas with underground mining stopes and coal piles (Fig. 7c).After the pollutants generated by the mining operation enter the soil, physical and mechanical absorption, retention, colloidal physicochemical adsorption, chemical precipitation, bioabsorption, etc. of the suspended pollutants through the soil continue to accumulate in the upper soil. When pollutants reach a certain maximum, they cause deterioration of the soil composition, cycle, properties, and functions and begin to accumulate in plants, which affects the normal growth and development of plants, decreases crop yield and quality, and ultimately affects human health. Metallic and non-metallic minerals have different effects on soil pollution. The pollution of soil by non-metallic minerals mainly occurs in coal mines and coal piles, and the buffer zone is centred on the coal mines and coal piles. Coal production activities can cause heavy metals in coal piles to enter the soil and cause pollution. Due to different types of heavy metals, the range of soil contamination is different47. Combining the non-metallic mineral industrial squares and coal mine-based non-metallic minerals around the heavy metal soil pollution range, the buffers are graded at 30, 200, and 1000 m (Fig. 7d1)43. The metal mines in Gongyi are mainly aluminium ore and iron ore. Referencing the spread range of heavy metal pollution in the soil of aluminium ore and iron ore mines48,49,50,51, the buffers are graded at 50, 100, 300, and 500 m (Fig. 7d2).The four new elements in the improved evaluation system are normalized from 0 to 1 during the calculation.Intelligent evaluation modelArtificial neural networks, decision trees, and SVMs were calculated using IBM SPSS Modeler software to find an intelligent model suitable for environmental assessment of the mine in the study area. Then, several models with high evaluation accuracy were selected. The SVM, CART, and C5.0 models were chosen for further comparison. The sampling points were selected randomly; 700 sampling points were selected from the area away from the mining area; 100 sampling points were selected from the mining area after random sampling by mine type, and these points were used as training samples. Non-mining evaluation scores were based on the national environmental protection standard, while the mining area scores were based on field investigation. In the field investigation, preliminary scoring of the sampling points was conducted according to mining type, mining intensity, air quality, and surrounding environment. Then, a photo of the field was taken at every sampling point in the mine, and experts were invited to further score the area according to the photo. This score is the relative score obtained by referencing the national environmental protection standard.The index layers of the training samples were used as input, and the scores were used as the output to train the machine learning models. The trained models were applied to the entire study area, and all points except the training sample points were used for verification. After further comparison with SVM, CART, and C5.0, the evaluation accuracy rates of the three methods in the mining area and non-mining area were obtained. In the non-mining area, the model evaluation results of various land use types were compared with the national environmental protection standards. The accuracy in various land use types is shown in Table 5. In the mining area, the model evaluation score is compared with the score from the experts, and the obtained accuracy table is shown in Table 6.Table 5 Accuracy of each algorithm in various land use types in non-mining areas.Full size tableTable 6 Accuracy of each algorithm in the mining area.Full size tableIn non-mining areas, the accuracy of the SVM model is significantly better than that of the other two methods. However, in the mining area, the accuracy of the CART model is higher. Therefore, the SVM model was used to evaluate the area away from the mine, and the CART model was used to evaluate the mining area. The evaluation results of these two models were combined to obtain the evaluation map of the entire study area. More

1.Pan, Y. et al. Science 333, 988–993 (2011).CAS 
Article 

Google Scholar 
2.Vanhaven, H. et al. (eds) Making Boreal Forests Work for People and Nature (IUFRO, 2012); https://www.iufro.org3.Stokland, J. N. For. Ecol. Manage. 488, 119017 (2021).Article 

Google Scholar 
4.Snäll et al. Nat. Sustain. https://doi.org/10.1038/s41893-021-00764-w (2021).5.Assessment Report on Land Degradation and Restoration (IPBES, 2018).6.Felipe-Lucia, M. R. et al. Nat. Commun. 9, 4839 (2018).Article 

Google Scholar 
7.Gamfeldt, L. et al. Nat. Commun. 4, 1340 (2013).Article 

Google Scholar 
8.Holland, R. A. et al. Biodivers. Conserv. 20, 3285–3294 (2011).Article 

Google Scholar 
9.National Forest Inventory: National Forest Monitoring (FAO, 2021); http://www.fao.org/national-forest-monitoring/areas-of-work/nfi/en/10.Simons, N. K. et al. For. Ecosyst. 8, 5 (2021).Article 

Google Scholar 
11.Sweden’s Forest Industry in Brief (FIS, 2021).12.Triviño, M. et al. J. Appl. Ecol. 54, 61–70 (2017).Article 

Google Scholar  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