Ecology

Aissi A, Beghami Y, Heuertz M (2019) Le chêne faginé (Quercus faginea, Fagaceae) en Algérie: potentiel germinatif et variabilité morphologique des glands et des semis. Plant Ecol Evol 152:437–449Article 

Google Scholar 
Aissi A, Beghami Y, Lepais O, Véla E (2021) Morphological and taxonomic analysis of Quercus faginea (Fagaceae) complex in Algeria. Botany 99:99–113Article 

Google Scholar 
Alberto F, Niort J, Derory J, Lepais O, Vitalis R, Galop D et al. (2010) Population differentiation of sessile oak at the altitudinal front of migration in the French Pyrenees. Mol Ecol 19:2626–2639CAS 
PubMed 
Article 

Google Scholar 
Allendorf FW, Luikart G (2007) Conservation and the genetics of populations. Blackwell PubDe Barba M, Miquel C, Lobréaux S, Quenette PY, Swenson JE, Taberlet P (2016) High-throughput microsatellite genotyping in ecology: Improved accuracy, efficiency, standardization and success with low-quantity and degraded DNA. Mol Ecol Resour 17:492–507PubMed 
Article 
CAS 

Google Scholar 
Barthe S, Gugerli F, Barkley NA, Maggia L, Cardi C, Scotti I (2012) Always look on both sides: phylogenetic information conveyed by simple sequence repeat allele sequences. PLoS One 7:e40699CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Beaumont MA, Zhang W, Balding DJ (2002) Approximate Bayesian computation in population genetics. Genetics 162:2025–35PubMed 
PubMed Central 
Article 

Google Scholar 
Belkhir K, Borsa P, Chikhi L, Raufaste N, Bonhomme F (2004) GENETIX 4.05, logiciel sous Windows TM pour la génétique des populations. Laboratoire Génome, Populations, Interactions, CNRS UMR 5171, Université de Montpellier II, Montpellier (France)
Google Scholar 
Bradbury IR, Wringe BF, Watson B, Paterson I, Horne J, Beiko R et al. (2018) Genotyping-by-sequencing of genome-wide microsatellite loci reveals fine-scale harvest composition in a coastal Atlantic salmon fishery. Evol Appl 11:918–930CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Buschbom J, Yanbaev Y, Degen B (2011) Efficient long-distance gene flow into an isolated relict oak stand. J Hered 102:464–472PubMed 
Article 

Google Scholar 
Chapuis M, Raynal L, Plantamp C, Meynard CN, Blondin L, Marin J et al. (2020) A young age of subspecific divergence in the desert locust inferred by ABC Random Forest. Mol Ecol 29:4542–4558PubMed 
Article 

Google Scholar 
Cornuet J-M, Ravigné V, Estoup A (2010) Inference on population history and model checking using DNA sequence and microsatellite data with the software DIYABC (v1.0). BMC Bioinforma 11:401Article 
CAS 

Google Scholar 
Crow JF, Aoki K (1984) Group selection for a polygenic behavioral trait: Estimating the degree of population subdivision. Proc Natl Acad Sci USA 81:6073–6077CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Curto M, Winter S, Seiter A, Schmid L, Scheicher K, Barthel LMF et al. (2019) Application of a SSR-GBS marker system on investigation of European Hedgehog species and their hybrid zone dynamics. Ecol Evol 9:2814–2832PubMed 
PubMed Central 
Article 

Google Scholar 
Darby BJ, Erickson SF, Hervey SD, Ellis-Felege SN (2016) Digital fragment analysis of short tandem repeats by high-throughput amplicon sequencing. Ecol Evol 6:4502–4512PubMed 
PubMed Central 
Article 

Google Scholar 
Dickey AM, Hall PM, Shatters RG, Mckenzie CL (2013) Evolution and homoplasy at the Bem6 microsatellite locus in three sweetpotato whitefly (Bemisia tabaci) cryptic species. BMC Res Notes 6:249PubMed 
PubMed Central 
Article 

Google Scholar 
Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR (2013) NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol Ecol Resour 14:209–214PubMed 
Article 

Google Scholar 
Durand J, Bodenes C, Chancerel E, Frigerio JM, Vendramin G, Sebastiani F et al. (2010) A fast and cost-effective approach to develop and map EST-SSR markers: oak as a case study. BMC Genomics 11:570PubMed 
PubMed Central 
Article 

Google Scholar 
Earl DA, vonHoldt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour 4:359–361Article 

Google Scholar 
Estoup A, Jarne P, Cornuet J-M (2002) Homoplasy and mutation model at microsatellite loci and their consequences for population genetics analysis. Mol Ecol 11:1591–1604CAS 
PubMed 
Article 

Google Scholar 
Estoup A, Raynal L, Verdu P, Marin J-M (2018) Model choice using Approximate Bayesian Computation and Random Forests: analyses based on model grouping to make inferences about the genetic history of Pygmy human populations. J la Société Fr Stat 159:167–190
Google Scholar 
Excoffier L, Dupanloup I, Huerta-Sánchez E, Sousa VC, Foll M (2013) Robust demographic inference from genomic and SNP data. PLoS Genet 9:e1003905PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Excoffier L, Lischer HEL (2010) Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour 10:564–567PubMed 
Article 

Google Scholar 
Falush D, Stephens M, Pritchard JK (2003) Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164:1567–87CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Feliner GN (2014) Patterns and processes in plant phylogeography in the Mediterranean Basin. A review. Perspect Plant Ecol, Evol Syst 16:265–278Article 

Google Scholar 
Flagel L, Brandvain Y, Schrider DR (2019) The unreasonable effectiveness of convolutional neural networks in population genetic inference. Mol Biol Evol 36:220–238CAS 
PubMed 
Article 

Google Scholar 
Foll M, Gaggiotti O (2008) A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective. Genetics 180:977–993PubMed 
PubMed Central 
Article 

Google Scholar 
Gaggiotti OE, Chao A, Peres-Neto P, Chiu CH, Edwards C, Fortin MJ et al. (2018) Diversity from genes to ecosystems: A unifying framework to study variation across biological metrics and scales. Evol Appl 11:1176–1193PubMed 
PubMed Central 
Article 

Google Scholar 
García Murillo P., Harvey-Brown Y (2017) Quercus canariensis. In: The IUCN Red List of Threatened Species,, p e.T78809256A80570536GBIF Secratariat (2021a) Quercus faginea Lam. In: GBIF Backbone Taxonomy. Checklist dataset https://doi.org/10.15468/39omeiaccessed via GBIF.org on 2022-05-10GBIF Secratariat (2021b). Quercus canariensis Willd. In: GBIF Backbone Taxonomy. Checklist dataset https://doi.org/10.15468/39omeiaccessed via GBIF.org on 2022-05-10Gómez A, Lunt DH (2007) Refugia within refugia: Patterns of phylogeographic concordance in the Iberian peninsula. In: Weiss S, Ferrand N (eds) Phylogeography of Southern European Refugia. Springer Netherlands, Dordrecht, p 155–188Chapter 

Google Scholar 
Gorener V, Harvey-Brown Y, Barstow M (2017) Quercus canariensis. IUCN red List Threat species e.T7880925Goudet J (1995) FSTAT (Version 1.2): A computer program to calculate F-statistics. J Hered 86:485–486Article 

Google Scholar 
Hampe A, Petit RJ (2005) Conserving biodiversity under climate change: the rear edge matters. Ecol Lett 8:461–7PubMed 
Article 

Google Scholar 
Hardy OJ, Charbonnel N, Fréville H, Heuertz M (2003) Microsatellite allele sizes: a simple test to assess their significance on genetic differentiation. Genetics 163:1467–82CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hardy OJ, Vekemans X (2002) SPAGeDi: a versatile computer program to analyse spatial genetic structure at the individual or population levels. Mol Ecol Notes 2:618–620Article 
CAS 

Google Scholar 
Harvey-Brown Y, García Murillo PG, Buira A (2017) Quercus faginea. IUCN Red List Threat Species: e.T78916251A80570540.Henriques R, von der Heyden S, Matthee CA (2016) When homoplasy mimics hybridization: a case study of Cape hakes (Merluccius capensis and M. paradoxus). PeerJ 4:e1827PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Hewitt GM (1999) Post-glacial re-colonization of European biota. Biol J Linn Soc 68:87–112Article 

Google Scholar 
Hey J, Won YJ, Sivasundar A, Nielsen R, Markert JA (2004) Using nuclear haplotypes with microsatellites to study gene flow between recently separated Cichlid species. Mol Ecol 13:909–919CAS 
PubMed 
Article 

Google Scholar 
Hoban S, Bruford M, D’Urban Jackson J, Lopes-Fernandes M, Heuertz M, Hohenlohe PA et al. (2020) Genetic diversity targets and indicators in the CBD post-2020 Global Biodiversity Framework must be improved. Biol Conserv 248:108654Article 

Google Scholar 
Hoogenboom J, de Knijff P, Laros JFJ, de Leeuw RH, van der Gaag KJ, Sijen T (2016) FDSTools: A software package for analysis of massively parallel sequencing data with the ability to recognise and correct STR stutter and other PCR or sequencing noise. Forensic Sci Int Genet 27:27–40PubMed 
Article 
CAS 

Google Scholar 
Hothorn T, Hornik K, van de Wiel MA, Zeileis A (2008) Implementing a class of permutation tests: the coin package. J Stat Softw 28:1–23Article 

Google Scholar 
Jamieson IG, Allendorf FW (2012) How does the 50/500 rule apply to MVPs? Trends Ecol Evol 27:578–584PubMed 
Article 

Google Scholar 
Jerome D, Vasquez F (2018) Quercus faginea. IUCN Red List Threat Species e.T7891625Kalinowski ST (2005) HP-RARE 1.0: A computer program for performing rarefaction on measures of allelic richness. Mol Ecol Notes 5:187–189CAS 
Article 

Google Scholar 
Kampfer S, Lexer C, Glössl J, Steinkellner H (1998) Characterization of (GA)n microsatellite loci from Quercus robur. Hereditas 129:183–186CAS 
Article 

Google Scholar 
Kivelä M, Arnaud-Haond S, Saramäki J (2015) EDENetworks: A user-friendly software to build and analyse networks in biogeography, ecology and population genetics. Mol Ecol Resour 15:117–122PubMed 
Article 

Google Scholar 
Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I (2015) Clumpak: A program for identifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour 15:1179–1191CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Layton KKS, Dempson B, Snelgrove PVR, Duffy SJ, Messmer AM, Paterson IG et al. (2020) Resolving fine‐scale population structure and fishery exploitation using sequenced microsatellites in a northern fish. Evol Appl: eva.12922.Lepais O, Chancerel E, Boury C, Salin F, Manicki A, Taillebois L et al. (2020) Fast sequence-based microsatellite genotyping development workflow. PeerJ 8:e9085PubMed 
PubMed Central 
Article 

Google Scholar 
Lepais O, Leger V, Gerber S (2006) Short note: high throughput microsatellite genotyping in oak species. Silvae Genet 55:238Article 

Google Scholar 
Lepais O, Muller SD, Ben Saad-Limam S, Benslama M, Rhazi L, Belouahem-Abed D et al. (2013) High genetic diversity and distinctiveness of rear-edge climate relicts maintained by ancient tetraploidisation for Alnus glutinosa. PLoS One 8:e75029CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Leroy T, Roux C, Villate L, Bodénès C, Romiguier J, Paiva JAP et al. (2017) Extensive recent secondary contacts between four European white oak species. N. Phytol 214:865–878CAS 
Article 

Google Scholar 
Lye GC, Lepais O, Goulson D (2011) Reconstructing demographic events from population genetic data: the introduction of bumblebees to New Zealand. Mol Ecol 20:2888–900CAS 
PubMed 
Article 

Google Scholar 
Magri D, Fineschi S, Bellarosa R, Buonamici A, Sebastiani F, Schirone B et al. (2007) The distribution of Quercus suber chloroplast haplotypes matches the palaeogeographical history of the western Mediterranean. Mol Ecol 16:5259–66CAS 
PubMed 
Article 

Google Scholar 
Marin J, Pudlo P, Estoup A, Robert C (2018) Likelihood-free model choice. In: Sisson S A, Fan Y, Beaumont M (eds) Handbook of Approximate Bayesian Computation, CRC Press, pp. 153.Médail F, Diadema K (2009) Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J Biogeogr 36:1333–1345Article 

Google Scholar 
Moracho E, Moreno G, Jordano P, Hampe A (2016) Unusually limited pollen dispersal and connectivity of Pedunculate oak (Quercus robur) refugial populations at the species’ southern range margin. Mol Ecol 25:3319–3331CAS 
PubMed 
Article 

Google Scholar 
Mountain JL, Knight A, Jobin M, Gignoux C, Miller A, Lin AA et al. (2002) SNPSTRs: Empirically derived, rapidly typed, autosomal haplotypes for inference of population history and mutational processes. Genome Res 12:1766–1772CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Muir G, Lowe AJ, Fleming CC, Vogl C (2004) High nuclear genetic diversity, high levels of outcrossing and low differentiation among remnant populations of Quercus petraea at the margin of its range in Ireland. Ann Bot 93:691–697CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Neophytou C, Gärtner SM, Vargas-Gaete R, Michiels H-G (2015) Genetic variation of Central European oaks: shaped by evolutionary factors and human intervention? Tree Genet Genomes 11:1–15Article 

Google Scholar 
Payseur BA, Cutter AD (2006) Integrating patterns of polymorphism at SNPs and STRs. Trends Genet 22:424–429CAS 
PubMed 
Article 

Google Scholar 
Petit RJ, Brewer S, Bordács S, Burg K, Cheddadi R, Coart E et al. (2002) Identification of refugia and post-glacial colonisation routes of European white oaks based on chloroplast DNA and fossil pollen evidence. Ecol Manag 156:49–74Article 

Google Scholar 
Petit RJ, Hampe A, Cheddadi R (2005) Climate changes and tree phylogeography in the Mediterranean. Taxon 54:877–885Article 

Google Scholar 
Press MO, Hall AN, Morton EA, Queitsch C (2019) Substitutions are boring: Some arguments about parallel mutations and high mutation rates. Trends Genet 35:253–264CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Press MO, Mccoy RC, Hall AN, Akey JM, Queitsch C (2018) Massive variation of short tandem repeats with functional consequences across strains of Arabidopsis thaliana. Genome Res 28:1169–1178CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155:945–59CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Pudlo P, Marin J-M, Estoup A, Cornuet J-M, Gautier M, Robert CP (2016) Reliable ABC model choice via random forests. Bioinformatics 32:859–866CAS 
PubMed 
Article 

Google Scholar 
Ramakrishnan U, Mountain JL (2004) Precision and accuracy of divergence time estimates from STR and SNPSTR variation. Mol Biol Evol 21:1960–1971CAS 
PubMed 
Article 

Google Scholar 
Raynal L, Marin J-M, Pudlo P, Ribatet M, Robert CP, Estoup A (2019) ABC random forests for Bayesian parameter inference. Bioinformatics 35:1720–1728CAS 
PubMed 
Article 

Google Scholar 
Rodríguez-Sánchez F, Hampe A, Jordano P, Arroyo J (2010) Past tree range dynamics in the Iberian Peninsula inferred through phylogeography and palaeodistribution modelling: A review. Rev Palaeobot Palynol 162:507–521Article 

Google Scholar 
Šarhanová P, Pfanzelt S, Brandt R, Himmelbach A, Blattner FR (2018) SSR-seq: Genotyping of microsatellites using next-generation sequencing reveals higher level of polymorphism as compared to traditional fragment size scoring. Ecol Evol 8:10817–10833PubMed 
PubMed Central 
Article 

Google Scholar 
Scotti-Saintagne C, Mariette S, Porth I, Goicoechea PG, Barreneche T, Bodénès C et al. (2004) Genome scanning for interspecific differentiation between two closely related oak species [Quercus robur L. and Q. petraea (Matt.) Liebl.]. Genetics 168:1615–26CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Vartia S, Villanueva-Cañas JL, Finarelli J, Farrell ED, Collins PC, Hughes GM et al. (2016) A novel method of microsatellite genotyping-by-sequencing using individual combinatorial barcoding. R Soc Open Sci 3:150565PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Viruel J, Haguenauer A, Juin M, Mirleau F, Bouteiller D, Boudagher-Kharrat M et al. (2018) Advances in genotyping microsatellite markers through sequencing and consequences of scoring methods for Ceratonia siliqua (Leguminosae). Appl Plant Sci 6:e01201PubMed 
PubMed Central 
Article 

Google Scholar 
Wang J (2016) Individual identification from genetic marker data: developments and accuracy comparisons of methods. Mol Ecol Resour 16:163–175CAS 
PubMed 
Article 

Google Scholar 
Waples RS, Do C (2010) Linkage disequilibrium estimates of contemporary Ne using highly variable genetic markers: A largely untapped resource for applied conservation and evolution. Evol Appl 3:244–262PubMed 
Article 

Google Scholar 
Xie KT, Wang G, Thompson AC, Wucherpfennig JI, Reimchen TE, MacColl ADC et al. (2019) DNA fragility in the parallel evolution of pelvic reduction in stickleback fish. Science (80-) 363:81–84CAS 
Article 

Google Scholar  More

In our simulations for the autotrophs, we varied two of the three trade-off properties (level of defence, plasticity costs and defence costs; see Fig. 1b) at a time and kept the third one constant. This results in three constellations reflecting three different trade-offs between these properties (Table 1):

parallel: trade-off between defence and plasticity costs;

