Philippa Kaur

Anderson, M. The evolution of eusociality. Annu. Rev. Ecol. Syst. 15, 165–189 (1984).Article 

Google Scholar 
Wilkinson, G. S. Reciprocal food sharing in the vampire bat. Nature 308, 181–184 (1984).Article 

Google Scholar 
Feistner, A. & Mcgrew, W. Food-sharing in primates: a critical review. Perspect. Primate Biol 3, (1989).Hoelzel, A. R. Killer whale predation on marine mammals at Punta Norte, Argentina; food sharing, provisioning and foraging strategy. Behav. Ecol. Sociobiol. 29, 197–204 (1991).Article 

Google Scholar 
Behmer, S. T. Animal behaviour: feeding the superorganism. Curr. Biol. 19, R366–R368 (2009).Article 
CAS 
PubMed 

Google Scholar 
Cassill, D. L. & Tschinkel, W. R. Information flow during social feeding in ant societies. in Information Processing in Social Insects (eds. Detrain, C., Deneubourg, J. L. & Pasteels, J. M.) 69–81 (Birkhäuser, 1999). https://doi.org/10.1007/978-3-0348-8739-7_4.Hunt, J. H. Trophallaxis and the evolution of eusocial Hymenoptera. in The Biology of Social Insects (CRC Press, 1982).Sorensen, A. A., Busch, T. M. & Vinson, S. B. Trophallaxis by temporal subcastes in the fire ant, Solenopsis invicta, in response to honey. Physiol. Entomol. 10, 105–111 (1985).Article 

Google Scholar 
Meurville, M.-P. & LeBoeuf, A. C. Trophallaxis: the functions and evolution of social fluid exchange in ant colonies (Hymenoptera: Formicidae). Myrmecol. News 31, 1–30 (2021).Bodner, L. et al. Nutrient utilization during male maturation and protein digestion in the Oriental hornet. Biology 11, 241 (2022).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Sorensen, A. A., Kamas, R. S. & Vinson, S. B. The influence of oral secretions from larvae on levels of proteinases in colony members of Solenopsis invicta Buren (Hymenoptera: Formicidae). J. Insect Physiol. 29, 163–168 (1983).Article 
CAS 

Google Scholar 
Erthal, M., Peres Silva, C. & Ian Samuels, R. Digestive enzymes in larvae of the leaf cutting ant, Acromyrmex subterraneus (Hymenoptera: Formicidae: Attini). J. Insect Physiol. 53, 1101–1111 (2007).Article 
CAS 
PubMed 

Google Scholar 
Went, F. W., Wheeler, J. & Wheeler, G. C. Feeding and digestion in some ants (Veromessor and Manica). BioScience 22, 82–88 (1972).Article 

Google Scholar 
Ishay, J. & Ikan, R. Food exchange between adults and larvae in Vespa orientalis F. Anim. Behav. 16, 298–303 (1968).Article 
CAS 
PubMed 

Google Scholar 
Hunt, J. H. The evolution of social wasps. (Oxford University Press, USA, 2007).Abe, T., Yoshiya, T., Hiromitsu, M. & Kawasaki, Y. Y. Comparative study of the composition of hornet larval saliva, its effect on behaviour and role of trophallaxis. Comp. Biochem. Physiol. Part C: Comp. Pharmacol. 99, 79–84 (1991).Article 

Google Scholar 
Ishay, J. & Ikan, R. Gluconeogenesis in the Oriental hornet Vespa orientalis F. Ecology 49, 169–171 (1968).Article 

Google Scholar 
Brock, R. E., Cini, A. & Sumner, S. Ecosystem services provided by aculeate wasps. Biol. Rev. 96, 1645–1675 (2021).Article 
PubMed 

Google Scholar 
Ueno, T. Flower-visiting by the invasive hornet Vespa velutina nigrithorax (Hymenoptera: Vespidae). Int. J. Chem., Environ. Biol. Sci. 3, 444–448 (2015).
Google Scholar 
Käfer, H., Kovac, H. & Stabentheiner, A. Respiration patterns of resting wasps (Vespula sp.). J. Insect Physiol. 59, 475–486 (2013).Article 
PubMed 
PubMed Central 

Google Scholar 
Bodner, L., Bouchebti, S. & Levin, E. Allocation and metabolism of naturally occurring dietary amino acids in the Oriental hornet. Insect Biochem. Mol. Biol. 139, 103675 (2021).Article 
CAS 
PubMed 

Google Scholar 
Levin, E., Lopez-Martinez, G., Fane, B. & Davidowitz, G. Hawkmoths use nectar sugar to reduce oxidative damage from flight. Science 355, 733–735 (2017).Article 
CAS 
PubMed 

Google Scholar 
Hunt, J. H., Baker, I. & Baker, H. G. Similarity of amino acids in nectar and larval saliva: the nutritional basis for trophallaxis in social wasps. Evolution 36, 1318–1322 (1982).Article 
CAS 
PubMed 

Google Scholar 
Hunt, J. H., Jeanne, R. L., Baker, I. & Grogan, D. E. Nutrient dynamics of a swarm-founding social wasp species, Polybia occidentalis (Hymenoptera: Vespidae). Ethology 75, 291–305 (1987).Article 

Google Scholar 
Cassill, D. L. & Tschinkel, W. R. Allocation of liquid food to larvae via trophallaxis in colonies of the fire ant, Solenopsis invicta. Anim. Behav. 50, 801–813 (1995).Article 

Google Scholar 
Buffin, A., Denis, D., Simaeys, G. V., Goldman, S. & Deneubourg, J.-L. Feeding and stocking up: radio-labelled food reveals exchange patterns in ants. PLOS ONE 4, e5919 (2009).Article 
PubMed 
PubMed Central 

Google Scholar 
Quque, M. et al. Hierarchical networks of food exchange in the black garden ant Lasius niger. Insect Sci. 28, 825–838 (2021).Article 
PubMed 

Google Scholar 
Cassill, D. L. & Tschinkel, W. R. A duration constant for worker-to-larva trophallaxis in fire ants. Ins. Soc. 43, 149–166 (1996).Article 

Google Scholar 
Cassill, D. L. & Tschinkel, W. R. Regulation of diet in the fire ant, Solenopsis invicta. J. Insect Behav. 12, 307–328 (1999).Article 

Google Scholar 
Wilson, E. O. & Eisner, T. Quantitative studies of liquid food transmission in ants. Ins. Soc. 4, 157–166 (1957).Article 

Google Scholar 
Markin, G. P. Food distribution within laboratory colonies of the argentine ant, Tridomyrmex humilis (Mayr). Ins. Soc. 17, 127–157 (1970).Article 

Google Scholar 
Howard, D. F. & Tschinkel, W. R. The flow of food in colonies of the fire ant, Solenopsis invicta: a multifactorial study. Physiol. Entomol. 6, 297–306 (1981).Article 

Google Scholar 
Suryanarayanan, S. & Jeanne, R. L. Antennal drumming, trophallaxis, and colony development in the social wasp Polistes fuscatus (Hymenoptera: Vespidae). Ethology 114, 1201–1209 (2008).Article 

Google Scholar 
Greenwald, E., Segre, E. & Feinerman, O. Ant trophallactic networks: simultaneous measurement of interaction patterns and food dissemination. Sci. Rep. 5, 12496 (2015).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Baltiansky, L., Sarafian-Tamam, E., Greenwald, E. & Feinerman, O. Dual-fluorescence imaging and automated trophallaxis detection for studying multi-nutrient regulation in superorganisms. Methods Ecol. Evol. 12, 1441–1457 (2021).Article 

Google Scholar 
Feldhaar, H. et al. Stable isotopes: past and future in exposing secrets of ant nutrition (Hymenoptera: Formicidae). Myrmecol. N. 13, 3–13 (2010).
Google Scholar 
Bouchebti, S., Bodner, L., Bergman, M., Magory Cohen, T. & Levin, E. The effects of dietary proline, β-alanine, and γ-aminobutyric acid (GABA) on the nest construction behavior in the Oriental hornet (Vespa orientalis). Sci. Rep. 12, 7449 (2022).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Motro, M., Motro, U., Ishay, J. S. & Kugler, J. Some social and dietary prerequisites of oocyte development in Vespa orientalis L. workers. Ins. Soc. 26, 155–164 (1979).Article 

Google Scholar 
Levin, E., McCue, M. D. & Davidowitz, G. More than just sugar: allocation of nectar amino acids and fatty acids in a Lepidopteran. Proc. R. Soc. B: Biol. Sci. 284, 20162126 (2017).Article 

Google Scholar 
Wright, G. A., Nicolson, S. W. & Shafir, S. Nutritional physiology and ecology of honey bees. Annu. Rev. Entomol. 63, 327–344 (2018).Article 
CAS 
PubMed 

Google Scholar 
Helm, B. R. et al. The geometric framework for nutrition reveals interactions between protein and carbohydrate during larval growth in honey bees. Biol. Open 6, 872–880 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
Paoli, P. P. et al. Nutritional balance of essential amino acids and carbohydrates of the adult worker honeybee depends on age. Amino Acids 46, 1449–1458 (2014).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Stabler, D., Paoli, P. P., Nicolson, S. W. & Wright, G. A. Nutrient balancing of the adult worker bumblebee (Bombus terrestris) depends on the dietary source of essential amino acids. J. Exp. Biol. 218, 793–802 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Arganda, S. et al. Parsing the life-shortening effects of dietary protein: effects of individual amino acids. Proc. R. Soc. B 284, 20162052 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Csata, E. & Dussutour, A. Nutrient regulation in ants (Hymenoptera: Formicidae): a review. Myrmecol. N. 29, 111–124 (2019).
Google Scholar 
Gottsberger, G., Schrauwen, J. & Linskens, H. F. Amino acids and sugars in nectar, and their putative evolutionary significance. Pl. Syst. Evol. 145, 55–77 (1984).Article 
CAS 

Google Scholar 
Ozimek, L. et al. Nutritive value of protein extracted from honey bees. J. Food Sci. 50, 1327–1329 (1985).Article 
CAS 

Google Scholar 
Nicolson, S. W. & Thornburg, R. W. Nectar chemistry. in Nectaries and Nectar (eds. Nicolson, S. W., Nepi, M. & Pacini, E.) 215–264 (Springer Netherlands, 2007). https://doi.org/10.1007/978-1-4020-5937-7_5.Contrera, F. A. L., Imperatriz-Fonseca, V. L. & Koedam, D. Trophallaxis and reproductive conflicts in social bees. Insect Soc. 57, 125–132 (2010).Article 

