Ecology

Soil physicochemical propertiesTable 1 shows the basic chemical characteristics of the 81 soil samples. The pH values of most of the soil samples from the rape fields and the paddy fields ranged between 6.5 and 7.5, while the pH values of most of the soil samples from the wheat fields were less than 6.0. The relatively low pH values for soils from the wheat fields could be due to traditional farming practices adopted in the farms, including continuous cropping. In addition, acidic soil is conducive for wheat growth. Organic matter contents were generally low in all the soils, and soils from the rape fields had the lowest organic matter contents. Generally, the mean soil total N contents in the fields were in the order of rape fields (756 mg/kg)  Zn  > Cd, and the concentrations of metals in most soil samples collected from the wheat and paddy fields were in the order of Fe  > Mn  > Pb  > Cu  > Zn  > Cd. However, independent of the farmland use patterns, the concentrations of all the metals extracted using NH4OAC and NH4NO3 were much lower than those extracted with the other three types of extractants.Table 3 Extractable metal concentrations in the soils (mg/kg).Full size tableThe concentrations of 0.1 mol/L HCl-extracted heavy metals in soils from the rape fields at the R21, R23, and R32 sites, in the wheat fields at the W21 site, and in the paddy fields at the P1 and P2 sites were extremely low, especially in the cases of Fe and Pb. For example, 0.1 mol/L HCl-extracted Fe concentrations at the R21, R23, and R32 sites were 0.2 mg/kg, 0.6 mg/kg, and below the detection limit, respectively, whereas the range of 0.1 mol/L HCl-extracted Fe concentration in the other sites in the rape fields was 22.4–533 mg/kg. Therefore, 0.1 mol/L HCl-extracted heavy metal concentrations at the R21, R23, R32, W21, P1, and P2 sites were omitted from subsequent analyses to ensure homogeneity of variance. The basic pH values at the sites (pH values at R21, R23, R32, W21, P1, and P2 were 8.0, 8.0, 7.7, 8.0, 8.0, and 8.0, respectively) could explain this low HCl-extractability and reduced mobility35,36. Although HCl is considered a universal extractant, it may not be suitable under alkaline soil conditions.Heavy metal concentrations in rape, wheat and rice grainsFigure 4 illustrates the concentrations of six heavy metals in rape, wheat, and rice grains. Similar to the order of the soil metal concentrations, the Fe, Mn, Zn, and Cu concentrations in the grains were much higher than the Pb and Cd concentrations.Figure 4Concentrations of heavy metals ((a) Cu, (b) Zn, (c)Pb, (d) Cd, (e) Fe, and (f) Mn) in the grain of three different crops. (Rape grains (n = 36); Wheat grains (n = 25); Rice grains (n = 20), on dry weight basis).Full size imageThe results are similar to those reported in previous studies with regard to heavy metal concentration trends in rice grains16,37. The reason could be Fe, Mn, Zn, and Cu are all essential for crop growth as micronutrients, leading to the higher levels in soils. Among the studied crops, rice grains accumulated relatively lower amounts of Fe, Mn, and Zn than rape and wheat grains, where Pb and Cd concentrations exhibited opposite trends. The trends are consistent with the findings of Liu et al.38 and Du et al.1, and indicate that rice has a stronger Cd uptake capacity from soil. Williams et al.39 also reported higher rice Cd concentrations than wheat, barley, and maize Cd concentrations. The results could also be attributed to water management and soil oxidation–reduction status in soil-rice systems. Paddy fields are irrigated considerably more than wheat and rape fields. In addition, the water in the study area was contaminated by Cd, Cu, As, and Zn27. Previous studies have also reported that water management practices influence Cd uptake by rice and its bioavailability in soils19,20,40. In the cases of Zn and Pb, according to Feng et al.23, rice grains accumulated lower amounts of Zn than wheat grains, whereas rice grains accumulated higher Pb amounts than wheat grains. Although soil Cd concentrations were generally high in the rape fields, the concentrations of Cd in rape grains were lower than those in the wheat and rice grains.On the other hand, Cu concentrations in grains in all the three crops were below the maximum allowable Cu levels in food (10 mg/kg (GB 15199-94)). Similarly, Zn concentrations in rice grains were below the maximum allowable Zn levels in food (50 mg/kg (GB 13106-91)), while 25% of the rape grain samples and 20% of the wheat grain samples exceeded the threshold value. Although total Pb concentrations in soils were generally low in the three farmland use patterns, Pb concentrations in 70% of the rice grain samples exceeded the maximum allowable Pb levels in food (0.2 mg/kg (GB 2762-2005)). Only four rape and two wheat grain samples exceeded the Pb threshold values, respectively. The varying Pb trends are potentially linked to physical contamination from direct atmospheric deposition8, differences in physiological activities among the crops, and the fruit structures of the studied crops23. Similar to the soil Cd contamination, 96% and 60% of wheat and rice grain samples, respectively exceeded the maximum allowable Cd levels in food (0.1 mg/kg and 0.2 mg/kg for wheat and rice, respectively (GB 2762-2005)). Furthermore, Cd concentrations in approximately 10% of the rice grain samples exceeded 1.0 mg/kg. Conversely, only 14% of rape grain samples exceeded the maximum allowable Cd levels in food (Cd: 0.1 mg/kg for rape (GB 2762-2005)), and the maximum Cd concentration in rape grains was 0.18 mg/kg. According to the results, rape grains were generally safe for consumption whereas wheat and rice grains posed health threats in the study area.Soil to grain bioaccumulation factorsBAF values have been used widely to evaluate the capacity of crop grains to accumulate metals from soil30,41. Similar to a previous study on food crops37, Fe and Pb had the lowest BAF values (Fig. 5). Generally, metal accumulation in crop grains did not increase considerably with an increase in total concentrations of metals in soil. Heavy metal accumulation could have been regulated by crops, so that only low amounts were accumulated into grains. In addition, there were significant differences in some BAF values of the same metal across different crop species (Fig. 5). Different crop species have different accumulation capacities for the same metal42. Overall, the average BAF values of Cu (0.22), Zn (0.37), and Mn (0.14) in wheat grains were significantly higher than those in rape and rice grains. The average BAF value of Pb (0.005) in rice grains was significantly higher than those in rape and wheat grains, and the average BAF values of Cu (0.07) and Cd (0.06) in rape grains were significantly lower than those in wheat and rice grains. The results indicated that rape grains have lower heavy metal accumulation capacity than wheat and rice grains, except in the cases of Zn and Fe. However, the finding is not consistent with the results of a previous study43, which report that grasses have lower accumulation capacity than dicotyledonous plants. Nevertheless, as mentioned above, numerous factors could influence the accumulation capacity of metals in crop grains.Figure 5Bioaccumulation factors (BAF) of heavy metals ((a) Cu, (b) Zn, (c) Pb, (d) Cd, (e) Fe, and (f) Mn) from soil to the grains of three crop species. The error bars indicate the standard deviation. Different letters on bars indicate significant difference (p  More

Field, C. B., Behrenfeld, M. J., Randerson, J. T. & Falkowski, P. Primary production of the biosphere: integrating terrestrial and oceanic components. Science 281, 237 (1998).CAS 
PubMed 

Google Scholar 
Anbar, A. D. & Knoll, A. H. Proterozoic ocean chemistry and evolution: a bioinorganic bridge? Science 297, 1137–1142 (2002).CAS 
PubMed 

Google Scholar 
Saito, M. A., Sigman, D. M. & Morel, F. M. M. The bioinorganic chemistry of the ancient ocean: the co-evolution of cyanobacterial metal requirements and biogeochemical cycles at the Archean–Proterozoic boundary? Inorg. Chim. Acta 356, 308–318 (2003).CAS 

Google Scholar 
Morel, F. M. M., Lam, P. J. & Saito, M. A. Trace metal substitution in marine phytoplankton. Annu. Rev. Earth Planet Sci. 48, 491–517 (2020).CAS 

Google Scholar 
Morel, F. M. & Price, N. M. The biogeochemical cycles of trace metals in the oceans. Science 300, 944–947 (2003).CAS 
PubMed 

Google Scholar 
Twining, B. S. & Baines, S. B. The trace metal composition of marine phytoplankton. Annu. Rev. Mar. Sci. 5, 191–215 (2013).
Google Scholar 
Ho, T.-Y. et al. The elemental composition of some marine phytoplankton. J. Phycol. 39, 1145–1159 (2003).CAS 

Google Scholar 
Ellwood, M. J. Wintertime trace metal (Zn, Cu, Ni, Cd, Pb and Co) and nutrient distributions in the subantarctic zone between 40–52°S; 155–160°E. Mar. Chem. 112, 107–117 (2008).CAS 

Google Scholar 
Zhao, Y., Vance, D., Abouchami, W. & de Baar, H. J. W. Biogeochemical cycling of zinc and its isotopes in the Southern Ocean. Geochim. Cosmochim. Acta 125, 653–667 (2014).CAS 

Google Scholar 
John, S. G., Helgoe, J. & Townsend, E. Biogeochemical cycling of Zn and Cd and their stable isotopes in the Eastern Tropical South Pacific. Mar. Chem. 201, 256–262 (2018).CAS 

Google Scholar 
Middag, R., de Baar, H. J. W. & Bruland, K. W. The relationships between dissolved zinc and major nutrients phosphate and silicate along the GEOTRACES GA02 transect in the West Atlantic Ocean. Glob. Biogeochem. Cy. 33, 63–84 (2019).CAS 

Google Scholar 
Sunda, W. G. & Huntsman, S. A. Feedback interactions between zinc and phytoplankton in seawater. Limnol. Oceanogr. 37, 25–40 (1992).CAS 

Google Scholar 
Sunda, W. G. & Huntsman, S. A. Cobalt and zinc interreplacement in marine phytoplankton: biological and geochemical implications. Limnol. Oceanogr. 40, 1404–1417 (1995).CAS 

Google Scholar 
Vance, D. et al. Silicon and zinc biogeochemical cycles coupled through the Southern Ocean. Nat. Geosci. 10, 202 (2017).CAS 

Google Scholar 
Weber, T., John, S., Tagliabue, A. & DeVries, T. Biological uptake and reversible scavenging of zinc in the global ocean. Science 361, 72 (2018).CAS 
PubMed 

Google Scholar 
Roshan, S., DeVries, T., Wu, J. & Chen, G. The internal cycling of zinc in the ocean. Glob. Biogeochem. Cy. 32, 1833–1849 (2018).CAS 

Google Scholar 
Scott, C. et al. Bioavailability of zinc in marine systems through time. Nat. Geosci. 6, 125–128 (2012).
Google Scholar 
Mock, T. et al. Evolutionary genomics of the cold-adapted diatom Fragilariopsis cylindrus. Nature 541, 536–540 (2017).CAS 
PubMed 

Google Scholar 
Blaby-Haas, C. E. & Merchant, S. S. Comparative and functional algal genomics. Annu. Rev. Plant Biol. 70, 605–638 (2019).CAS 
PubMed 

Google Scholar 
Zhang, Z. H. et al. Adaptation to extreme Antarctic environments revealed by the genome of a sea ice green alga. Curr. Biol. 30, 3330–3341 (2020).CAS 
PubMed 

Google Scholar 
Clarke, A. et al. The Southern Ocean benthic fauna and climate change: a historical perspective. Philos. Trans. R. Soc. Lond. B 338, 299–309 (1992).
Google Scholar 
Klug, A. The discovery of zinc fingers and their applications in gene regulation and genome manipulation. Annu. Rev. Biochem. 79, 213–231 (2010).CAS 
PubMed 

Google Scholar 
Krishna, S. S., Majumdar, I. & Grishin, N. V. Structural classification of zinc fingers: survey and summary. Nucleic Acids Res. 31, 532–550 (2003).CAS 
PubMed 
PubMed Central 

Google Scholar 
Barlow, P. N. et al. Structure of the C3HC4 domain by 1H-nuclear magnetic resonance spectroscopy: a new structural class of zinc-finger. J. Mol. Biol. 237, 201–211 (1994).CAS 
PubMed 

Google Scholar 
Stephens, T. G. et al. Genomes of the dinoflagellate Polarella glacialis encode tandemly repeated single-exon genes with adaptive functions. BMC Biol. 18, 56 (2020).Aranda, M. et al. Genomes of coral dinoflagellate symbionts highlight evolutionary adaptations conducive to a symbiotic lifestyle. Sci. Rep. 6, 39734 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Liu, H. et al. Symbiodinium genomes reveal adaptive evolution of functions related to coral–dinoflagellate symbiosis. Commun. Biol. 1, 95 (2018).PubMed 
PubMed Central 

Google Scholar 
Shoguchi, E. et al. Draft assembly of the Symbiodinium minutum nuclear genome reveals dinoflagellate gene structure. Curr. Biol. 23, 1399–1408 (2013).CAS 
PubMed 

Google Scholar 
Shoguchi, E. et al. Two divergent Symbiodinium genomes reveal conservation of a gene cluster for sunscreen biosynthesis and recently lost genes. BMC Genomics 19, 458 (2018).PubMed 
PubMed Central 

Google Scholar 
Hoppe, C. J. M., Flintrop, C. M. & Rost, B. The Arctic picoeukaryote Micromonas pusilla benefits synergistically from warming and ocean acidification. Biogeosciences 15, 4353–4365 (2018).CAS 

Google Scholar 
Ferguson, R. E. et al. Housekeeping proteins: a preliminary study illustrating some limitations as useful references in protein expression studies. Proteomics 5, 566–571 (2005).CAS 
PubMed 

Google Scholar 
Aslam, S. N. et al. Identifying metabolic pathways for production of extracellular polymeric substances by the diatom Fragilariopsis cylindrus inhabiting sea ice. ISME J. 12, 1237–1251 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Valenzuela, J. J. et al. Ocean acidification conditions increase resilience of marine diatoms. Nat. Commun. 9, 2328 (2018).PubMed 
PubMed Central 

Google Scholar 
Martin, K. et al. The biogeographic differentiation of algal microbiomes in the upper ocean from pole to pole. Nat. Commun. 12, 5483 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
Mock, Thomas. Sea of Change: Eukaryotic Phytoplankton Communities in the Arctic Ocean. United States. https://doi.org/10.25585/1488054Duncan, A. et al. Metagenome-assembled genomes of phytoplankton communities across the Arctic Circle and Atlantic Oceans. Microbiome 10 https://doi.org/10.1186/s40168-022-01254-7 (2022).Persi, E., Wolf, Y. I. & Koonin, E. V. Positive and strongly relaxed purifying selection drive the evolution of repeats in proteins. Nat. Commun. 7, 13570 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Mock, T. & Gradinger, R. Determination of Arctic ice algal production with a new in situ incubation technique. Mar. Ecol. Prog. Ser. 177, 15–26 (1999).CAS 

Google Scholar 
Rühle, T., Hemschemeier, A., Melis, A. & Happe, T. A novel screening protocol for the isolation of hydrogen producing Chlamydomonas reinhardtii strains. BMC Plant Biol. 8, 107 (2008).PubMed 
PubMed Central 

Google Scholar 
Crawford, D. W. et al. Influence of zinc and iron enrichments on phytoplankton growth in the northeastern subarctic Pacific. Limnol. Oceanogr. 48, 1583–1600 (2003).CAS 

Google Scholar 
Provasoli, L. Media and prospects for the cultivation of marine algae. In Cultures and Collections of Algae. Proc. US-Japan Conference, Hakone, 12-15 September 1966 (eds Watanabe, A & Hattori, A.) 63–75 (Japanese Society of Plant Physiology, 1968).Ranallo-Benavidez, T. R., Jaron, K. S. & Schatz, M. C. GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes. Nat. Commun. 11, 1432 (2020).Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ye, C. X. et al. DBG2OLC: efficient assembly of large genomes using long erroneous reads of the third generation sequencing technologies. Sci. Rep. 6, 31900 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Chin, C. S. et al. Phased diploid genome assembly with single-molecule real-time sequencing. Nat. Methods 13, 1050–1054 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Qin, M. et al. LRScaf: improving draft genomes using long noisy reads. BMC Genomics 20, 955 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ruan, J. & Li, H. Fast and accurate long-read assembly with wtdbg2. Nat. Methods 17, 155–158 (2020).CAS 
PubMed 

Google Scholar 
Servant, N. et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biol. 16, 259 (2015).PubMed 
PubMed Central 

Google Scholar 
Ellinghaus, D., Kurtz, S. & Willhoeft, U. LTRharvest, an efficient and flexible software for de novo detection of LTR retrotransposons. BMC Bioinf. 9, 18 (2008).
Google Scholar 
Xu, Z. & Wang, H. LTR_FINDER: an efficient tool for the prediction of full-length LTR retrotransposons. Nucleic Acids Res. 35, W265–W268 (2007).PubMed 
PubMed Central 

Google Scholar 
Ou, S. J. & Jiang, N. LTR_retriever: a highly accurate and sensitive program for identification of long terminal repeat retrotransposons. Plant Physiol. 176, 1410–1422 (2018).CAS 
PubMed 

Google Scholar 
Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinf. 11, 119 (2010).
Google Scholar 
Stanke, M., Schöffmann, O., Morgenstern, B. & Waack, S. Gene prediction in eukaryotes with a generalized hidden Markov model that uses hints from external sources. BMC Bioinf. 7, 62 (2006).
Google Scholar 
Kim, D. et al. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14, R36 (2013).PubMed 
PubMed Central 

Google Scholar 
Trapnell, C. et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Haas, B. J. et al. Automated eukaryotic gene structure annotation using EVidenceModeler and the program to assemble spliced alignments. Genome Biol. 9, R7 (2008).PubMed 
PubMed Central 

Google Scholar 
Mitchell, A. L. et al. InterPro in 2019: improving coverage, classification and access to protein sequence annotations. Nucleic Acids Res. 47, D351–D360 (2018).PubMed Central 

Google Scholar 
Rice, P., Longden, I. & Bleasby, A. EMBOSS: The European Molecular Biology Open Software Suite. Trends Genet. 16, 276–277 (2000).CAS 
PubMed 

Google Scholar 
Kimura, M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol. 16, 111–120 (1980).CAS 
PubMed 

Google Scholar 
Emms, D. M. & Kelly, S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 20, 238 (2019).PubMed 
PubMed Central 

Google Scholar 
Zhang, Z. et al. KaKs_Calculator: calculating Ka and Ks through model selection and model averaging. Genom. Proteom. Bioinf. 4, 259–263 (2006).CAS 

Google Scholar 
Kumar, S., Stecher, G., Suleski, M. & Hedges, S. B. TimeTree: a resource for timelines, timetrees, and divergence times. Mol. Biol. Evol. 34, 1812–1819 (2017).CAS 
PubMed 

Google Scholar 
Yang, Z. H. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586–1591 (2007).CAS 
PubMed 

Google Scholar 
Wiśniewski, J. R., Hein, M. Y., Cox, J. & Mann, M. A ‘proteomic ruler’ for protein copy number and concentration estimation without spike-in standards. Mol. Cell Proteom. 13, 3497–3506 (2014).
Google Scholar 
Huntemann, M. et al. The standard operating procedure of the DOE-JGI Metagenome Annotation Pipeline (MGAP v.4). Stand. Genomic Sci. 10, 86 (2016).
Google Scholar 
Fu, L. et al. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28, 3150–3152 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Bushnell, B. BBMap: A Fast, Accurate, Splice-aware Aligner (Lawrence Berkeley National Laboratory, 2014).Löytynoja, A. Phylogeny-aware Alignment with PRANK: Multiple Sequence Alignment Methods (Humana Press, 2014).Castresana, J. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17, 540–552 (2000).CAS 
PubMed 

Google Scholar  More

Willcox, G., Fornite, S. & Herveux, L. Early Holocene cultivation before domestication in northern Syria. Veg. Hist. Archaeobot. 17, 313–325 (2008).Article 

Google Scholar 
Fuller, D. Q., Willcox, G. & Allaby, R. G. Cultivation and domestication had multiple origins: arguments against the core area hypothesis for the origins of agriculture in the Near East. World Archaeol. 43, 628–652 (2011).Article 

Google Scholar 
Ibáñez, J. J., Anderson, P. C., González-Urquijo, J. & Gibaja, J. Cereal cultivation and domestication as shown by microtexture analysis of sickle gloss through confocal microscopy. J. Archaeol. Sci. 73, 62–81 (2016).Article 

Google Scholar 
Weiss, E., Kislev, M. E. & Hartmann, A. Autonomous cultivation before domestication. Science 312, 1608–1610 (2006).CAS 
PubMed 
Article 

Google Scholar 
Willcox, G. Measuring grain size and identifying Near Eastern cereal domestication: evidence from the Euphrates Valley. J. Archaeol. Sci. 31, 145–150 (2004).Article 

Google Scholar 
White, C. E. & Makarewicz, C. A. Harvesting practices and early Neolithic barley cultivation at el-Hemmeh, Jordan. Veg. Hist. Archaeobot. 21, 85–94 (2012).Article 

Google Scholar 
Colledge, S., Conolly, J., Finlayson, B. & Kuijt, I. New insights on plant domestication, production intensification, and food storage: the archaeobotanical evidence from PPNA Dhra‘. Levant 50, 14–31 (2018).Article 

Google Scholar 
Kuijt, I. & Finlayson, B. Evidence for food storage and predomestication granaries 11,000 years ago in the Jordan Valley. Proc. Natl Acad. Sci. USA 106, 10966–10970 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Willcox, G. & Stordeur, D. Large-scale cereal processing before domestication during the tenth millennium cal bc in northern Syria. Antiquity 86, 99–114 (2012).Article 

Google Scholar 
Colledge, S. in The Origins of Agriculture and Crop Domestication (eds Damania, A. B. et al.) 121–131 (ICARDA, 1998).Hillman, G. C., Hedges, R., Moore, A. M. T., Colledge, S. & Pettitt, P. New evidence of Lateglacial cereal cultivation at Abu Hureyra on the Euphrates. Holocene 11, 383–393 (2001).Article 

Google Scholar 
Willcox, G. Searching for the origins of arable weeds in the Near East. Veg. Hist. Archaeobot. 21, 163–167 (2012).Article 

Google Scholar 
Snir, A. et al. The origin of cultivation and proto-weeds, long before neolithic farming. PLoS ONE 10, e0131422 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Harris, D. R. & Fuller, D. Q. in Encyclopedia of Global Archaeology (ed. Smith, C.) 104–113 (Springer, 2014).Grime, J. P., Hodgson, J. G. & Hunt, R. Comparative Plant Ecology: A Functional Approach to Common British Species (Springer, 2014).Harlan, J. R., de Wet, J. M. J. & Price, E. G. Comparative evolution of cereals. Evolution 27, 311–325 (1973).PubMed 
Article 

Google Scholar 
Fuller, D. Q. Contrasting patterns in crop domestication and domestication rates: recent archaeobotanical insights from the Old World. Ann. Bot. 100, 903–924 (2007).PubMed 
PubMed Central 
Article 

Google Scholar 
Asouti, E. in Neolithic Corporate Identities. Studies in Early Near Eastern Production, Subsistence and Environment 20 (eds Benz, M. et al.) 21–53 (Ex oriente, 2017).Harris, D. R. in Foraging and Farming: the Evolution of Plant Exploitation (eds Harris, D. R. & Hillman, G.) 11–26 (Unwin Hyman, 1989).Smith, B. D. Low-level food production. J. Archaeol. Res. 9, 1–43 (2001).Article 

Google Scholar 
Rindos, D. The Origins of Agriculture: an Evolutionary Perspective (Academic, 1984).Weide, A. Towards a socio-economic model for southwest Asian cereal domestication. Agronomy 11, 2432 (2021).Article 

Google Scholar 
Hillman, G. C. & Davies, M. S. Measured domestication rates in wild wheats and barley under primitive cultivation, and their archaeological implications. J. World Prehist. 4, 157–222 (1990).Article 

Google Scholar 
Kislev, M. E., Hartmann, A. & Weiss, E. Impetus for sowing and the beginning of agriculture: ground collecting of wild cereals. Proc. Natl Acad. Sci. USA 101, 2692–2695 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Weide, A. et al. The association of arable weeds with modern wild cereal habitats: implications for reconstructing the origins of plant cultivation in the Levant. Environ. Archaeol. https://doi.org/10.1080/14614103.2021.1882715 (2021).Zohary, M. The segetal plant communities of Palestine. Vegetatio 2, 387–411 (1950).Article 

Google Scholar 
Abbo, S., Lev-Yadun, S. & Gopher, A. Plant domestication and crop evolution in the Near East: on events and processes. Crit. Rev. Plant Sci. 31, 241–257 (2012).Article 

Google Scholar 
Wood, D. & Lenné, J. M. A natural adaptive syndrome as a model for the origins of cereal agriculture. Proc. R. Soc. Lond. B 285, 20180277 (2018).
Google Scholar 
Bogaard, A., Palmer, C., Jones, G., Charles, M. & Hodgson, J. G. A FIBS approach to the use of weed ecology for the archaeobotanical recognition of crop rotation regimes. J. Archaeol. Sci. 26, 1211–1224 (1999).Article 

Google Scholar 
Jones, G., Bogaard, A., Charles, M. & Hodgson, J. G. Distinguishing the effects of agricultural practices relating to fertility and disturbance: a functional ecological approach in archaeobotany. J. Archaeol. Sci. 27, 1073–1084 (2000).Article 

Google Scholar 
Díaz, S. et al. The global spectrum of plant form and function. Nature 529, 167–171 (2016).PubMed 
Article 
CAS 

Google Scholar 
Garnier, E., Navas, M.-L. & Grigulis, K. Plant Functional Diversity: Organism Traits, Community Structure, and Ecosystem Properties (Oxford Univ. Press, 2016).Bogaard, A. Neolithic Farming in Central Europe (Routledge, 2004).Bogaard, A. et al. From traditional farming in Morocco to early urban agroecology in northern Mesopotamia: combining present-day arable weed surveys and crop isotope analysis to reconstruct past agrosystems in (semi-)arid regions. Environ. Archaeol. 23, 303–322 (2018).Article 

Google Scholar 
Hamerow, H. et al. An integrated bioarchaeological approach to the medieval ‘agricultural revolution’: a case study from Stafford, England, c. ad 800–1200. Eur. J. Archaeol. 23, 585–609 (2020).Article 

Google Scholar 
Green, L., Charles, M. & Bogaard, A. Exploring the agroecology of Neolithic Çatalhöyük, Central Anatolia: an archaeobotanical approach to agricultural intensity based on functional ecological analysis of arable weed flora. Paléorient 44, 29–44 (2018).
Google Scholar 
Green, L. Assessing the Nature of Early Farming in Neolithic Western Asia: A Functional Ecological Approach to Emerging Arable Weeds. Univ. of Oxford (2017).Atran, S. Hamula organisation and masha’a tenure in Palestine. Man 21, 271–295 (1986).Article 

Google Scholar 
Palmer, C. ‘Following the plough’: the agricultural environment of northern Jordan. Levant 30, 129–165 (1998).Article 

Google Scholar 
Håkansson, S. in Biology and Ecology of Weeds (eds Holzner, W. & Numata, M.) 123–135 (Springer Netherlands, 1982).Charles, M., Bogaard, A., Jones, G., Hodgson, J. & Halstead, P. Towards the archaeobotanical identification of intensive cereal cultivation: present-day ecological investigation in the mountains of Asturias, northwest Spain. Veg. Hist. Archaeobot. 11, 133–142 (2002).Article 

Google Scholar 
Hartmann-Shenkman, A., Kislev, M. E., Galili, E., Melamed, Y. & Weiss, E. Invading a new niche: obligatory weeds at Neolithic Atlit-Yam, Israel. Veg. Hist. Archaeobot. 24, 9–18 (2015).Article 

Google Scholar 
Kuijt, I. in The Neolithic Demographic Transition and its Consequences (eds Bocquet-Appel, J.-P. & Bar-Yosef, O.) 287–313 (Springer Netherlands, 2008).Bogaard, A. et al. Private pantries and celebrated surplus: storing and sharing food at Neolithic Çatalhöyük, Central Anatolia. Antiquity 83, 649–668 (2009).Article 

Google Scholar 
Jones, G. et al. The origins of agriculture: intentions and consequences. J. Archaeol. Sci. 125, 105290 (2021).Article 

Google Scholar 
Weiss, E., Kislev, M. E., Simchoni, O., Nadel, D. & Tschauner, H. Plant-food preparation area on an Upper Paleolithic brush hut floor at Ohalo II, Israel. J. Archaeol. Sci. 35, 2400–2414 (2008).Article 

Google Scholar 
Kluyver, T. A., Charles, M., Jones, G., Rees, M. & Osborne, C. P. Did greater burial depth increase the seed size of domesticated legumes? J. Exp. Bot. 64, 4101–4108 (2013).CAS 
PubMed 
Article 

Google Scholar 
Preece, C., Jones, G., Rees, M. & Osborne, C. P. Fertile Crescent crop progenitors gained a competitive advantage from large seedlings. Ecol. Evol. 11, 3300–3312 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Halstead, P. Two Oxen Ahead: Pre-mechanized Farming in the Mediterranean (Wiley, 2014).Anderson, P. C. in The Origins of Agriculture and Crop Domestication (eds Damania, A. B. et al.) 145–159 (ICARDA, 1998).Mercuri, A. M., Fornaciari, R., Gallinaro, M., Vanin, S. & di Lernia, S. Plant behaviour from human imprints and the cultivation of wild cereals in Holocene Sahara. Nat. Plants 4, 71–81 (2018).PubMed 
Article 

Google Scholar 
Spengler, R. N. & Mueller, N. G. Grazing animals drove domestication of grain crops. Nat. Plants 5, 656–662 (2019).PubMed 
Article 

Google Scholar 
Smith, B. D. General patterns of niche construction and the management of ‘wild’ plant and animal resources by small-scale pre-industrial societies. Phil. Trans. R. Soc. Lond. B 366, 836–848 (2011).Article 

Google Scholar 
Bogaard, A. et al. Reconsidering domestication from a process archaeology perspective. World Archaeol. https://doi.org/10.1080/00438243.2021.1954990 (2021).Coqueugniot, E. in Espace Naturel, Espace Habité En Syrie Du Nord (10e–2e millénaires av. J.-C.) (eds M. Fortin & O. Aurenche) 109–114 (Maison de l’Orient et de la Méditerranée, 1998).Douché, C. Émergence et développement des sociétés agricoles au Néolithique acéramique (Xe-VIIIe millénaires av. n. ère) étude archéobotanique de Dja’de El-Mughara et Tell Aswad, Syrie. PhD thesis (Archaeological Mission of Dja’de el Mughara, 2018).Noy, T. Gilgal I: a pre-pottery Neolithic site, Israel. The 1985–1987 seasons. Paléorient 15, 11–18 (1989).Article 

Google Scholar 
Bar-Yosef, O. & Gopher, A. in An Early Neolithic Village in the Jordan Valley (eds Bar-Yosef, O. & Gopher, A.) 41–69 (Harvard Univ., 1997).Wright, K. I. The social origins of cooking and dining in early villages of western Asia. Proc. Prehist. Soc. 66, 89–121 (2000).Article 

Google Scholar 
Finlayson, B. Egalitarian societies and the earliest Neolithic of southwest Asia. Prehist. Archaeol. J. Interdiscip. Stud. 3, 27–43 (2020).Article 

Google Scholar 
Bowles, S. & Choi, J.-K. The Neolithic agricultural revolution and the origins of private property. J. Polit. Econ. 127, 2186–2228 (2019).Article 

Google Scholar 
Kuijt, I. The Neolithic refrigerator on a Friday night: how many people are coming to dinner and just what should I do with the slimy veggies in the back of the fridge? Environ. Archaeol. 20, 321–336 (2015).Article 

Google Scholar 
Danin, A. Flora and vegetation of Israel and adjacent areas. Zoogeogr. Isr. 30, 251–276 (1988).
Google Scholar 
Noy-Meir, I., Gutman, M. & Kaplan, Y. Responses of Mediterranean grassland plants to grazing and protection. J. Ecol. 77, 290–310 (1989).Article 

Google Scholar 
Noy-Meir, I. The effect of grazing on the abundance of wild wheat, barley and oat in Israel. Biol. Conserv. 51, 299–310 (1990).Article 

Google Scholar 
Jones, G., Bogaard, A., Halstead, P., Charles, M. & Smith, H. Identifying the intensity of crop husbandry practices on the basis of weed floras. Annu. Br. Sch. Athens 94, 167–189 (1999).Article 

Google Scholar 
Sternberg, M., Gutman, M., Perevolotsky, A., Ungar, E. D. & Kigel, J. Vegetation response to grazing management in a Mediterranean herbaceous community: a functional group approach. J. Appl. Ecol. 37, 224–237 (2000).Article 

Google Scholar 
Sternberg, M. et al. Testing the limits of resistance: a 19-year study of Mediterranean grassland response to grazing regimes. Glob. Change Biol. 21, 1939–1950 (2015).Article 

Google Scholar 
Calev, A. et al. High-intensity thinning treatments in mature Pinus halepensis plantations experiencing prolonged drought. Eur. J. For. Res. 135, 551–563 (2016).Article 

Google Scholar 
Osem, Y., Perevolotsky, A. & Kigel, J. Grazing effect on diversity of annual plant communities in a semi-arid rangeland: interactions with small-scale spatial and temporal variation in primary productivity. J. Ecol. 90, 936–946 (2002).Article 

Google Scholar 
Temper, L. Creating facts on the ground: agriculture in Israel and Palestine (1882–2000). Hist. Agrar. 48, 75–110 (2009).
Google Scholar 
Dan, J., Yaalon, D., Koyumdjisky, H. & Raz, Z. The soil association map of Israel (1:1,000,000). Isr. J. Earth Sci. 21, 29–49 (1970).
Google Scholar 
Sans, F. X. & Masalles, R. M. Phenological patterns in an arable land weed community related to disturbance. Weed Res. 35, 321–332 (1995).Article 

Google Scholar 
Zohary, M. & Feinbrun-Dothan, N. Flora Palaestina Vol. 1–4 (Israel Academy of Sciences and Humanities, 1966).Davis, P. Flora of Turkey and the East Aegean Islands Vol. 1–10 (Edinburgh Univ. Press, 1965).Mortimer, A. M. in Weed Control Handbook: Principles (eds Hance, R. J. & Holly, K.) 1–42 (Blackwell, 1990).Douché, C. & Willcox, G. New archaeobotanical data from the Early Neolithic sites of Dja’de el-Mughara and Tell Aswad (Syria): a comparison between the northern and the southern Levant. Paléorient 44, 45–58 (2018).
Google Scholar 
Jones, G. The application of present-day cereal processing studies to charred archaeobotanical remains. Circaea 6, 91–96 (1990).
Google Scholar 
Bogaard, A., Jones, G. & Charles, M. The impact of crop processing on the reconstruction of crop sowing time and cultivation intensity from archaeobotanical weed evidence. Veg. Hist. Archaeobot. 14, 505–509 (2005).Article 

Google Scholar 
Bogaard, A. et al. in Humans and Landscapes of Çatalhöyük: Reports from the 2000–2008 Seasons (ed. Hodder, I.) 93–128 (Cotsen Institute of Archaeology/British Institute at Ankara, 2013).Filipović, D. Early Farming in Central Anatolia: an Archaeobotanical Study of Crop Husbandry, Animal Diet and Land Use at Neolithic Çatalhöyük (British Archaeological Reports, 2014).Helmer, D. et al. in New Methods and the First Steps of Mammal Domestication (eds Vigne, J.-D. et al.) 86–95 (Oxbow Books, 2005).Charles, M. Fodder from dung: the recognition and interpretation of dung-derived plant material from archaeological sites. Environ. Archaeol. 1, 111–122 (1998).Article 

Google Scholar 
Kislev, M. E. in An Early Neolithic Village in the Jordan Valley (eds Ofer Bar-Yosef & Avi Gopher) 209–236 (Harvard Univ., 1997).Kislev, M. E. et al. in Gilgal: Early Neolithic Occupations in the Lower Jordan Valley. The Excavations of Tamar Noy (eds Bar-Yosef, O. et al.) 251–257 (Oxbow Books, 2010).Snir, A., Nadel, D. & Weiss, E. Plant-food preparation on two consecutive floors at Upper Paleolithic Ohalo II, Israel. J. Archaeol. Sci. 53, 61–71 (2015).Article 

Google Scholar 
Jones, G., Charles, M., Bogaard, A. & Hodgson, J. Crops and weeds: the role of weed functional ecology in the identification of crop husbandry methods. J. Archaeol. Sci. 37, 70–77 (2010).Article 

Google Scholar 
Šmilauer, P. & Lepš, J. Multivariate Analysis of Ecological Data Using CANOCO 5 (Cambridge Univ. Press, 2014).Galili, E. et al. Atlit-Yam: a Prehistoric site on the sea floor off the Israeli coast. J. Field Archaeol. 20, 133–157 (1993).
Google Scholar 
Brenet, M., Sanchez-Priego, J. & Ibáñez-Estévez, J. J. in Préhistoire et Approche Expérimentale (eds Bourguignon, L. et al.) 121–164 (Monique Mergoil, 2001).Bar-Yosef, O., Gopher, A., Goring-Morris, A. N. & Kozlowski, S. K. in Gilgal: Early Neolithic Occupations in the Lower Jordan Valley. The Excavations of Tamar Noy (eds Bar-Yosef, O. et al.) 11–26 (Oxbow Books, 2010). More

Lattice data geoprocessing and temporal extentWe latticed the data49 using a worldwide grid composed of 18,874 hexagonal 7774 km2 units, built using Discrete Global for R (https://github.com/r-barnes/dggridR)50. All the information we processed on yellow fever cases, on urban and sylvatic vectors presences, and on zoogeographic, spatial and environmental variables (see details on this information below) was aggregated at this spatial resolution. We used zonal statistics to calculate average variable values using ArcMAP 10.7.The temporal extent for our analysis was divided into three periods: the late 20th century (1970–2000), the early 21st century (2001–2017), and the period 2018–2020. Predictions estimated by the late 20th century models were validated using cases reported in the early 21st century, and predictions from the early 21st century models were validated with records from 2018‒2020. Although the limit between periods at the turn of the century is arbitrary, it reflects: 1) Distributional changes in the ranges of the Ae. aegypti and Ae. albopictus vectors51; 2) after 1999, the yellow fever genotype I has spread outside the endemic regions, and the genotype I modern-lineage has caused all major yellow fever outbreaks detected in non-endemic regions of South America since 200013; 3) the maximum potential of globalization was realised at the beginning of the 21st century with the opening of international borders, the widespread access to the Internet and to cell phones, and the generalization of online travel booking and of low-cost flights34. The end of the second period, 2017, was chosen in order to include three years with occurrence of yellow fever cases in south-western Brazil (and two since its occurrence in Angola and the DRC), while retaining three later years for predictive testing purposes (details on this testing are given below).Yellow fever datasetsWe used georeferenced cases of yellow fever in humans for a period of 51 years (from 1970 to 2020). This study period starts immediately after the suspension of the use of DDT due to to the appearance of resistance of Ae. aegypti in the late 1960s in several countries, after 50 years of eradication efforts10. We took from Shearer et al.6 the distribution of yellow fever cases for the period 1970–2016. We extracted additional cases for the period 1970–2020 from various sources (Supplementary data 1), including ProMED-mail: Program of International society for infectious diseases; World Health Organization (WHO): Yellow fever outbreak weekly situation reports, Rapport de situation fievre jaune en RD Congo and Weekly epidemiological record; Health Ministry of different countries: Epidemiological Bulletins of yellow fever in Brazil, Peru, Colombia, and Paraguay; Pan American Health Organization (PAHO): Epidemiological Update Yellow Fever; European Centre for Disease Prevention and Control (ECDC): Communicable disease threats report and Rapid risk assessment report; Nigeria Centre for Disease Control (NCDC): Situation report, yellow fever outbreak in Nigeria and Global Infectious Disease and Epidemiology Online Network (GIDEON). The reported cases were complemented with publications available since 2016 with geo-referenced information on case location (Supplementary data 1). In addition, information was also sought on cases reported in French and Portuguese from local news reports in Africa.We only represented in the hexagonal lattice the reported cases of yellow fever that had a precise location or that were referred to administrative unit was smaller than or of similar size to the hexagons. This dataset was transformed into a binary variable per study period representing the presence (n = 218 hexagons in the late 20th century; 493 hexagons in the early 21st century, see Supplementary data 2) or absence (n = 18,656 hexagons in the late 20th century; 18,381 hexagons in the early 21st century), hereafter the distribution of reported cases of yellow fever.Vector datasetThe georeferenced presences of vectors involved in the urban cycle of yellow fever (i.e., the mosquitoes Ae. aegypti and Ae. albopictus) were taken from “The global compendium of the Ae. aegypti and Ae. Albopictus occurrence”26 for the period 1970–2014. We complemented these records with georeferenced data scientifically validated for the period 2014–2017, taken from VectorBase (https://www.vectorbase.org/) and Mosquito Alert (http://www.mosquitoalert.com/). We included both species because, although Ae. Aegypti is the main vector of yellow fever, Ae. albopictus can also transmit the yellow fever virus to humans4,52.In addition, we included georeferenced occurrence data of sylvatic vectors (Haemagogus janthinomys, H. leucocelaenus and Sabethes chloropterus in South America; Ae. africanus and Ae. vittatus in Africa), which were obtained from Vectormap (vectormap.si.edu) and Gbif (https://gbif.org).We represented in the hexagonal lattice the reported occurrence of mosquitoes that had a precise location or were located in administrative smaller than or of similar size to the hexagons. With this information, we built binary variables representing the presence or absence of each mosquito species in each hexagon. For species involved in the urban cycle, we built two binary variables per species: one for the late 20th century, and another for the early 21st century. For species involved in the sylvatic cycle, we merged the data of late 20th century and early 21st century in order to build a binary variable per species, due the scarcity of data and under the assumption that their distributions have been stable during the four last decades53,54,55.Zoogeographic, spatial and environmental variablesWe built zoogeographic variables based on chorotypes, or types of distribution ranges, of all non-human primate species, as all are potentially vulnerable to yellow fever56. A chorotype is a distribution pattern shared by a group of species57. For obtaining these zoogeographic variables, we proceeded in 4 steps: (1) Distribution maps of non-human primates were obtained from the IUCN for South-America and Africa; (2) the species ranges were classified hierarchically using the classification algorithm UPGMA according to the Baroni-Urbani & Buser´s similarity index58; (3) we evaluated the statistical significance of all clusters obtained as a result of the classification using RMacoqui 1.0 software (http://rmacoqui.r-forge.r-project.org/)59; (4) in each hexagon, the number of species belonging to each chorotype was quantified. We generated a zoogeographic model based on the non-human primates chorotypes by running a forward-backward stepwise logistic regression using presence/absence of yellow fever cases and the number of species of each chorotype as dependent and predictor variables, respectively. This procedure was made for two periods: late 20th century and early 21st century. Henceforth, only the selected chorotype variables were considered in the baseline disease favourability models explained below.We built a yellow fever spatial variable for each continent (South-America and Africa), which were calculated through the trend surface approach, by performing a backward-stepwise logistic regression of the distribution of yellow fever cases on a ensemble of variables defined for polynomial combinations of longitude (X) and latitude (Y) up to the third degree: X, Y, XY, X2, Y2, X2Y, XY2, X3, and Y3. We transformed probability values derived from logistic regression into spatial favourability values by applying the Favourability Function60,61, using the following equation:$$F=frac{P}{1-P}Big/left(frac{{n}_{1}}{{n}_{0}}+frac{P}{1-P}right)$$
(1)
where P is the spatial probability of occurrence of at least a case of yellow fever at each hexagon, and n1 and n0 are the numbers of hexagons with presence and absence of yellow fever cases, respectively. We built a different spatial variable for each continent and time period.We used environmental variables related to the following factors: climate, human activity, topography, hydrography, biome, ecosystem type, and forest loss. For details about the source and description of the environmental variables selected, see Supplementary Table 3.Pathogeographical approach to transmission risk modellingOur objectives were to construct a global yellow fever transmission risk map, and to identify areas where primates contribute to increased risk, using the methodology previously used to analyse the worldwide dynamic biogeography of zoonotic and anthroponotic dengue34 (see flowchart in Fig. 1 and Supplementary Methods). We produced a transmission model focused on the late 20th century and another for the early 21st century.The risk of transmission was assessed by combining a first model describing areas favourable to the presence of yellow fever, i.e., the “baseline disease model”; and another model describing areas favourable to the presence of mosquitoes known to act as vectors, i.e., the “vector model”. For this combination, we used the fuzzy intersection62, i.e., the risk of transmission at each hexagon was valued at the minimum between favourability in the baseline disease model and favourability in the vector model.In this way, we considered that the vectors are a limiting factor, and that the risk of transmission derives from the degree to which the environmental conditions are simultaneously favourable for the presence of vectors and for disease cases to occur63. In order to analyze the spatio-temporal dynamic of yellow fever, we made comparable models for the late 20th century and the early 21st century, using predictor variables that are available for both periods. Later, we made a 21st-century enhanced model that optimized the predictive capacity of availabe information in the search for current risk areas. For this purpose, we included, in the variable set, predictors that are only accessible for the 21st century (e.g., high-resolution population density, livestock, irrigation, infrastructures, intact forest, and GlobCover land cover classes; see Supplementary Table 3).Baseline disease modelsThe baseline disease model in the late 20th century was expressed in terms of favourability values, using the Eq. (1) (see above). This time, P was calculated through a multivariable forward-backward stepwise logistic regression of the 20th-century yellow fever presences/absences on a set of zoogeographic, environmental and spatial variables. This was made in two blocks: 1) a stepwise selection of environmental and spatial variables; 2) a later stepwise addition of chorotypes whose presence contribute to improve significantly the likelihood of the model based only on the first block. Variables for each block were preselected using RAO´s score tests (which estimated the significance of its association to the distribution of yellow fever cases), and Benjamini and Hochberg´s (1995) false discovery rate (FDR) to control for Type I errors, which could pass due to the number of variables analysed. We also avoided excesive multicollinearity by preventing that variables with Spearman correlation values >0.8 were included in the same model. In case this happened, only the variable with the most significant RAO´s score-test value was retained, and the multivariable model was re-run. The parameters in the models were estimated using a gradient ascent machine learning algorithm, and the significance of these paremeters was assessed using the test of Wald. The goodness of fit of the models was established using the test of Hosmer and Lemeshow, which checks the significance of the difference between expected and observed values, so that non significant differences mean that the fit is good. We used IBM-SPSS Statistics 24 software package to perform the models and all the associated tests.We subsequently updated the baseline disease model to explain the distribution of yellow fever cases in the early 21st century. Compared to the procedure described for the 20th-century model, we included a new block before the two ones mentioned above. Thus, the methodological sequence was as follows: (1) forcing the input, as a predictor variable, of the logit of the late 20th century baseline disease model (the logit being the linear combination of variables in the 20th-century model); (2) making a later stepwise selection of spatial and environmental variables; and (3) a stepwise addition of chorotypes that contribute to improving the model’s likelihood. In this way, we took into account that the current spread of yellow fever is influenced by the inertia of previous situations. This is equivalent to assuming that there is temporal autocorrelation (i.e., disease cases in the early 21st century are more probable to occur in areas where they already occurred in the late 20th century). In the 21st-century model, the variables entering in blocks (2) and (3) represent the drivers potentially favouring the spread34. The preselection of variables for blocks (2) and (3) and the control for excessive multicollinearity between environmental variables were made as explained for the late 20th-century model.Vector modelsWe produced a favourabuility model for each vector species for the late 20th century and for the early 21st century separately. We built multivariable favourability models for urban vectors using the distribution of each urban mosquito species in the late 20th century and the spatial and environmental variables for the late 20th century, following the same procedure used for block (1) in the 20th-century baseline disease model. We also updated each urban vector model for the early 21st century as in the baseline disease model, using the procedure described for blocks (1) and (2).A single model, referred to both the late 20th and the early 21st centuries, was made for sylvatic vectors, for the reasons explained above. Finally, we built up the vector models for the late 20th century and for the early 21st century by joining all individual vector models of each period using the fuzzy union64 (i.e., considering for each hexagon the maximum value shown by any of the species models). This criterion was taken into account because, if the pathogen were present, the occurrence of a single vector species would involve potential for yellow fever transmission.Model fit assessment and validation of its predictive capacityFavourability models were assessed according to their classification and discrimination capacities respect to the training data set (i.e., to the observations used for model training). The classification capacity was based on two classification thresholds: F = 0.5, which represents the neutral favourability, and F = 0.2, below which the risk of transmission was considered to be low61. Six classification assessment indices were used65: (1) sensitivity (i.e., proportion of presences correctly classified in favourable hexagons), (2) specificity (i.e., proportion of absences correctly classified in unfavourable hexagons), (3) CCR (i.e., proportion of presences and absences correctly classified in favourable hexagons respectively), (4) TSS (that is sensitivity + specifity – 1), (5) underprediction rate (i.e., proportion of favourable areas that are recorded to have presences), and (6) overprediction rate (i.e., proportion of favourable areas that are not recorded to have presences). The discrimination capacity was assessed using the area under the receiver operating characteristic (ROC) curve (AUC)66.The validation of the predictive capacity of the late 20th century disease and transmission-risk models was done by evaluating, with the same indices used above, classification and discrimination capacities with respect to the cases of the period 2001‒2020. The predictive capacity of the models for the early 21st century was validated with respect to the yellow fever cases reported in the period 2018‒2020.Relative importance of the zoogeographical factorWe estimated the pure contribution of non-human primates to the baseline disease model, i.e., how much of the variation in favourability for yellow fever cases was explained exclusively by the zoogeographical factor, by performing a variation partitioning analysis67. This implied the use of the zoogeographic model and a spatio-environmental model constructed with the environmental and spatial variables that entered the baseline disease model. This approach also allowed us to calculate how much is the variation of the baseline disease model attributable simultaneously to the zoogeographical and other factors. We built maps identifying the zones where the non-human primates could increase yellow fever cases in humans, that is, where the presence of primates could favour the occurrence of yellow fever regardless of correlations with other factors. To map these areas we identified the hexagons that fulfilled these conditions: 1) favourability values for the baseline disease model were ≥ 0.2; and 2) the difference between the favourability values provided by the baseline disease model and the spatio-environmental model was positive and ≥ 0.01.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