crossing: trade-off between defence costs and plasticity costs;

angle: trade-off between defence and defence costs.

Table 1 Description of the three constellations parallel, crossing, and angle defining the position of the four phenotypes in the trait space of defence and growth rate.Full size tableIn all three constellations, the autotrophic species B spanned the entire defence range, i.e. it had a completely undefended phenotype Bu and a maximally defended phenotype Bd. A either had a more limited defence range (in constellations parallel and angle) or spanned the entire range as well (in constellation crossing), representing three distinct ways that the trade-off between defence, growth rate, and plasticity range may play out. For each constellation, we varied the maximum switching rate χmax over 5 orders of magnitude to investigate the effect of plasticity (Table 1, middle row). This parameter determines how rapidly a species can switch between phenotypes (see “Methods”, “Exchange rates”); higher values indicate faster adaptation. These results were also compared with a non-plastic baseline scenario where both phenotypes of each species are presented but χmax = 0 (Table 1, upper row), as well as a rigid scenario where the species have only a single phenotype (Table 1, bottom row). All parameters and their values can be found in Supplementary Table S1.In the following, we give a detailed description of the results for constellation parallel, where the autotroph species A and B have the same defence costs resulting in parallel trade-off lines between defence and growth rate, while varying the level of defence for A and varying the plasticity costs for B (Table 1, left column). We start with examining patterns for the phenotype biomasses, coexistence and community stability in the non-plastic baseline scenario “parallel 0”, and then compare the corresponding scenarios with a low exchange rate (“parallel 0.01”) and a high exchange rate (“parallel 1”). We next discuss the other two constellations (crossing and angle, Table 1) more briefly. Finally, we generalize across all scenarios and focus on the coexistence, the degree of maladaptive switching, and the consumer and total autotroph biomasses.Non-plastic baseline dynamics: scenario parallel 0
In this scenario, four single phenotypes unconnected by exchange compete with each other. Thus, species coexistence here depends entirely on phenotype coexistence: the trade-offs have to be such that for each species, at least one phenotype is a good enough competitor to survive. Which phenotypes survive depend on the two trade-off parameters, defence of the defended phenotype of species A (dAu) and plasticity costs for species B (pcB), which thus determine whether coexistence is possible.The defence costs were kept constant at an intermediate value of 0.3 for both species, resulting in parallel trade-off lines (Table 1, scenario “parallel 0”). The undefended phenotype of A, Au, is a growth-specialist with the highest growth rate of all phenotypes. The defended phenotype of the same species, Ad, has a defence between 0 and 0.9 and a relatively high growth rate, and can be viewed as a generalist. Species B has variable plasticity costs that lower the growth rate of both phenotypes. The defended phenotype of species B, Bd, has the lowest growth rate of all phenotypes but is very well-defended, and thus a defence-specialist. Its undefended phenotype, Bu, is as undefended as Au but has a lower growth rate; it is thus always an inferior competitor and inevitably goes extinct (Fig. 2c).Figure 2Biomasses, coexistence and trait space for scenario parallel 0. Biomasses of the four autotrophic phenotypes (a–d), their coexistence patterns (e), the consumer biomass (f) and the autotrophs’ trait values (g–j) (higher biomasses are shown by darker colours). Lines in (a–f) separate the regions I–III of different coexistence patterns. Note that in (a–f), the y-axis is reversed to show increasing fitness along all axes. An exemplary trait combination for every region is shown in (g–j); larger symbols indicate the surviving phenotypes. Shaded areas in (e) depict oscillating systems (quarter-lag predator–prey cycles in dense shading, antiphase cycles in loose shading).Full size imageAs Bu never survives, coexistence of the autotroph species requires the survival of defence-specialist Bd. Bd can only survive if Ad is not too defended, because Ad has a higher growth rate than Bd and will outcompete Bd in the “defended” niche otherwise (region Ib; Fig. 2d,h). A second criterion is that the plasticity costs for B must not be too high, because then the benefits of the defence of Bd no longer outweigh the costs, and it will go extinct even if there are no other highly defended phenotypes around (region Ia; Fig. 2d,g). In the regions where Bd goes extinct, species coexistence is not possible (Fig. 2e). The generalist Ad either survives by itself (region Ia in Fig. 2b,g) if its defence is low to intermediate, or together with the growth-specialist Au if its defence is high (region Ib in Fig. 2a,b,h). In the regions II and III where Bd survives, it never survives on its own, but always together with one of the phenotypes of A. It coexists with the growth-specialist Au if the plasticity costs are very low (region II in Fig. 2,i), and together with Ad if they are low to intermediate (region III in Fig. 2,j). These two regions do support species coexistence (Fig. 2e).In three of the four regions (Ib, II and III in Fig. 2f), consumer biomass is low, because the final community always contains a well-defended phenotype (Ad in region Ib, and Bd in regions II and III); the overall level of defence of the community is relatively high in these regions (Supplementary Figure S1). Conversely, consumer biomass is relatively high in region Ia, because the only surviving autotroph phenotype is relatively fast-growing and fairly undefended (Fig. 2f,g). The regions where a well-defended phenotype survives often show antiphase cycles (Ib, II and III in Fig. 2e). These cycles do not occur in the region where only Ad survives (Ia in Fig. 2e); but regular quarter-lag predator–prey cycles can be found here if Ad is almost entirely undefended.While the community defence (i.e. mean defence of the autotroph community) depends strongly on the coexisting phenotypes, the community growth rate is roughly constant because over the entire trait space, at least one phenotype with a high growth rate always survives (Supplementary Figure S1). The standing variance of the community defence was high when two phenotypes coexist as they occupy different niches along the defence axis (Fig. 2h–j). In contrast, the variance of the community growth rate was very low and almost constant across all regions.Effect of phenotypic plasticityEven a little bit of plasticity in the scenario parallel 0.01 (χmax = 0.01) can change the above patterns for coexistence, stability, and average consumer biomass (Fig. 3a–d). While the autotrophs are intuitively expected to benefit from being plastic, the effect of plasticity on consumer biomass always turned out to be positive (Fig. 3a). This may be explained by the fact that switching was always, on average, maladaptive (Fig. 3c,d), measured by the adaptation index Φ (see Eqs. (11–13) in “Methods”). This index combines information on the net “flow” of individuals due to switching (i.e. whether more undefended individuals switch to defended or vice versa) with the fitness difference between the two phenotypes, and thus measures whether overall, more individuals switch from a low-fitness to a high-fitness phenotype (adaptive) or the reverse (maladaptive). This index can approach zero, but is always negative at equilibrium (see Appendix B), indicating maladaptive switching.Figure 3Consumer biomass, autotroph coexistence and maladaptive switching for the scenarios parallel 0.01 (a–d) and parallel 1 (e–h). Consumer biomass (a,e), the autotroph coexistence patterns (b,f), and the autotrophs’ maladaptive switching Φ (c,d,g,h) (higher biomasses or more intensive maladaptive switching are shown by darker colours). Lines separate the regions I–III of different autotroph coexistence. The y-axis is reversed to follow the pattern of increasing fitness. Grey areas in (c,d,g,h) depict areas where the species was extinct. Shaded areas in b and f depict oscillating systems (quarter-lag predator–prey cycles in dense shading, antiphase cycles in loose shading).Full size imageThe most striking effect of plasticity was on coexistence, which was affected both positively and negatively by plasticity in different regions of the parameter space (Fig. 3b, Supplementary Figure S4a–d). A negative effect on coexistence is seen in region II, where the autotroph species previously coexisted (Fig. 2e), while with plasticity, B outcompeted A (Fig. 3b). Without plasticity, coexistence was possible in this region because Au and Bd survived; importantly, Au outcompeted Bu due to its higher growth rate, even though the difference between their growth rates is very small in this region (Fig. 2i). Plasticity reverses the competitive exclusion pattern between the two undefended phenotypes: Bu receives a constant flow of biomass from the well-defended Bd, which compensates for its slightly lower growth rate and allows it to outcompete Au. Thus, coexistence is reduced as a direct consequence of maladaptive switching.Plasticity can also promote coexistence, as the coexistence region now extends into former region Ib where the generalist Ad is highly defended (Fig. 2b, Supplementary Figure S1a). This is also an effect of maladaptive switching, though in this case the effect is indirect, mediated through the effect of plasticity on consumer biomass. Without plasticity, coexistence was impossible in region Ib because Bd was always outcompeted by Ad: even though the latter had a slightly lower level of defence, this was outweighed by its higher growth rate, making Ad the superior competitor over Bd. However, plasticity changes this because maladaptive switching increases the consumer biomass, which in turn alters the cost/ benefit balance of defence: Bd derives a stronger benefit from its high level of defence, which now outweighs the cost and allows it to survive. Coexistence through this mechanism is not possible when the plasticity costs for B are too high or when Ad is too well-defended, explaining the narrowing of the coexistence “tail” for high defence of Ad (Fig. 3b).While the patterns of coexistence changed when allowing for plasticity, the patterns in the trait values were nearly indistinguishable from the previous scenario (Supplementary Figure S1, S2). Finally, plasticity had a strong impact on the community dynamics, as most of the antiphase cycles were stabilized (Ib, II, III in Fig. 3b). Their area decreased sharply as these cycles were characterized by asynchronous dynamics between the two prey phenotypes, which were reduced by plasticity. In contrast, the area of the quarter-lag predator–prey cycles remained unaffected by plasticity.All the above patterns were found to a far stronger degree with a higher amount of plasticity (χmax = 1; Fig. 3e–h, Supplementary Figure S4e–h). Consumer biomass increased strongly everywhere (cf. Fig. 3a,e), reflecting the strong increase in the degree of maladaptive switching (cf. Fig. 3c,d,g,h). The higher exchange rates led to more synchronization between the phenotypes, extinguishing the antiphase cycles completely (Fig. 3f). It also decreased the biomass of both defended phenotypes (cf. Supplementary Figure S4b,d,f,h). This in turn led to a lower community defence and a higher community growth rate (Supplementary Figure S3) both contributing to a higher consumer biomass. Finally, there was a sharp decrease in the coexistence region for high plasticity (Fig. 3e). Region II, where B outcompetes A through maladaptive switching, doubled in size due to the much higher degree of maladaptive switching (Fig. 3g,h). Region I, where A outcompetes B, now also increased, when the level of defence of Ad is relatively low (Fig. 3e). This is again an indirect effect of maladaptive switching causing a strong increase in consumer biomass, affecting the cost/ benefit balance of defence: while Bd derives a strong benefit from its high level of defence, Bu is completely undefended, and is at an extra disadvantage because of its low growth rate. Thus, while Bd would have been able to survive by itself, the high exchange rate causes a strong source-sink dynamic that drives B extinct.Effect of plasticity in constellations crossing and angle
In constellation crossing the trade-off lines of both species cross in the trait space, as the level of defence is the same for both defended phenotypes; species B has a lower growth rate for its undefended phenotype than species A due to plasticity costs, while its defence costs are low and thus the growth rate of its defended phenotype is higher than for species A (Table 1, Supplementary Figure S5). Without plasticity the crossing trade-off lines lead to coexistence of both species in all simulations as Au and Bd were always the only survivors, mostly showing antiphase oscillations (Supplementary Figure S5).Allowing for phenotypic plasticity has the same results as were observed for constellation parallel: consumer biomass sharply increases (Fig. 4a,e); antiphase cycles are dampened or absent; and the area of coexistence decreases (Fig. 4b,f). All these changes are more pronounced for higher exchange rates (cf. Fig. 4a,b,e,f). Again, the biomass of the defended phenotypes decreased for high exchange rates (Supplementary Figure S6). Switching was always maladaptive for high exchange rates (Fig. 4g,h), and mostly maladaptive for low exchange rates (Fig. 4c,d). As was seen for constellation parallel, maladaptive switching was the reason for the decrease in coexistence. B can outcompete A when B has low plasticity costs. Bd has a much higher growth rate than Ad, while the undefended phenotypes have similar growth rates. The direction of competitive exclusion between Au and Bu is thus easily reversed by Bd donating biomass to the sink Bu, allowing B to occupy both niches and outcompete A (region II in Fig. 4b,f). The same mechanism happens in reverse for high plasticity and defence costs of B: the differences in growth rate for the undefended phenotypes are high, while the defended phenotypes have very similar growth rates. Au can support Ad, and A outcompetes B (region III in Fig. 4b,f).Figure 4Coexistence and maladaptive switching for scenario crossing 0.01 (a-d) and crossing 1 (e–h). Consumer biomass (a,e), the autotroph coexistence patterns (b,f), and the autotrophs’ maladaptive switching Φ (c,d,g,h) (higher biomasses or more intensive maladaptive switching are shown by darker colours). Lines separate the regions I–III of different autotroph coexistence. The x- and y-axis are reversed to follow the pattern of increasing fitness. Shaded areas in (b) depict antiphase cycles. Grey areas in (c,d,g,h) depict areas where the species was extinct. Shaded grey areas depict areas without simulations (cf. “Methods”). Note that (c,d,g,h) have each a different colour scale.Full size imageIn constellation angle there are no plasticity costs, and thus the undefended phenotypes Au and Bu have identical growth rates. The defended phenotypes take the same places in trait space as in the parallel constellation: Ad is a generalist, with a lower level of defence and a relatively high growth rate due to low defence costs, whereas Bd is a defence-specialist with a high level of defence but a low growth rate. This leads to the trade-off lines forming an angle (see Table 1). Without phenotypic plasticity, the coexistence patterns are the same as in constellation parallel, except that no competitive exclusion occurs between the undefended phenotypes; instead, they (neutrally) coexist in regions Ib, II and III (Supplementary Figure S7; cf. Fig. 2).With plasticity, neutral coexistence vanished: the defended phenotype that survived (Ad in region Ib, Bd in region III) could support the undefended phenotype of its own species, driving the other species extinct (Fig. 5b,f). As in the other constellations, the area of coexistence and the biomasses of the defended phenotypes decreased and antiphase cycles vanished with increasing χmax (Fig. 5b,f, Supplementary Figure S8), while maladaptive switching and the consumer biomass increased (Fig. 5).Figure 5Coexistence and maladaptive switching for scenario angle 0.01 (a–d) and angle 1 (e–h). Consumer biomass (a,e), the autotroph coexistence patterns (b,f), and the autotrophs’ maladaptive switching Φ (c,d,g,h) (higher biomasses or more intensive maladaptive switching are shown by darker colours). Lines separate the regions I–III of different autotroph coexistence. The y-axis is reversed to follow the pattern of increasing fitness. Shaded areas in (b) depict antiphase cycles. Grey areas in (c,d,g,h) depict areas where the species was extinct. Note that (c,d,g,h) have each a different colour scale.Full size imageGeneral resultsAs plasticity had very similar effects across all three constellations, we here generalize our results: we compare the three constellations for exchange rates over 5 orders of magnitude, as well as the non-plastic scenario and the rigid scenario (Table 1). That is, all simulations from one scenario (e.g. parallel 0) were summarized into one bar respective point in Fig. 6.Figure 6General patterns for coexistence, maladapative switching and biomasses. Share of surviving species in percent (A, B, coexistence or neutral coexistence) (a–c), median absolute value of maladaptive switching Φ (d–f) and median of total autotroph biomass (A + B), median consumer biomass C and share of defended phenotypes ((Ad + Bd)/(A + B)) (g–i) for the three constellations and increasing maximum exchange rates χmax. χmax = 0 denotes the non-plastic scenarios; *denotes the rigid scenarios. Maladaptive switching and the share of defended phenotypes do not apply for the rigid scenarios.Full size imageFor all constellations, the fraction of simulation runs leading to coexistence was highest in the non-plastic scenario and decreased with increasing χmax (Fig. 6a–c). In constellation parallel the share of coexistence for increasing χmax continuously decreased from 51 to 3% (Fig. 6a). In crossing, the share decreased from full to no coexistence (Fig. 6b). In angle, the share of coexistence was 88% in the non-plastic scenario when taking also neutral coexistence into account (Fig. 6c). Its share decreased to 9% for a χmax of 10 and increased again to 25% for the rigid scenario. Maladaptive switching increased for both species and all constellations for increasing χmax (Fig. 6d–f). The increased plasticity led to a lower total autotroph biomass and a lower share of defended phenotypes (Fig. 6g–i), which resulted in higher consumer biomass (Fig. 6g–i).Interestingly, and counterintuitively, the above patterns show that increasing the speed of plasticity (by increasing χmax) makes the system behave more like the rigid system. The coexistence patterns in scenarios with high χmax approach those of the rigid scenarios in two of the constellations (Fig. 6a,b). Similarly, the total autotroph and consumer biomasses approach the ones in the rigid scenarios (Fig. 6g–i). Thus, we found the higher χmax make the autotrophs not more adaptive, but behave more like non-adaptive species. More