Google Scholar 
Carter, G. G. & Wilkinson, G. S. Food sharing in vampire bats: reciprocal help predicts donations more than relatedness or harassment. Proc. R. Soc. B: Biol. Sci. 280, 20122573 (2013).Article 

Google Scholar 
Nalepa, C. A. Origin of termite eusociality: trophallaxis integrates the social, nutritional, and microbial environments. Ecol. Entomol. 40, 323–335 (2015).Article 

Google Scholar 
Werenkraut, V., Arbetman, M. P. & Fergnani, P. N. The Oriental hornet (Vespa orientalis L.): a threat to the Americas? Neotrop. Entomol. 51, 330–338 (2022).Article 
PubMed 

Google Scholar 
Darchen, R. Biologie de Vespa orientalis. Les premiers stades de développement. Ins. Soc. 11, 141–157 (1964).Article 

Google Scholar 
Van Itterbeeck, J. et al. Rearing techniques for hornets with emphasis on Vespa velutina (Hymenoptera: Vespidae): A review. J. Asia-Pac. Entomol. 24, 103–117 (2021).Article 

Google Scholar 
R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. (2020).Bouchebti, S., Bodner,L. & Levin, E. Continuous exchange of nectar nutrients in an Oriental hornet colony- Dataset [Data set]. Zenodo. https://doi.org/10.5281/zenodo.7135100 (2022). More

Data collectionIn phase 1 (2017), we first circulated a comprehensive multi-language questionnaire and associated guidelines on protected area personnel numbers to major national protected area agencies, focusing on the 50 countries listed in the WDPA as having the most protected areas. The questionnaire requested information on personnel numbers, type of employers and management levels (from executive to skilled practical workers). Protected area personnel were defined as those spending at least 50% of their work time on protected area-related tasks. The questionnaire also requested information about job titles used for personnel equivalent to rangers. This phase produced usable data for 28 countries/territories.In phase 2 (2018 onwards), we conducted online searches for published data on protected area personnel numbers in the countries/territories not included in the questionnaire survey or where questionnaire responses were incomplete or unclear. The resulting information came from official organizational reports (10 countries/territories), published external studies, project documents and journal papers (35 countries/territories) and websites of protected area organizations or individual sites (9 countries/territories).In phase 3 (2018–2021), we directly requested personal contacts to locate or supply information from official sources both for the remaining countries/territories and to improve or verify data from phases 1 and 2. The minimum data requested were the overall number of protected area personnel, the number of those personnel that could be categorized as rangers, the terrestrial area of protected areas managed by the listed personnel and the source of the information. This phase contributed usable data for 68 countries and territories. Data for a further 17 countries/territories were assembled from multiple sources.The final dataset covered 176 countries/territories: 167 surveyed countries/territories and a further 9 countries/territories that have no WDPA-listed protected areas (Supplementary Table 1), with contributions from more than 150 individuals.Initial data processingTo assess and, where necessary, improve the reliability of data obtained in a wide range of formats and levels of detail and from multiple sources, we scored the data for each country/territory from 0 to 5 for each of four criteria—detail, accuracy, source and age of the data—with a maximum score of 20 (Supplementary Table 1 and Supplementary Fig. 1). For all low-scoring records (a score of less than 15), we sought more-reliable sources in later phases of the study, rejecting any final scores of less than 10.On reviewing the data, we excluded from the analysis protected areas identified in the WDPA as predominantly or entirely marine, Antarctica and countries/territories categorized in the WDPA as polar (Greenland, French Southern Territories, Bouvet Island, Heard Island and McDonald Islands, South Georgia and the South Sandwich Islands). These large, remote and/or largely uninhabited areas are likely to have quite different management models and scales of staffing from terrestrial protected areas (although marine protected areas are also widely understaffed11). For example, in 2012 the 972,000 km2 of Northeast Greenland Protected Area (categorized by the WDPA as polar) was only periodically visited by six two-person teams of naval personnel47, and the 2008 management plan of the 1.51 million km2 Papahānaumokuākea Marine National Monument (Hawai’i, USA) specifies just nine personnel, working in conjunction with several other agencies48. Data for one country were supplied by officials on the agreement that the country was not specifically identified in publications (the country is given the three-letter code ZZZ in relevant tables and figures).Because the format, completeness and level of detail of the data varied widely, from comprehensive personnel lists to single figures, we restricted our raw dataset to six variables that could be consistently extracted from data obtained for each country/territory:

1.

Total number of non-ranger personnel (if known)

2.

Total number of rangers (if known)

3.

Total number of protected area personnel (either the sum of 1 and 2 or provided as an undifferentiated total)

4.

Terrestrial area of protected areas covered by surveyed personnel (km2)

5.

Total terrestrial area of protected areas of the country/territory (km2)

6.

Year of the data