General observation and methods for ploidy analysis on aphid bacteriome cellsConsistent with previous observations9,21,22,40, the bacteriome of viviparous aphids consisted of two types of cells: bacteriocytes and sheath cells (Fig. 2). Bacteriocytes contained Buchnera cells and were much larger than sheath cells. Sheath cells exhibited a flattened morphology and surrounded the bacteriocytes. Both cell types possessed a single nucleus. Bacteriocytes had a single prominent nucleolus, which was not stained using DAPI, but using “Nucleolus Bright Red” staining (Fig. 2). Most sheath cells also had a single nucleolus, yet a small number had two. “Nucleolus Bright Red” also stained the peripheral region of Buchnera, probably because of the richness of RNA around Buchnera cells.Figure 2Morphology of bacteriocytes and sheath cells from each morph of aphids visualized using DAPI/Phalloidin/Nucleolus Bright Red staining. DNA and F-actin were stained by DAPI (gray or blue) and Phalloidin (green), respectively. The nucleolus, which is the site of ribosome biogenesis, was visualized by Nucleolus Bright Red (red). This dye binds RNA electrostatically, therefore the cytoplasm of bacteriocytes and Buchnera cells were also stained. Bacteriocytes (white arrows) had single prominent nucleolus, and the cell sizes were much larger than sheath cells (white arrowheads) in all aphid morphs.Full size imageTo determine the most suitable methods for ploidy analysis of aphid bacteriocytes, three types of methods, flow cytometry, Feulgen densitometry, and fluorometry were compared. First, flow cytometry successfully detected the nuclei of bacteriome cells and heads, and distinct peaks were present (Fig. S3). There were several peaks, which can be categorized as ploidy classes based on head peaks, assuming that the smallest peaks correspond to a diploid population. We recognized peaks up to 256C (256-ploidy) cells but could not distinguish cell types (i.e., bacteriocytes or sheath cells) in this method due to a lack of cytological information. Note that “C” means haploid genome size, for example, 2C = diploid and 8C = octoploid. Second, Feulgen densitometry also showed several ploidy levels of up to 128C (Fig. S4) in bacteriocytes. Sheath cells mainly consisted of 16-32C cells. However, we found that many cells were lost during the experimental procedures, probably due to the repeated washing processes and the long incubation time.We found the third method, image-based fluorometry for isolated nuclei, the best for quantitative ploidy analysis of aphid bacteriocytes (Fig. 3). Fluorometry showed distinct peaks of integrated fluorescence intensity, and they could be categorized as each ploidy class based on the intensity of the smallest peak in head cells (diploid population). The results were consistent with other methods; ploidy levels were 32C-256C in bacteriocytes and 16C-32C in sheath cells. In this analysis, the nucleolus size was used to discriminate between cell types. During cytological observation, we obtained the size distribution of the nucleolus, and it was revealed that the nucleolus of bacteriocytes was always larger than that of sheath cells (Fig. S5). Based on the results, we determined the threshold of the size of the nucleolus. More specifically, in viviparous females, nuclei that have nucleoli larger than 20 μm2 were categorized into bacteriocytes. Note that the peaks of sheath cells were not distinct or reliable for categorizing their ploidy class; therefore, we showed results focusing on bacteriocytes in the following sections.Figure 3Ploidy analysis of aphid bacteriocytes using DAPI-fluorometry. A representative result from the analysis of adult viviparous females is presented. An image of DAPI-stained nuclei was also shown (the blue channel was extracted). Isolated nuclei of bacteriome cells were stained using DAPI, image-captured with a CCD camera, and their integrated fluorescence intensity was measured using ImageJ software. Nuclei were categorized into “bacteriocytes” or “sheath cells,” based on the size distribution of nucleolus (see “Materials and Methods”). Relative ploidy levels were calculated based on the data from head cells which are mainly diploid. Bacteriocytes of adult viviparous aphids consisted of 16C-256C cells, and 64–128 cells were dominant, while sheath cells exhibited lower ploidy levels (mainly 16C). “C” means haploid genome size, for example, 2C = diploid and 8C = octoploid.Full size imageCellular features of bacteriome cells in viviparous and oviparous females, and malesThe cellular features were generally consistent among young adults (within 5 days of adult eclosion) of three morphs, viviparous and oviparous females, and males (Fig. 2). Nevertheless, Buchnera-absence zones in the cytoplasm of bacteriocytes, which are considered to be degeneration of Buchnera45, and bacteriocytes degeneration46 were both observed more frequently in male bacteriocytes than in females (Fig. 2). The cell size of bacteriocytes was significantly different among morphs (LM with type II test, F = 286.15, df = 2, p  More

Swedish environment minister Annika Strandhäll before the start of the Stockholm +50 Climate Summit. Few world leaders will be attending.Credit: Fredrik Persson/TT News Agency/AFP/Getty

Sustainability is now a household term, but it wasn’t always so.Fifty years ago, the United Nations held its Conference on the Human Environment in Stockholm. This landmark event gave the concept of sustainable development its first international recognition. Sweden and the UN are marking the occasion this week with Stockholm+50, an international meeting that serves as both commemoration and call to action.The world is deep in planetary and human crises, with the UN’s Sustainable Development Goals off track and multilateral agreements on climate change and biodiversity behind schedule. Governments need to integrate sustainability into economic planning — and listen to researchers, who are ready with evidence-based arguments and tools to help them do so.Fifty years ago, the time was ripe for an environmental agenda to enter the world stage. Optimistic ideas of economic growth as a driver of progress, propelled by the Industrial Revolution, needed to accommodate concerns over damage to the natural environment. Books such as Rachel Carson’s Silent Spring (1962) — which raised awareness about harms caused by pesticides — brought scientific information about environmental risks into the mainstream.In March 1972, a team of researchers and policymakers sounded another alarm in The Limits to Growth, one of the first reports to forecast catastrophic consequences if humans kept exploiting Earth’s limited supply of natural resources. The conference in Stockholm followed a few months later, steered to success by its secretary-general, Canadian industrialist Maurice Strong. That set crucial institutions in motion, starting with the establishment of the UN Environment Programme (UNEP), based in Nairobi — the first UN body to be headquartered in a developing country. UNEP went on to facilitate a new international law — the 1987 Montreal Protocol to phase out ozone-depleting substances — and co-founded the Intergovernmental Panel on Climate Change (IPCC). It assisted in establishing the first action plans for sustainable development through landmark international agreements on biodiversity, climate and desertification.But there were mistakes and missed opportunities. The establishment of multiple agencies and policy instruments created a disjointed governance system. Newly created environment ministers wielded little power. In national budgets, environmental protection was siloed away from economic development and social concerns. For a long time, action on climate change remained unfocused. And the economic drivers of environmental change were overlooked.And so, 50 years after that momentous conference, the world remains in crisis. With impending climate and biodiversity crises, the warnings issued by visionaries now hit even closer.Stockholm+50 promises “clear and concrete recommendations and messages for action at all levels”. More than 90 ministers are expected to attend, but only 10 heads of government. That’s a missed opportunity for high-level action. World leaders are needed because their presence signals that sustainability remains at the top of their agendas.Awareness of the need to embed sustainability into policymaking has broken into the mainstream, although much of it is still talk. City governments around the world are implementing ambitious climate action plans through the C40 Cities network. Some companies, too, are adopting sustainability principles, from reporting (and reducing) their carbon footprints to ensuring that investments, as far as possible, do not harm the environment.But this urgency has not ascended to heads of state and government. With a handful of exceptions — such as Finland, Iceland, New Zealand, Scotland and Wales — most nations seem unwilling to systemically integrate their economic, environmental and social policymaking.Doing so is not only good for the environment; it is also sound economics and good for well-being. The food and energy crisis driving poverty and diminishing living standards around the world might have been triggered by the shocks of a pandemic and war on Ukraine — but it is driven just as much by the depletion of natural resources.Ahead of the 1972 conference, 2,200 environmental scientists signed a letter — called the Menton Message — to then UN secretary-general U Thant. The signatories had a sense that the world was moving towards multiple crises. They urged “massive research into the problems that threaten the survival of mankind”, such as hunger, wars, environmental degradation and natural-resource depletion. The UN system went on to play a big part in building the body of knowledge that has shown why sustainability is necessary, and in creating the policy architecture to make it happen. But to do the Stockholm vision justice, there must be bolder action from heads of government and from the UN system. The planned creation of a board of science advisers to UN secretary-general António Guterres needs to be accelerated. Once established, the board must find a way to bring joined-up action on sustainability closer to world leaders.Researchers can now join a successor to the Menton Message that has been organized by the International Science Council, the global science network Future Earth and the Stockholm Environment Institute. In an open letter addressed to world citizens, the authors write: “After 50 years, pro-environmental action seems like one step forward and two back. The world produces more food than needed, yet many people still go hungry. We continue to subsidize and invest in fossil fuels, even though renewable energy is increasingly cost-effective. We extract resources where the price is lowest, often in direct disregard of local rights and values.”World leaders must listen to the research community, and accept the evidence and narrative offered to help them to navigate meaningful change. Environmental sustainability does not impede prosperity and well-being — in fact, it is vital to them. People in power need to sit up and take notice. More

Di Giacomo, G., Hadrich, J., Hutchison, W. D., Peterson, H. & Rogers, M. Economic impact of spotted wing drosophila (Diptera: Drosophilidae) yield loss on minnesota raspberry farms: A grower survey. J. Integr. Pest Manag. https://doi.org/10.1093/jipm/pmz006 (2019).Article 

Google Scholar 
Farnsworth, D. et al. Economic analysis of revenue losses and control costs associated with the spotted wing drosophila, Drosophila suzukii (Matsumura), in the California raspberry industry. Pest. Manag. Sci. 73, 1083–1090. https://doi.org/10.1002/ps.4497 (2017).CAS 
Article 
PubMed 

Google Scholar 
Beers, E. H., Van Steenwyk, R. A., Shearer, P. W., Coates, W. W. & Grant, J. A. Developing Drosophila suzukii management programs for sweet cherry in the western United States. Pest Manag. Sci. 67, 1386–1395. https://doi.org/10.1002/ps.2279 (2011).CAS 
Article 
PubMed 

Google Scholar 
Tait, G. et al. Drosophila suzukii (Diptera: Drosophilidae): A decade of research towards a sustainable integrated pest management program. J. Econ. Entomol. 114, 1950–1974. https://doi.org/10.1093/jee/toab158 (2021).Article 
PubMed 

Google Scholar 
Daane, K. M. et al. First exploration of parasitoids of Drosophila suzukii in South Korea as potential classical biological agents. J. Pest Sci. 89, 823–835. https://doi.org/10.1007/s10340-016-0740-0 (2016).ADS 
Article 

Google Scholar 
Abram, P. K. et al. New records of Leptopilina, Ganaspis, and Asobara species associated with Drosophila suzukii in North America, including detections of L. japonica and G. brasiliensis. J. Hymenoptera Res. 78, 1–17. https://doi.org/10.3897/jhr.78.55026 (2020).Article 

Google Scholar 
Chabert, S., Allemand, R., Poyet, M., Eslin, P. & Gibert, P. Ability of European parasitoids (Hymenoptera) to control a new invasive Asiatic pest, Drosophila suzukii. Biol. Control 63, 40–47. https://doi.org/10.1016/j.biocontrol.2012.05.005 (2012).Article 

Google Scholar 
Gonzalez-Cabrera, J., Moreno-Carrillo, G., Sanchez-Gonzalez, J. A., Mendoza-Ceballos, M. Y. & Arredondo-Bernal, H. C. Single and Combined Release of Trichopria drosophilae (Hymenoptera: Diapriidae) to Control Drosophila suzukii (Diptera: Drosophilidae). Neotrop. Entomol. 48, 949–956. https://doi.org/10.1007/s13744-019-00707-3 (2019).CAS 
Article 
PubMed 

Google Scholar 
Rossi Stacconi, M. V., Grassi, A., Ioriatti, C. & Anfora, G. Augmentative releases of Trichopria drosophilae for the suppression of early season Drosophila suzukii populations. Biocontrol 64, 9–19. https://doi.org/10.1007/s10526-018-09914-0 (2018).CAS 
Article 

Google Scholar 
Poyet, M. et al. The wide potential trophic niche of the asiatic fruit fly Drosophila suzukii: The key of its invasion success in temperate Europe?. PLoS ONE 10, e0142785. https://doi.org/10.1371/journal.pone.0142785 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Mazzetto, F. et al. Drosophila parasitoids in northern Italy and their potential to attack the exotic pest Drosophila suzukii. J. Pest Sci. 89, 837–850. https://doi.org/10.1007/s10340-016-0746-7 (2016).Article 

Google Scholar 
Wang, X. G., Kacar, G., Biondi, A. & Daane, K. M. Foraging efficiency and outcomes of interactions of two pupal parasitoids attacking the invasive spotted wing drosophila. Biol. Control 96, 64–71. https://doi.org/10.1016/j.biocontrol.2016.02.004 (2016).Article 

Google Scholar 
Rossi Stacconi, M. V. et al. Host location and dispersal ability of the cosmopolitan parasitoid Trichopria drosophilae released to control the invasive spotted wing Drosophila. Biol. Control 117, 188–196. https://doi.org/10.1016/j.biocontrol.2017.11.013 (2018).Article 

Google Scholar 
Esteban-Santiago, J. M., Rodríguez-Leyva, E., Lomeli-Flores, J. R. & González-Cabrera, J. Demographic parameters of Trichopria drosophilae in three host species. Entomol. Exp. Appl. 169, 330–337. https://doi.org/10.1111/eea.13026 (2021).CAS 
Article 

Google Scholar 
Häussling, B. J. M., Lienenluke, J. & Stokl, J. The preference of Trichopria drosophilae for pupae of Drosophila suzukii is independent of host size. Sci. Rep. 11, 995. https://doi.org/10.1038/s41598-020-80355-5 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Wang, X. G., Kacar, G., Biondi, A. & Daane, K. M. Life-history and host preference of Trichopria drosophilae, a pupal parasitoid of spotted wing drosophila. Biocontrol 61, 387–397. https://doi.org/10.1007/s10526-016-9720-9 (2016).CAS 
Article 

Google Scholar 
Woltz, J. M. & Lee, J. C. Pupation behavior and larval and pupal biocontrol of Drosophila suzukii in the field. Biol. Control 110, 62–69. https://doi.org/10.1016/j.biocontrol.2017.04.007 (2017).Article 

Google Scholar 
Ballman, E. S., Collins, J. A. & Drummond, F. A. Pupation behavior and predation on Drosophila suzukii (Diptera: Drosophilidae) Pupae in Maine wild blueberry fields. J. Econ. Entomol. 110, 2308–2317. https://doi.org/10.1093/jee/tox233 (2017).Article 
PubMed 

Google Scholar 
Guillén, L., Aluja, M. N., Equihua, M. & Sivinski, J. Performance of two fruit fly (Diptera: Tephritidae) pupal parasitoids (Coptera haywardi [Hymenoptera: Diapriidae] and Pachycrepoideus vindemiae [Hymenoptera: Pteromalidae]) under different environmental soil conditions. Biol. Control 23, 219–227. https://doi.org/10.1006/bcon.2001.1011 (2002).Article 

Google Scholar 
Yi, C. et al. Life history and host preference of Trichopria drosophilae from Southern China, one of the effective pupal parasitoids on the Drosophila species. Insects 11, 103. https://doi.org/10.3390/insects11020103 (2020).Article 
PubMed Central 

Google Scholar 
BoychevaWoltering, S., Romeis, J. & Collatz, J. Influence of the rearing host on biological parameters of Trichopria drosophilae, a potential biological control agent of Drosophila suzukii. Insects. https://doi.org/10.3390/insects10060183 (2019).Article 

Google Scholar 
Otto, M. & Mackauer, M. The developmental strategy of an idiobiont ectoparasitoid, Dendrocerus carpenteri: Influence of variations in host quality on offspring growth and fitness. Oecologia 117, 353–364. https://doi.org/10.1007/s004420050668 (1998).ADS 
Article 
PubMed 

Google Scholar 
Bates, D., Machler, M., Bolker, B. M. & Walker, S. C. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48. https://doi.org/10.18637/jss.v067.i01 (2015).Article 

Google Scholar 
R: A Language and Environment for Statistical Computing (Vienna, Austria, 2008).Johnson, S. N. & Gregory, P. J. Chemically-mediated host-plant location and selection by root-feeding insects. Physiol. Entomol. 31, 1–13. https://doi.org/10.1111/j.1365-3032.2005.00487.x (2006).CAS 
Article 

Google Scholar 
Bezerra Da Silva, C. S., Park, K. R., Blood, R. A. & Walton, V. M. Intraspecific competition affects the pupation behavior of spotted-wing drosophila (Drosophila suzukii). Sci. Rep. 9, 7775. https://doi.org/10.1038/s41598-019-44248-6 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Renkema, J. M. & Devkota, S. Pupation depth of spotted wing drosophila (Drosophila suzukii) and effects of field sanitation in Florida strawberries. Viii Int. Strawberry Symp. 1156, 849–855. https://doi.org/10.17660/ActaHortic.2017.1156.125 (2017).Article 

Google Scholar 
Tsitsipis, J. A. & Papanicolaou, E. P. Pupation depth in artificially reared olive fruits-flies Dacus-oleae (Diptera, Tephritidae), as affected by several physical characteristics of the substrates. Annales De Zoologie Ecologie Animale 11, 31–40 (1979).
Google Scholar 
Dimou, I., Koutsikopoulos, C., Economopoulos, A. P. & Lykakis, J. Depth of pupation of the wild olive fruit fly, Bactrocera (Dacus) oleae (Gmel.) (Dipt., Tephritidae), as affected by soil abiotic factors. J. Appl. Entomol. 127, 12–17. https://doi.org/10.1046/j.1439-0418.2003.00686.x (2003).Article 

Google Scholar 
de Belle, J. S., Hilliker, A. J. & Sokolowski, M. B. Genetic localization of foraging (for): A major gene for larval behavior in Drosophila melanogaster. Genetics 123, 157–163. https://doi.org/10.1093/genetics/123.1.157 (1989).Article 
PubMed 

Google Scholar 
Sokolowski, M. B. et al. Ecological genetics and behaviour of Drosophila melanogaster larvae in nature. Anim. Behav. 34, 403–408. https://doi.org/10.1016/S0003-3472(86)80109-9 (1986).Article 

Google Scholar 
Rodriguez, L., Sokolowski, M. B. & Shore, J. S. Habitat selection by Drosophila melanogaster larvae. J. Evol. Biol. 5, 61–70. https://doi.org/10.1046/j.1420-9101.1992.5010061.x (1992).Article 

Google Scholar 
McIntosh, H., Atucha, A., Townsend, P. A., Hills, W. B. & Guédot, C. Plastic mulches reduce adult and larval populations of Drosophila suzukii in fall-bearing raspberry. J. Pest. Sci. 95, 525–536. https://doi.org/10.1007/s10340-021-01456-2 (2021).Article 

Google Scholar 
Ballman, E. & Drummond, F. Larval movement of spotted wing drosophila, Drosophila suzukii (Matsumura) (Diptera: Drosophilidae). J. Kansas Entomol. Soc. 92, 412–421. https://doi.org/10.2317/0022-8567-92.1.412 (2019).Article 

Google Scholar  More

In this study, the five main characteristics of green spaces that were measured were area, perimeter, perimeter-area ratio, leaf area index, and canopy density. The structure of parameter between them is shown in Table 3.Table 3 Parameter structure of the cooling and humidification effect based on the spatial characteristics of green spaces.Full size tableCorrelation between various spatial characteristics and cooling and humidifying intensity in green spacesSmall-size green spacesFigures 4 and 6 shows the results of linear regressions between spatial characteristics and the cooling effect in small-size green spaces. There were relatively weak correlations between area, perimeter, perimeter-area ratio, leaf area index and cooling intensity, and a strong correlation between canopy density and cooling intensity. Small-size green space has the weakest positive correlation between perimeter-area ratio and cooling intensity (R2 = 0.11), and its canopy density and cooling intensity have the strongest positive correlation (R2 = 0.64). Meanwhile, small-size green space has weakest negative correlation between perimeter and humidifying intensity (R2 = 0.17), and its leaf area index and humidifying intensity have significant positive correlation (R2 = 0.42). Figures 4a and 5a show that for every 1 ha increase in area of small-size green spaces, the cooling intensity increased by 1.026 °C, and the humidifying intensity decreased by 1.56%. Figures 4b and 5b show that for every 100 m increase in perimeter, the cooling intensity decreases by 1.06 °C, and the humidifying intensity decreased by 1.19%. Figures 4c and 5c show that for every 0.01 increase in the perimeter-area ratio, the cooling intensity increases by 1.12 °C, and the humidifying intensity increased by 1.46%. Figures 4d and 5d show that for every 0.1 increase in the leaf area index, the cooling intensity increases by 1.11 °C, and the humidifying intensity increased by 1.12%. Figures 4e and 5e show that each 0.01 increase in the canopy density, the cooling intensity increases by 1.60 °C, and each 0.1 increase in canopy density, the humidifying intensity increased by 1.15% (Fig. 6).
Figure 4Linear regressions between spatial characteristics and cooling intensity of small-size green spaces.Full size imageFigure 5Linear regressions of spatial characteristics and humidifying intensity of small-size green spaces.Full size imageFigure 6The correlation between the spatial characteristics of small-size green spaces and the intensity of cooling and humidifying (GA means green area; GP means green perimeter; GPAR means green perimeter-area ratio; LAI means leaf area index; CD means canopy density).Full size imageMedium-size green spacesFigures 7 and 9 shows the linear regressions between spatial characteristics and cooling intensity in medium-size green spaces. There was an extremely significant positive correlation between area and cooling intensity, an insignificant positive correlation between the leaf area index and cooling intensity, and a relatively weak negative correlation between the other three characteristics and cooling intensity. Medium-size green space has the weakest negative correlation between canopy density and cooling intensity (R2 = 0.12), and its green area and cooling intensity have the strongest positive correlation (R2 = 0.83). Meanwhile, medium-size green space has weakest negative correlation between perimeter-area ratio and humidifying intensity (R2 = 0.41), and its area and humidifying intensity have most significant positive correlation (R2 = 0.81). Figures 7a and 8a show that for every 1 ha increase in area of medium-size green spaces, the cooling intensity increased by 1.19 °C, and the humidifying intensity increased by 1.24%. Figures 7b and 8b show that for every 100 m increase in perimeter, the cooling intensity decreases by 1.02 °C, and the humidifying intensity increased by 1.17%. Figures 7c and 8c show that for every 0.01 increase in the perimeter-area ratio, the cooling intensity decreases by 1.29 °C, and the humidifying intensity decreased by 2.40%. Figures 7d and 8d show that for every 0.1 increase in the leaf area index, the cooling intensity increases by 1.37 °C, and the humidifying intensity decreased by 1.92%. Figures 7e and 8e show that each 0.01 increase in the canopy density, increases the cooling intensity decreases by 1.23 °C, and the humidifying intensity decreased by 6.48% (Fig. 9).Figure 7Linear regressions between spatial characteristics and cooling intensity of medium-size green spaces.Full size imageFigure 8Linear regressions of spatial characteristics and humidifying intensity of medium-size green spaces.Full size imageFigure 9The correlation between the spatial characteristics of medium-size green spaces and the intensity of cooling and humidifying (GA means green area; GP means green perimeter; GPAR means green perimeter-area ratio; LAI means leaf area index; CD means canopy density).Full size imageLarge-size green spacesFigures 10 and 12 shows the linear regressions between spatial characteristics and cooling intensity in large-size green spaces. There was an insignificant correlation between area and cooling intensity, a weak correlation between canopy density and cooling intensity, and a significant correlation between perimeter, perimeter-area ratio and the leaf area index and cooling intensity. Medium-size green space has the weakest negative correlation between green area and cooling intensity (R2 = 0.35), and its leaf area index and cooling intensity have the strongest positive correlation (R2 = 0.92). Meanwhile, medium-size green space has weakest negative correlation between perimeter-area ratio and humidifying intensity (R2 = 0.11), and its leaf area index and humidifying intensity have most significant positive correlation (R2 = 0.39). Figures 10a and 11a show that for every 1 ha increase in area of large-size green spaces, the cooling intensity decreased by 1.02 °C, and the humidifying intensity decreased by 1.22%. Figures 10b and 11b show that for every 100 m increase in perimeter, the cooling intensity decreases by 1.05 °C, and the humidifying intensity decreased by 1.34%. Figures 10c and 11c show that for every 0.005 increase in the perimeter-area ratio, the cooling intensity decreases by 1.43 °C, and each 0.01 increase in perimeter-area ratio, the humidifying intensity decreased by 1.27%. Figures 10d and 11d show that for every 0.1 increase in the leaf area index, the cooling intensity increases by 2.41 °C, and the humidifying intensity increased by 1.37%. Figures 10e and 11e show that each 0.1 increase in the canopy density, the cooling intensity increased by 3.69 °C, and the humidifying intensity decreased by 2.84% (Fig. 12).Figure 10Linear regressions of spatial characteristics and cooling intensity of large-size green spaces.Full size imageFigure 11Linear regressions of spatial characteristics and humidifying intensity of large-size green spaces.Full size imageFigure 12The correlation between the spatial characteristics of large-size green spaces and the intensity of cooling and humidifying (GA means green area; GP means green perimeter; GPAR means green perimeter-area ratio; LAI means leaf area index; CD means canopy density).Full size imageQuantitative analysis of the microclimatic effects of different types of green spacesQuantitative analysis of the effects of different types of green space on cooling intensityFigure 13 shows the linear regressions between the different types of green spaces and cooling intensity. There were negative correlations between green spaces a short, medium, and long distance from a water body and cooling intensity in small-size green spaces, medium-size green spaces and large-size green spaces. The negative correlation between the distance to a water body and cooling intensity in medium-size green spaces was most significant (R2 = 0.985). The greater the distance to a water body, the lower the cooling intensity. For medium-size green spaces, for every 1/4 increase in the distance ratio, the cooling intensity decreased by 0.81 °C. For small-size green spaces, for every 1/4 increase in the distance ratio, the cooling intensity decreased by 1.04 °C. For large-size green spaces, for every 1/4 increase in the distance ratio, the cooling intensity decreased by 1.36 °C. For small-, medium-, and large-size green spaces, there was a positive correlation between canopy density and cooling intensity. There was a most significant positive correlation between canopy density and cooling intensity in large-size green spaces (R2 = 0.941). The greater the canopy density, the greater the cooling intensity. For large green spaces, for every 0.5 increase in canopy density, the cooling intensity increased by 0.16 °C. For small-size green spaces, for every 0.5 increase in canopy density, the cooling effect increased by 0.15 °C. For medium-size green spaces, for every 0.5 increase in canopy density, the cooling intensity increased by 0.16 °C.Figure 13Linear regressions between the distance from different types of green spaces to water areas, canopy density and cooling intensity.Full size imageQuantitative analysis of the effects of different types of green space on humidifying intensityFigure 14 shows the linear regression between the distance of a green space from a water body, canopy density and humidifying intensity. There was a negative correlation between the distance to a water body and humidifying intensity in small, medium, and large green spaces. The negative correlation between the distance to a water body and humidifying intensity in small green spaces was most significant (R2 = 0.996). The longer the distance, the lower the humidifying intensity. For small green spaces, for every 1/4 in-crease in the distance ratio, the humidifying intensity decreased by 4.23%. For medium-size green spaces, for every 1/4 increase in the distance ratio, the humidifying intensity decreased by 3.02%. For large-size green spaces, for every 1/4 increase in the distance ratio, the humidifying intensity de-creased by 6.14%. For small, medium, and large green spaces, there was a positive correlation between canopy density and humidifying intensity. The positive correlation between canopy density and humidifying intensity in medium-size green spaces was extremely significant (R2 = 0.925). The greater the canopy density, the greater the humidifying intensity. For medium-size green spaces, for every 0.5 increase in canopy density, the humidifying intensity increased by 3.29%. For small-size green spaces, for every 0.5 increase in canopy density, the humidifying intensity increased by 3.17%. For large-size green spaces, for every 0.5 increase in canopy density, the humidifying intensity increased by 4.06% (Fig. 15).
Figure 14Linear regressions between the distance from different types of green space to water area, canopy density and humidifying intensity.Full size imageFigure 15Correlation of different green space types with water distance, canopy density and cooling and humidifying intensity.Full size imageEffect of shape and area of water bodies on microclimatic effects based on numerical simulationBanded waterWe constructed a numerical simulation model to explore the effects of a simulated increase in water body area on cooling and humidification. Figure 16 shows the simulated distribution characteristics of temperature and relative humidity after a 5% and 10% increase in water area at 14:00 when temperatures were high. The results suggest that between 7:00 and 10:00, with a 5% and 10% increase in water area, the air temperature was basically the same and the cooling effect was insignificant. However, between 12:00 and 19:00 and particularly in the hours between 13:00 and 16:00 when temperatures were highest, a 5% increase in water area produced a significant cooling effect, with a daily average value of 0.05 °C and a maximum value of 0.09 °C. A 10% increase in water area produced an extremely significant cooling effect, with a daily average value of 0.07 °C and a maximum value of 0.14 °C. From 11:00 to 19:00, a 5% increase in water area produced a significant humidifying effect, with a daily average value of 0.08% and a maximum value of 0.17%. A 10% increase produced an extremely significant humidifying effect, with a daily average value of 0.13% and a maximum value of 0.26% (See supplementary file).Figure 16Distribution characteristics of cooling and humidifying effects of simulated increase of banded water area at 14:00. (a) original cooling effect of banded water in the sample area; (b) cooling effect of 5% increase in water area; (c) cooling effect of 10% increase in water area; (d) original humidifying effect of banded water in the sample area; (e) humidifying effect of 5% increase in water area; (f) humidifying effect of 10% increase of water area.Full size imageMassive waterFigure 17 shows the simulated distribution characteristics of the cooling and humidifying effects after a 5% and 10% increase in the water area at 14:00 when temperatures were high. Between 8:00 and 19:00, a 5% and 10% increase in water area produced a significant cooling effect. At 19:00, the numerical simulation result was abnormal when the water area increased by 5% and 10%; at 13:00, the numerical simulation result was also ab-normal when the water area increased by 10%. After excluding the abnormal simulated data, a 5% increase in water area produced a cooling effect, with a daily average value of 0.06 °C and a maximum value of 0.10 °C. A 10% increase in water area produced an extremely significant cooling effect, with a daily average value of 0.10 °C and a maximum value of 0.18 °C. Between 11:00 and 19:00, a 5% increase in water area produced a significant humidifying effect, with a daily average value of 0.05% and a maximum value of 0.13%. A 10% increase in water area produced an extremely significant humidifying effect, with a daily average value of 0.13% and a maximum value of 0.27% (See supplementary file).Figure 17Distribution characteristics of cooling and humidifying effects of simulated increase of massive water area at 14:00. (a) original cooling effect of massive water in the sample area; (b) cooling effect of 5% increase in water area; (c) cooling effect of 10% increase in water area; (d) original humidifying effect of massive water in the sample area; (e) humidifying effect of 5% increase in water area; (f) humidifying effect of 10% increase of water area.Full size imageAnnular waterFigure 18 shows the simulated distribution characteristics of the cooling and humidifying effects after a 5% and 10% increase in the area of the annular water body at 14:00 when temperatures were high. Between 7:00 and 19:00, a 5% and 10% increase in water area produced a significant cooling effect. Between 11:00 and 16:00 when temperatures were high, a 5% increase in water area produced a cooling effect, with a daily average value of 0.06 °C and a maximum value of 0.14 °C°C and a 10% increase in water area produced an extremely significant cooling effect, with a daily average value of 0.13 °C and a maximum value of 0.28 °C. Between 7:00 and 19:00, a 5% and 10% increase in water area produced significant humidifying effects. Between 11:00 and 16:00 when temperatures were high, a 5% increase in water area produced an extremely significant humidifying effect, with a daily average value of 0.17% and a maximum value of 0.39% and a 10% increase in water area produced an extremely significant humidifying effect with a daily average value of 0.38% and a maximum value of 0.81% (See supplementary file).Figure 18Distribution characteristics of cooling and humidifying effects of simulated increase of annular water area at 14:00. (a) original cooling effect of annular water in the sample area; (b) cooling effect of 5% increase in water area; (c) cooling effect of 10% increase in water area; (d) original humidifying effect of annular water in the sample area; (e) humidifying effect of 5% increase in water area; (f) humidifying effect of 10% increase of water area.Full size image More