Insects and plantsLarvae of fall armyworm (FAW, Spodoptera frugiperda) were cultured at the Department of Entomology at the Max Planck Institute for Chemical Ecology, and reared on a semi-artificial diet based on pinto bean59, and maintained under controlled light and temperature conditions (12:12 h light/dark, 21 °C).Feeding experiments3rd–4th instar FAW larvae were utilized for all experiments. Insects were starved overnight prior to feeding experiments. The following day insects were fed with a semi-artificial, pinto bean-based diet or put on maize leaves in small plastic cups and allowed to feed on the respective diets for a day. Insects were dissected in cold phosphate buffered saline (PBS, pH = 7.4) to harvest larval tissues (guts, Malphigian tubules, fat bodies, cuticle), which were stored at − 80 °C until further use. For droplet feeding, 12.5 mM DIMBOA was prepared by dissolving the compound in DMSO. This DIMBOA solution was further diluted in 10% aqueous sucrose solution. The larvae were stimulated with forceps to encourage regurgitation, and 2 μL DIMBOA-sucrose solution was administered directly to the larval mouthparts. Insects were then fed on semi-artificial diet for up to 6 h; following which gut tissue was dissected using cold phosphate buffer and the tissue samples were stored at − 80 °C until further use.Insect cell culturesSpodoptera frugiperda Sf9 cells and Trichoplusia ni Hi5 cells were cultured in Sf-900 II serum-free medium (Gibco) and ExpressFive serum-free medium (Gibco), respectively. Adherent cultures were maintained at 27 °C, and sub-cultured every 3–4 days.Cell treatmentsInsect cells were seeded in 6 well culture plates (Corning) and left at 27 °C overnight. For transcript stability tests, a fresh cycloheximide (CHX) stock (50 mg/mL) was prepared in ethanol and added to the cultured cells at a concentration of 50 µg/mL. Incubations with CHX were performed up to 6 h. For testing substrate specificity, cells were then treated with the following compounds for 1 h—DIMBOA (25–100 μM), indole (50–100 μM), quercetin (50–100 μM), and esculetin (50–100 μM). All the stocks were prepared in DMSO and cells treated with the corresponding volume of pure DMSO served as a control. The range of concentrations used for the substrates was based on previous work38.RNA extraction, reverse transcription and real time-PCR analysisTissue samples from the larvae were homogenized and total RNA extracted using the innuPREP RNA Mini Kit (Analytik Jena). Cell cultures used for RNA extraction were obtained during sub-culturing at full confluency, and centrifuged at 500×g for 5 min. The culture medium was discarded, and the fresh pellets were directly used for RNA extraction. RNA concentrations were measured with the NanoDrop 2000 UV–Vis Spectrophotometer (Thermo Scientific). First strand cDNA was synthesized from 1 μg total RNA using SuperScript III Reverse Transcriptase and OligodT primers from Invitrogen. Sequences were successfully amplified using Phusion High Fidelity DNA Polymerase (New England Biolabs) (PCR protocol: 30 s at 98 °C; 35 cycles of 10 s at 98 °C, 20 s at 55 °C, 45 s at 72 °C; and 5 min at 72 °C). The PCR products were purified with a PCR cleanup kit (Qiagen) and cloned into pCR-Blunt II-TOPO vector (Life Technologies) and transformed into NEB cells (Life Technologies), which were plated on selective LB agar medium containing 100 μg/mL ampicillin and incubated overnight at 37 °C. Positive colonies were identified by PCR using vector-specific M13 primers. Positive clones were confirmed by sequencing. Real time PCR analyses were carried out using Brilliant III SYBR Master Mix, employing SYBR Green chemistry. Relative quantification of the transcript levels was done using the 2−∆∆Ct method60. SfRPL10 was used as reference gene for all analyses. The primer pairs used for distinguishing between the variants are listed in Supplementary Table 1. As the expression of full-length and variants of SfUGT33F28 differed according to the strains, tissues, and treatments being analyzed, variant expression is reported as ratios relative to the canonical transcript to facilitate comparisons.Preparation of minigenes for alternative splicing studiesGenomic DNA was isolated from S. frugiperda larvae using the cetyl trimethyl ammonium bromide (CTAB) protocol61. DNA concentration was measured with the NanoDrop 2000 UV–Vis Spectrophotometer (Thermo Scientific). The minigene was amplified using Phusion High Fidelity DNA Polymerase (New England Biolabs) (PCR protocol: 30 s at 98 °C; 35 cycles of 10 s at 98 °C, 30 s at 55–60 °C, 1 min 30 s at 72 °C; and 10 min at 72 °C), cloned into a pCR-Blunt II-TOPO vector (Life Technologies) and sequenced using M13 primers. The confirmed sequence was eventually cloned into a pIB/V5-His-TOPOvector (Life Technologies) and transformed into NEB cells (Life Technologies). Positive colonies were identified by colony PCR using vector-specific OpIE2 primers, sub-cultured overnight at 37 °C in liquid LB medium containing 100 μg/mL ampicillin and used for plasmid DNA purification with the NucleoSpin Plasmid kit (Macherey-Nagel). Concentration and purity of the obtained construct was assessed by the NanoDrop 2000 UV–Vis Spectrophotometer (Thermo Scientific) and the correct orientation of the PCR products was confirmed by DNA sequencing.Nuclear protein isolationNuclear proteins were isolated from insect cells62 using the protocol originally described with few modifications. Cells grown to concentrations of up to 1 × 106 cells/well were harvested and washed with PBS (pH 7.4). The extracts were centrifuged at 12,000×g for 10 min and pellets were re-suspended in 400 μL cell lysis buffer (10 mM HEPES, pH 7.5, 10 mM KCl, 0.1 mM EDTA pH 8.0, 1 mM DTT, 0.5% Nonidet-40 and 10 μL protease inhibitor cocktail). Cells were allowed to swell on ice for 20 min with intermittent mixing. Suspensions were vortexed to disrupt the cell membranes and then centrifuged at 12,000×g for 10 min at 4 °C. Pelleted nuclei were washed thrice with cell lysis buffer, re-suspended in 50 μL nuclear extraction buffer (20 mM HEPES pH 7.5, 400 mM KCl, 1 mM EDTA pH 8.0, 1 mM DTT, 10% glycerol and protease inhibitor) and incubated on ice for 30 min. Nuclear fractions were collected by centrifugation at 12,000g for 15 min at 4 °C. Protein concentrations were measured by Bradford and extracts were stored at − 80 °C until further use.Electrophoretic mobility shift assay (EMSA)EMSA was performed using the LightShift Chemiluminescent EMSA kit (Thermo Scientific) following the manufacturer’s instructions. Genomic DNA fragments of 20–25 bp corresponding to the 5′ flanking region of UGT33F28 exon 1 (with and without AhR-ARNT motif deletion) were synthesized with covalently linked biotin (Sigma). The DNA probes used in the experiment are listed in Supplementary Table 6. EMSA was performed in 20 µL reactions containing 20 fmol biotinylated DNA probe with 3.5–4 µg nuclear protein extracted from insect cells, according to manufacturer’s instructions. A reaction comprising the above along with the excess of unlabeled canonical DNA probe (200 molar excess) was further employed as a control. The reaction was assembled at room temperature and incubated for 30 min. The reactions were separated on a 5% TBE gel in 0.5X TBE at 100 V for 60 min. The samples were then transferred to a positively charged nylon membrane (Hybond N+, Amersham Bioscience) using semi-dry transfer at 15 V for 30 min. The membrane was cross-linked for 1 min using the auto cross-link function on the UV cross-linker (Stratagene). The biotinylated DNA–protein complex was detected by the streptavidin–horseradish peroxidase conjugated antibody provided in the kit. The membrane was washed and incubated with the chemiluminescence substrate for 5 min and the signals were developed by exposing the membrane to an X-ray film for 1 min.Streptavidin affinity purificationStreptavidin agarose (Sigma-Aldrich) was employed for protein purification. Briefly, 50–100 μL of agarose was packed into a 1.5 mL Eppendorf tube for each sample. The agarose was allowed to settle with a short centrifugation (500×g, 5 min) and the supernatant was discarded. The agarose was washed 4–5 times with binding buffer (PBS containing 1 mM EDTA, 1 mM DTT, 4 µg poly dI. dC as non-specific competitor DNA and protease inhibitor). Simultaneously, the binding reaction with the nuclear protein fraction and the DNA probe was assembled as described above. A 100 μg amount of total nuclear protein was incubated with 4 μg of biotinylated DNA probe at room temperature for 20 min. The reaction was loaded onto the streptavidin column equilibrated with the binding buffer and incubated for another 1 h at room temperature with gentle shaking. Subsequently, the agarose was washed 4–5 times with the binding buffer. After the final wash, the supernatant was aspirated and 10 μL was left above the beads. For protein separation, 20–30 μL pf the SDS loading buffer was added onto the agarose, boiled at 95 °C for 5 min and the sample thus obtained was utilized for electrophoresis.Deletion mutagenesisFor deletion mutagenesis, a pair of primers flanking the sequence to be deleted (non-overlapping) was designed. The pCR-Blunt II-TOPO vector (Life Technologies) clone for the SfUGT33F28 exon 1–2 minigene was utilized as a template. Sequence was successfully amplified using Phusion High Fidelity DNA Polymerase (New England Biolabs) (PCR protocol: 30 s at 98 °C; 20 cycles of 10 s at 98 °C, 30 s at 55–60 °C, 4 min at 72 °C; and 10 min at 72 °C). A DpnI digest was performed to remove the background DNA, followed by ligation and transformation into fresh cells. The sequence of the mutant TOPO clone was then confirmed and utilized as a template for cloning into pIB/V5-His-TOPO vector (Life Technologies) for transfection into insect cells.Cloning and heterologous expression of SfUGTsSequences were amplified from S. frugiperda gut cDNA samples using Phusion High Fidelity DNA Polymerase (New England Biolabs) (PCR protocol: 30 s at 98 °C; 35 cycles of 10 s at 98 °C, 20 s at 55–60 °C, 45 s at 72 °C; and 5 min at 72 °C). The resulting amplified products were purified with a PCR cleanup kit (Qiagen) and incubated with GoTaq DNA polymerase (Promega) for 15 min at 72 °C in order to add A overhangs. The products were cloned into the pIB/V5-His-TOPO vector (Life Technologies) and transformed into NEB cells (Life Technologies), which were plated on selective LB agar medium containing 100 μg/mL ampicillin and incubated overnight at 37 °C. Positive colonies were identified by PCR using vector-specific OpIE2 primers, sub-cultured overnight at 37 °C in liquid LB medium containing 100 μg/mL ampicillin and used for plasmid DNA purification with the NucleoSpin Plasmid kit (Macherey-Nagel). Concentration and purity of the obtained constructs were assessed by NanoDrop 2000 UV–Vis Spectrophotometer (Thermo Scientific) and the correct orientation of the PCR products was confirmed by DNA sequencing.Insect cell transfectionFor transfection, Sf9 cells and Hi5 cells were sub-cultured at full confluency in a 6-well plate in a 1:3 dilution and left overnight to adhere to the flask surface. The medium was replaced, and transfections were carried out using FuGENE HD Transfection Reagent (Promega) in a 1:3 plasmid/lipid ratio (1.7 μg plasmid and 5.0 μL lipid for 3 mL medium). Cells were incubated for 48–72 h at 27 °C and re-suspended in fresh medium containing 50 μg/mL blasticidin for 2 weeks. Stable cell cultures were subsequently maintained at 10 μg/mL blasticidin.Cell lysate preparationCells were obtained from cultures 2 weeks post transfection growing stably on 50 μg/mL blasticidin. A 1 mL quantity of cells was harvested for each construct and re-suspended into 100 µL buffer. Protein concentrations were measured using the Bradford reagent, and 1–2 μg of the cell lysate was used for enzyme assays.Microsome preparationFor microsome extraction, confluent, stably transfected cells from five T-75 flasks (10 mL culture) per recombinant plasmid were harvested by scraping the cells off the bottom using a sterile cell scraper (Sarstedt AG, Nuembrecht, Germany). The obtained cell suspensions were combined into a 50 mL falcon tube and centrifuged at 1000×g for 15 min at 4 °C (AvantiTM J-20 XP Centrifuge, Beckman Coulter, Krefeld, Germany). The supernatant was discarded, the cells were washed twice with ice-cold PBS buffer (pH 7.4) and centrifuged at 1000×g for 15 min. The resulting cell pellet was re-suspended in 10 mL hypotonic buffer (20 mM Tris, 5 mM EDTA, 1 mM DTT, 20% glycerol, pH 7.5), containing 0.1% BenzonaseR nuclease and 100 μL Protease Inhibitor Cocktail (Serva) followed by incubation on ice for 30 min. After cell lysis, the cells were homogenized by 20–30 strokes in a Potter–Elvehjem tissue grinder (Kontes Glass Co., Vineland, USA) and were subsequently mixed with an equal volume of sucrose buffer (20 mM Tris, 5 mM EDTA, 1 mM DTT, 500 mM sucrose, 20% glycerol, pH 7.5). The homogenate was centrifuged at 1200×g and 4 °C for 10 min (AvantiTM J-20 XP Centrifuge, Beckman Coulter), and the supernatant was transferred into Beckman polycarbonate ultracentrifugation bottles (25 × 89 mm) (Beckman Coulter) and centrifuged at 100,000×g and 4 °C for 1.5 h in a fixed angle Type 70 Ti rotor (OptimaTM L-90K Ultracentrifuge, Beckman Coulter). After ultracentrifugation, the clear supernatant, containing the cytosolic fraction, was aliquoted into 1.5 mL Eppendorf tubes. The pellet, containing the microsomal fractions, was re-suspended in 1 mL of phosphate buffer (100 mM K2HPO4, pH 7.0), containing 10 μL Protease Inhibitor Cocktail (Serva) and stored at − 80 °C until further use. Typically, 5–10 μg of the microsome fraction so obtained was utilized for the enzyme assays.Cross-linking assaysCross-linking assays were performed using dimethyl suberimidate (DMS) as the cross-linking agent. A fresh stock of DMS (5 mg/mL) was prepared in 0.2 M triethanolamine (pH 8.0) at the start of each assay. DMS was added to a final concentration of 2.5 mg/mL to insect cell microsomes with gentle shaking up to 3 h, and samples were subsequently stored at − 20 °C until further use. All protein samples were electrophoresed using a 12% Mini-PROTEAN tris glycine gel, blotted onto PVDF membrane using wet transfer at 70 V for 30–45 min, followed by detection using the V5-HRP conjugate.V5-based affinity purificationAnti-V5 agarose affinity gel (Sigma-Aldrich) was employed for protein purification. Briefly, 50–75 μL of the agarose was packed into a 1.5 mL Eppendorf tube for each sample. The agarose was allowed to settle with a short centrifugation and the supernatant was discarded. The agarose was washed 4–5 times with PBS (pH 7.4). Samples to be purified were incubated with 5% digitonin on ice for 20 min and subject to centrifugation at 16,000×g for 30 min. Clarified cell lysate or microsomal extract was added onto the resin (up to 200 μL, volume adjusted by addition of PBS) and incubated for 1.5 h on a shaker. Subsequently, the agarose was washed 4–5 times with PBS. After the final wash, the supernatant was aspirated and 10 μL was left above the beads. This fraction was used for both protein electrophoresis and enzyme assays (separate purifications). For SDS-PAGE, 20–30 μL pf the SDS loading buffer was added onto the agarose, boiled at 95 °C for 5 min and sample thus obtained was utilized for electrophoresis.LC–MS/MS peptide analysisProtein bands of Coomassie Brilliant blue R250 stained gels were cut from the gel matrix and tryptic digestion was carried out63. For LC–MS/MS analysis of the resulting peptides, samples were reconstituted in 20 μL aqueous 1% formic acid, and 1 μL was injected onto an UPLC M-class system (Waters, Manchester, UK) coupled to a Synapt G2-si mass spectrometer (Waters, Manchester, UK). Samples were first pre-concentrated and desalted using a Symmetry C18 trap column (100 Å, 180 µm × 20 mm, 5 µm particle size) at a flow rate of 15 µL/min (0.1% aqueous formic acid). Peptides were eluted onto a ACQUITY UPLC HSS T3 analytical column (100 Å, 75 µm × 200 mm, 1.8 µm particle size) at a flow rate of 350 nL/min with the following gradient: 3–15% over 3 min, 15–20% B over 7 min, 20–40% B over 30 min, 40–50% B over 5 min, 50–70% B over 5 min, 70–95% B over 3 min, isocratic at 95% B for 1 min, and a return to 1% B over 1 min. Phases A and B were composed of 0.1% formic acid and 100% acetonitrile in 0.1% formic acid, respectively). The analytical column was re-equilibrated for 10 min prior to the next injection. The eluted peptides were transferred into the mass spectrometer operated in V-mode with a resolving power of at least 20,000 full width at half height FWHM. All analyses were performed in a positive ESI mode. A 100 fmol/μL sample of human Glu-Fibrinopeptide B in 0.1% formic acid/acetonitrile (1:1 v/v) was infused at a flow rate of 1 μL/min through the reference sprayer every 45 s to compensate for mass shifts in MS and MS/MS fragmentation mode. Data were acquired using data-dependent acquisition (DDA). The acquisition cycle for DDA analysis consisted of a survey scan covering the range of m/z 400–1800 Da followed by MS/MS fragmentation of the ten most intense precursor ions collected at 0.5 s intervals in the range of 50–2000 m/z. Dynamic exclusion was applied to minimize multiple fragmentations for the same precursor ions. MS data were collected using MassLynx v4.1 software (Waters, Manchester, UK).Data processing and protein identificationDDA raw data were processed and searched against a sub-database containing common contaminants (human keratins and trypsin) using ProteinLynx Global Server (PLGS) version 2.5.2 (Waters, Manchester, UK). Spectra remaining unmatched by database searching were interpreted de novo to yield peptide sequences and subjected to homology-based searching using the MS BLAST program64 installed on a local server. MS BLAST searches were performed against a Spodoptera frugiperda database obtained by in silico translation of the S. frugiperda transcriptome37 and against arthropoda database (NCBI). PKL-files of MS/MS spectra were generated and searched against Spodoptera frugiperda database combined with NCBI nr (downloaded on May 24, 2020) using MASCOT software version 2.6.2. The following searching parameters were applied: fixed precursor ion mass tolerance of 15 ppm for the survey peptide, fragment ion mass tolerance of 0.1 Da, 1 missed cleavage, fixed carbamidomethylation of cysteines and possible oxidation of methionine.Enzymatic assaysFor UGT assays, samples from insect cell cultures (transient or stable) were prepared in phosphate buffer (pH 7.0, 100 mM). Typical enzyme reactions included 5–10 µg cell microsomal extracts, 2 μL of 12.5 mM DIMBOA in DMSO (25 nmol), 4 μL of 12.5 mM UDP-glucose in water (50 nmol), and phosphate buffer (pH 7.0, 100 mM) to give an assay volume of 50 μL. Controls containing either boiled enzymatic preparation, or only the protein suspension and buffer were included. After incubation at 30 °C for 60 min, the enzyme reactions were interrupted by adding 50 μL of 1:1 (v:v) methanol/formic acid solution. For enzyme assays involving resin purified microsomal extracts, equal amounts of extracts were employed for resin purification and the enzyme assay (buffer + substrate) was pipetted directly onto the resin. Post incubation, samples were centrifuged, supernatant was collected, and reaction was stopped by addition of methanol/formic acid solution. Assays were centrifuged at 5000g for 5 min and the obtained supernatant was collected and analyzed by LC–MS/MS.Chromatographic methodsFor all analytical chromatography procedures, formic acid (0.05%) in water and acetonitrile were used as mobile phases A and B, respectively, and the column temperature was maintained at 25 °C. Analyses of enzymatic assays and plant samples used a Zorbax Eclipse XDB-C18 column (50 × 4.6 mm, 1.8 μm, Agilent Technologies) with a flow rate of 1.1 mL/min and with the following elution profile: 0–0.5 min, 95% A; 0.5–6 min, 95–67.5% A; 6.02–7 min, 100% B; 7.1–9.5 min, 95% A. LC–MS/MS analyses were performed on an Agilent 1200 HPLC system (Agilent Technologies) coupled to an API 6500 tandem spectrometer (AB Sciex) equipped with a turbospray ion source operating in negative ionization mode. Multiple reaction monitoring (MRM) was used to monitor analyte parent ion to product ion conversion with parameters from the literature for DIMBOA65 and DIMBOA-Glc16. Analyst (version 1.6.3, Applied Biosystems) software was used for data acquisition and processing.Statistical analysisAll statistical analyses were carried out using SigmaPlot 12.0 and R studio (version 3.6.3). Data were tested for homogeneity of variance and normality and were appropriately transformed to meet these criteria where required. The specific statistical method used for each data set is described in the figure legends. More