We used the WDPA, official publications and websites to determine (or verify) the area of terrestrial protected areas covered by the personnel listed for each country/territory, using WDPA data if there were discrepancies. Total national terrestrial protected area coverage was taken from the WDPA, with the exception of Turkey, where the area officially reported to the WDPA is significantly less than the nationally published area.The raw data from the survey are shown in Supplementary Table 1.Candidate predictorsTo predict the number of rangers and non-rangers in countries and territories for which we had no data (Statistical analysis), we collected information on the following set of variables, hereafter referred to as candidate predictors:Location dataThe WGS84 latitude and longitude of the centroid of the largest land mass associated with each country/ territory (to obtain the polygons defining the land masses, we used the R package rnaturalearth version 0.1.0; https://github.com/ropensci/rnaturalearth)2020 data from the World Bank (https://data.worldbank.org/indicator)

Area of the country/territory

Population density: the mid-year population divided by land area

Gross domestic product (GDP) in US dollars

GDP per capita in US dollars (GDP divided by mid-year population)

Growth rate of GDP

The proportion of rural inhabitants

The proportion of unemployed inhabitants

The forested proportion of the country/territory

2020 data for each country/territory from the WDPA (https://www.protectedplanet.net/)

The total terrestrial area of WDPA-listed protected areas

The proportion of the terrestrial area of all IUCN-categorized protected areas (Categories I–VI) that falls within protected areas in Category I or II

The proportion of the terrestrial area of all IUCN-categorized protected areas (Categories I–VI) that falls within protected areas in Categories I–IV

2020 data from the Yale Center for Environmental Law and Policy Environmental Performance Index (https://epi.yale.edu/)

Environmental Performance Index (EPI): a composite index using 32 performance indicators across 11 categories

Ecosystem Vitality Index (EVI): an indicator of how well countries preserve, protect and enhance ecosystems and the services they provide

Species Protection Index (SPI): an indicator of the species-level ecological representativeness of each country’s/territory’s protected area network

Not all this information was available for all countries/territories. Most of the missing data were for small territories that account for only a very small proportion of the total area of protected areas worldwide (Supplementary Table 2c).Statistical analysisOur primary objective was to estimate the total number of all personnel engaged in managing all the world’s WDPA-listed terrestrial protected areas and the number categorized as rangers. Our raw data collection yielded full, partial or no information on total personnel and ranger numbers for each country/territory (Supplementary Table 1 shows the completeness of all the data collected). Our first task, therefore, was (1) to impute the information for unsurveyed protected areas on the basis of information from surveyed protected areas within the same countries/territories and (2) to predict those numbers for countries/territories where no information was available on overall personnel numbers and/or ranger numbers on the basis of relationships we could establish between available information and candidate predictors in other countries/territories (Supplementary Table 7). A brief description of these two approaches follows, and full details on the analysis are provided in Supplementary Information.Data imputationFor countries/territories where we had obtained information about numbers of personnel and/or rangers for only some protected areas, our strategy was to populate the unsurveyed protected areas in proportion to the densities of personnel or rangers from the surveyed protected areas of the same countries/territories. For example, for Spain we obtained evidence that there are 619 rangers responsible for protected areas covering 44,328 km2, out of a national total protected area system covering 142,573 km2. To impute the number of rangers for the remaining 98,245 km2, we used the density of rangers in the surveyed area (one ranger per 44,328/619 = 71.6 km2) and applied that to the unsurveyed area, giving a total of 1,991 rangers (619 + (98,245/71.6)). This imputation assumes that unsurveyed areas are staffed at the same density as surveyed areas, whereas in reality the relative densities are likely to vary in unknown ways within different countries/territories. To study the sensitivity of our results to the assumed proportion, we repeated our analysis using the following proportions of the observed densities: 0, 0.25, 0.50, 0.75 and 1.00. This provided a range of personnel numbers from a minimum (based on a proportion of 0) to a presumed maximum (based on a proportion of 1.00). From the data obtained, it was not possible to calculate the actual proportions, but based on the experience of the practitioners in the author team, the unsurveyed areas are highly unlikely to be staffed at higher densities than surveyed areas and, on average, are very likely to be staffed at lower densities. After all, most survey respondents were national or subnational agencies responsible for protected areas subject to stronger formal requirements for protection and management and therefore likely to have larger workforces. Unsurveyed protected areas are more likely to be managed by local entities, with fewer resources, less-stringent management obligations and therefore fewer personnel. The range of proportions we considered to populate unsurveyed areas should therefore yield predictions encompassing the actual (unknown) numbers of rangers and non-rangers with a conservative margin of error. In the main text, we have reported the results of imputation assuming a proportion of 1, which is probably the most optimistic assessment of the current workforce in protected areas within the proportions of the observed densities considered. Results using lower proportions are shown in Extended Data Fig. 2 and Supplementary Tables 4 and 5.Data predictionOur imputation approach was not possible for countries/territories where (1) zero ranger or personnel data had been obtained and (2) specific data had not been obtained that allowed imputation either for rangers or for total personnel (where only total personnel numbers or only ranger numbers had been obtained). To predict the missing information, we used two different statistical approaches: linear mixed models (LMMs)49 and a general implementation of random forests, which we term RF/ETs because it encompasses both random forests sensu stricto (RFs)50 and a variant called extremely randomized trees (ETs)51. LMMs and RFs have been extensively discussed and reviewed in the literature49,52,53. We adopted these approaches because both have proved successful in producing accurate predictions for a wide range of applications and because both are well suited to our data since they both produce predictions from a set of predictors and allow for the consideration of spatial effects54,55. Furthermore, comparing predictions generated through very different methods informs us about the robustness of our results with respect to key statistical assumptions. LMMs come from the ‘data modelling culture’56 and belong to parametric statistics; RF/ETs come from the ‘algorithmic modelling culture’ and belong to non-parametric statistics.We followed the same workflow for both statistical approaches, comprising eight steps: (1) general data preparation; (2) preparation of initial training datasets; (3) selection of predictor variables and of the method used for handling spatial autocorrelation; (4) preparation of final training datasets; (5) fine tuning; (6) final training; (7) preparation of datasets for predictions and simulations; and (8) predictions and simulations (see Supplementary Information for details).Both approaches yielded very similar results with our data. We chose to present the LMM results in the main text, but we provide and compare the results obtained by both approaches in Supplementary Information.SoftwareWe performed all the data analyses using the free open-source statistical software R version 4.157. We used the R package spaMM version 3.9.13 to implement LMMs58 and the R package ranger version 0.13.1 to implement RF/ETs59. To reformat and plot the data, we used the Tidyverse suite of packages60. Details are provided in an R package we specifically developed so that findings presented in this paper can readily be reproduced (see Code availability). Using a workstation with an AMD Ryzen Threadripper 3990 × 64-core processor and 256 GB of RAM, our complete workflow ran in ~3,000 CPU hours.Estimation of required numbers and densities of personnelTo estimate the numbers of personnel and rangers required for effective management of existing protected areas, we referred to ref. 25. This estimates that the minimum budget needed to adequately manage the existing protected area system is US$67.6 billion per year and that current annual expenditure is US$24.3 billion. From these figures, we can calculate that resources invested in the current global system of protected areas are approximately 36% of what is required. We consulted data from https://ourworldindata.org to determine that the proportion of global public expenditure on employee compensation has remained between 21.01% and 23.33% in the years from 2006 to 2019. We obtained these figures from the ‘Government Spending’ section of the site, consulting the chart ‘Share of employee compensation in public spending, 2002 to 2019’ and selecting data for ‘World’. On the basis of this broadly constant proportion and the assumption that total employee compensation is an indicator of total employee numbers, we inferred that current numbers of protected area employees are also around 36% of what is required. We therefore multiplied our estimations of personnel and ranger numbers by 1/0.36 and recalculated the densities on this basis (current requirement = 1/0.36 × current estimate).To estimate staffing requirements for 30% global coverage of protected areas—the global target intended to be reached by 2030—we used the mean personnel and ranger densities calculated as being required at present to ‘populate’ a global area of terrestrial protected areas if increased from the percentage at the time of our study (15.7%) to 30% (current requirement × (0.300/0.157)).Economic calculationsWe based our calculations on published data from 202025, which estimate that expanding the protected areas to 30% would generate higher overall output (revenues) than non-expansion (an extra US$64–454 billion per year by 2050). This figure is only an indicative, partial estimate, generated for the purposes of comparison and to illustrate the substantial return on investment that protected area staff investments imply. Using these figures and our estimates of personnel requirements to ensure effective management of 30% coverage, we calculated the range of sums that each additional protected area staff member has the potential to generate (Supplementary Table 8). For clarity, we rounded these figures to the nearest hundred US dollars in the main text.Our estimates of the gross value added per worker in forestry and agriculture (sectors responsible for similar proportions of the world as protected areas) are included to provide a point of comparison for the figures showing the economic benefit generated per protected area personnel member (see the preceding). The data for the gross annual value of world agricultural production (US$3,550,231,736,000) and the number of workers employed in agriculture (343,527,711) come from the Food and Agriculture Organization of the United Nations30, providing an average gross value of annual agricultural production per worker of US$10,335. We adjusted these 2018 data to 2020 price levels using a deflator based on the US consumer price index (CPI) from the World Economic Outlook database61 (Supplementary Table 9). This ensures that all the economic value data we present are directly comparable for protected area, agricultural and forestry workers. We calculated the gross value of forest production per worker on the basis of direct contribution of forestry of more than US$539 billion to world GDP in 201162 and total forest-sector employment of 11.881 million full-time-equivalent jobs in 201032. These were the most up-to-date global estimates we could locate from credible sources that presented comparable estimates of forest-sector employment and contribution to GDP. This gives an average gross value of forest production per worker of US$45,367 per year. We used the same method as for agriculture to bring these figures to 2020 price levels (Supplementary Table 9). These figures are rounded to the nearest hundred US dollars in the main text. More

Manjeri, G., Muhamad, R. & Tan, S. G. Oryctes rhinoceros beetles, an oil palm pest in Malaysia. Annu. Res. Rev. Biol. 4, 3429–3439 (2014).Article 

Google Scholar 
Allou, K., Morin, J. P., Kouassi, P., Nklo, F. H. & Rochat, D. Oryctes monoceros trapping with synthetic pheromone and palm material in Ivory Coast. J. Chem. Ecol. 32, 1743–1754 (2006).Article 
CAS 
PubMed 

Google Scholar 
Alibert, H. Study on the insect pests of oil palm in Dahomey. Rev. Botan. Appl. 18, 745–773 (1936).
Google Scholar 
Catley, A. The coconut rhinoceros beetle Oryctes rhinoceros (L) [Coleoptera: Scarabaeidae: Dynastinae]. PANS Pest Articles News Summar. 15, 18–30 (1969).Article 

Google Scholar 
Fauzana, H., Sutikno, A. & Salbiah, D. Population fluctuations Oryctes rhinoceros L. beetle in plant oil palm (Elaeis guineensis Jacq.) given mulching oil palm empty bunch. Cropsaver Int. J. Trop. Insect Sci. 1, 42–47 (2018).
Google Scholar 
Paudel, S., Mansfield, S., Villamizar, L. F., Jackson, T. A. & Marshall, S. D. Can biological control overcome the threat from newly invasive coconut rhinoceros beetle populations (Coleoptera: Scarabaeidae)? A review. Ann. Entomol. Soc. Am. 114, 247–256 (2021).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Molet, T. In CPHST Pest Datasheet for Oryctes rhinoceros. USDA-APHIS-PPQCPHST. Revised July 2014 (2013).Hinckley, A. D. Ecology of the coconut rhinoceros beetle, Oryctes rhinoceros (L.) (Coleoptera: Dynastidae). Biotropica 1973, 111–116 (1973).Article 

Google Scholar 
Sitepu, D., Kharie, S., Waroka, JS & Motulo, HFJ. Methods for the production and use of Marhizium anisopliae against Oryctes rhinoceros. In Integrated Coconut Pest Control Project—Annual report of Coconut Research Institute—Manado, North Sulawesi, Indonesia 104–111 (1988).Philippe, R. & Dery, S. K. Coconut research and development. CORD 20, 43–51 (2004).
Google Scholar 
Purrini, K. Baculovirus oryctes release into Oryctes monoceros population in Tanzania, with special reference to the interaction of virus isolates used in our laboratory infection experiments. J. Invertebr. Pathol. 53, 285–300 (1989).Article 

Google Scholar 
Ukeh, D. A., Usua, E. J. & Umoetok, S. B. A. Notes on the biology of Oryctes monoceros (OLIV.) A pest of palms in Nigeria. World J. Agric. Res. 2, 33–36 (2003).
Google Scholar 
Dry, F. W. Notes on the coconut beetle (Oryctes monoceros, Ol.) in Kenya Colony. Bull. Entomol. Res. 13, 103–107 (1922).Article 

Google Scholar 
Bedford, G. O. Biology, ecology, and control of palm rhinoceros beetles. Annu. Rev. Entomol. 25, 309–339 (1980).Article 

Google Scholar 
Khoo, K. C., Yusoff, M. N. M. & Lee, T. W. Pulp and paper of oil palm trunk. In Research Pamphlet No.107: Oil Palm Stem Utilisation, Kuala Lumpur, Malaysia, FRIM 51–65 (1991).Giblin-Davis, R. M. Borers of palms. In Insects on Palms (eds Moore, D. et al.) (CABI Publishing, Wallingford, 2001).
Google Scholar 
Drumoni, A. & Ponchel, Y. Première capture au Yémen d’ Oryctes (Rykanoryctes) monoceros (Olivier, 1789) et confirmation de la présence de cette espèce africaine dans la Péninsule Arabique (Coleoptera, Dynastidae). Entomol. Afr. 15, 25–29 (2010).
Google Scholar 
Lever, R. J. A. W. Pests of the Coconut Palm (Food and Agriculture Organization of the United Nations, Rome, 1969).Moore, A. Rhinoceros beetle pest found in Guam and Saipan. In Pest Alert, Suva, Fiji: Plant Protection Service, Secretariat of the Pacific Community (2007).Zhang, K., Yao, L., Meng, J. & Tao, J. Maxent modeling for predicting the potential geographical distribution of two peony species under climate change. Sci. Total Environ. Sci. 634, 1326–1334 (2018).Article 
ADS 
CAS 

Google Scholar 
Ding, F., Fu, J., Jiang, D., Hao, M. & Lin, G. Mapping the spatial distribution of Aedes aegypti and Aedes albopictus. Acta Trop. 178, 155–162 (2018).Article 
PubMed 

Google Scholar 
Valencia-Rodríguez, D., Jiménez-Segura, L., Rogéliz, C. A. & Parra, J. L. Ecological niche modeling as an effective tool to predict the distribution of freshwater organisms: The case of the Sabaleta Brycon henni (Eigenmann, 1913). PLoS ONE 16, e0247876 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Escobar, L. E., Qiao, H., Cabello, J. & Peterson, A. T. Ecological niche modeling re-examined: A case study with the Darwin’s fox. Ecol. Evol. 8, 4757–4770 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Warren, D. L. & Seifert, S. N. Ecological niche modeling in Maxent: The importance of model complexity and the performance of model selection criteria. Ecol. Appl. 21, 335–342 (2011).Article 
PubMed 

Google Scholar 
Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259 (2006).Article 

Google Scholar 
Phillips, S. J. Transferability, sample selection bias and background data in presence-only modelling: A response to Peterson et al. (2007). Ecography 31, 272–278 (2008).Article 

Google Scholar 
Elith, J. et al. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 17, 43–57 (2011).Article 

Google Scholar 
Phillips, S. J. & Dudík, M. Modeling of species distributions with MaxEnt: New extensions and a comprehensive evaluation. Ecography 31, 161–175 (2008).Article 

Google Scholar 
Arnold, J. D., Brewer, S. C. & Dennison, P. E. Modeling climate-fire connections within the Great basin and Upper Colorado River Basin. Fire Ecol. 10, 64–75 (2014).Article 

Google Scholar 
Phillips, J. S. & Elith, J. On estimating probability of presence from use-availability or presence-background data. Ecology 94, 1409–1419 (2013).Article 
PubMed 

Google Scholar 
Santana, P. A. Jr., Kumar, L., Da Silva, R. S., Pereira, J. L. & Picanço, M. C. Assessing the impact of climate change on the worldwide distribution of Dalbulus maidis (DeLong) using MaxEnt. Pest. Manag. Sci. 75, 2706–2715 (2019).Article 
CAS 
PubMed 

Google Scholar 
Li, et al. Predicting the current and future distributions of Brontispa longissima (Coleoptera: Chrysomelidae) under climate change in China. Glob. Ecol. Conserv. 25, e01444 (2021).Article 

Google Scholar 
Li, T. et al. Direct and indirect effects of environmental factors, spatial constraints, and functional traits on shaping the plant diversity of montane forests. Ecol. Evol. 10, 557–568 (2020).Article 
PubMed 

Google Scholar 
Namgung, H., Kim, M. J., Baek, S., Lee, J. H. & Kim, H. Predicting potential current distribution of Lycorma delicatula (Hemiptera: Fulgoridae) using MaxEnt model in South Korea. J. Asia Pac. Entomol. 23, 291–297 (2020).Article 

Google Scholar 
Ji, W., Gao, G. & Wei, J. Potential global distribution of Daktulosphaira vitifoliae under climate change based on MaxEnt. Insects. 12, 347 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Ji, W., Han, K., Lu, Y. & Wei, J. Predicting the potential distribution of the vine mealybug, Planococcus ficus under climate change by MaxEnt. J. Crop. Prot. 137, 105268 (2020).Article 

Google Scholar 
Sharma, HC & Prabhakar, CS. Impact of climate change on pest management and food security. In Integrated Pest Management 23–36 (Academic Press, Cambridge, 2014).Skendžić, S., Zovko, M., Živković, I. P., Lešic, V. & Lemić, D. The impact of climate change on agricultural insect pests. Insects. 12, 440 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Ward, N. L. & Masters, G. J. Linking climate change and species invasion: An illustration using insect herbivores. Glob. Change Biol. 13, 1605–1615 (2007).Article 
ADS 

Google Scholar 
De Queiroz, D. L., Burckhardt, D. & Majer, J. Integrated pest management of eucalypt psyllids (Insecta, Hemiptera, Psylloidea). In Integrated pest management and pest control-current and future tactics. INTECH 2012, 385–412 (2012).
Google Scholar 
Hochberg, M. E. & Waage, J. K. A model for the biological control of Oryctes rhinoceros (Coleoptera: Scarabaeidae) by means of pathogens. J. Appl. Ecol. 28, 514–531 (1991).Article 

Google Scholar 
Liu, Y. et al. MaxEnt modelling for predicting the potential distribution of a near threatened rosewood species (Dalbergia cultrata Graham ex Benth). Ecol. Eng. 141, 105612 (2019).Article 

Google Scholar 
Wang, R. et al. Predictions of potential geographical distribution of Diaphorina citri (Kuwayama) in China under climate change scenarios. Sci. Rep. 10, 1–9 (2020).CAS 

Google Scholar 
Wood, B. J. Studies on the effect of ground vegetation on infestations of Oryctes rhinoceros (L.) (Col., Dynastidae) in young oil palm replantings in Malaysia. Bull Entomol. Res. 59, 85–96 (1969).Article 

Google Scholar 
Mittal, I. C. Survey of scarabaeid (Coleoptera) fauna of Himachal Pradesh (India). J. Entomol. Res. 24, 259–269 (2000).
Google Scholar 
Zheng, C., Jiang, D., Ding, F., Fu, J. & Hao, M. Spatiotemporal patterns and risk factors for scrub typhus from 2007 to 2017 in southern China. Clin. Infect. Dis. 69, 1205–1211 (2019).Article 
PubMed 

Google Scholar 
Chen, S., Ding, F., Hao, M. & Jiang, D. Mapping the potential global distribution of red imported fire ant (Solenopsis invicta Buren) based on a machine learning method. Sustainability. 12, 10182 (2020).Article 

Google Scholar 
Ding, F. et al. Infection and risk factors of human and avian influenza in pigs in south China. Prev. Vet. Med. 190, 105317 (2021).Article 
PubMed 

Google Scholar 
Jiang, D. et al. Spatiotemporal patterns and spatial risk factors for Visceral leishmaniasis from 2007 to 2017 in Western and Central China: A modelling analysis. Sci. Total Environ Sci. 764, 144275 (2021).Article 
ADS 
CAS 

Google Scholar 
Méndez-Rojas, D. M., Cultid-Medina, C. & Escobar, F. Influence of land use change on rove beetle diversity: A systematic review and global meta-analysis of a mega-diverse insect group. Ecol. Indic. 122, 107239 (2021).Article 

Google Scholar 
Oke, T. R. City size and the urban heat island. Atmos. Environ. 7, 769–779 (1973).Article 
ADS 

Google Scholar 
Briere, J. F., Pracros, P., Le Roux, A. Y. & Pierre, J. S. A novel rate model of temperature-dependent development for arthropods. Environ. Entomol. 28, 22–29 (1999).Article 

Google Scholar 
Zeng, Y., Low, B. W. & Yeo, D. C. Novel methods to select environmental variables in MaxEnt: A case study using invasive crayfish. Eco. Model. 341, 5–13 (2016).Article 

Google Scholar 
Fand, B. B. et al. Invasion risk of the South American tomato pinworm Tuta absoluta (Meyrick) (Lepidoptera: Gelechiidae) in India: Predictions based on MaxEnt ecological niche modelling. Int. J. Trop. Insect Sci. 40, 1–11 (2020).Article 

Google Scholar 
Li, W. J. et al. Potential distribution prediction of natural Pseudotsuga sinensis forest in Guizhou based on Maxent model. J. For. Res. 48, 47–52 (2019).
Google Scholar 
McIntyre, S., Rangel, E. F., Ready, P. D. & Carvalho, B. M. Species-specific ecological niche modelling predicts different range contractions for Lutzomyia intermedia and a related vector of Leishmania braziliensis following climate change in South America. Parasit. Vectors 10, 1–15 (2017).Article 

Google Scholar 
Hao, M. et al. Global potential distribution of Oryctes rhinoceros, as predicted by boosted regression tree model. Glob. Ecol. Conserv. 37, e02175 (2022).Article 

Google Scholar 
Aidoo, O. F. et al. The impact of climate change on potential invasion risk of Oryctes monoceros worldwide. Front. Ecol. Evol. 10, 633 (2022).Article 

Google Scholar 
Aidoo, O. F. et al. Lethal yellowing disease: Insights from predicting potential distribution under different climate change scenarios. J. Plant Dis. Prot. 2021, 1–13 (2021).
Google Scholar 
Ruheili, A. M. A., Boluwade, A. & Subhi, A. M. A. Assessing the Impact of Climate Change on the Distribution of Lime (16srii-B) and Alfalfa (16srii-D) Phytoplasma Disease Using MaxEnt. Plants. 10, 460 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Wang, R. et al. Predicting the potential distribution of the Asian citrus psyllid, Diaphorina citri (Kuwayama), in China using the MaxEnt model. PeerJ 7, e7323 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
He, S. T. & Jing, P. F. Prediction of potential distribution areas of Salvia bowleyana Dunn. in China based on MaxEnt and suitability analysis. J Anhui Agri. Sci. 8, 2311–2314 (2014).
Google Scholar 
Chahouki, M. A. Z. & Sahragard, H. P. Maxent modelling for distribution of plant species habitats of rangelands (Iran). Pol. J. Ecol. 64, 453–467 (2016).
Google Scholar 
Shabani, F., Kumar, L. & Ahmadi, M. Assessing accuracy methods of species distribution models: AUC, specificity, sensitivity and the true skill statistic. Glob. Int. J. Hum. Soc. Sci. 18, 6–18 (2018).
Google Scholar 
Baloch, M. N., Fan, J., Haseeb, M. & Zhang, R. Mapping potential distribution of Spodoptera frugiperda (Lepidoptera: Noctuidae) in central Asia. Insects. 11, 172 (2020).Article 
PubMed Central 

Google Scholar 
Wang, N., Li, Z., Wu, J., Rajotte, E. G., Wan, F & Wang, Z. The potential geographical distribution of Bactrocera dorsalis (Diptera: Tephrididae) in China based on emergence rate model and ArcGIS. In International Conference on Computer and Computing Technologies in Agriculture 399–411. (Springer, Boston, 2008).Manrique, V., Cuda, J. P., Overholt, W. A. & Diaz, R. Temperature-dependent development and potential distribution of Episimus utilis (Lepidoptera: Tortricidae), a candidate biological control agent of Brazilian peppertree (Sapindales: Anacardiaceae) in Florida. Environ. Entomol. 37, 862–870 (2008).Article 
PubMed 

Google Scholar 
Das, D. K., Singh, J. & Vennila, S. Emerging crop pest scenario under the impact of climate change–a brief review. AgroPhysics. 11, 13–20 (2011).CAS 

Google Scholar 
Porter, J. H., Parry, M. L. & Carter, T. R. The potential effects of climatic change on agricultural insect pests. Agric. For. Meteorol. 57, 221–240 (1991).Article 
ADS 

Google Scholar 
Trenberth, K. E. Climate change caused by human activities is happening and it already has major consequences. J. Energy Nat. Resour. Law. 36, 463–481 (2018).Article 

Google Scholar 
Xu, D., Zhuo, Z., Li, X. & Wang, R. Distribution and invasion risk assessment of Oryctes rhinoceros (L.) in China under changing climate. J. Appl. Entomol. 146, 385–395 (2022).Article 

Google Scholar 
Sushil, K. & Mukhtar, A. Effect of temperature and humidity on biology of rhinoceros beetle, Oryctes rhinoceros Linn. on oil palm. J. Appl. Anim. Res. 18, 108–112 (2007).
Google Scholar 
Sabidin, N. N. E. The effect of climate change to the population of rhinoceros beetle (Oryctes rhinoceros) at selected oil palm plantation. In Bachelor of Science Thesis Dissertation. Universiti Teknologi MARA. https://ir.uitm.edu.my/id/eprint/22754. (2018).Yadav, R. & Chang, N. T. Effects of temperature on the development and population growth of the melon thrips, Thrips palmi, on eggplant, Solanum melongena. J. Insect Sci. 14, 78 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Ju, R. T., Wang, F. & Li, B. Effects of temperature on the development and population growth of the sycamore lace bug, Corythucha ciliata. J. Insect Sci. 11, 1–12 (2011).Article 

Google Scholar 
Zheng, F. S., Du, Y. Z., Wang, Z. J. & Xu, J. J. Effect of temperature on the demography of Galerucella birmanica (Coleoptera: Chrysomelidae). Insect Sci. 15, 375–380 (2008).Article 

Google Scholar 
Azrag, A. G. et al. Modelling the effect of temperature on the biology and demographic parameters of the African coffee white stem borer, Monochamus leuconotus (Pascoe) (Coleoptera: Cerambycidae). J. Therm. Biol. 89, 102534 (2020).Article 
CAS 
PubMed 

Google Scholar 
Aidoo, O. F. et al. The African citrus triozid Trioza erytreae Del Guercio (Hemiptera: Triozidae): Temporal dynamics and susceptibility to entomopathogenic fungi in East Africa. Int. J. Trop. Insect Sci. 41, 563–573 (2021).Article 

Google Scholar 
Leonard, A. et al. Predicting the current and future distribution of the edible long-horned grasshopper Ruspolia differens (Serville) using temperature-dependent phenology models. J. Therm. Biol. 95, 102786 (2021).Article 
PubMed 

Google Scholar 
Roy, B. A. et al. Increasing forest loss worldwide from invasive pests requires new trade regulations. Front. Ecol. Environ. 12, 457–465 (2014).Article 

Google Scholar 
Shabani, F., Kumar, L. & Ahmadi, M. A comparison of absolute performance of different correlative and mechanistic species distribution models in an independent area. Ecol. Evol. 6, 5973–5986 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Cianci, D., Hartemink, N. & Ibáñez-Justicia, A. Modelling the potential spatial distribution of mosquito species using three different techniques. Int. J. Health Geogr. 14, 10 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Zelazny, B. & Alfiler, A. Oryctes rhinoceros (Coleoptera: Scarabaeidae) larva abundance and mortality factors in the Philippines. Environ. Entomol. 15, 84–87 (1986).Article 

Google Scholar 
Wood, B.J. Studies on the effect of ground vegetation on infestations of Oryctes rhinoceros (L.)(Col., Dynastidae) in young oil palm replantings in Malaysia. Bull. Entomol. Res. 59, 85–96 (1969). More

Retrieval of E. coli plasmid sequencesAll E. coli sequences were downloaded from the NCBI FTP server in May 2020. To establish an initial collection of plasmids, only complete genomes with an associated plasmid were retained. All genomes were verified for belonging to the species E. coli using kmerfinder (https://cge.cbs.dtu.dk/services/KmerFinder/). Sequence type (ST) was determined via multi-locus sequence typing (MLST) based on the 7-gene Achtman scheme using pubMLST (https:/github.com/tseemann/mlst). Only genomes with exact matches were assigned for each ST and used for subsequent analysis. To ensure our sequences were sufficiently representative of E. coli pathogens expected in nature, a systematic literature search (see description below and Fig. S1) was conducted to establish an expected distribution of STs (Table S1). This information was used to update our initial collection to match the top 4 most prevalent STs (131, 11, 73, and 95). Specifically, to identify supplementary plasmid sequences, genome accession IDs were chosen from EnteroBase based on the following criteria: the strain was matched to the correct ST and had a high-quality genome sequence (based on N50  > 20,000 and the number of contigs  0.1, 2-tailed student t test). For the second method, all kanR plasmids were used, and instead changed the hosts such that DH5αPro cells were in competition with DH5αPro containing a spontaneous rifampicin-resistant mutant (rifR). Any rifR strain was quantified on rifampicin-containing plates, and the second strain was quantified by rifampicin CFU minus CFU obtained on blank plates. We established that rifR exhibited no fitness defects by (1) growth rates between the wild-type (WT) strain (W) and rifR (M) (Fig. S5D), and (2) directly competing the two control strains (Fig. S5E). In both cases, results were statistically indistinguishable (p  > 0.1, two-tailed student t test). KanR/cmR and WT/rifR experiments were each conducted in LB or M9CAG, respectively. In all cases, experiments were repeated with at least three independent biological replicates.Time-kill measurements in the presence of carbenicillinAll strains were grown as previously described. Time-kill experiments entailed hourly measurements of CFU in presence of carbenicillin at either 3.75 μg/mL (3x IC50) or 5 μg/mL (4x IC50) over a span of 2 or 3 h, including time 0. Specifically, overnight cultures were first diluted 1:100 into LB media containing 1 mM IPTG and 50 μg/mL kanamycin and sub-cultured for two hours in a 37 °C incubator with shaking at 250 rpm. Following this, cell density was adjusted as necessary to achieve a starting OD600 of ~0.15 in all cases. Adjusted subcultures were then aliquoted into a 96-well plate and the appropriate carbenicillin treatments were added directly to the well. Plates were sealed with a paper film and placed in a 37 °C incubator with shaking at 250 rpm. Initial collection for time=0 was acquired before carbenicillin treatment. Thereafter, 10 μL of culture was removed from the well every hour, 10-fold serial dilutions were performed and 10 μL was plated on blank LB agar with three technical replicates at each time point. Colonies were counted after plates were grown for 16 h in a 37 °C incubator to determine CFU. This procedure utilized 14 strains of DH5αPro transformed with kanR plasmids of interest – ctrl, katG, lpxM, yfbR, aroH, pld, fdtC, agp, eptC, arcA, argF, mmuM, ahr, and fabG. CFUs were averaged for all technical replicates, and experiments were conducted with at least three independent biological replicates.Oxygen consumption rateOxygen consumption rates (OCR) were obtained with the Resipher device from Lucid Scientific. The selected strains were grown overnight as previously described. Overnight cultures were resuspended in M9CAG media with 1 mM IPTG and 50 μg/mL kanamycin, and placed in 25 °C for one hour to initiate gene expression. Following this, cells were diluted 10x into M9CAG media containing kanamycin and IPTG, and 100 μL was aliquoted per well into a 96-well microtiter plate according to the manufacturer’s instructions. Plates were placed at 30 °C to minimize growth, and oxygen concentration (μM) was measured immediately thereafter. 24 wells were measured consisting of 6 technical replicates for each strain. Given the clear well-well variability (Fig. S8B, C), data shown are for one biological replicate. However, qualitative trends were consistently reproduced in multiple independent experiments.StatisticsIn all cases where t tests and ANOVA’s were used, data was first verified to be normally distributed using Kolmogorov test for normality. Otherwise, Mann-Whitney U-tests were conducted. For panels with multiple tests, Bonferroni correction was used to adjust the p values. To determine whether any metabolic category was significantly dependent on incompatibility groups, we implemented logistic regressions in MATLAB with the function fitglm. Random forest classification was used to establish the relative importance of prevalent metabolic genes and gene categories predicting the presence of antibiotic resistance genes. Chi-square tests were conducted to determine significant co-occurrence of individual antibiotic resistant and metabolism genes. Dissociative relationships were distinguished by the odds ratios from the chi-square tests. To investigate whether the strong associations and disassociations were driven by evolutionary constraints, or simply artifacts of a common ancestor, we re-ran our statistical analysis using Coinfinder [29] to take in our gene presence-absence data, along with the genome phylogeny, and compute the Bonferroni-corrected statistical likelihood of coincidence (either associations or dissociations), thereby accounting for evolutionary relatedness. More

Dorning, M. & Cipollini, D. Leaf and root extracts of the invasive shrub, Lonicera maackii, inhibit seed germination of three herbs with no autotoxic effects. Plant Ecol. 184, 287–296 (2006).Article 

Google Scholar 
Greer, M. J., Wilson, G. W., Hickman, K. R. & Wilson, S. M. Experimental evidence that invasive grasses use allelopathic biochemicals as a potential mechanism for invasion: Chemical warfare in nature. Plant Soil 385, 165–179 (2014).Article 
CAS 

Google Scholar 
Möhler, H., Diekötter, T., Herrmann, J. D. & Donath, T. W. Allelopathic vs. autotoxic potential of a grassland weed-evidence from a seed germination experiment. Plant Ecol. Divers. 11, 539–549 (2018).Article 

Google Scholar 
Callaway, R. M. & Aschehoug, E. T. Invasive plants versus their new and old neighbors: A mechanism for exotic invasion. Science 290, 521–523 (2000).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Niu, H. B., Liu, W. X., Wan, F. H. & Liu, B. An invasive aster (Ageratina adenophora) invades and dominates forest understories in China: Altered soil microbial communities facilitate the invader and inhibit natives. Plant Soil 294, 73–85 (2007).Article 
CAS 

Google Scholar 
Wardle, D. A., Karban, R. & Callaway, R. M. The ecosystem and evolutionary contexts of allelopathy. Trends Ecol. Evol. 26, 655–662 (2011).Article 
PubMed 

Google Scholar 
Meiners, S. J., Kong, C. H., Ladwig, L. M., Pisula, N. L. & Lang, K. A. Developing an ecological context for allelopathy. Plant Ecol. 213, 1221–1227 (2012).Article 

Google Scholar 
Liebhold, A. M., Brockerhoff, E. G., Kalisz, S., Nunez, M. A. & Wardle, D. A. Biological invasions in forest ecosystems. Biol. Invasions 19, 3437–3458 (2017).Article 

Google Scholar 
Liao, H. X. et al. Soil microbes regulate forest succession in a subtropical ecosystem in China: Evidence from a mesocosm experiment. Plant Soil 430, 277–289 (2018).Article 
CAS 

Google Scholar 
Wardle, D. A., Nilsson, M. C., Gallet, C. & Zackrisson, O. An ecosystem-level perspective of allelopathy. Biol. Rev. 73, 305–319 (2010).Article 

Google Scholar 
Hierro, J. L. & Callaway, R. M. Allelopathy and exotic plant invasion. Plant Soil 256, 29–39 (2003).Article 
CAS 

Google Scholar 
Uddin, M. N., Robinson, R. W., Buultjens, A., Harun, M. A. & Shampa, S. H. Role of allelopathy of Phragmites australis in its invasion processes. J. Exp. Mar. Biol. Ecol. 486, 237–244 (2017).Article 

Google Scholar 
Thiébaut, G., Tarayre, M. & Rodríguez-Pérez, H. Allelopathic effects of native versus invasive plants on one major invader. Front. Plant Sci. 2, 854 (2019).Article 

Google Scholar 
Smith, M., Cecchi, L., Skjøth, C. A., Karrer, G. & Šikoparijae, B. Common ragweed: A threat to environmental health in Europe. Environ. Int. 61, 115–126 (2013).Article 
CAS 
PubMed 

Google Scholar 
Montagnani, C., Gentili, R., Smith, M., Guarino, M. F. & Citterio, S. The worldwide spread, success, and impact of ragweed (Ambrosia spp.). Crit. Rev. Plant Sci. 36, 1–40 (2017).Article 

Google Scholar 
Zeng, K., Zhu, Y. Q. & Liu, J. X. Research progress on ragweed (Ambrosia). Acta Prataculturae Sin. 19, 212–219 (2010).
Google Scholar 
Jacobs, R. L. et al. Responses to ragweed pollen in a pollen challenge chamber versus seasonal exposure identify allergic rhinoconjunctivitis endotypes. J. Allergy Clin. Immun. 130, 122-127.e8 (2012).Article 
PubMed 

Google Scholar 
Lake, R. I. et al. Climate change and future pollen allergy in Europe. Environ. Health Perspect. 125, 385–391 (2017).Article 
PubMed 

Google Scholar 
Wang, J. J., Zhao, B. Y., Li, M. T. & Li, R. Ecological invasion plant-bitter weed (Ambrosia artemisiifolia) and integrated control strategy. Pratacultural Sci. 023, 71–75 (2006).CAS 

Google Scholar 
Deng, Z. Z., Bai, J. D., Zhao, C. Y. & Li, J. S. Advance in invasion mechanisms of Ambrosia artemisiifolia. Pratacultural Sci. 32, 54–63 (2015).
Google Scholar 
Dong, H. G. et al. Diffusion and intrusion features of Ambrosia artemisiifolia and Ambrosia trifida in Yili River Valley. J. Arid Land Resour. Environ. 31, 175–180 (2017).
Google Scholar 
Vink, J. P. et al. Glyphosate-resistant giant ragweed (Ambrosia trifida) control in dicamba-tolerant soybean. Weed Technol. 26, 422–428 (2012).Article 
CAS 

Google Scholar 
Simard, M. J. & Benoit, D. L. Effect of repetitive mowing on common ragweed (Ambrosia artemisiifolia L.) pollen and seed production. Ann. Agric. Environ. Med. 18, 55–62 (2011).PubMed 

Google Scholar 
Goplen, J. J. et al. Seedbank depletion and emergence patterns of giant ragweed (Ambrosia trifida) in Minnesota cropping systems. Weed Sci. 65, 52–60 (2017).Article 

Google Scholar 
Jurik, T. W. Population distributions of plant size and light environment of giant ragweed (Ambrosia trifida L.) at three densities. Oecologia 87, 539–550 (1991).Article 
ADS 
PubMed 

Google Scholar 
Patracchini, C., Vidotto, F. & Ferrero, A. Common ragweed (Ambrosia artemisiifolia) growth as affected by plant density and clipping. Weed Technol. 25, 268–276 (2011).Article 

Google Scholar 
Kazinczi, G. Ragweed seed bank in the soils of arable fields of Transdanubia, Hungary. Hung. Weed Res. Technol. 19(1), 21–36 (2018).
Google Scholar 
Essl, F. et al. Biological flora of the British Isles: Ambrosia artemisiifolia. J. Ecol. 103, 1069–1098 (2015).Article 

Google Scholar 
Goplen, J. J. Giant Ragweed (Ambrosia trifida) Seed Bank Dynamics and Management. (Master’s dissertation, University of Minnesota.) Retrieved from https://hdl.handle.net11299174767 (2015).Yoda, K. Self-thinning in overcrowded pure stands under cultivated and natural conditions. J. Biol. 14, 107–129 (1963).
Google Scholar 
Friedman, J. & Waller, G. R. Allelopathy and autotoxicity. Trends Biochem. Sci. 10, 47–50 (1985).Article 
CAS 

Google Scholar 
Weller, D. E. The interspecific size-density relationship among crowded plant stands and its implications for the −3/2 power rule of self-thinning. Am. Nat. 133, 20–41 (1989).Article 

Google Scholar 
Deng, J. et al. Autotoxicity of phthalate esters in tobacco root exudates: Effects on seed germination and seedling growth. Pedosphere 27, 1073–1082 (2017).Article 
CAS 

Google Scholar 
Sudatti, D. B., Duarte, H. M., Soares, A. R., Salgado, L. T. & Pereira, R. C. New ecological role of seaweed secondary metabolites as autotoxic and allelopathic. Front. Plant Sci. 11, 347 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Singh, H. P., Batish, D. & Kohil, R. Autotoxicity: Concepts, organisms, and ecological significance. Plant Sci. 18, 757–772 (1999).CAS 

Google Scholar 
Chon, S. U. et al. Effects of alfalfa leaf extracts and phenolic allelochemicals on early seedling growth and root morphology of alfalfa and barnyard grass. Crop Prot. 21, 1077–1082 (2002).Article 
CAS 

Google Scholar 
Chen, B. M., D’Antonio, C. M., Molinari, N. & Peng, S. L. Mechanisms of influence of invasive grass litter on germination and growth of coexisting species in California. Biol. Invasions 20, 1881–1897 (2018).Article 

Google Scholar 
Chen, L. C., Wang, S. L., Wang, P. & Kong, C. H. Autoinhibition and soil allelochemical (cyclic dipeptide) levels in replanted Chinese fir (Cunninghamia lanceolata) plantations. Plant Soil 374, 793–801 (2014).Article 
CAS 

Google Scholar 
Perry, L. G. et al. Retracted: Dual role for an allelochemical: catechin from Centaurea maculosa root exudates regulates conspecific seedling establishment. J. Ecol. 93, 1126–1135 (2005).Article 
CAS 

Google Scholar 
Yu, J. Q., Ye, S. F., Zhang, M. F. & Hu, W. H. Effects of root exudates and aqueous root extracts of cucumber (Cucumis sativus) and allelochemicals, on photosynthesis and antioxidant enzymes in cucumber. Biochem. Syst. Ecol. 31, 129–139 (2003).Article 
CAS 

Google Scholar 
Kong, C. H., Wang, P. & Xu, X. H. Allelopathic interference of Ambrosia trifida with wheat (Triticum aestivum). Agric. Ecosyst. Environ. 119, 416–420 (2007).Article 
CAS 

Google Scholar 
Béres, I., Kazinczi, G. & Narwal, S. S. Allellopathic plants. 4. Common ragweed (Ambrosia elatior L. syn. A. artemisiifolia). Allelopathy J. 9, 27–34 (2002).
Google Scholar 
Bauer, J. T., Shannon, S. M., Stoops, R. E. & Reynolds, H. L. Context dependency of the allelopathic effects of Lonicera maackii on seed germination. Plant Ecol. 213, 1907–1916 (2012).Article 

Google Scholar 
Renne, I. J., Sinn, B. T., Shook, G. W., Sedlacko, D. M. & Hierro, J. L. Eavesdropping in plants: Delayed germination via biochemical recognition. J. Ecol. 102, 86–94 (2014).Article 

Google Scholar 
Loydi, A., Donath, T. W., Eckstein, R. L. & Otte, A. Non-native species litter reduces germination and growth of resident forbs and grasses: Allelopathic, osmotic or mechanical effects?. Biol. Invasions 17, 581–595 (2014).Article 

Google Scholar 
Bais, H. P., Weir, T. L., Perry, L. G., Gilroy, S. & Vivanco, J. M. The role of root exudates in rhizosphere interactions with plants and other organisms. Annu. Rev. Plant Biol. 57, 233–266 (2006).Article 
CAS 
PubMed 

Google Scholar 
Bonea, D., Bonciu, E., Niculescu, M. & Olaru, A. L. The allelopathic, cytotoxic and genotoxic effect of Ambrosia artemisiifolia on the germination and root meristems of Zea mays. Caryologia 71, 24–28 (2017).Article 

Google Scholar 
Dadkhah, A. Allelopathic effect of sugar beet (Beta vulgaris) and eucalyptus (Eucalyptus camaldulensis) on seed germination and growth of Portulaca oleracea. Russ. Agric. Sci. 39, 117–123 (2013).Article 

Google Scholar 
Zheng, L. & Feng, Y. L. Allelopathic effects of Eupatorium adenophorum Spreng on. seed germination and seedling growth in ten herbaceous species. Acta Ecol. Sin. 25, 2782–2787 (2005).CAS 

Google Scholar 
Brückner, D. J. The allelopathic effect of ragweed (Ambrosia artemisiifolia L.) on the germination of cultivated plants. Novenytermeles 47, 635–644 (1998).
Google Scholar 
Qin, R. M. et al. The evolution of increased competitive ability, innate competitive advantages, and novel biochemical weapons act in concert for a tropical invader. New Phytol. 197, 979–988 (2012).Article 
PubMed 

Google Scholar 
Zheng, Y. L. et al. Integrating novel chemical weapons and evolutionarily increased competitive ability in success of a tropical invader. New Phytol. 205, 1350–1359 (2015).Article 
PubMed 

Google Scholar 
Kaushal, R., Verma, K. S. & Singh, K. N. Effect of Grewia optiva and Populus deltoides leachatesv on field crops. Allelopathy J. 11, 229–234 (2003).
Google Scholar 
Kumari, A. & Kohli, R. Autotoxicity of ragweed parthenium (Parthenium hysterophorus). Weed Sci. 35, 629–632 (1987).Article 

Google Scholar 
Einhellig, F. A. Allelopathy: Current status and future goals. In Allelopathy: Organisms, processes and applications (ed. Inderjit Dakshini, K. M. M.) 1–24 (Am Chem. Soc, Washington, 1995).
Google Scholar 
Hadack, F. Secondary metabolites as plant traits: Current assessment and future perspectives. Crit. Rev. Plant Sci. 21, 273–322 (2002).Article 

Google Scholar 
Rice, E. L. Biological Control of Weeds and Plant Diseases (Oklahomka Press, 1995).
Google Scholar 
Choi, B. et al. Common ragweed-derived phenolic compounds and their effects on germination and seedling growth of weed species. Weed Turfgrass Sci. 30, 396–404 (2010).
Google Scholar 
Friedman, J. & Waller, G. R. Seeds as allelopathic agents. Chem. Ecol. 9, 1107–1117 (1983).Article 
CAS 

Google Scholar 
Canals, R. M., Emeterio, L. S. & Peralta, J. Autotoxicity in Lolium rigidum: Analyzing the role of chemically mediated interactions in annual plant populations. J. Theor. Biol. 235, 402–407 (2005).Article 
ADS 
CAS 
PubMed 

Google Scholar 
San Emeterio, L., Damgaard, C. & Canals, R. M. Modelling the combined effect of chemical interference and resource competition on the individual growth of two herbaceous populations. Plant Soil 292, 95–103 (2007).Article 
CAS 

Google Scholar 
Dickerson, C. T. Studies on the germination, growth, development and control of Common Ragweed (Ambrosia artemisiifolia L.). PhD thesis, Cornell University, Ann Arbor (1968).Nuutinen, V. & Butt, K. R. Homing ability widens the sphere of influence of the earthworm Lumbricus terrestris L. Soil Biol. Biochem. 37, 805–807 (2005).Article 
CAS 

Google Scholar 
Favaretto, A., Scheffer-basso, S. M. & Perez, N. B. Autotoxicity in tough lovegrass (Eragrostis plana). Planta Daninha 35(35), e017164046 (2017).
Google Scholar 
Sinkkonen, A. Modelling the effect of autotoxicity on density-dependent phytotoxicity. J. Theor. Biol. 244, 218–227 (2007).Article 
ADS 
MathSciNet 
CAS 
PubMed 
MATH 

Google Scholar 
Zhang, S. S., Shi, F. Q., Yang, W. Z., Xiang, Z. Y. & Duan, Z. L. Autotoxicity as a cause for natural regeneration failure in Nyssa yunnanensis and its implications for conservation. Isr. J. Plant Sci. 62, 187–197 (2015).Article 

Google Scholar 
Liu, Y. et al. Relationship between seed germination and invasion of Ambrosia artemisiifolia and A. trifida at different positions. Acta Ecol. Sin. 39, 9079–9088 (2019).

Google Scholar  More

Insect rearingThe cotton mealybugs used in this study were originally collected from Rose of Sharon, Hibiscus syriacus L. (Malvales: Malvaceae) in Jinhua, Zhejiang Province, China, in June 2016. They were maintained on fresh tomato plants (cv. Hezuo-903, Shanghai Changzhong Seeds Industry Co., Ltd, China) in a climatically controlled chamber maintained at 27 ± 1 °C, 75% relative humidity (RH), and a photoperiod of 14:10 (L:D). For detailed insect rearing and tomato cultivation methods see ref. 56.Scanning electron microscopy (SEM) of P. solenopsis waxSEM was used to observe changes in wax on the body surface of adult P. solenopsis females according to the methods of Huang et al.57. Briefly, collected insects were taped onto a stub and dried in an ion sputter (Hatachi, Tokyo, Japan) under a vacuum. After gold sputtering, the samples were observed using a TM-1000 SEM (Hatachi, Tokyo, Japan). Photos were scanned from the dorsal part of the third thoracic segment. Thirty insects were used for both RNAi-treated and control groups.Chemical composition analysis of mealybug waxA small soft brush was used to collect wax filaments from the body surface of P. solenopsis females. Prior to use, the brush was washed successively by 70% ethanol, sterile water, and 1× sterile phosphate-buffered saline (PBS, pH 7.4). The wax was collected into a clean chromatography vial for the following experiments. Two vials of wax, each collected from 1000 adult females, were dissolved in 1 ml of methanol and 1 ml of n-hexane, respectively. The vials were stirred gently for 3 min, kept at room temperature for 30 min, and then put into an S06H ultrasonic vibrator (Zealway, Xiamen, China) for 30 min to dissolve the wax sufficiently. The samples were analyzed on a TRACE 1310 (Thermo Scientific, Waltham, USA) gas chromatograph (GC) equipped with an ISQ single quadrupole MS and interfaced with the Chromeleon 7.2 data analysis system (Thermo Scientific, Waltham, USA), with a constant flow of helium at 1 ml/min. For each sample, a splitless injection of 1.0 μl was respectively made into a polar TG-WaxMS (Thermo Scientific, Waltham, USA) and a nonpolar TG-5MS (Thermo Scientific, Waltham, USA) 30 m × 0.25 mm × 0.25 μm capillary column. The temperature program for polar column samples was as follows: 40 °C for 2 min, then 5 °C/min to 240 °C, hold 10 min; the program for nonpolar column samples was: 40 °C for 2 min, then 5 °C/min to 300 °C, hold 5 min. Injector and detector temperatures were, respectively, set at 250 and 230 °C for polar column samples, and at 300 and 300 °C for nonpolar column samples. Mass detection for all samples was run under an EI mode with a 70 eV ionization potential and an effective m/z range of 35–450 at a scan rate of 5 scan/s. Chemical compounds were identified by mapping against the NIST database. The relative content of each compound was calculated by peak area which was determined using the Agilent MassHunter system.RNA extraction and RT-qPCRTotal RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA) following the manufacturer’s instructions, and RNA quality was accessed using agarose gel electrophoresis and a Biodrop μLite. 800 ng of total RNA was used for cDNA synthesis using the HiScript III RT SuperMixfor qPCR (+gDNA wiper) (Vazyme Biotech Co., Ltd., Nanjing, China), according to the manufacturer’s instructions. Quantitative RT-PCR (RT-qPCR) was conducted using an AriaMx real-time PCR system (Agilent Technologies, USA), using a 20 μl reaction containing 2 μl of 10-fold diluted cDNA, 0.8 μl of each primer, and 10 μl ChamQ SYBR Color qPCR Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China). The RT-qPCR thermocycling protocol was 95 °C for 30 s, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s. The PsActin gene was used as an internal control. At least three biological replicates were used for each experiment. Quantitative variations were evaluated using the relative quantitative method (2−ΔΔCt)58.Transcriptome analysis of integumentary and non-integumentary tissuesTo obtain the integument and other tissues, adult P. solenopsis females were dissected in 1× sterile PBS (pH 7.4) on a sterile Petri dish. Dissected fresh tissues were directly used or frozen in liquid nitrogen and stored at −80 °C for follow-up experiments. We sequenced the transcriptomes of integumentary and non-integumentary tissues (all other tissues without integument) dissected from 150 adult females, with each sample being repeated in triplicate. mRNAs were purified from total RNA via oligo (dT) magnetic beads, and the fragmented mRNAs were then reverse transcribed into cDNA using random primers. Constructed pair-end libraries were sequenced using an Illumina HiSeq X Ten platform in Novogene (Beijing, China). After quality control, the clean RNA-Seq data of the six libraries were aligned with the P. solenopsis genome (http://v2.insect-genome.com/Organism/624) using HISTAT259. Then featureCounts60 and DESeq261 were used for the differential expression analysis of genes. The threshold for differentially expressed genes (DEGs) was defined by log2fold ≥ 1 or ≤−1 and a padj-value  More

Basic characteristics of river sedimentsA considerable variation was found in the distribution of clay (81 to 48.4 g kg−1), silt (145 to 656 g kg−1), and sand (38 to 821 g kg−1) particles among sediment materials. The associated coefficient of variations (CV) was 57, 59.5, and 41%, respectively. Statistical data related to the physicochemical properties of sediments and their main elements are reported in Table 2. The variations in particle size distribution located sediment material in seven textural classes ranging from loamy sand to silty clay. The high variability in particle size distribution suggests that different sets of geogenic and anthropogenic processes are enacted in the development and distribution of sediments in the rivers. The pH and CCE ranged from 7.4 to 8.2 and 31 to 251 g kg−1, respectively, indicating the dominancy of alkaline-calcareous condition. None of the sediment samples exhibited salinity conditions (EC  > 4 dS m−1) with EC in the range of 0.3 to 1.4 dS m−1. A relatively low range of OM was found in all samples ranging from 7 to 61 g kg−1 with a mean value of 19 g kg−1. This range of OM coincides with the corresponding values in regional soils47. Except for pH, other sediments properties demonstrated above 35% of CV illustrating a wide range of variability in sediments’ physicochemical properties across the study rivers.Table 2 Summary statistics of sediment properties.Full size tableThe highest concentration among major elements was observed in SiO2, varying between 37.5 and 55.2%, with a mean percentage of 44.9%. This element followed in magnitude by Al2O3 (8.9–15.9%), CaO (5–14.3%), Fe2O3 (4.8–10%), MgO (2.4–17.2%), K2O (1.2–3.1%), Na2O (0.68–2.7%), SO3 (0.01–4.8% g kg−1) (Table 2). Considering the semi-arid climatic condition of the study region, higher levels of SiO2 and lower levels of Al2O3 may indicate that the silicate minerals forming the sediments of the area have not been subjected to severe weathering processes. Likewise, the Na2/K2O ratio was greater than 1 in the majority of sediment samples, implying an enrichment of potassium feldspar and the relatively intense weathering of Na-bearing minerals in the region48,49. The CIA value was in the range of 64.9 to 85.7% with a mean percentage of 72.9%, representing a moderate chemical weathering intensity of lithological materials (65%  Pb  > Cu  > Cd which varied largely among the sampling points. The level of Zn, Cu, Cd, Pb, and Ni varied in the ranges of 32.6–87.5, 14.2–33.3, 0.42–4.8, 14.5–69.5, and 20.1–183.5 mg kg-1, respectively, for winter, and 35.3–92.5, 15.6–35.1, 0.47–5.1, 15.5–73.1, 23.2–188.3 mg kg−1 for summer. The obtained ranges are comparable with data found in previous studies in Asia4,54,55,54.Figure 2The comparison of the mean concentration of Zn, Cu, Cd, Pb, and Ni elements in the study rivers’ sediments during summer and winter. Different letters show significant differences in metal content among rivers pooled over seasons at P  More

This study quantitatively assessed waterbody loss due to urban expansion of large Chinese cities. We first extracted multi-temporal urban boundaries to determine the expansion of cities of over one million in population from 1990 to 2018. The monthly surface-water dataset was then used to identify surface waterbodies in the study period. Depending on the ratio of surface waterbody area to urban area, cities were further divided into three categories (i.e. water-abundant, water-medium, water-deficient). Finally, we quantified the rate of waterbody loss and evaluated the spatial and temporal variation of waterbody loss as a function of urban expansion and according to city type.GUB datasetThe Global Urban Boundary (GUB) dataset (http://data.ess.tsinghua.edu.cn) was used to determine urban expansion. GUB provides data on built-up areas over 30 years, with a spatial resolution of 30 m. In the GUB dataset, nonurban areas (such as green space and water space) surrounded by artificial impervious areas are filled within the urban boundary and removed by the algorithm, which is consistent with global mapping methods. The continuous urban boundary was demarcated by morphological image processing methods, which have an overall accuracy of over 90%. In this dataset, extensive water and forests are excluded, and the impervious surface within the urban boundaries accounts for about 60% of the total surface area47. Compared with urban boundaries obtained from night-time light, GUB better separates urban areas from surrounding nonurban areas.Monthly waterbody datasetWe selected the JRC Monthly Water History V1.3 dataset(https://global-surface-water.appspot.com/), which is available from the Google Earth Engine, as the basis for representing surface waterbodies48. This data collection, which was produced by using images from the Landsat series, contains 442 images of global monthly waterbody area from March 1984 to December 2020. In this dataset, the validation confirmed that fewer than 1% of waterbodies were incorrectly detected, and fewer than 5% of waterbodies were missed altogether. We chose this dataset due to the long-term spatial distribution of waterbodies and due to mountain shadows and urban-constructions masking, which reflects the real changes in waterbodies.Theoretical backgroundIt is well known that cities have high concentrations of population and resources and expand spatially during development. There are many different perspectives on the size of cities, and studies have mostly used urban density and population to characterize them. However, because it is challenging to standardize data sources and quality, there is no unified quantitative standard49. Urban construction has concentrated human activity and brought about changes in land types. Cities are also identified as physical spaces, which can be defined as the built environment50,51. The built environment, which includes structures like buildings, roads, and other artificial constructions, is sometimes referred to as a non-natural environment52.Rural is the antithesis of urban. As large cities have spread outward in developing nations like Asia, a transitional fringe has been created by the gradual blurring of the line separating urban and rural areas53. According to McGee, good locations, easy access, and sizable agricultural land all contribute to the development potential of large cities. Thus, between urban and rural areas, there are transitional areas of active spatial morphological change known as desakota33,54. The peri-urban areas, like desakota, are gradually developed and incorporated into original built-up urban areas in urbanization. The original landscape, which included agricultural land, vegetation, and waterbodies, gradually changed into an urban land use type, i.e. impervious surface, and thus the city continues to expand outwards. Waterbody, an essential ecological element, has been heavily developed or filled in during urbanization, which may present dangerous ecological risks. In this paper, we identified the urban boundaries based on physical space to explore the encroachment activities on waterbodies during the urbanization of large cities. We determined whether existing waterbodies were transformed into urban waterbodies or encroached upon and whether waterbodies were increased in the expansion of urban boundaries, thus proposing strategies for protecting waterbodies in the future.Extracting the extent of large Chinese cities from GUB datasetTo characterize urban expansion, GUB data are selected as the original data for urban boundary selection. The Chinese administrative scale of municipalities is not exclusively urban, but also includes rural areas. In our study, cities were defined as municipal districts excluding the vast countryside within the administrative boundaries of prefecture-level cities. We identified urban areas based on the physical boundaries from the perspective of remote sensing, which can precisely track urban expansion51.In this work, we selected 159 cities with a population of over one million in 2018 based on the average annual population of urban districts from the 2019 China City Statistical Yearbook (Fig. S1). Taiwan, Hong Kong, and Macau are omitted. According to statistics, China had 160 cities with populations exceeding one million in 2018. However, due to the lack of data for the built-up area in 1990, Guang’an was not included in the study. We thus obtained 159 cities from the GUB dataset. Due to numerous fragmented patches within the administrative boundary, the population identified the main urban areas, and max patch areas were comprehensively based on the urban boundaries. Through manual detection and adjustment of the map, we determined that the location of the extracted urban area was consistent with that of the municipal government, and the boundary was extracted for each period. We took the growth area as the expansion area, with the original area being the city at the onset of each period (Fig. S3).We used the average annual urban growth (AUG) rate to characterize the rate of urban expansion, as is widely done to evaluate urban expansion55,56. It is calculated as$${text{AUG}} = left[ {frac{{Land_{t1} }}{{Land_{t0} }}^{{frac{1}{t1 – t0}}} – 1} right] times 100% ,$$
where (Land_{t0}) and (Land_{t1}) represent the urban land area at time t0 and t1, where t0 and t1 are the start and end of the given study period.Identification of urban waterbodiesUrban waterbodies contain all the components of urban flow networks above the ground and include natural waterbodies such as lakes, rivers, streams, and wetlands and artificial waterbodies such as parks and ponds48. We identified all waterbodies existing within the urban boundary as urban waterbody. Considering urban expansion, urban waterbodies vary as urban boundary shift at different stages. Our study explored how the original waterbodies changed under urban expansion, including whether they were kept as urban waterbodies or encroached upon. Considering the dryness or wetness of each year, we used the data for 3 years (36 months) around each period (1990, 1995, 2000, 2005, 2010, 2015, and 2018) to describe the waterbody. Not all waterbodies could be detected for each month of the year; for example, freezing may prevent waterbodies from being detected. To cover seasonal and permanent waterbodies, we used the waterbody frequency index (WFI), which is calculated as the fraction of waterbody months within the 3 years to identify stable waterbodies pixel by pixel57. The spatial distribution of each waterbody was then mapped comprehensively for each period. By comparing the extracted waterbody with the long-time-series high-resolution remote-sensing images from Google Earth, we found that the extracted waterbodies fit the actual waterbody distribution quite well (Fig. S2):$$WFIleft( i right) = frac{WMleft( i right)}{{DMleft( i right)}}$$
where WFI(i) is the water occurrence for pixel i in the images before and after the given year, and i is the pixel number for the study area. WM(i) is the number of months during which the waterbody is detected in i pixel over the 3 years. DM(i) is the number of months during which the data are available in pixel i. If the waterbody frequency index of a pixel is greater than 25%, this pixel is considered as a waterbody; otherwise, it is not.City classification based on surface waterbodyCities with over one million in population may not be short of waterbodies, but significant differences remain in surface waterbody abundance. Due to large differences in city size, it is inappropriate to use waterbody area as a criterion. Considering the influence of urban expansion, we ranked 159 cities according to the indicator of waterbody fraction (WF), namely the fraction of the original surface water within the urban boundary in 2018. Waterbodies not impacted by urbanization were taken as the original surface waterbody, which used the average surface waterbody from 1985 to 1991 as baseline. We used the natural break method to divide cities into abundant, moderate, and deficient levels (referred to as Type I, Type II, and Type III, respectively) and evaluate the abundance of waterbodies in cities. Based on the waterbody fraction (WF) value, which is calculated as follows:$${text{WF}} = frac{{Water_{origin} }}{{Land_{2018} }}$$
where WF is used to judge the urban waterbody abundance in cities. (Water_{1990}) is the origin surface waterbody area (used the year in 1985–1991) in the urban boundary of 2018, (Land_{2018}) the urban land area in the urban boundary of 2018.Temporal characteristic of waterbody loss and gainTo understand the spatial–temporal features of surface waterbodies, we used five normalized indicators to compare waterbody variations between cities during urban expansion from the overall perspective and from the city perspective.The variation in original natural waterbodies reflects the intensity of the natural resource development in urban expansion. We summarized the reduction and preservation of original waterbodies in urban expansion areas with a population of over one million to represent the encroachment of urban expansion on waterbodies:$$WL = frac{{sum NWL_{t0_t1} left( i right)}}{{sum W_{t0} left( i right)}} times 100%$$$$WP = frac{{sum (W_{t0} left( i right) – NWL_{t0_t1} left( i right))}}{{sum W_{t0} left( i right)}} times 100%$$
where i labels the city within the 159 cities, WL and WP are the fractions of waterbody loss and preservation in urban expansion areas of all cities, (NWL_{t0_t1}) is the net waterbody loss during period t0–t1 (, and W_{t0}) is the natural waterbody in the urban expansion area at time t0.To estimate the net waterbody loss caused by urban expansion at various stages, we used the standardized indicator, annual average net waterbody loss rate (ANWL), to compare waterbody loss speeds over time. This indicator is independent of the difference in waterbody abundance and can be compared over time. Waterbody loss is one part of the impact of urbanization; the other is waterbody gain. We used the same method to evaluate the annual average net waterbody gain rate (ANWG). The formulas are$$A{text{NWL}} = frac{{NWL_{t0_t1} }}{{W_{t0} left( {t1 – t0} right)}} times 100%$$$$ANWG = frac{{NWG_{t0 – t1} }}{{W_{t0} left( {t1 – t0} right)}} times 100%$$
where NWL and NWG are the net waterbody loss and gain, respectively, and the other abbreviations are the same as above.Considering the direct impact of urban expansion, we used a normalized indicator, the average net waterbody loss velocity of urban expansion ((AWLV)), which refers to the amount of waterbody encroachment per unit urban expansion area. It quantifies the time-heterogeneity of waterbody loss due to urban expansion and is calculated as follows:$$AWLV = frac{{NWL_{t0_t1} }}{{Land_{t1} – Land_{t0} }}$$We calculated these indicators for the six expansion periods (1990–1995, 1995–2000, 2000–2005, 2005–2010, 2010–2015, and 2015–2018) (Fig. 3). In the study, if the waterbody pixel count is zero at the onset of the period, the indicator for the period is abnormal and thus excluded. More