The conceptual flow chart of the process is provided in Fig. 1. We used seven reanalysis SM data (Table 2) masked with soil temperature (ST) and soil freeze/thaw status to calculate water table depth, i.e. the input of TOPMODEL, given the obvious disagreements between the input datasets. The diagnostic algorithms based on TOPMODEL were used following Stocker et al. (ref. 20) and Xi et al. (ref. 25), where the optimized parameters were calibrated with long-term maximum wetland areas from four observation-based wetland datasets (Table 1). Details about these datasets and computational processing are shown as follows.Fig. 1Diagram of workflow for parameter calibration and the simulation of global wetland dynamics.Full size imageTable 2 Key characteristics of seven global soil moisture reanalysis data used in this study.Full size tableReanalysis soil moisture datasetsSeven long-term reanalysis SM datasets used in this study include NCEP-DOE (National Centers for Environmental Prediction-the Department of Energy)26, MERRA-Land (the Modern-Era Retrospective Analysis for Research and Applications)27, MERRA-2 (ref. 28), GLDAS-Noah v2.0 (the Global Land Data Assimilation System)29, GLDAS-Noah v2.1 (ref. 29), ERA5 (European Environment Agency)30,31, and ERA5-Land30,31. Key characteristics of the seven SM data are listed in Table 2. The datasets differ by their spatial and temporal resolutions, the time-period they cover, as well as the definition of the soil layers. More details are provided for each dataset below.

NCEP-DOE
NCEP-DOE is an updated version of the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) Reanalysis 1 project, which uses a state-of-the-art analysis/forecast system to perform data assimilation with past data from 1948 to the present32. NCEP-DOE features the newer physics and observed SM forcing and also eliminates several previous errors, such as oceanic albedo and snowmelt term during the entire period, and snow cover analysis error from 1974 to 1994 (ref. 26). With a spatial resolution of about 210 km, there are two vertical soil layers in NCEP-DOE for both SM and ST: 0–0.1 and 0.1–2 m.

MERRA-Land and MERRA-2
MERRA-Land soil moisture is generated by driving the Goddard Earth Observing System model version 5.7.2 (GEOS-5.7.2) with meteorological forcing from the MERRA reanalysis product27. The precipitation forcing in MERRA-Land merges MERRA precipitation with a gauge-based data product from the National Oceanic and Atmospheric Administration (NOAA) Climate Prediction Center, and the Catchment land surface model used in MERRA-Land is updated to the “Fortuna-2.5” version27. MERRA-2 intends to replace the original MERRA reanalysis and ingests important new data types28. The Catchment model in MERRA-2 has been updated with rainfall interception and snow model parameters of MERRA-Land, and the precipitation correction is a refined version of MERRA-Land. For MERRA-Land and MERRA-2, there is only one layer for SM from the surface to the bedrock, with “depth-to-rock” depending on local conditions. ST is computed on six vertical soil layers: 0–0.10, 0.10–0.29, 0.29–0.68, 0.68–1.44, 1.44–2.95, and 2.95–12.95 m.

ERA5 and ERA5-Land
ERA5 is the fifth generation ECMWF (European Centre for Medium-Range Weather Forecasts) reanalysis of global climate and weather, replacing ERA-Interim30,31. Based on a decade of developments in model dynamics and data assimilation, there is a significantly enhanced horizontal resolution (31 km), temporal resolution (hourly) and uncertainty estimation. ERA5 covers 1979–2020 and continues to be updated in near-real-time. ERA5-Land is produced with a finer horizontal resolution of 9 km by running the land component of the ERA5 climate reanalysis but without data assimilation. By March of 2021, the ERA5-Land outputs are only available since 1981. SM and ST are computed on four vertical soil layers (0–0.07, 0.07–0.28, 0.28–1, and 1–2.89 m) for both ERA5 and ERA5-Land.

GLDAS-Noah v2.0 and GLDAS-Noah v2.1

GLDAS is a global, moderate-resolution (0.25° × 0.25°) offline terrestrial modeling system developed by NASA Goddard Space Flight Center (GSFC) and the NOAA National Centers for Environmental Prediction29, thus similar to ERA5. To produce optimal fields of land surface variables in near-real-time, it incorporates satellite- and ground-based observations. GLDAS-Noah drives the Noah land surface model and has two components: one forced with the Princeton meteorological forcing data (i.e. GLDAS-Noah v2.0) and the other forced with a combination of model and observation (i.e. GLDAS-Noah v2.1). GLDAS-Noah v2.0 covers the period 1948–2014, while GLDAS-Noah v2.1 is available from 2000 to the present. There are four vertical layers in the Noah land surface model for both ST and SM: 0–0.1, 0.1–0.4, 0.4–1, and 1–2 m.Observation-based wetland/flooded area dataIn terms of large uncertainties in current wetland datasets (Table 1) we selected four widely used and available satellite/satellite-based wetland/flooded area data including GIEMS-2 (ref. 14), RFW (the Regularly Flooded Wetland map)10, WAD2M (a global dataset of Wetland Area and Dynamics for Methane Modeling)33, and G2017 (the pantropical wetland extent from an expert system model)9 for parameter calibration. Among them, GIEMS-2 and WAD2M include monthly wetland dynamics, while RFW and G2017 are static. The comparison of the four wetland datasets is shown in Supplementary Fig. 1; details on each data are provided below.

GIEMS-2
The GIEMS-1 is the first global estimate of monthly inundated areas, derived from passive microwave land surface emissivity34. With a 0.25° × 0.25° resolution, GIEMS-1 documents a mean annual maximum inundated area of 9.5 Mkm2 for 1993–2007 (including open water, wetlands, and rice paddies, but excluding large lakes), which shows good agreement with existing independent, static inventories as well as regional high-resolution synthetic aperture radar observations34. Based on similar retrieval principles with GIEMS-1, GIEMS-2 is developed to less depend on ancillary data with an updated microwave emissivity, and correct a known overestimation over low vegetated areas from GIEMS-1 (ref. 14). The period is extended to 1992–2015 for GIEMS-2 and can be updated with the availability of observations. Globally, the mean annual maximum and long-term maximum inundated extent after removing the rice paddies using the Monthly Irrigated and Rainfed Crop Areas dataset (MIRCA2000)35 for the period 1992–2015 are 6.7 and 10.9 million km2 (hereafter Mkm2; sum of mean annual maximum or long-term maximum inundated extent for each grid cell) respectively. The rice paddies are removed here as they are not natural wetlands and cannot be simulated with TOPMODEL.

RFW
RFW is a static, high-resolution map (15 arc-sec) of regularly flooded wetlands, developed by overlapping flooded areas (permanent wetlands and flooded vegetation classes) for 2008–2012 from the ESA-CCI land cover map36, mean annual maximum inundated areas (including wetlands, rivers, small lakes, and irrigated rice) for 1993–2004 from GIEMS-D15 global inundation extent (downscaled using GIEMS-1)37, and long-term maximum surface water areas for 1984–2015 from JRC global surface water bodies product13. The large permanent lakes and reservoirs are distinguished using the HydroLAKES database38. Globally, RFW covers 9.7% of the land surface area (~13.0 Mkm2) including wetlands, river channels, deltas, and flooded lake margins, but excluding large lakes10. Due to the mean annual maximum or long-term maximum inundation/surface water extent for 1984–2016 from the three wetland data is used, we treated RFW as the long-term maximum wetland extent in this study. Besides, given that GIEMS-D15 includes artificial rice paddies, we removed them with MIRCA2000 from RFW (~11.9 Mkm2 after removing rice paddies).

WAD2M
WAD2M dataset used in this study is an improved version of the SWAMPS v3.2 from Jensen et al. (ref. 15), covering the years 2000 to 2018. With a spatial resolution of 25 km × 25 km, this data was used as input wetland area data of phase 2 of the Global Methane Budget33. Given that the initial SWAMPS failed to detect wetlands lacking surface inundation and to differentiate between lakes, wetlands, and other surface water bodies, Zhang et al. (ref. 33) modified it using a series of independent static wetland distribution data7,9,39,40,41 in an attempt to include missing wetlands under dense canopies. Besides, they removed inland waters (lakes, rivers, and ponds) and rice agriculture with JRC and MIRCA2000, respectively. Globally, the mean annual maximum and long-term maximum wetland extent for the period 2000–2018 estimated by WAD2M are 8.1 Mkm2 and 13.2 Mkm2 (sum of mean annual maximum or long-term maximum inundated extent for each grid cell) respectively.

G2017

G2017 (ref. 9) is a static, pantropical wetland and peatland extent map (covering 60°S–40°N) at 232 m × 232 m resolution, derived from a hybrid expert model system. With three biophysical indices related to wetland and peat formation (long-term water supply exceeding atmospheric water demand, annually or seasonally waterlogged soils, and geomorphological position where water is supplied and retained), G2017 identifies not only permanently and seasonally wetland areas, but also soil wetness and topographic conditions that favor waterlogging in the absence of flooding for the end of the 20th century. Given the broad coverage of different types of wetlands, we also treated this map as long-term maximum wetland areas. This ‘pantropical’ data (60°S to 40°N) offers the advantage to include non-flooded wetland areas that are missed in satellite-based wetland products. However, note that not all detected wetlands or peatlands in G2017 have been observed. Rice agriculture was also removed with MIRCA2000 from G2017. The resulting wetland and peatland area for 60°S–40°N is 4.0 Mkm2.The TOPMODEL-based diagnostic modelTOPMODEL as improved by Stocker et al. (ref. 20) and Xi et al. (ref. 25) was used to calculate the inundated fraction from WTD at grid-scale in this study. Based on the assumptions that the local hydraulic gradient is approximated by the local topographic slope and the water table variations can be assimilated to a succession of steady states with uniform recharge, the classical TOPMODEL establishes an analytical relationship between the soil moisture deficit and the distributions of local topographic index within a catchment. At grid-scale, the analytical relationship can be represented as:$$CT{I}_{i}-overline{CT{I}_{x}}=mathrm{-M}left({{Gamma }}_{i}-overline{{{Gamma }}_{x}}right)$$
(1)
where CTI indicates the topographic index, defined as the log of the ratio of contributing area to the local slope. We used the CTI data at 500 m × 500 m resolution from Marthews et al. (ref. 22), where lakes, reservoirs, mountain glaciers, and ice caps are removed using the Global Lakes and Wetlands Database7. The (overline{CT{I}_{x}}) indicates the average of CTIi of all sub-grids (index i) within the grid cell x. M indicates a tunable parameter that describes the exponential decrease of soil transmissivity with depth21. Γi is the water table of the pixel i and (overline{{{Gamma }}_{{x}}}) is the mean water table of the grid x. When Γi is at the soil surface (i.e. Γi = 0), the threshold (CT{I}_{x}^{* }) above which all pixels are flooded for the grid x is derived:$$CT{I}_{x}^{* }=overline{CT{I}_{x}}+{rm{M}}cdot overline{{{Gamma }}_{x}}$$
(2)
The wetlands area is defined as the flooded areas (i.e. Γ ≤ 0), the flooded fraction in the grid x (fx) being the percentage of pixels with CTIi larger than a threshold (CT{I}_{x}^{* }):$${f}_{x}=frac{1}{{A}_{x}}{sum }_{i}{A}_{i}^{* }$$with$${A}_{i}^{* }=left{begin{array}{c}{A}_{i},if,CT{I}_{i}ge CT{I}_{x}^{* }\ 0,if,CT{I}_{i} < CT{I}_{x}^{* }end{array}right.$$ (3) To reduce the computational costs from the high-resolution CTI data for predicting long time series of wetland area, we used the asymmetric sigmoid function from Stocker et al. (ref. 20) to fit the “empirical” relationship (widehat{{Psi }}) between (widehat{f}) and Γ:$${{rm{psi }}}_{x}left({{Gamma }}_{x}right)={left(1+{v}_{x}cdot {e}^{-{k}_{x}left({{Gamma }}_{x}-{q}_{x}right)}right)}^{-1/{v}_{x}}$$ (4) where vx, kx, qx are three parameters of the function. Given a value of parameter M, the three parameters can be derived with a sequence of Γx spanning a plausible range of values (−1 m to 2 m) and corresponding fx from the initial TOPMODEL approach (Eq. (3)). Thus, the wetlands in our study are defined as the flooded area simulated by TOPMODEL. As for the range of parameter M, Stocker et al. (ref. 20) used a global uniform value for M (M = 8) after testing simulated wetland fraction for a range of M (7, 8, 9). Nevertheless, given that distinct topography, soil types, and other intrinsic characteristics in different regions, we considered M as a tunable, spatially heterogeneous, and grid-specific parameter, with a range of 1–15 following Xi et al. (ref. 25). Thus, for each grid cell x there are 15 choices for M, and then 15 sets of (vx, kx, qx). The optimized parameter combination of (vx, kx, qx) is determined by selecting minimum root-mean-square-error (RMSE) between simulated inundated fractions and observations:$$RMSE=sqrt{frac{{sum }_{i=1}^{n}{left({O}_{i}-{P}_{i}right)}^{2}}{n}}$$ (5) where Oi and Pi are observed and simulated wetland fraction, respectively. n represents the time-series length for wetland extent. For simulations calibrated with RFW and G2017, the RMSE was computed with the long-term maximum (hereafter called MAX) monthly wetland area because the two data sets are static and only record the MAX wetland extent. While for simulations calibrated with GIEMS-2 and WAD2M which include temporal variations of wetland area, we calibrated the parameters with all months, mean seasonal cycle, yearly maximum, and MAX wetland area, but only showed the optimal simulations calibrated with MAX wetland area in this work to keep consistency with RFW and G2017. Besides, to provide more choices for users, we combined all of the four wetland datasets (i.e. the union of long-term maximum wetland extent) to generate a new wetland map (hereafter called MAX_all), and then used the MAX_all to calibrate the parameters to produce seven sets of global wetland extent products with seven soil moisture datasets. The simulations calibrated with yearly maximum wetland area from GIEMS-2 and WAD2M and long-term maximum wetland area from MAX_all are also provided in our resulting products.Finally, to avoid unrealistically high wetland fraction output from the function, the simulated maximum wetland fraction fx is constrained by the observed MAX wetland area with a parameter ({f}_{x}^{max}) (Eq. (6)), which is different from Stocker et al. (ref. 20). The determination of ({f}_{x}^{max}) is analyzed in the supplemental material in detail (Supplementary Text 1). Once the value of (vx, kx, qx) are determined, the wetland fraction fx can be directly derived from the monthly water table Γx according to Eqs. (4) and (6).$${f}_{x}=minleft({{Psi }}_{x}left({{Gamma }}_{x}right),{f}_{x}^{max}right)$$ (6) Calculation of water table depthWater table depth is not computed by land surface models, given their coarse soil vertical discretization. We thus used the saturation deficit of soil moisture (θSD) as a surrogate of water table depth, θSD being defined as an index consisting of saturated volumetric water content and the “actual” soil depth modified by soil freeze/thaw status:$${theta }_{SD}={z}_{{l}_{0}}-{sum }_{l=1}^{{l}_{0}}{theta }_{l}cdot frac{Delta {z}_{l}}{{theta }_{S}}$$ (7) Subscript l represents the lth soil layer, l0 is the number of layers above the first frozen soil layer counted from the top (l = 1 at the soil surface), θl is the monthly volumetric water content in the lth soil layer (m3 m−3), (Delta {z}_{l}) is the thickness of the lth soil layer (m), θS is the saturated volumetric water content (in m3 m−3 units, uniform over depth).As formulated in Eq. (7), ({z}_{{l}_{0}}) is the thickness of all soil layers (or depth to bedrock) when there is no frozen soil layer. If there exists at least one frozen layer, ({z}_{{l}_{0}}) is set to the depth of the uppermost frozen soil layer. We excluded the frozen soil layers here given that some important wetland processes such as methane production and transport are insignificant when the soils are frozen. In high latitudes, the presence of frozen soil layers may lead to an overestimation of the wetland fraction due to relatively large θSD values even if there is little liquid soil water above the uppermost frozen soil layer. Hence, we used monthly soil temperature (ST) at 70 cm, the Global Record of Daily Landscape Freeze/Thaw Status data42, and the Köppen climate classification system43 to refine the frozen mask. When the monthly mean ST at 70 cm is below 0 °C, or soil freezing days are more than 5 in a month, or the grid is classified as the Hot desert (BWh) in the Köppen climate classification system, the wetland fraction for the grid is set to zero. However, it should be noted that the algorithm using the ST at 70 cm could omit some unfrozen soil layers above 70 cm, which could lead to bias in estimation of methane emissions from these unfrozen layers. We provided the global wetland maps in our resulting products, but the potential uncertainties in wetland estimation due to the omitted unfrozen layers should be considered, particularly at high latitudes. We used seven reanalysis SM products to compute θSD to provide the uncertainty in SM input (Table 2). All data are re-interpolated to 0.25° × 0.25° resolution.Evaluation against wetland calibration data and independent satellite productsAlthough we calibrated parameters of the TOPMODEL-based diagnostic model with the observation-based wetland data, to what extent the simulations can reproduce the spatial patterns and temporal dynamics of the calibration wetland data must be evaluated. For spatial patterns, we calculated the RMSE of wetland area between our simulations and corresponding wetland calibration data following Eq. (5), and evaluated the spatial patterns of simulated wetland extent in two wetland hotspots including Amazon basin and Western Siberia lowlands with three independent global/regional water products. For Amazon basin, we used the global surface water dataset from JRC13 (optical satellite images) and the wetland map produced using mosaics of Japanese Earth Resources Satellite (JERS-1) L-band SAR imagery from Hess et al. (ref. 44, hereafter H2015). For West Siberian lowlands, we used JRC and the Boreal–Arctic Wetland and Lake Dataset (BAWLD, only covers the north of ~55°N) produced using an expert assessment and extrapolated using random forest modelling from climate, topography, soils, permafrost conditions, vegetation, wetlands, and surface water extents and dynamics45. For temporal dynamics, since we only used the static wetland area (long-term maximum) from all of the four observation-based wetland products to calibrate parameters, the simulated temporal dynamics can be evaluated with the two dynamic wetland products (GIEMS-2 and WAD2M). Besides, we also used the terrestrial water storage (TWS) from the Gravity Recovery and Climate Experiment (GRACE), which retrieves relative change in TWS from the monthly anomalies of the Earth’s gravity field for 2003–2016 measured by the twin GRACE satellites46,47 to evaluate the simulated temporal dynamics. More

Liu, H. et al. Optimal nitrogen input for higher efficiency and lower environmental impacts of winter wheat production in China. Agr. Ecosyst. Environ. 224, 1–11 (2016).Article 

Google Scholar 
Guang-hao, L., Gui-gen, C., Wei-ping, L. & Da-lei, L. Differences of yield and nitrogen use efficiency under different applications of slow-release fertilizer in spring maize. J. Integr. Agric. 20(2), 554–564 (2021).Article 

Google Scholar 
Kumar, V. V. Role of Rhizospheric Microbes in Soil 377–398 (Springer, 2018).Book 

Google Scholar 
Ullah, A. et al. Factors affecting the adoption of organic farming in Peshawar-Pakistan. Agric. Sci. 6(06), 587–593 (2015).
Google Scholar 
Cui, Z. et al. Pursuing sustainable productivity with millions of smallholder farmers. Nature 555(7696), 363–366 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Zhang, X. et al. Managing nitrogen for sustainable development. Nature 528(7580), 51–59 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Alzaidi, A. A., Baig, M. B. & Elhag, E. A. An investigation into the farmers ’ attitudes towards organic farming in Riyadh Region–Kingdom of Saudi Arabia. Bulg. J. Agric. Sci. 19(3), 426–431 (2013).
Google Scholar 
Zhihui, W. et al. Combined applications of nitrogen and phosphorus fertilizers with manure increase maize yield and nutrient uptake via stimulating root growth in a long-term experiment. Pedosphere 26(1), 62–73 (2016).Article 
CAS 

Google Scholar 
Guang-hao, L., Gui-gen, C., Wei-ping, L. & Da-lei, L. Differences of yield and nitrogen use efficiency under different applications of slow release fertilizer in spring maize. J. Integr. Agric. 20(2), 554–564 (2020).
Google Scholar 
Zant, W. Is organic fertilizer going to be helpful in bringing a green revolution to sub-Saharan Africa? Economic explorations for Malawi agriculture (Working Paper). International House Hold Survey Network (2010).Barman, M., Paul, S., Choudhury, A. G., Roy, P. & Sen, J. Biofertilizer as prospective input for sustainable agriculture in India. Int. J. Curr. Microbiol. App. Sci. 6(11), 1177–1186 (2017).Article 

Google Scholar 
Kalhapure, A. H., Shete, B. T. & Dhonde, M. B. Integrated nutrient management in maize (Zea Mays L.) for increasing production with sustainability. Int. J. Agric. Food Sci. Technol. 4(3), 2249–3050 (2013).
Google Scholar 
Nazli, R. I., Kuşvuran, A., Inal, I., Demirbaş, A. & Tansi, V. Effects of different organic materials on forage yield and quality of silage maize (Zea mays L.). Turk. J. Agric. For. 38(1), 23–31 (2014).CAS 
Article 

Google Scholar 
Niu, Z. et al. Total factor productivity growth in china’s corn farming: an application of generalized productivity indicator. J. Bus. Econ. Manag. 22(5), 1189–1208 (2021).Article 

Google Scholar 
van Wesenbeeck, C. F. A., Keyzer, M. A., van Veen, W. C. M. & Qiu, H. Can China’s overuse of fertilizer be reduced without threatening food security and farm incomes?. Agric. Syst. 190, 103093 (2021).Article 

Google Scholar 
Ji, Y., Liu, H. & Shi, Y. Will China’s fertilizer use continue to decline? Evidence from LMDI analysis based on crops, regions and fertilizer types. PLoS ONE 15, e0237234 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Jiao, X. et al. Grain production versus resource and environmental costs: towards increasing sustainability of nutrient use in China. J. Exp. Bot. 67(17), 4935–4949 (2016).CAS 
PubMed 
Article 

Google Scholar 
Sher, A. et al. Response of maize grown under high plant density; performance, issues and management: a critical review. Adv. Crop Sci. Technol. 5(3), 1–8 (2017).Article 

Google Scholar 
De-yang, S. H. I. et al. Increased plant density and reduced N rate lead to more grain yield and higher resource utilization in summer maize. J. Integr. Agric. 15(11), 2515–2528 (2016).Article 

Google Scholar 
Du, X., Wang, Z., Lei, W. & Kong, L. Increased planting density combined with reduced nitrogen rate to achieve high yield in maize. Sci. Rep. 11(1), 1–12 (2021).CAS 
Article 

Google Scholar 
Li, T., Zhang, W., Yin, J., Chadwick, D., Norse, D., Lu, Y., Liu, X., Chen, X., Zhang, F., Powlson, D., & Dou, Z. Enhanced-efficiency fertilizers are not a panacea for resolving the nitrogen problem (2017).Adu-gyamfi, R. et al. One-time fertilizer briquettes application for maize production in savanna agroecologies of Ghana. Soil Fertil. Crop Prod. 111(6), 3339–3350 (2019).CAS 

Google Scholar 
Jiang, C. et al. Optimal nitrogen application rates of one-time root zone fertilization and the effect of reducing nitrogen application on summer maize. Sustainability 11, 2979 (2019).Article 

Google Scholar 
Jiang, C. et al. One-time root-zone N fertilization increases maize yield, NUE and reduces soil N losses in lime concretion black soil. Sci. Rep. 8(1), 1–10 (2018).ADS 

Google Scholar 
Li, G., Zhao, B., Dong, S., Liu, P. & Vyn, T. J. Impact of controlled release urea on maize yield and nitrogen use efficiency under different water conditions. PLoS ONE 12(7), 1–16 (2017).
Google Scholar 
Sikora, J. et al. Assessment of the efficiency of nitrogen slow-release fertilizers in integrated production of carrot depending on fertilization strategy. Sustainability (Switzerland) 12(5), 1–10 (2020).
Google Scholar 
Tian, C. et al. Effects of a controlled-release fertilizer on yield, nutrient uptake, and fertilizer usage efficiency in early ripening rapeseed (Brassica napus L.). J. Zhejian Univ. Sci. B (Biomed. Biotechnol.) 17(14), 775–786 (2016).CAS 
Article 

Google Scholar 
Tong, D. & Xu, R. Effects of urea and ( NH4)2SO4 on nitrification and acidification of Ultisols from Southern China. J. Environ. Sci. 24(4), 682–689 (2012).CAS 
Article 

Google Scholar 
El-rokiek, K. G., Ahmed, S. A. & Abd-elsamad, E. E. H. Effect of adding urea or ammonium sulphate on some herbicides efficiency in controlling weeds in onion plants. J. Am. Sci. 6(11), 536–543 (2010).
Google Scholar 
FAO. Guidelines for soil description. Enhanced Recovery After Surgery, (2006).Landon, J. Booker Tropical Soil manual: A Handbook for Soil Survey and Agriculture Land Evaluation in the Tropics and Subtropics (2013).Zhao, R. F. et al. Fertilization and nitrogen balance in a wheat-maize rotation system in North China. Agron. J. 98(4), 938–945 (2006).CAS 
Article 

Google Scholar 
Huang, S. et al. Estimation of nitrogen supply for summer maize production through a long-term field trial in china. Agronomy 11(7), 1358 (2021).CAS 
Article 

Google Scholar 
Dong, Y. J. et al. Effects of new coated release fertilizer on the growth of maize. J. Soil Sci. Plant Nutr. 16(3), 637–649 (2016).CAS 

Google Scholar 
Ngosong, C., Bongkisheri, V., Tanyi, C. B., Nanganoa, L. T. & Tening, A. S. Optimizing nitrogen fertilization regimes for sustainable maize (Zea mays L.) production on the volcanic soils of Buea Cameroon. Adv. Agric. 2019, 1–8 (2019).
Google Scholar 
Su, W., Ahmad, S., Ahmad, I. & Han, Q. Nitrogen fertilization affects maize grain yield through regulating nitrogen uptake, radiation and water use efficiency, photosynthesis and root distribution. PeerJ 8, 1–21 (2020).CAS 

Google Scholar 
Sainju, U. M, Ghimire, R., & Pradhan, G.P. Nitrogen Fertilization I: Impact on Crop, Soil, and Environment. IntechOpen https://doi.org/10.5772/intechopen.86028 (2020).Sha, Z. et al. Effect of N stabilizers on fertilizer-N fate in the soil-crop system: a meta- analysis. Agr. Ecosyst. Environ. 2020, 290 (2019).
Google Scholar 
Chen, K. & Vyn, T. J. Post-silking factor consequences for N efficiency changes over 38 years of commercial maize hybrids. Front. Plant Sci. https://doi.org/10.3389/fpls.2017.01737 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Jia, X. P. et al. Farmer’s adoption of improved nitrogen management strategies in maize production in China: an experimental knowledge training. J. Integr. Agric. 12(2), 364–373 (2013).Article 

Google Scholar 
Amanullah,. Rate and timing of nitrogen application influence partial factor productivity and agronomic NUE of maize (Zea mays L.) planted at low and high densities on calcareous soil in northwest Pakistan. J. Plant Nutr. 39(5), 683–690 (2016).CAS 
Article 

Google Scholar 
Draman, A., Almas, L. K. Partial factor productivity, agronomic efficiency, and economic analyses of maize in wheat-maize cropping system in Pakistan. Southern Agricultural Economics Association Annual Meetings, 2009 (January 2009).Yan, P. et al. Interaction between plant density and nitrogen management strategy in improving maize grain yield and nitrogen use efficiency on the North China Plain. Agric. Sci. 154, 978–988 (2016).Article 

Google Scholar 
Oenema, O. Nitrogen use efficiency (NUE) an indicator for the utilization of nitrogen in food systems. EU Nitrogen Expert Panel, January 2017, 1–4 (2015).Venterea, R. T., Coulter, J. A. & Dolan, M. S. Evaluation of intensive “4R” strategies for decreasing nitrous oxide emissions and nitrogen surplus in rainfed corn. J. Environ. Qual. 45(4), 1186–1195 (2016).CAS 
PubMed 
Article 

Google Scholar 
Zhang, C., Ju, X., Powlson, D., Oenema, O. & Smith, P. Nitrogen surplus benchmarks for controlling N pollution in the main cropping systems of China. Environ. Sci. Technol. 53(12), 6678–6687 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Fernández, C., Koop, G. & Steel, M. F. J. Multiple-output production with undesirable outputs multiple-output production with undesirable outputs : an application to nitrogen surplus in agriculture. J. Am. Stat. Assoc. 97(458), 432–442 (2013).MATH 
Article 

Google Scholar 
Børsting, C. F., Kristensen, T., Misciattelli, L., Hvelplund, T. & Weisbjerg, M. R. Reducing nitrogen surplus from dairy farms. Effects of feeding and management. Livest. Prod. Sci. 83(2–3), 165–178 (2003).Article 

Google Scholar 
Liang, K. et al. Reducing nitrogen surplus and environmental losses by optimized nitrogen and water management in double rice cropping system of South China. Agric. Ecosyst. Environ. 286, 106680 (2019).CAS 
Article 

Google Scholar 
Klages, S. et al. Nitrogen surplus-a unified indicator for water pollution in Europe?. Water (Switzerland) 12(4), 1197 (2020).CAS 

Google Scholar 
Muratoglu, A. Grey water footprint of agricultural production: an assessment based on nitrogen surplus and high-resolution leaching runoff fractions in Turkey. Sci. Total Environ. 742, 140553 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Niemiec, M. & Komorowska, M. The use of slow-release fertilizers as a part of optimization of celeriac production technology. Agric. Eng. 22(2), 59–68 (2018).
Google Scholar 
Ranum, P., Peña-Rosas, J. P. & Garcia-Casal, M. N. Global maize production, utilization, and consumption. Ann. N. Y. Acad. Sci. 1312(1), 105–112 (2014).ADS 
PubMed 
Article 

Google Scholar 
HLPE. Biofules and food security. High Level Panel of Experts on Food Security and Nutrition of the Committee on World Food Security, Rome (2013).Karp, A., Beale, M. H., Beaudoin, F. & Eastmond, P. J. Growing innovations for the bioeconomy. Nat. Plants https://doi.org/10.1038/nplants.2015.193 (2015).Article 
PubMed 

Google Scholar 
Chavarria, H., Trigo, E., Villarreal, F., Elverdin, P., & Piñeiro, V. Policy brief bioeconomy: a sustainable development strategy task force 10 sustainable energy, water, and food systems. T20, Saudi Arabia (2020). More

Here we have used existing surveillance to detect an emerging wildlife disease and appraise its impact by combining traditional host and vector screening with utilisation of national datasets generated by citizen scientists. Following the detection of USUV in the UK in 20207, whilst national surveillance identified no further cases of USUV infection in wild birds that year, we discovered a significant cluster of blackbird DIRs and an overlapping regional reduction in reported blackbird observations, possibly indicating disease-mediated population decline. Our investigation also identified mosquito vectors at the index site that were positive for USUV RNA, suggesting that ongoing virus transmission was likely.The most prevalent and notable histological changes in the blackbirds and house sparrow with confirmed USUV infection were those in the liver and spleen, consisting of necrosis and lymphohistiocytic inflammation along with moderate to abundant virus antigen labelling. Whilst neurotropism resulting in brain necrosis and lymphohistiocytic inflammation has been reported in studies which examined large numbers of wild blackbirds with USUV infection in continental Europe4,12, we found minimal evidence of neural lesions in the five wild birds examined in this study. Although, histopathological changes in other tissues were generally non-specific, immunolabelling demonstrated widespread virus antigen distribution in both bird species, which is similar to reports of USUV infection elsewhere4,13,14. Immunolabelling was disproportionately greater in the brain and heart in contrast to the minimal or absent histological changes observed in these organs: similar contrasting results of histological and immunohistochemical examinations of USUV-infected wild birds have previously been reported12. Although only brain and kidney samples were examined using USUV RT-PCR, our findings, together with earlier reports4,14, demonstrate that viral antigen can be detected in abundance in the heart and liver, suggesting that these organs could be useful for molecular diagnostic sampling.A differential for necrotising lesions in European passerines, and a comorbidity detected in blackbirds with USUV infection, is Plasmodium spp. infection4,8,15. DNA of the same Plasmodium spp. as detected in the tissues of USUV-positive blackbirds from the ZSL London Zoo site in 2020 was identified in Cx. pipiens s.l. that fed on blackbird hosts at this site previously in 2015, supporting endemic avian haemoparasite infection of this wild bird species at this location. In contrast to the results reported from USUV-positive blackbirds in the Netherlands4, no exo-erythrocytic stages of haemoprotozoa indicative of avian malaria were observed histologically in the two UK blackbirds positive for Plasmodium DNA. Since histological examination has limited sensitivity, in situ hybridisation could be used to further appraise the clinical significance of this co-infection in the future16.Zoological collections are ideally placed to form part of wildlife disease surveillance networks and have already contributed to flavivirus detection in mainland Europe10,13,17,18. The collection grounds at ZSL London Zoo are well monitored for evidence of morbidity or mortality in synanthropic wildlife; this unusually high level of vigilance is considered the likely explanation for detection of USUV at such a location. Recent import of infected captive birds can be excluded as a potential route of USUV introduction as the COVID-19 pandemic had led to suspension of animal movements into the zoological collection. Following USUV detection in synanthropic wildlife, preemptive management practices were employed to safeguard the health of captive animals (Supplementary Materials 1); there was no evidence of USUV-associated disease in the collection animals.The majority of mosquitoes trapped in 2020 were primarily ornithophagic Cx. pipiens s.l., a known vector for USUV1 and a common species in temperate urban habitats. This mosquito species was also the most frequently detected at the ZSL London Zoo site in 2015, during historical trapping sessions19 and at two zoological collections in northern England20. Bloodmeal analyses from mosquitoes at ZSL London Zoo in 2015 and 2020 demonstrate that this species feeds on both wild and collection birds, as would be expected for a generalist ornithophagic mosquito21. In addition, targeted mosquito surveillance in 2020 confirmed circulating USUV in multiple Cx. pipiens s.l. pools at the index site over a three-week period subsequent to the detection of USUV-associated wild bird mortality. This further demonstrates that local mosquito trapping combined with PCR screening is useful as part of an integrated surveillance programme22 and provides evidence that native vectors in the UK may facilitate the onward transmission of USUV to susceptible hosts following an emergence event.Wild bird flavivirus surveillance in Great Britain integrates submissions from three schemes, each with a different taxonomic focus. These convenience samples inevitably lead to skews in species coverage. Although a common garden bird, the number of blackbirds tested for USUV was modest at 2–8 per annum over the period 2012–2019. A communication programme to raise awareness of blackbirds as a sentinel species for USUV, involving a range of stakeholder communities (e.g. non-governmental organisations, wildlife rehabilitators and veterinary surgeons) could help to increase the volume of submissions and, by extension, the ability to rapidly identify the occurrence of USUV in this species. The potential value of target species as sentinels within wild bird surveillance networks has been highlighted for other pathogens, e.g. highly pathogenic H5N1 avian influenza and West Nile virus23,24. In addition to this passive surveillance focused on disease detection in avian hosts, active targeted serosurveys could be conducted to identify cryptic exposure of subclinically affected birds in the future. Given the logistical challenges around active serosurveys in wild birds, screening of archived samples from captive birds in the zoological collection may provide a means to further appraise the extent of USUV circulation, as has previously been undertaken at other collections in mainland Europe25,26.Local reductions of blackbird populations have been reported following USUV outbreaks in mainland Europe27,28,29, but numbers recorded by the BBS have been stable in the UK and Greater London since 2011 when USUV incursion would be predicted most likely to have occurred on the basis of spatio-temporal patterns of spread in mainland Europe3 until the latest data are available from 2019 (Supplementary Figure 5). Whilst our index site detection of USUV is unlikely to represent the incursion event, and earlier sporadic or localised USUV incidents prior to 2020 may have occurred7, based on historical blackbird population trends it seems plausible that the existing surveillance system enabled rapid detection of this emerging infectious disease.Significant clustering of blackbird DIRs was observed in the Greater London, South East and East of England regions in 2020. These results should be interpreted with care given the potential for biases with these opportunistic data and the absence of confirmed aetiology for the DIRs, however, these findings are consistent with a regional increase in blackbird morbidity and mortality in summer 2020 around the USUV index site. Consequently, it is likely that further blackbirds, in addition to those recovered for PME, were infected and died with USUV. Whilst no evidence of an increase in generalised ill health or neurological disease category blackbird DIRs was found in 2020, particular attention should be paid to early detection of clusters of DIRs of these categories as a potential signal of USUV occurrence in the future.One indicator, the dead bird ringed recovery dataset, did not support increased scale of blackbird mortality in Greater London; however, the dataset is small and vulnerable to variation in observer bias (e.g. related to COVID-19 induced lockdown and travel restrictions). In contrast, using the GBW dataset, we identified a substantial seasonal decline in the blackbird weekly reporting rate which was associated with a concomitant reduction in weekly count in gardens, but not in ecologically similar control species, which was contemporaneous with the period of detected USUV activity in Greater London. These population trends are consistent with a hypothesis of disease-mediated decline. Alternative explanations, such as variation in climate, food availability or bird movement need consideration and are discussed next.Exploration of climate data indicates that, whilst the spring and early summer of 2020 was noteworthy with a high daily temperature average and low rainfall, at the time of USUV detection and the decline in the blackbird reporting rate, these parameters were within historical ranges (Supplementary Table 7). Consequently, while the climate may have been permissive for USUV transmission, there is no evidence to support variation in the weather alone as an explanation for the seasonal pattern of blackbird reporting rate decline; nor were declines observed in the robin or starling data, the control species with similar soil invertebrate diet and therefore similar vulnerability to summer drought. Blackbird, robin and starling populations in the UK are partially migratory; however, birds from mainland Europe do not migrate to overwinter in England until mid-October (i.e. after the decline in blackbird reporting rate occurred): consequently international bird movement does not offer an explanation for the observed regional blackbird decline. During the late summer season, short-distance movement from garden to non-garden habitats typically occurs, during the period of moult; however, the extent of the decline in blackbird reporting rate in gardens that occurred in Greater London in 2020 markedly exceeds that of the historical trend (2011–2019 inclusive; BTO unpubl. data). In summary, despite the fact that surveillance did not confirm further cases of wild bird USUV infection in 2020, and whilst it is not possible to ascribe causality, or exclude the chance that other factors may have contributed to the observed population trend, it remains possible that large-scale blackbird mortality due to USUV occurred in Greater London in summer 2020.Our study and others30 illustrate the need to integrate disease surveillance and long-term population monitoring schemes to evaluate disease impact, and to use control species to explore potential confounding drivers of population change (e.g. climate, food availability). Since GBW reporting rates are generated online in real-time, and nationwide, they offer a tool to rapidly detect changes in species presence (i.e. reporting rate) or flock size in gardens (i.e. weekly maximum count) that can be used to strategically enhance surveillance effort for disease detection. As wild bird ring recovery reports are also submitted online, there is also the potential to develop a complementary system that monitors for trends in occurrence of dead birds that might signal a disease outbreak. The BBS survey provides the most robust available data on population trends to appraise disease impact, however there is a delay of some months until data from this scheme become available. Since repeated incursions have occurred in mainland Europe following first detection1,17,31, it is likely that USUV will emerge in the UK again, either through overwintering or repeat incursion(s). Integrated disease surveillance in combination with bird population monitoring using the various available datasets, as we have capitalised on here, is required to assess whether USUV re-occurs, or becomes endemic, in UK wild birds and to identify any associated population impacts.By combining a range of professional and citizen science datasets our study approach facilitates the rapid detection of an emerging disease in free-living wildlife and enables insights into its incipient impact. We believe this multidisciplinary approach presents a framework for the early detection of disease outbreaks and incursion, thus helping to safeguard animal and public health. Such early warning systems could facilitate prompt mitigation action, for example targeted biosecurity measures and enhanced vigilance by medical and veterinary authorities. In addition, there is opportunity to further develop collaboration with ornithologists through active surveillance of wild birds, as was recently employed to detect West Nile virus in a migratory bird in the Netherlands32. Whilst population monitoring schemes are most developed for wild birds, lessons learned may be applied for the surveillance of diseases affecting other taxa. More

Blosser, E. M. et al. Environmental drivers of seasonal patterns of host utilization by Culiseta melanura (Diptera: Culicidae) in Florida. J. Med. Entomol. 54, 1365–1374 (2017).PubMed 
PubMed Central 

Google Scholar 
Burkett-Cadena, N. D., Hassan, H. K., Eubanks, M. D., Cupp, E. W. & Unnasch, T. R. Winter severity predicts the timing of host shifts in the mosquito Culex erraticus. Biol. Lett. 8, 567–569 (2012).PubMed 
PubMed Central 

Google Scholar 
Gürtler, R. E., Cecere, M. C., Vazquez, D. P., Chuit, R. & Cohen, J. E. Host-feeding patterns of domiciliary Triatoma infestans (Hemiptera: Reduviidae) in northwest Argentina: Seasonal and instar variation. J. Med. Entomol. 33, 15–26 (1996).PubMed 

Google Scholar 
Rabinovich, J. E. et al. Ecological patterns of blood-feeding by kissing-bugs (Hemiptera: Reduviidae: Triatominae). Mem. Inst. Oswaldo Cruz 106, 479–494 (2011).PubMed 

Google Scholar 
Kent, R. J. Molecular methods for arthropod bloodmeal identification and applications to ecological and vector-borne disease studies. Mol. Ecol. Resour. 9, 4–18 (2009).CAS 
PubMed 

Google Scholar 
Cecere, M. C. et al. Host-feeding sources and infection with Trypanosoma cruzi of Triatoma infestans and Triatoma eratyrusiformis (Hemiptera: Reduviidae) from the Calchaqui Valleys in Northwestern Argentina. J. Med. Entomol. 53, 666–673 (2016).CAS 
PubMed 

Google Scholar 
Logue, K. et al. Unbiased characterization of Anopheles mosquito blood meals by targeted high-throughput sequencing. PLoS Negl. Trop. Dis. 10, e0004512 (2016).PubMed 
PubMed Central 

Google Scholar 
Keller, J. I., Schmidt, J. O., Schmoker, A. M., Ballif, B. A. & Stevens, L. Protein mass spectrometry extends temporal blood meal detection over polymerase chain reaction in mouse-fed Chagas disease vectors. Mem. Inst. Oswaldo Cruz 113, e180160 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Borland, E. M. & Kading, R. C. Modernizing the toolkit for arthropod bloodmeal identification. Insects 12, 1–27 (2021).
Google Scholar 
Clark, K., Karsch-Mizrachi, I., Lipman, D. J., Ostell, J. & Sayers, E. W. GenBank. Nucleic Acids Res. 44, D67–D72 (2016).CAS 
PubMed 

Google Scholar 
Hamer, S. A. et al. Comparison of DNA and carbon and nitrogen stable isotope-based techniques for identification of prior vertebrate hosts of ticks. J. Med. Entomol. 52, 1043–1049 (2015).CAS 
PubMed 

Google Scholar 
Scott, M. C., Harmon, J. R., Tsao, J. I., Jones, C. J. & Hickling, G. J. Reverse line blot probe design and polymerase chain reaction optimization for bloodmeal analysis of ticks from the eastern United States. J. Med. Entomol. 49, 697–709 (2012).CAS 
PubMed 

Google Scholar 
Arias-Giraldo, L. M. et al. Identification of blood-feeding sources in Panstrongylus, Psammolestes, Rhodnius and Triatoma using amplicon-based next-generation sequencing. Parasit. Vectors 13, 1–14 (2020).
Google Scholar 
Kieran, T. J. et al. Blood meal source characterization using Illumina sequencing in the Chagas Disease vector Rhodnius pallescens (Hemiptera: Reduviidae) in Panamá. J. Med. Entomol. https://doi.org/10.1093/jme/tjx170 (2017).Article 
PubMed 

Google Scholar 
Dumonteil, E. et al. Detailed ecological associations of triatomines revealed by metabarcoding and next-generation sequencing: Implications for triatomine behavior and Trypanosoma cruzi transmission cycles. Sci. Rep. 8, 1–13 (2018).CAS 

Google Scholar 
Estrada-Franco, J. G. et al. Vertebrate-Aedes aegypti and Culex quinquefasciatus (Diptera)-arbovirus transmission networks: Non-human feeding revealed by meta-barcoding and nextgeneration sequencing. PLoS Negl. Trop. Dis. 14, 1–22 (2020).
Google Scholar 
Campana, M. G. et al. Simultaneous identification of host, ectoparasite and pathogen DNA via in-solution capture. Mol. Ecol. Resour. 16, 1224–1239 (2016).CAS 
PubMed 

Google Scholar 
Klotz, S. A. et al. Free-roaming kissing bugs, vectors of Chagas disease, feed often on humans in the Southwest. Am. J. Med. 127, 421–426 (2014).PubMed 
PubMed Central 

Google Scholar 
Waleckx, E., Suarez, J., Richards, B. & Dorn, P. L. Triatoma sanguisuga blood meals and potential for Chagas Disease, Louisiana, USA. Emerg. Infect. Dis. 20, 2141–2143 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
Kjos, S. A. et al. Identification of bloodmeal sources and Trypanosoma cruzi infection in triatomine bugs (Hemiptera: Reduviidae) from residential settings in Texas, the United States. J. Med. Entomol. 50, 1126–1139 (2013).PubMed 

Google Scholar 
Gürtler, R. E., Cohen, J. E., Cecere, M. C. & Chuit, R. Shifting host choices of the vector of Chagas Disease, Triatoma Infestans, in relation to the availability of host in houses in North-West Argentina. J. Appl. Ecol. 34, 699–715 (1997).
Google Scholar 
Minuzzi-Souza, T. et al. Molecular bloodmeal analyses reveal that Trypanosoma cruzi-infected, native triatomine bugs often feed on humans in houses in central Brazil. Med. Vet. Entomol. 32, 504–508 (2018).CAS 
PubMed 

Google Scholar 
Lent, H. & Wygodzinsky, P. W. Revision of the Triatominae (Hemiptera, Reduviidae), and their significance as vectors of Chagas’ Disease. Bull. Am. Museum Nat. Hist. 163, 123–520 (1979).
Google Scholar 
World Health Organization. Chagas disease in Latin America: An epidemiological update based on 2010 estimates. Wkly. Epidemiol. Rec. 6, 33–44 (2015).
Google Scholar 
Dorn, P. L. et al. Autochthonous transmission of Trypanosoma cruzi, Louisiana. Emerg. Infect. Dis. 13, 13–15 (2007).
Google Scholar 
Cantey, P. T. et al. The United States Trypanosoma cruzi infection study: Evidence for vector-borne transmission of the parasite that causes Chagas disease among United States blood donors. Transfusion 52, 1922–1930 (2012).PubMed 

Google Scholar 
Garcia, M. N. et al. Molecular identification and genotyping of Trypanosoma cruzi DNA in autochthonous Chagas disease patients from Texas, USA. Infect. Genet. Evol. 49, 151–156 (2017).CAS 
PubMed 

Google Scholar 
Barr, S., Gossett, K. A. & Klei, T. R. Clinical, clinicopathologic, and parasitologic observations of trypanosomiasis in dogs infected with North American Trypanosoma cruzi isolates. Am. J. Vet. Res. 52, 954–960 (1991).CAS 
PubMed 

Google Scholar 
Meyers, A. C., Meinders, M. & Hamer, S. A. Widespread Trypanosoma cruzi infection in government working dogs along the Texas-Mexico border: Discordant serology, parasite genotyping and associated vectors. PLoS Negl. Trop. Dis. 11, 1–19 (2017).
Google Scholar 
Meyers, A. C., Edwards, E. E., Sanders, J. P., Saunders, A. B. & Hamer, S. A. Fatal Chagas myocarditis in government working dogs in the southern United States: Cross-reactivity and differential diagnoses in five cases across six months. Vet. Parasitol. Reg. Stud. Rep. 24, 1–7 (2021).
Google Scholar 
Hodo, C. L. & Hamer, S. A. Toward an ecological framework for assessing reservoirs of vector-borne pathogens: Wildlife reservoirs of Trypanosoma cruzi across the southern United States. ILAR J. 58, 379–392 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
Guarneri, A. A., Pereira, M. H. & Diotaiuti, L. Influence of the blood meal source on the development of Triatoma infestans, Triatoma brasiliensis, Triatoma sordida, and Triatoma pseudomaculata (Heteroptera, Reduviidae). J. Med. Entomol. 37, 373–379 (2000).CAS 
PubMed 

Google Scholar 
Pippin, W. F. The biology and vector capability of Triatoma sanguisuga texana Usinger and Triatoma gerstaeckeri (Stål) compared with Rhodnius prolixus (Stål) (Hemiptera: Triatominae). J. Med. Entomol. 7, 30–45 (1970).CAS 
PubMed 

Google Scholar 
Bern, C., Kjos, S., Yabsley, M. J. & Montgomery, S. P. Trypanosoma cruzi and Chagas’ disease in the United States. Clin. Microbiol. Rev. 24, 655–681 (2011).PubMed 
PubMed Central 

Google Scholar 
Kjos, S. et al. Distribution and characterization of canine Chagas disease in Texas. Vet. Parasitol. 152, 249–256 (2008).CAS 
PubMed 

Google Scholar 
Tenney, T. D., Curtis-Robles, R., Snowden, K. F. & Hamer, S. A. Shelter dogs as sentinels for Trypanosoma cruzi transmission across Texas. Emerg. Infect. Dis. 20, 1323–1326 (2014).PubMed 
PubMed Central 

Google Scholar 
Curtis-Robles, R., Wozniak, E. J., Auckland, L. D., Hamer, G. L. & Hamer, S. A. Combining public health education and disease ecology research: Using citizen science to assess Chagas disease entomological risk in Texas. PLoS Negl. Trop. Dis. 9, e0004235 (2015).PubMed 
PubMed Central 

Google Scholar 
Curtis-Robles, R., Hamer, S. A., Lane, S., Levy, M. Z. & Hamer, G. L. Bionomics and spatial distribution of triatomine vectors of Trypanosoma cruzi in Texas and other southern states, USA. Am. J. Trop. Med. Hyg. 98, 113–121 (2018).PubMed 

Google Scholar 
Curtis-Robles, R., Aukland, L. D., Snowden, K. F., Hamer, G. L. & Hamer, S. A. Analysis of over 1500 triatomine vectors from across the US, predominantly Texas, for Trypanosoma cruzi infection and discrete typing units. Infect. Genet. Evol. 58, 171–180 (2018).PubMed 

Google Scholar 
Hodo, C. L., Wilkerson, G. K., Birkner, E. C., Gray, S. B. & Hamer, S. A. Trypanosoma cruzi transmission among captive nonhuman primates, wildlife, and vectors. EcoHealth 15, 426–436 (2018).PubMed 
PubMed Central 

Google Scholar 
Duffy, T. et al. Analytical performance of a multiplex Real-Time PCR assay using TaqMan probes for quantification of Trypanosoma cruzi satellite DNA in blood samples. PLoS Negl. Trop. Dis. 7, e2000 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
Piron, M. et al. Development of a real-time PCR assay for Trypanosoma cruzi detection in blood samples. Acta Trop. 103, 195–200 (2007).CAS 
PubMed 

Google Scholar 
Curtis-Robles, R. et al. Parasitic interactions among Trypanosoma cruzi, triatomine vectors, domestic animals, and wildlife in Big Bend National Park along the Texas-Mexico border. Acta Trop. 188, 225–233 (2018).PubMed 

Google Scholar 
Cupp, E. W. et al. Identification of reptilian and amphibian blood meals from mosquitoes in an eastern equine encephalomyelitis virus focus in central Alabama. Am. J. Trop. Med. Hyg. 71, 272–276 (2004).PubMed 

Google Scholar 
Medeiros, M. C. I., Ricklefs, R. E., Brawn, J. D. & Hamer, G. L. Plasmodium prevalence across avian host species is positively associated with exposure to mosquito vectors. Parasitology 142, 1612–1620 (2015).CAS 
PubMed 

Google Scholar 
Hamer, G. L. et al. Host selection by Culex pipiens mosquitoes and West Nile virus amplification. Am. J. Trop. Med. Hyg. 80, 268–278 (2009).PubMed 

Google Scholar 
Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).CAS 
PubMed 

Google Scholar 
Hathaway, N. J., Parobek, C. M., Juliano, J. J. & Bailey, J. A. SeekDeep: Single-base resolution de novo clustering for amplicon deep sequencing. Nucleic Acids Res. 46, e21 (2018).CAS 
PubMed 

Google Scholar 
Zeledón, R. et al. An Appraisal of the Status of Chagas Disease in the United States (Elsevier Inc., Amsterdam, 2012).
Google Scholar 
Gorchakov, R. et al. Trypanosoma cruzi infection prevalence and bloodmeal analysis in triatomine vectors of Chagas disease from rural peridomestic locations in Texas, 2013–2014. J. Med. Entomol. 53, 911–918 (2016).CAS 
PubMed 

Google Scholar 
Stevens, L. et al. Vector blood meals and Chagas Disease transmission potential, United States. Emerg. Infect. Dis. 18, 646–650 (2012).PubMed 
PubMed Central 

Google Scholar 
Polonio, R., López-Domínguez, J., Herrera, C. & Dumonteil, E. Molecular ecology of Triatoma dimidiata in southern Belize reveals risk for human infection and the local differentiation of Trypanosoma cruzi parasites. Int. J. Infect. Dis. 108, 320–329 (2021).CAS 
PubMed 

Google Scholar 
Sasaki, H., Rosales, R. & Tabaru, Y. Host feeding profiles of Rhodnius prolixus and Triatoma dimidiata in Guatemala (Hemiptera: Reduviidae: Triatominae). Med. Entomol. Zool. 54, 283–289 (2003).
Google Scholar 
Villalobos, G., Martínez-Hernández, F., de la Torre, P., Laclette, J. P. & Espinoza, B. Entomological indices, feeding sources, and molecular identification of Triatoma phyllosoma (Hemiptera: Reduviidae) one of the main vectors of Chagas disease in the Istmo de Tehuantepec, Oaxaca, Mexico. Am. J. Trop. Med. Hyg. 85, 490–497 (2011).PubMed 
PubMed Central 

Google Scholar 
Mota, J. et al. Identification of blood meal source and infection with Trypanosoma cruzi of Chagas disease vectors using a multiplex cytochrome b polymerase chain reaction assay. Vector Borne Zoonotic Dis. 7, 617–627 (2007).PubMed 

Google Scholar 
Pizarro, J. C. & Stevens, L. A new method for forensic DNA analysis of the blood meal in Chagas disease vectors demonstrated using Triatoma infestans from Chuquisaca, Bolivia. PLoS ONE 3, e3585 (2008).ADS 
PubMed 
PubMed Central 

Google Scholar 
Abad-Franch, F. & Gurgel-Gonçalves, R. The ecology and natural history of wild triatominae in the Americas. In Triatominae—The Biology of Chagas Disease Vectors (eds Guarneri, A. & Lorenzo, M.) 387–445 (Springer Nature Switzerland AG, 2021). https://doi.org/10.1007/978-3-030-64548-9_16.Chapter 

Google Scholar 
Busselman, R. E. & Hamer, S. A. Chagas disease ecology in the United States: Recent advances in understanding Trypanosoma cruzi transmission among triatomines, wildlife, and domestic animals and a quantitative synthesis of vector-host interactions. Annu. Rev. Anim. Biosci. 10, 325–348 (2022).PubMed 

Google Scholar 
Minuzzi-Souza, T. T. C. et al. Vector-borne transmission of Trypanosoma cruzi among captive Neotropical primates in a Brazilian zoo. Parasit. Vectors 9, 1–7. https://doi.org/10.1186/s13071-016-1334-7 (2016).CAS 
Article 

Google Scholar 
Reis, F. C. et al. Trypanosomatid infections in captive wild mammals and potential vectors at the Brasilia Zoo, Federal District, Brazil. Vet. Med. Sci. 6, 248–256. https://doi.org/10.1002/vms3.216 (2020).CAS 
Article 
PubMed 

Google Scholar 
Martínez-Hernández, F., Oria-Martínez, B., Rendón-Franco, E., Villalobos, G. & Muñoz-García, C. I. Trypanosoma cruzi, beyond the dogma of non-infection in birds. Infect. Genet. Evol. https://doi.org/10.1016/j.meegid.2022.105239 (2022).Article 
PubMed 

Google Scholar 
Botto-Mahan, C. et al. Lizards as silent hosts of Trypanosoma cruzi. Emerg. Infect. Dis. 28, 1250–1253. https://doi.org/10.3201/eid2806.220079 (2022).Article 
PubMed 
PubMed Central 

Google Scholar  More

Allen, R. G. et al. Crop evapotranspiration-guidelines for computing crop water requirements-FAO irrigation and drainage paper 56. FAO, Rome 300, D05109 http://www.fao.org/docrep/X0490E/X0490E00.htm (1998).
Google Scholar 
Trenberth, K. E., Fasullo, J. T. & Kiehl, J. Earth’s global energy budget. Bulletin of the American Meteorological Society 90, 311–324, https://doi.org/10.1175/2008BAMS2634.1 (2009).ADS 
Article 

Google Scholar 
Wang, K. & Dickinson, R. E. A review of global terrestrial evapotranspiration: Observation, modeling, climatology, and climatic variability. Reviews of Geophysics 50, https://doi.org/10.1029/2011RG000373 (2012).Liu, W. Evaluating remotely sensed monthly evapotranspiration against water balance estimates at basin scale in the Tibetan Plateau. Hydrology Research 49, 1977–1990, https://doi.org/10.2166/nh.2018.008 (2018).Article 

Google Scholar 
Valipour, M., Bateni, S. M., Gholami Sefidkouhi, M. A., Raeini-Sarjaz, M. & Singh, V. P. Complexity of forces driving trend of reference evapotranspiration and signals of climate change. Atmosphere 11, 1081, https://doi.org/10.1007/s41748-021-00252-3 (2020).ADS 
Article 

Google Scholar 
Anwar, S. A., Mamadou, O., Diallo, I. & Sylla, M. B. On the influence of vegetation cover changes and vegetation-runoff systems on the simulated summer potential evapotranspiration of Tropical Africa using RegCM4. Earth Systems and Environment 5, 883–897, https://doi.org/10.3390/atmos11101081 (2021).ADS 
Article 

Google Scholar 
Tomas-Burguera, M., Vicente-Serrano, S. M., Beguera, S., Reig, F. & Latorre, B. Reference crop evapotranspiration database in Spain (1961–2014). Earth System Science Data 11, 1917–1930, https://doi.org/10.5194/essd-11-1917-2019 (2019).ADS 
Article 

Google Scholar 
Córdova, M., Carrillo-Rojas, G., Crespo, P., Wilcox, B. & Célleri, R. Evaluation of the Penman-Monteith (FAO 56 PM) Method for Calculating Reference Evapotranspiration Using Limited Data. Mountain Research and Development 35, 230–239, https://doi.org/10.1659/MRD-JOURNAL-D-14-0024.1 (2015).Article 

Google Scholar 
Valle Júnior, L. C. G. d. et al. Evaluation of FAO-56 procedures for estimating reference evapotranspiration using missing climatic data for a Brazilian tropical savanna. Water 13, 1763, https://doi.org/10.3390/w13131763 (2021).Article 

Google Scholar 
Djaman, K., Irmak, S. & Futakuchi, K. Daily reference evapotranspiration estimation under limited data in Eastern Africa. Journal of Irrigation and Drainage Engineering 143, 06016015, https://doi.org/10.1061/(ASCE)IR.1943-4774.0001154 (2017).Article 

Google Scholar 
Čadro, S., Uzunović, M., Žurovec, J. & Žurovec, O. Validation and calibration of various reference evapotranspiration alternative methods under the climate conditions of Bosnia and Herzegovina. International Soil and Water Conservation Research 5, 309–324, https://doi.org/10.1016/j.iswcr.2017.07.002 (2017).Article 

Google Scholar 
Vicente-Serrano, S. M. et al. Recent changes in monthly surface air temperature over Peru, 1964–2014. International Journal of Climatology 38, 283–306, https://doi.org/10.1002/joc.5176 (2018).ADS 
Article 

Google Scholar 
Huerta, A., Aybar, C. & Lavado-Casimiro, W. PISCO temperatura versión 1.1 (PISCOt v1. 1). Lima, Peru: National Meteorology and Hydrology Service of Peru (SENAMHI) https://iridl.ldeo.columbia.edu/SOURCES/.SENAMHI/.HSR/.PISCO/.Temp/ (2018).Lavado Casimiro, W. S., Labat, D., Guyot, J. L. & Ardoin-Bardin, S. Assessment of climate change impacts on the hydrology of the Peruvian Amazon–Andes basin. Hydrological Processes 25, 3721–3734, https://doi.org/10.1002/hyp.8097 (2011).ADS 
Article 

Google Scholar 
Rau, P. et al. Assessing multidecadal runoff (1970–2010) using regional hydrological modelling under data and water scarcity conditions in Peruvian Pacific catchments. Hydrological Processes 33, 20–35, https://doi.org/10.1002/hyp.13318 (2019).ADS 
Article 

Google Scholar 
Olsson, T. et al. Downscaling climate projections for the Peruvian coastal Chancay-Huaral basin to support river discharge modeling with WEAP. Journal of Hydrology: Regional Studies 13, 26–42, https://doi.org/10.1016/j.ejrh.2017.05.011 (2017).Article 

Google Scholar 
Lavado-Casimiro, W., Lhomme, J. P., Labat, D. & Loup, J. Estimating reference evapotranspiration (FAO 56 Penman Monteith) with limited climatic data in the Peruvian Amazon-Andes basin. Revista Peruana Geo-Atmosferica 4, 31–43 (2015).
Google Scholar 
Hargreaves, G. H. & Samani, Z. A. Reference crop evapotranspiration from temperature. Applied engineering in agriculture 1, 96–99, https://doi.org/10.13031/2013.26773 (1985).Article 

Google Scholar 
Laqui, W. et al. Can artificial neural networks estimate potential evapotranspiration in Peruvian highlands? Modeling Earth Systems and Environment 5, 1911–1924, https://doi.org/10.1007/s40808-019-00647-2 (2019).Article 

Google Scholar 
Baigorria, G. A., Villegas, E. B., Trebejo, I., Carlos, J. F. & Quiroz, R. Atmospheric transmissivity: distribution and empirical estimation around the central Andes. International Journal of Climatology 24, 1121–1136, https://doi.org/10.1002/joc.1060 (2004).ADS 
Article 

Google Scholar 
Huerta, A. PISCO potential evapotranspiration, https://iridl.ldeo.columbia.edu/SOURCES/.SENAMHI/.HSR/.PISCO/.PET/.Oudin, L., Michel, C. & Anctil, F. Which potential evapotranspiration input for a lumped rainfall-runoff model?: Part 1—can rainfall-runoff models effectively handle detailed potential evapotranspiration inputs? Journal of Hydrology 303, 275–289, https://doi.org/10.1016/j.jhydrol.2004.08.025 (2005).ADS 
Article 

Google Scholar 
Xiang, K., Li, Y., Horton, R. & Feng, H. Similarity and difference of potential evapotranspiration and reference crop evapotranspiration – a review. Agricultural Water Management 232, 106043, https://doi.org/10.1016/j.agwat.2020.106043 (2020).Article 

Google Scholar 
Harris, I., Osborn, T. J., Jones, P. & Lister, D. Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. Scientific data 7, 1–18, https://doi.org/10.1038/s41597-020-0453-3 (2020).Article 

Google Scholar 
Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Scientific data 5, 1–12, https://doi.org/10.1038/sdata.2017.191 (2018).Article 

Google Scholar 
Kalnay, E. et al. The NCEP/NCAR 40-year reanalysis project. Bulletin of the American Meteorological Society 77, 437–472, 10.1175/1520-0477(1996)077  2.0.CO;2 (1996).Hersbach, H. et al. The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society 146, 1999–2049, https://doi.org/10.1002/qj.3803 (2020).ADS 
Article 

Google Scholar 
Muñoz Sabater, J. et al. ERA5-Land: a state-of-the-art global reanalysis dataset for land applications. Earth System Science Data 13, 4349–4383, https://doi.org/10.5194/essd-13-4349-2021 (2021).ADS 
Article 

Google Scholar 
Aybar, C. et al. Construction of a high-resolution gridded rainfall dataset for Peru from 1981 to the present day. Hydrological Sciences Journal 65, 770–785, https://doi.org/10.1080/02626667.2019.1649411 (2020).Article 

Google Scholar 
Llauca, H., Lavado-Casimiro, W., Montesinos, C., Santini, W. & Rau, P. PISCO_HyM_GR2M: A model of monthly water balance in Peru (1981–2020). Water 13, 1048, https://doi.org/10.3390/w13081048 (2021).Article 

Google Scholar 
Gubler, S. et al. The influence of station density on climate data homogenization. International Journal of Climatology 37, 4670–4683, https://doi.org/10.1002/joc.5114 (2017).ADS 
Article 

Google Scholar 
Hunziker, S. et al. Identifying, attributing, and overcoming common data quality issues of manned station observations. International Journal of Climatology 37, 4131–4145, https://doi.org/10.1002/joc.5037 (2017).ADS 
Article 

Google Scholar 
Hunziker, S. et al. Effects of undetected data quality issues on climatological analyses. Climate of the Past 14, 1–20, https://doi.org/10.5194/cp-14-1-2018 (2018).ADS 
Article 

Google Scholar 
Paredes, P., Pereira, L., Almorox, J. & Darouich, H. Reference grass evapotranspiration with reduced data sets: Parameterization of the FAO Penman-Monteith temperature approach and the Hargeaves-Samani equation using local climatic variables. Agricultural Water Management 240, 106210, https://doi.org/10.1016/j.agwat.2020.106210 (2020).Article 

Google Scholar 
Irmak, S., Kabenge, I., Skaggs, K. E. & Mutiibwa, D. Trend and magnitude of changes in climate variables and reference evapotranspiration over 116-yr period in the Platte River Basin, central Nebraska–USA. Journal of Hydrology 420, 228–244, https://doi.org/10.1016/j.jhydrol.2011.12.006 (2012).ADS 
Article 

Google Scholar 
Tomas-Burguera, M., Vicente-Serrano, S. M., Grimalt, M. & Beguera, S. Accuracy of reference evapotranspiration (ETo) estimates under data scarcity scenarios in the Iberian Peninsula. Agricultural water management 182, 103–116, https://doi.org/10.1016/j.agwat.2016.12.013 (2017).Article 

Google Scholar 
Mardikis, M., Kalivas, D. & Kollias, V. Comparison of interpolation methods for the prediction of reference evapotranspiration—an application in Greece. Water Resources Management: An International Journal, Published for the European Water Resources Association (EWRA) 19, 251–278, https://doi.org/10.1007/s11269-005-3179-2 (2005).Article 

Google Scholar 
McVicar, T. R. et al. Spatially distributing monthly reference evapotranspiration and pan evaporation considering topographic influences. Journal of Hydrology 338, 196–220, https://doi.org/10.1016/j.jhydrol.2007.02.018 (2007).ADS 
Article 

Google Scholar 
Robinson, E. L., Blyth, E. M., Clark, D. B., Finch, J. & Rudd, A. C. Trends in atmospheric evaporative demand in Great Britain using high-resolution meteorological data. Hydrology and Earth System Sciences 21, 1189–1224, https://doi.org/10.5194/hess-21-1189-2017 (2017).ADS 
Article 

Google Scholar 
Mohammadi, B. & Moazenzadeh, R. Performance analysis of daily global solar radiation models in Peru by regression analysis. Atmosphere 12, https://doi.org/10.3390/atmos12030389 (2021).Vicente-Serrano, S. M., Beguería, S., López-Moreno, J. I., García-Vera, M. A. & Stepanek, P. A complete daily precipitation database for northeast Spain: reconstruction, quality control, and homogeneity. International Journal of Climatology 30, 1146–1163, https://doi.org/10.1002/joc.1850 (2010).ADS 
Article 

Google Scholar 
Lanzante, J. R. Resistant, robust and non-parametric techniques for the analysis of climate data: Theory and examples, including applications to historical radiosonde station data. International Journal of Climatology: A Journal of the Royal Meteorological Society 16, 1197–1226, 10.1002/(SICI)1097-0088(199611)16:11 < 1197::AID-JOC89 > 3.0.CO;2-L (1996).ADS 
Article 

Google Scholar 
Wood, W. H., Marshall, S. J., Whitehead, T. L. & Fargey, S. E. Daily temperature records from a mesonet in the foothills of the Canadian Rocky Mountains, 2005–2010. Earth System Science Data 10, 595–607, https://doi.org/10.5194/essd-10-595-2018 (2018).ADS 
Article 

Google Scholar 
Guentchev, G., Barsugli, J. J. & Eischeid, J. Homogeneity of gridded precipitation datasets for the Colorado River basin. Journal of Applied Meteorology and Climatology 49, 2404–2415, https://doi.org/10.1175/2010JAMC2484.1 (2010).ADS 
Article 

Google Scholar 
Oyler, J. W., Ballantyne, A., Jencso, K., Sweet, M. & Running, S. W. Creating a topoclimatic daily air temperature dataset for the conterminous United States using homogenized station data and remotely sensed land skin temperature. International Journal of Climatology 35, 2258–2279, https://doi.org/10.1002/joc.4127 (2015).ADS 
Article 

Google Scholar 
McAfee, S. A., McCabe, G. J., Gray, S. T. & Pederson, G. T. Changing station coverage impacts temperature trends in the Upper Colorado River basin. International Journal of Climatology 39, 1517–1538, https://doi.org/10.1002/joc.5898 (2019).ADS 
Article 

Google Scholar 
Beguería, S., Vicente-Serrano, S. M., Tomás-Burguera, M. & Maneta, M. Bias in the variance of gridded data sets leads to misleading conclusions about changes in climate variability. International Journal of Climatology 36, 3413–3422, https://doi.org/10.1002/joc.4561 (2016).ADS 
Article 

Google Scholar 
Thevakaran, A. & Sonnadara, D. Estimating missing daily temperature extremes in Jaffna, Sri Lanka. Theoretical and applied climatology 132, 145–152, https://doi.org/10.1007/s00704-017-2082-0 (2018).ADS 
Article 

Google Scholar 
Hubbard, K. Spatial variability of daily weather variables in the high plains of the USA. Agricultural and Forest Meteorology 68, 29–41, https://doi.org/10.1016/0168-1923(94)90067-1 (1994).ADS 
Article 

Google Scholar 
Camargo, M. B. & Hubbard, K. G. Spatial and temporal variability of daily weather variables in sub-humid and semi-arid areas of the United States high plains. Agricultural and forest meteorology 93, 141–148, https://doi.org/10.1016/S0168-1923(98)00122-1 (1999).ADS 
Article 

Google Scholar 
Brugnara, Y., Good, E., Squintu, A. A., van der Schrier, G. & Brönnimann, S. The EUSTACE global land station daily air temperature dataset. Geoscience Data Journal 6, 189–204, https://doi.org/10.5285/7925ded722d743fa8259a93acc7073f2 (2019).ADS 
Article 

Google Scholar 
Gonzalez-Hidalgo, J. C., Peña-Angulo, D., Brunetti, M. & Cortesi, N. MOTEDAS: a new monthly temperature database for mainland Spain and the trend in temperature (1951–2010). International Journal of Climatology 35, 4444–4463, https://doi.org/10.1002/joc.4298 (2015).ADS 
Article 

Google Scholar 
Gudmundsson, L., Bremnes, J. B., Haugen, J. E. & Engen-Skaugen, T. Technical note: Downscaling RCM precipitation to the station scale using statistical transformations – a comparison of methods. Hydrology and Earth System Sciences 16, 3383–3390, https://doi.org/10.5194/hess-16-3383-2012 (2012).ADS 
Article 

Google Scholar 
Stanley, T., Kirschbaum, D. B., Huffman, G. J. & Adler, R. F. Approximating long-term statistics early in the global precipitation measurement era. Earth Interactions 21, 1–10, https://doi.org/10.1175/EI-D-16-0025.1 (2017).Article 

Google Scholar 
Haimberger, L. Homogenization of radiosonde temperature time series using innovation statistics. Journal of Climate 20, 1377–1403, https://doi.org/10.1175/JCLI4050.1 (2007).ADS 
Article 

Google Scholar 
Menne, M. J. & Williams, C. N. Homogenization of temperature series via pairwise comparisons. Journal of Climate 22, 1700–1717, https://doi.org/10.1175/2008JCLI2263.1 (2009).ADS 
Article 

Google Scholar 
Vincent, L. A., Zhang, X., Bonsal, B. R. & Hogg, W. D. Homogenization of daily temperatures over Canada. Journal of Climate 15, 1322–1334, 10.1175/1520-0442(2002)015 < 1322:HODTOC > 2.0.CO;2 (2002).ADS 
Article 

Google Scholar 
Jin, M. & Dickinson, R. E. Land surface skin temperature climatology: Benefitting from the strengths of satellite observations. Environmental Research Letters 5, 044004, https://doi.org/10.1088/1748-9326/5/4/044004 (2010).ADS 
Article 

Google Scholar 
Wan, Z., Hook, S. & Hulley, G. MOD11A2 MODIS/Terra land surface temperature/emissivity 8-day l3 global 1 km SIN grid v006. NASA EOSDIS Land Processes DAAC 10, https://doi.org/10.5067/MODIS/MOD11A2.006 (2015).Wilson, A. M. & Jetz, W. Remotely sensed high-resolution global cloud dynamics for predicting ecosystem and biodiversity distributions. PLoS biology 14, e1002415, https://doi.org/10.1371/journal.pbio.1002415 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Danielson, J. J. & Gesch, D. B. Global multi-resolution terrain elevation data 2010 (GMTED2010) (US Department of the Interior, US Geological Survey Washington, DC, USA, 2011).Holden, Z. A., Abatzoglou, J. T., Luce, C. H. & Baggett, L. S. Empirical downscaling of daily minimum air temperature at very fine resolutions in complex terrain. Agricultural and Forest Meteorology 151, 1066–1073, https://doi.org/10.1016/j.agrformet.2011.03.011 (2011).ADS 
Article 

Google Scholar 
Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International journal of climatology 37, 4302–4315, https://doi.org/10.1002/joc.5086 (2017).ADS 
Article 

Google Scholar 
Willmott, C. J. & Robeson, S. M. Climatologically aided interpolation (CAI) of terrestrial air temperature. International Journal of Climatology 15, 221–229, https://doi.org/10.1002/joc.3370150207 (1995).ADS 
Article 

Google Scholar 
Hunter, R. D. & Meentemeyer, R. K. Climatologically aided mapping of daily precipitation and temperature. Journal of Applied Meteorology 44, 1501–1510, https://doi.org/10.1175/JAM2295.1 (2005).ADS 
Article 

Google Scholar 
Parmentier, B. et al. Using multi-timescale methods and satellite-derived land surface temperature for the interpolation of daily maximum air temperature in Oregon. International Journal of Climatology 35, 3862–3878, https://doi.org/10.1002/joc.4251 (2015).ADS 
Article 

Google Scholar 
Hengl, T., Heuvelink, G. B. & Rossiter, D. G. About regression-kriging: From equations to case studies. Computers & Geosciences 33, 1301–1315, https://doi.org/10.1016/j.cageo.2007.05.001. Spatial Analysis (2007).Sun, X.-L., Yang, Q., Wang, H.-L. & Wu, Y.-J. Can regression determination, nugget-to-sill ratio and sampling spacing determine relative performance of regression kriging over ordinary kriging? CATENA 181, 104092, https://doi.org/10.1016/j.catena.2019.104092 (2019).Article 

Google Scholar 
Trajkovic, S. & Gocic, M. Evaluation of three wind speed approaches in temperature-based ET 0 equations: a case study in Serbia. Arabian Journal of Geosciences 14, 1–8, https://doi.org/10.1007/s12517-020-06331-5 (2021).Article 

Google Scholar 
Fotheringham, A. S., Brunsdon, C. & Charlton, M. Geographically weighted regression: the analysis of spatially varying relationships (John Wiley & Sons, 2003).Comber, A. et al. A route map for successful applications of geographically weighted regression. Geographical Analysis https://doi.org/10.1111/gean.12316 (2021).Li, X., Zhou, Y., Asrar, G. R. & Zhu, Z. Developing a 1 km resolution daily air temperature dataset for urban and surrounding areas in the conterminous United States. Remote Sensing of Environment 215, 74–84, https://doi.org/10.1016/j.rse.2018.05.034 (2018).ADS 
Article 

Google Scholar 
Wu, J., Zhong, B., Tian, S., Yang, A. & Wu, J. Downscaling of urban land surface temperature based on multi-factor geographically weighted regression. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 12, 2897–2911, https://doi.org/10.1109/JSTARS.2019.2919936 (2019).ADS 
Article 

Google Scholar 
Zhang, G., Zhu, S., Zhang, N., Zhang, G. & Xu, Y. Downscaling hourly air temperature of WRF simulations over complex topography: A case study of Chongli District in Hebei Province, China. Journal of Geophysical Research: Atmospheres 127, e2021JD035542, https://doi.org/10.1029/2021JD035542 (2022).ADS 
Article 

Google Scholar 
Callañaupa Gutierrez, S. et al. Seasonal variability of daily evapotranspiration and energy fluxes in the Central Andes of Peru using eddy covariance techniques and empirical methods. Atmospheric Research 261, 105760, https://doi.org/10.1016/j.atmosres.2021.105760 (2021).Article 

Google Scholar 
Zotarelli, L., Dukes, M. D., Romero, C. C., Migliaccio, K. W. & Morgan, K. T. Step by step calculation of the Penman-Monteith Evapotranspiration (FAO-56 Method). Institute of Food and Agricultural Sciences. University of Florida https://edis.ifas.ufl.edu/pdf/AE/AE45900.pdf (2010).Huerta, A. PISCOeo_pm, a reference evapotranspiration gridded database based on FAO Penman-Monteith in Peru. figshare https://doi.org/10.6084/m9.figshare.c.5633182.v3 (2021).Willmott, C. J., Robeson, S. M. & Matsuura, K. A refined index of model performance. International Journal of climatology 32, 2088–2094, https://doi.org/10.1002/joc.2419 (2012).ADS 
Article 

Google Scholar 
Legates, D. R. & McCabe, G. J. A refined index of model performance: a rejoinder. International Journal of Climatology 33, 1053–1056, https://doi.org/10.1002/joc.3487 (2013).ADS 
Article 

Google Scholar 
Durre, I., Menne, M. J., Gleason, B. E., Houston, T. G. & Vose, R. S. Comprehensive automated quality assurance of daily surface observations. Journal of Applied Meteorology and Climatology 49, 1615–1633, https://doi.org/10.1175/2010JAMC2375.1 (2010).ADS 
Article 

Google Scholar 
Huerta, A. & Lavado-Casimiro, W. Atlas de Zonas Áridas del Perú: una evaluación presente y futura (Servicio Nacional de Meteorología e Hidrología del Perú, Lima, Perú, 2021).Singer, M. B. et al. Hourly potential evapotranspiration at 0.1° resolution for the global land surface from 1981-present. Scientific Data 8, 1–13, https://doi.org/10.1038/s41597-021-01003-9 (2021).Article 

Google Scholar 
Pelosi, A. & Chirico, G. Regional assessment of daily reference evapotranspiration: Can ground observations be replaced by blending ERA5-Land meteorological reanalysis and CM-SAF satellite-based radiation data? Agricultural Water Management 258, 107169, https://doi.org/10.1016/j.agwat.2021.107169 (2021).Article 

Google Scholar 
McCuen, R. H. A sensitivity and error analysis cf procedures used for estimating evaporation. JAWRA Journal of the American Water Resources Association 10, 486–497, https://doi.org/10.1111/j.1752-1688.1974.tb00590.x (1974).ADS 
Article 

Google Scholar 
Coleman, G. & DeCoursey, D. G. Sensitivity and model variance analysis applied to some evaporation and evapotranspiration models. Water Resources Research 12, 873–879, https://doi.org/10.1029/WR012i005p00873 (1976).ADS 
Article 

Google Scholar 
Beven, K. A sensitivity analysis of the Penman-Monteith actual evapotranspiration estimates. Journal of Hydrology 44, 169–190, https://doi.org/10.1016/0022-1694(79)90130-6 (1979).ADS 
Article 

Google Scholar 
Hupet, F. & Vanclooster, M. Effect of the sampling frequency of meteorological variables on the estimation of the reference evapotranspiration. Journal of Hydrology 243, 192–204, https://doi.org/10.1016/S0022-1694(00)00413-3 (2001).ADS 
Article 

Google Scholar 
Ali, M. et al. Sensitivity of Penman–Monteith estimates of reference evapotranspiration to errors in input climatic data. Journal of Agrometeorology 11, 1–8, https://journal.agrimetassociation.org/index.php/jam/article/view/1214 (2009).
Google Scholar 
Field, C. B., Jackson, R. B. & Mooney, H. A. Stomatal responses to increased CO2: implications from the plant to the global scale. Plant, Cell & Environment 18, 1214–1225, https://doi.org/10.1111/j.1365-3040.1995.tb00630.x (1995).Article 

Google Scholar 
Roderick, M. L., Greve, P. & Farquhar, G. D. On the assessment of aridity with changes in atmospheric CO2. Water Resources Research 51, 5450–5463, https://doi.org/10.1002/2015WR017031 (2015).ADS 
CAS 
Article 

Google Scholar 
Swann, A. L., Hoffman, F. M., Koven, C. D. & Randerson, J. T. Plant responses to increasing CO2 reduce estimates of climate impacts on drought severity. Proceedings of the National Academy of Sciences 113, 10019–10024, https://doi.org/10.1073/pnas.1604581113 (2016).ADS 
CAS 
Article 

Google Scholar 
Milly, P. C. & Dunne, K. A. Potential evapotranspiration and continental drying. Nature Climate Change 6, 946–949, https://doi.org/10.1038/nclimate3046 (2016).ADS 
Article 

Google Scholar 
Greve, P., Roderick, M. L. & Seneviratne, S. I. Simulated changes in aridity from the last glacial maximum to 4xCO2. Environmental Research Letters 12, 114021, https://doi.org/10.1088/1748-9326/aa89a3 (2017).ADS 
CAS 
Article 

Google Scholar 
Scheff, J. Drought indices, drought impacts, CO2, and warming: a historical and geologic perspective. Current Climate Change Reports 4, 202–209, https://doi.org/10.1007/s40641-018-0094-1 (2018).Article 

Google Scholar 
Swann, A. L. Plants and drought in a changing climate. Current Climate Change Reports 4, 192–201, https://doi.org/10.1007/s40641-018-0097-y (2018).Article 

Google Scholar 
Lemordant, L., Gentine, P., Swann, A. S., Cook, B. I. & Scheff, J. Critical impact of vegetation physiology on the continental hydrologic cycle in response to increasing CO2. Proceedings of the National Academy of Sciences 115, 4093–4098, https://doi.org/10.1073/pnas.1720712115 (2018).ADS 
CAS 
Article 

Google Scholar 
Yang, Y., Roderick, M. L., Zhang, S., McVicar, T. R. & Donohue, R. J. Hydrologic implications of vegetation response to elevated CO2 in climate projections. Nature Climate Change 9, 44–48, https://doi.org/10.1038/s41558-018-0361-0 (2019).ADS 
CAS 
Article 

Google Scholar 
Greve, P., Roderick, M., Ukkola, A. & Wada, Y. The aridity index under global warming. Environmental Research Letters 14, 124006, https://doi.org/10.1088/1748-9326/ab5046 (2019).ADS 
CAS 
Article 

Google Scholar 
Huerta, A. & Gutierrez, L. PISCOeo_pm. figshare https://doi.org/10.6084/m9.figshare.19686738.v1 (2022). More