Ecology

1.Eisner, R., Seabrook, L. M. & McAlpine, C. A. Are changes in global oil production influencing the rate of deforestation and biodiversity loss?. Biol. Conserv. 196, 147–155. https://doi.org/10.1016/j.biocon.2016.02.017 (2016).Article 

Google Scholar 
2.Fahrig, L. Effects of habitat fragmentation on biodiversity. Annu. Rev. Ecol. Evol. Syst. 34, 487–515. https://doi.org/10.1146/132419 (2003).Article 

Google Scholar 
3.Pfeifer, M. et al. Creation of forest edges has a global impact on forest vertebrates. Nature 551, 187–191. https://doi.org/10.1038/nature24457 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
4.Tilman, D., May, R., Lehman, C. & Nowak, M. Habitat destruction and the extinction debt. Nature 371, 65–66. https://doi.org/10.1038/371065a0 (1994).ADS 
Article 

Google Scholar 
5.Fisher, J. T. & Burton, C. A. Wildlife winners and losers in an oil sands landscape. Front Ecol. Environ. https://doi.org/10.1002/fee.1807 (2018).Article 

Google Scholar 
6.Heim, N., Fisher, J. T., Volpe, J., Clevenger, A. P. & Paczkowski, J. Carnivore community response to anthropogenic landscape change: species-specificity foils generalizations. Landscape Ecol. 34, 2493–2507. https://doi.org/10.1007/s10980-019-00882-z (2019).Article 

Google Scholar 
7.Pereira, H. M., Navarro, L. & Martins, I. Global biodiversity change: the bad, the good, and the unknown. Annu. Rev. Environ. Resour. https://doi.org/10.1146/annurev-environ-042911-093511 (2012).Article 

Google Scholar 
8.Northrup, J. M., Anderson, C. R. Jr. & Wittemyer, G. Quantifying spatial habitat loss from hydrocarbon development through assessing habitat selection patterns of mule deer. Glob Change Biol. 21, 3961–3970. https://doi.org/10.1111/gcb.13037 (2015).ADS 
Article 

Google Scholar 
9.Holbrook, S. J. & Schmitt, R. J. The combined effects of predation risk and food reward on patch selection. Ecology 69, 125–134. https://doi.org/10.2307/1943167 (1988).Article 

Google Scholar 
10.Moody, A. L., Houston, A. I. & McNamara, J. M. Ideal free distributions under predation risk. Behav. Ecol. Sociobiol. 38, 131–143 (1996).Article 

Google Scholar 
11.Dietz, H. & Edwards, P. J. Recognition that causal processes change during plant invasion helps explain conflicts in evidence. Ecology 87, 1359–1367 (2006).Article 

Google Scholar 
12.Hobbs, R. J. & Huenneke, L. F. Disturbance, diversity, and invasion: implications for conservation. Conserv. Biol. 6, 324–337 (1992).Article 

Google Scholar 
13.Van der Graaf, S., Stahl, J., Klimkowska, A. & Drent, J. P. B. Surfing on a green wave—How plant growth drives spring migration in the Barnacle Goose Branta leucopsis. Ardea -Wageningen- 94, 567 (2006).
Google Scholar 
14.Parker, I. M. et al. Impact: toward a framework for understanding the ecological effects of invaders. Biol. Invasions 1, 3–19. https://doi.org/10.1023/A:1010034312781 (1999).Article 

Google Scholar 
15.Pimentel, D., Zuniga, R. & Morrison, D. Update on the environmental and economic costs associated with alien-invasive species in the United States. Ecol. Econ. 52, 273–288. https://doi.org/10.1016/j.ecolecon.2004.10.002 (2005).Article 

Google Scholar 
16.Shackelford, N. et al. Primed for change: developing ecological restoration for the 21st Century. Restor. Ecol. 21, 297–304. https://doi.org/10.1111/rec.12012 (2013).Article 

Google Scholar 
17.Pickell, P. D., Pickell, P. D., Andison, D. W., Coops, N. C. & Gergel, S. E. The spatial patterns of anthropogenic disturbance in the western Canadian boreal forest following oil and gas development. Can. J. For. Res. 45, 732–743. https://doi.org/10.1139/cjfr-2014-0546 (2015).Article 

Google Scholar 
18.Fisher, J. T. & Wilkinson, L. The response of mammals to forest fire and timber harvest in the North American boreal forest. Mammal Rev. 35, 51–81 (2005).Article 

Google Scholar 
19.Wittische, J., Heckbert, S., James, P. M. A., Burton, A. C. & Fisher, J. T. Community-level modelling of boreal forest mammal distribution in an oil sands landscape. Sci. Total Environ. 755, 142500. https://doi.org/10.1016/j.scitotenv.2020.142500 (2021).ADS 
CAS 
Article 
PubMed 

Google Scholar 
20.Hewitt, D. G. Biology and management of white-tailed deer (CRC Press, Boca Raton, 2011).Book 

Google Scholar 
21.McCabe, R. E. & McCabe, T. R. in White tailed deer: ecology and management Ch. Chapter 2, 19–72 (Stackpole, A Wildlife Management Institute Book, 1984).22.Webb, R. The range of white-tailed deer in Alberta (Alberta Fish and Wildlife Division Edmonton, Alberta, 1967).
Google Scholar 
23.Dawe, K. L. & Boutin, S. Climate change is the primary driver of white-tailed deer (Odocoileus virginianus) range expansion at the northern extent of its range; land use is secondary. Ecol. Evol. 6, 6435–6451. https://doi.org/10.1002/ece3.2316 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
24.DeCesare, N. J., Hebblewhite, M., Robinson, H. S. & Musiani, M. Endangered, apparently: the role of apparent competition in endangered species conservation. Anim. Conserv. 13, 353–362. https://doi.org/10.1111/j.1469-1795.2009.00328.x (2010).Article 

Google Scholar 
25.Latham, A. D. M., Latham, M. C., McCutchen, N. A. & Boutin, S. Invading white-tailed deer change wolf-caribou dynamics in northeastern Alberta. J. Wildl. Manag. 75, 204–212. https://doi.org/10.1002/jwmg.28 (2011).Article 

Google Scholar 
26.Latham, A. D. M., Latham, M. C., Boyce, M. C. & Boutin, S. Movement responses by wolves to industrial linear features and their effect on woodland caribou in northeastern Alberta. Ecol. Appl. 21, 11 (2011).Article 

Google Scholar 
27.Fisher, J. T., Burton, A. C., Nolan, L. & Roy, L. Influences of landscape change and winter severity on invasive ungulate persistence in the Nearctic boreal forest. Sci. Rep. 10, 8742. https://doi.org/10.1038/s41598-020-65385-3 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
28.Dabros, A., Pyper, M. & Castilla, G. Seismic lines in the boreal and arctic ecosystems of North America: environmental impacts, challenges, and opportunities. Environ. Rev. 26, 214–229. https://doi.org/10.1139/er-2017-0080 (2018).Article 

Google Scholar 
29.Dickie, M., Serrouya, R., McNay, R. S., Boutin, S. & du Toit, J. Faster and farther: wolf movement on linear features and implications for hunting behaviour. J. Appl. Ecol. 54, 253–263. https://doi.org/10.1111/1365-2664.12732 (2017).Article 

Google Scholar 
30.Finnegan, L., MacNearney, D. & Pigeon, K. E. Divergent patterns of understory forage growth after seismic line exploration: implications for caribou habitat restoration. For. Ecol. Manag. 409, 634–652. https://doi.org/10.1016/j.foreco.2017.12.010 (2018).Article 

Google Scholar 
31.Prokopenko, C. M., Boyce, M. S., Avgar, T. & Tulloch, A. Characterizing wildlife behavioural responses to roads using integrated step selection analysis. J. Appl. Ecol. 54, 470–479. https://doi.org/10.1111/1365-2664.12768 (2017).Article 

Google Scholar 
32.Waring, G. H., Griffis, J. L. & Vaughn, M. E. White-tailed deer roadside behavior, wildlife warning reflectors, and highway mortality. Appl. Anim. Behav. Sci. 29, 215–223. https://doi.org/10.1016/0168-1591(91)90249-W (1991).Article 

Google Scholar 
33.Bowman, J., Ray, J. C., Magoun, A. J., Johnson, D. S. & Dawson, F. N. Roads, logging, and the large-mammal community of an eastern Canadian boreal forest. Can. J. Zool. 88, 454–467. https://doi.org/10.1139/z10-019 (2010).Article 

Google Scholar 
34.Munro, K. G., Bowman, J. & Fahrig, L. Effect of paved road density on abundance of white-tailed deer. Wildl. Res. 39, 478. https://doi.org/10.1071/wr11152 (2012).Article 

Google Scholar 
35.Fisher, J. T. & Burton, A. C. Spatial structure of reproductive success infers mechanisms of ungulate invasion in Nearctic boreal landscapes. Ecol. Evol. 11, 900–911. https://doi.org/10.1002/ece3.7103 (2021).Article 
PubMed 

Google Scholar 
36.Kie, J. G. Optimal foraging and risk of predation effects on behavior and social structure in ungulates. J. Mammal. 80, 1114–1129 (1999).Article 

Google Scholar 
37.Brown, J. S., Laundré, J. W. & Gurung, M. The ecology of fear: optimal foraging, game theory, and trophic interactions. J. Mammal. 80, 385–399. https://doi.org/10.2307/1383287 (1999).Article 

Google Scholar 
38.Kittle, A. M., Fryxell, J. M., Desy, G. E. & Hamr, J. The scale-dependent impact of wolf predation risk on resource selection by three sympatric ungulates. Oecologia 157, 163–175. https://doi.org/10.1007/s00442-008-1051-9 (2008).ADS 
Article 
PubMed 

Google Scholar 
39.Moen, A. N. Energy conservation by white-tailed deer in the winter. Ecology 57, 192–198. https://doi.org/10.2307/1936411 (1976).Article 

Google Scholar 
40.Schmidt, K. Winter ecology of nonmigratory Alpine red deer. Oecologia 95, 226–233. https://doi.org/10.1007/BF00323494 (1993).ADS 
Article 
PubMed 

Google Scholar 
41.Kilgo, J. C., Ray, H. S., Vukovich, M., Goode, M. J. & Ruth, C. Predation by coyotes on white-tailed deer neonates in South Carolina. J. Wildl. Manag. https://doi.org/10.1002/jwmg.393 (2012).Article 

Google Scholar 
42.Laurent, M., Dickie, M., Becker, M., Serrouya, R. & Boutin, S. Evaluating the mechanisms of landscape change on white-tailed deer populations. J. Wildl. Manag. 85, 340–353. https://doi.org/10.1002/jwmg.21979 (2020).Article 

Google Scholar 
43.Schneider, R. R., Hauer, G., Adamowicz, W. L. & Boutin, S. Triage for conserving populations of threatened species: the case of woodland caribou in Alberta. Biol. Conserv. 143, 1603–1611. https://doi.org/10.1016/j.biocon.2010.04.002 (2010).Article 

Google Scholar 
44.Kilkenny, C., Browne Wj Fau – Cuthill, I. C., Cuthill Ic Fau – Emerson, M., Emerson M Fau – Altman, D. G. & Altman, D. G. Improving bioscience research reporting: the ARRIVE guidelines for reporting animal research. PLoS biol. 8(6), e1000412 (2010).45.DelGiudice, G. D., Mangipane, B. A., Sampson, B. A. & Kochanny, C. O. Chemical immobilization, body temperature, and post-release mortality of white-tailed deer captured by clover trap and net-gun. Wildl. Soc. Bull. (1973-2006) 29, 1147–1157 (2001).
Google Scholar 
46.Droge, E., Creel, S., Becker, M. S. & M’Soka, J. Risky times and risky places interact to affect prey behaviour. Nat. Ecol. Evol. 1, 1123–1128. https://doi.org/10.1038/s41559-017-0220-9 (2017).Article 
PubMed 

Google Scholar 
47.Kunkel, K. E. & Mech, L. D. Wolf and bear predation on white-tailed deer fawns in northeastern Minnesota. Can. J. Zool. 72, 1557–1565 (1994).Article 

Google Scholar 
48.Latham, A., Latham, M., Knopff, K., Hebblewhite, M. & Boutin, S. Wolves, white-tailed deer, and beaver: Implications of seasonal prey switching for woodland caribou declines. Ecography https://doi.org/10.1111/j.1600-0587.2013.00035.x (2013).Article 

Google Scholar 
49.Alberta Environment and Sustainable Resource Development. Alberta Vegetation Index. Accessed October 2016. https://geodiscover.alberta.ca/50.Manly, B., McDonald, L., Thomas, D., McDonald, T. & Erickson, W.Resource selection by animals: statistical design and analysis for field studies. Vol. 63, pp. 1-10 (Springer Science & Business Media, 2007).51.Boyce, M. S., Vernier, P. R., Nielsen, S. E. & Schmiegelow, F. K. A. Evaluating resource selection functions. Ecol. Model. 157, 281–300. https://doi.org/10.1016/S0304-3800(02)00200-4 (2002).Article 

Google Scholar 
52.Hijmans, R. & van Etten, J. Raster: Geographic data analysis and modeling. CRAN R package 2 (2016).53.R: A language and environment for statistical computing. (Vienna, Austria, 2013).54.Zuur, A., Hilbe, J. & Ieno, E. A Beginner’s Guide to GLM and GLMM with R: a frequentist and Bayesian perspective for ecologists. (Highland Statistics, 2013).55.Gillies, C. S. et al. Application of random effects to the study of resource selection by animals. J. Anim. Ecol. 75, 887–898. https://doi.org/10.1111/j.1365-2656.2006.01106.x (2006).Article 
PubMed 

Google Scholar 
56.Craney, T. A. & Surles, J. G. Model-dependent variance inflation factor cutoff values. Qual. Eng. 14, 391–403. https://doi.org/10.1081/QEN-120001878 (2002).Article 

Google Scholar 
57.Akaike, H. Information theory and an extension of the maximum likelihood principle. Selected papers of hirotugu akaike 199–213 (Springer, New York, 1998).Book 

Google Scholar 
58.Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding AIC and BIC in model selection. Sociol. Methods Res. 33, 261–304. https://doi.org/10.1177/0049124104268644 (2004).MathSciNet 
Article 

Google Scholar 
59.Boulanger, Y. et al. Climate change impacts on forest landscapes along the Canadian southern boreal forest transition zone. Landscape Ecol. 32, 1415–1431. https://doi.org/10.1007/s10980-016-0421-7 (2017).Article 

Google Scholar 
60.Sulla-Menashe, D., Woodcock, C. E. & Friedl, M. A. Canadian boreal forest greening and browning trends: an analysis of biogeographic patterns and the relative roles of disturbance versus climate drivers. Environ. Res. Lett. 13, 014007. https://doi.org/10.1088/1748-9326/aa9b88 (2018).ADS 
Article 

Google Scholar 
61.St-Pierre, F., Drapeau, P. & St-Laurent, M.-H. Drivers of vegetation regrowth on logging roads in the boreal forest: Implications for restoration of woodland caribou habitat. For. Ecol. Manag. 482, 118846. https://doi.org/10.1016/j.foreco.2020.118846 (2021).Article 

Google Scholar 
62.Berger, J. Fear, human shields and the redistribution of prey and predators in protected areas. Biol. Let. 3, 620–623. https://doi.org/10.1098/rsbl.2007.0415 (2007).Article 

Google Scholar 
63.Heyes, A., Leach, A. & Mason, C. F. The economics of Canadian oil sands. Rev. Environ. Econ. Policy 12, 242–263. https://doi.org/10.1093/reep/rey006 (2018).Article 

Google Scholar 
64.Komers, P. E. & Stanojevic, Z. Rates of disturbance vary by data resolution: implications for conservation schedules using the Alberta boreal forest as a case study. Global Change Biol. 19, 2916–2928 (2013).ADS 
CAS 
Article 

Google Scholar 
65.Hebblewhite, M. & Merrill, E. H. Trade-offs between predation risk and forage differ between migrant strategies in a migratory ungulate. Ecology 90, 3445–3454. https://doi.org/10.1890/08-2090.1 (2009).Article 
PubMed 

Google Scholar 
66.Mech, D. L. & Boitani, L. Wolves: behavior, ecology, and conservation Vol. 57 (University of Chicago Press, Chicago, 2004).
Google Scholar 
67.Creel, S., Winnie, J. A., Christianson, D. & Liley, S. Time and space in general models of antipredator response: tests with wolves and elk. Anim. Behav. 76, 1139–1146. https://doi.org/10.1016/j.anbehav.2008.07.006 (2008).Article 

Google Scholar 
68.Steenweg, R. et al. Scaling-up camera traps: monitoring the planet’s biodiversity with networks of remote sensors. Front. Ecol. Environ. 15, 26–34. https://doi.org/10.1002/fee.1448 (2017).Article 

Google Scholar 
69.Hebblewhite, M. Billion dollar boreal woodland caribou and the biodiversity impacts of the global oil and gas industry. Biol. Cons. 206, 102–111. https://doi.org/10.1016/j.biocon.2016.12.014 (2017).Article 

Google Scholar 
70.Côté, S. D., Rooney, T. P., Tremblay, J.-P., Dussault, C. & Waller, D. M. Ecological impacts of deer overabundance. Annu. Rev. Ecol. Evol. Syst. 35, 113–147 (2004).Article 

Google Scholar 
71.McCullough, D. R. Evaluation of night spotlighting as a deer study technique. J. Wildl. Manag. 46, 963–973. https://doi.org/10.2307/3808229 (1982).Article 

Google Scholar 
72.Preston, T., Wildhaber, M., Green, N., Albers, J. & Debenedetto, G. Enumerating white-tailed deer using unmanned aerial vehicles. Wildlife Soc. Bull. https://doi.org/10.1002/wsb.1149 (2021).Article 

Google Scholar 
73.Parks, A. E. Provincial woodland caribou range plan. 212 (Edmonton, Alberta, 2017).74.Tattersall, E. R., Burgar, J. M., Fisher, J. T. & Burton, A. C. Boreal predator co-occurrences reveal shared use of seismic lines in a working landscape. Ecol. Evol. 10, 1678–1691. https://doi.org/10.1002/ece3.6028 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
75.Diaz, S. et al. Pervasive human-driven decline of life on Earth points to the need for transformative change. Science (New York N.Y.) https://doi.org/10.1126/science.aax3100 (2019).Article 
PubMed Central 

Google Scholar 
76.Bayoumi, T. & Muhleisen, M. Energy, the exchange rate, and the economy: macroeconomic benefits of Canada’s oil sands production (International Monetary Fund, Washington, 2006).
Google Scholar 
77.Zhu, K., Song, Y. & Qin, C. Forest age improves understanding of the global carbon sink. Proc. Natl. Acad. Sci. 116, 3962. https://doi.org/10.1073/pnas.1900797116 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar  More

1.Zhang QQ, Ying GG, Pan CG, Liu YS, Zhao JL. Comprehensive evaluation of antibiotics emission and fate in the river basins of China: source analysis, multimedia modeling, and linkage to bacterial rResistance. Environ Sci Technol. 2015;49:6772–82.CAS 
Article 

Google Scholar 
2.Zhao WX, Wang B, Yu G. Antibiotic resistance genes in China: occurrence, risk, and correlation among different parameters. Environ Sci Pollut R. 2018;25:21467–82.CAS 
Article 

Google Scholar 
3.Han XM, Hu HW, Chen QL, Yang LY, Li HL, Zhu YG, et al. Antibiotic resistance genes and associated bacterial communities in agricultural soils amended with different sources of animal manures. Soil Biol Biochem. 2018;126:91–102.CAS 
Article 

Google Scholar 
4.Huerta B, Marti E, Gros M, López P, Pompêo M, Armengol J, et al. Exploring the links between antibiotic occurrence, antibiotic resistance, and bacterial communities in water supply reservoirs. Sci Total Environ. 2013;456:161–70.Article 

Google Scholar 
5.Martinez JL, Sánchez MB, Martínez-Solano L, Hernandez A, Garmendia L, Fajardo A, et al. Functional role of bacterial multidrug efflux pumps in microbial natural ecosystems. Fems Microbiol Rev. 2009;33:430–49.CAS 
Article 

Google Scholar 
6.Wright GD. The antibiotic resistome: the nexus of chemical and genetic diversity. Nat Rev Microbiol. 2007;5:175–86.CAS 
Article 

Google Scholar 
7.Meng F, Yang S, Wang X, Chen T, Wang X, Tang X, et al. Reclamation of Chinese herb residues using probiotics and evaluation of their beneficial effect on pathogen infection. J Infect Public Health. 2017;10:749–54.Article 

Google Scholar 
8.Zhou Y, Selvam A, Wong JWC. Chinese medicinal herbal residues as a bulking agent for food waste composting. Bioresour Technol. 2018;249:182–8.CAS 
Article 

Google Scholar 
9.Wu HW, Sun XQ, Liang BW, Chen JB, Zhou XF. Analysis of livestock and poultry manure pollution in China and its treatment and resource utilization. J Agro-Environ Sci. 2020;39:1168–76.
Google Scholar 
10.Chen J, Yu Z, Michel FC Jr., Wittum T, Morrison M. Development and application of real-time PCR assays for quantification of erm genes conferring resistance to macrolides-lincosamides-streptogramin B in livestock manure and manure management systems. Appl Environ Microbiol. 2007;73:4407–16.CAS 
Article 

Google Scholar 
11.Duan M, Gu J, Wang X, Li Y, Zhang S, Yin Y, et al. Effects of genetically modified cotton stalks on antibiotic resistance genes, intI1, and intI2 during pig manure composting. Ecotoxicol Environ Saf. 2018;147:637–42.CAS 
Article 

Google Scholar 
12.Cui E, Wu Y, Zuo Y, Chen H. Effect of different biochars on antibiotic resistance genes and bacterial community during chicken manure composting. Bioresour Technol. 2016;203:11–7.CAS 
Article 

Google Scholar 
13.Ma Y, Wilson CA, Novak JT, Riffat R, Aynur S, Murthy S, Pruden A. Effect of various sludge digestion conditions on sulfonamide, macrolide, and tetracycline Resistance Genes and Class I Integrons. Environ Sci Technol. 2011;45:7855–61.CAS 
Article 

Google Scholar 
14.Tien YC, Li B, Zhang T, Scott A, Murray R, Sabourin L, et al. Impact of dairy manure pre-application treatment on manure composition, soil dynamics of antibiotic resistance genes, and abundance of antibiotic-resistance genes on vegetables at harvest. Sci Total Environ. 2017;581-582:32–9.CAS 
Article 

Google Scholar 
15.Zhang L, Sun XY. Effects of waste lime and Chinese medicinal herbal residue amendments on physical, chemical, and microbial properties during green waste composting. Environ Sci Pollut Res. Int. 2018;25:31381–95.CAS 
Article 

Google Scholar 
16.Wang YQ, Wu XQ, Zhu TT, Ma QG, Chen HG. Study on utilization of solid slag compost of Chinese medicinal herbal. J Chin Medicinal Mater. 2008;31:1622–4.CAS 

Google Scholar 
17.Wu DL, Liu P, Luo YZ, Tian GM, Mahmood Q. Nitrogen transformations during co-composting of herbal residues, spent mushrooms, and sludge. J Zhejiang Univ Sci B. 2010;11:497–505.Article 

Google Scholar 
18.Ward T, Larson J, Meulemans J, Hillmann B, Lynch J, Sidiropoulos D, et al. BugBase predicts organism-level microbiome phenotypes. bioRxiv. 2017;133462.19.Chao A. Nonparametric estimation of the number of classes in a population. Scand J Stat. 1984;11:265–70.
Google Scholar 
20.Chao A, Yang MCK. Stopping rules and estimation for recapture debugging with unequal failure rates. Biometrika. 1993;80:193–201.Article 

Google Scholar 
21.Shannon CE. A mathematical theory of communication. Bell Syst Tech J. 1948;27:623–56.Article 

Google Scholar 
22.Simpson EH. Measurement of diversity. Nature 1949;163:688.Article 

Google Scholar 
23.Huang K, Xia H, Wu Y, Chen J, Cui G, Li F, et al. Effects of earthworms on the fate of tetracycline and fluoroquinolone resistance genes of sewage sludge during vermicomposting. Bioresour Technol. 2018;259:32–9.CAS 
Article 

Google Scholar 
24.Qian X, Sun W, Gu J, Wang XJ, Sun JJ, Yin YN, et al. Variable effects of oxytetracycline on antibiotic resistance gene abundance and the bacterial community during aerobic composting of cow manure. J Hazard Mater. 2016;315:61–9.CAS 
Article 

Google Scholar 
25.Zhang R, Gu J, Wang X, Li Y, Zhang K, Yin Y, Zhang X. Contributions of the microbial community and environmental variables to antibiotic resistance genes during co-composting with swine manure and cotton stalks. J Hazard Mater. 2018;358:82–91.CAS 
Article 

Google Scholar 
26.Wang H, Sangwan N, Li HY, Su JQ, Oyang WY, Zhang ZJ, et al. The antibiotic resistome of swine manure is significantly altered by association with the Musca domestica larvae gut microbiome. Isme J. 2017;11:100–11.Article 

Google Scholar 
27.Li J, Xin Z, Zhang Y, Chen J, Yan J, Li H, Hu H. Long-term manure application increased the levels of antibiotics and antibiotic resistance genes in a greenhouse soil. Appl Soil Ecol. 2017;121:193–200.Article 

Google Scholar 
28.Su JQ, Wei B, Ou-Yang WY, Huang FY, Zhao Y, Xu HJ, et al. Antibiotic resistome and its association with bacterial communities during sewage sludge composting. Environ Sci Technol. 2015;49:7356–63.CAS 
Article 

Google Scholar 
29.Li H, Duan M, Gu J, Zhang Y, Qian X, Ma J, et al. Effects of bamboo charcoal on antibiotic resistance genes during chicken manure composting. Ecotoxicol Environ Saf. 2017;140:1–6.Article 

Google Scholar 
30.Zhang J, Lin H, Ma J, Sun W, Yang Y, Zhang X. Compost-bulking agents reduce the reservoir of antibiotics and antibiotic resistance genes in manures by modifying bacterial microbiota. Sci Total Environ. 2019;649:396–404.CAS 
Article 

Google Scholar 
31.Ghosh S, Ramsden SJ, LaPara TM. The role of anaerobic digestion in controlling the release of tetracycline resistance genes and class 1 integrons from municipal wastewater treatment plants. Appl Microbiol Biotechnol. 2009;84:791–6.CAS 
Article 

Google Scholar 
32.Selvam A, Xu D, Zhao Z, Wong JW. Fate of tetracycline, sulfonamide and fluoroquinolone resistance genes and the changes in bacterial diversity during composting of swine manure. Bioresour Technol. 2012;126:383–90.CAS 
Article 

Google Scholar 
33.Antunes P, Machado J, Sousa JC, Peixe L. Dissemination of sulfonamide resistance genes (sul1, sul2, and sul3) in Portuguese Salmonella enterica strains and relation with integrons. Antimicrob Agents Chemother. 2005;49:836–9.CAS 
Article 

Google Scholar 
34.Zhu YG, Johnson TA, Su JQ, Qiao M, Guo GX, Stedtfeld RD, et al. Diverse and abundant antibiotic resistance genes in Chinese swine farms. Proc Natl Acad Sci USA. 2013;110:3435–40.CAS 
Article 

Google Scholar 
35.Chen Q, An X, Li H, Su J, Ma Y, Zhu YG. Long-term field application of sewage sludge increases the abundance of antibiotic resistance genes in soil. Environ Int. 2016;92-93:1–10.CAS 
Article 

Google Scholar  More

The GCBD is stored at figshare23. Below we describe 20 Tables (also see Fig. 3 schematic) that comprise the GCBD: (1) Site_Info_tbl, (2) Sample_Event_tbl, (3) R_Scripts_tbl, (4) Cover_tbl, (5) Bleaching_tbl, (6) Environmental_tbl, (7) Authors_LUT, (8) Bleaching_Level_LUT, (9) City_Town_Name_LUT, (10) Country_Name_LUT, (11) Data_Source_LUT, (12) Ecoregion_Name_LUT, (13) Exposure_LUT, (14) Ocean_Name_LUT, (15) Realm_Name_LUT, (16) State_Island_Province_Name_LUT, (17) Substrate_Type_LUT, (18) Relevant_Papers_tbl, (19) Severity_Code_LUT, and (20) Bleaching_Prevalence_Score_LUT, where LUT stands for look-up table.

1)

Site Information (Site_Info_tbl)
Latitude_Degrees: latitude coordinates in decimal degrees.
Longitude_Degrees: longitude coordinates in decimal degrees.
Ocean_Name: the ocean in which the sampling took place.
Realm_Name: identification of realm as defined by the Marine Ecoregions of the World (MEOW)12.
Ecoregion_Name: identification of the Ecoregions (150) as defined by Veron et al.13.
Country_Name: the country where sampling took place.
State_Island_Province_Name: the state, territory (e.g., Guam) or island group (e.g., Hawaiian Islands) where sampling took place.
City_Town_Name: the region, city, or nearest town, where sampling took place.
Site_Name: the accepted name of the site or the name given by the team that sampled the reef.
Distance_to_Shore: the distance (m) of the sampling site from the nearest land.
Exposure: a site was considered exposed if it had >20 km of fetch, if there were strong seasonal winds, or if the site faced the prevailing winds. Otherwise, the site was considered sheltered or ‘sometimes’. ‘Sometimes’ refers to a few sites with a >20 km fetch through a narrow geographic window, and therefore we considered that the site was potentially exposed during cyclone seasons. We left the category ‘sometimes’ in the database because those sites were not clearly exposed sites, nor were they clearly sheltered sites, and future researchers may be interested in temporary exposure.
Turbidity: kd490 with a 100-km buffer.
Cyclone_Frequency: number of cyclone events from 1964 to 2014.
Comments: comments of any issues with the site or additional information.

2)

Sample Event Information (Sample_Event_tbl)
Site_ID: site ID field from Site_Info_tbl.
Reef_ID: name of reef site that was adopted by sampling group (from ReefCheck).
Quadrat_No: quadrat number (from McClanahan et al.)20.
Date_Day: the date of the sampling event.
Date_Month: the month of sampling event.
Date_Year: the year of sampling event.
Depth: depth (m) of sampling site. Comments: comments of any issue or additional information of sampling event.

3)

R Code (R_Scripts_tbl)
Relevant_Papers_ID: relevant papers ID field from Relevant_Papers_tbl.
Project name: name of project associated with R code.
Paper_Title: title of paper where R code was published.
Code_Name: name of R code file.
Description: description of the R code.
Data_Source: data source ID field from Data_Source_LUT.
R_Code: attachment of R code file.
URL: hyperlink to R code or link to github.

4)

Coral Cover Information (Cover_tbl)
Sample_ID: sampled ID field from Sample_Event_tbl.
Substrate_Type: substrate type ID field from Substrate_LUT.
S1: Reef Check breaks down transects into four 20 m × 5 m segments, point data from segment one of transect.
S2: Reef Check breaks down transects into four 20 m × 5 m segments, point data from segment two of transect.
S3: Reef Check breaks down transects into four 20 m × 5 m segments, point data from segment three of transect.
S4: Reef Check breaks down transects into four 20 m × 5 m segments, point data from segment four of transect.
Perc_hardcoral: percent hard coral cover from McClanahan et al.20 data source.
Perc_macroalgae: percent macroalgae cover from McClanahan et al.20 data source.
Average_Ellipse_Transect: calculated percent hard coral cover per 10 m × 1 m transect using ellipse equation.
Average_Ellipse_Site: calculated percent hard coral cover per site using ellipse equation.
Comments: comments of any issue or additional information of sampling event

5)

Bleaching Information (Bleaching_tbl)
Sample_ID: sample ID field from Sample_Event_tbl.
Bleaching_Level: Reef Check data, coral population or coral colony.
S1: Reef Check breaks down transects into four 20 m × 5 m segments, percent bleaching from segment one of transect.
S2: Reef Check breaks down transects into four 20 m × 5 m segments, percent bleaching from segment two of transect.
S3: Reef Check breaks down transects into four 20 m × 5 m segments, percent bleaching from segment three of transect.
S4: Reef Check breaks down transects into four 20 m × 5 m segments, percent bleaching from segment four of transect.
Percent_Bleaching_RC_Old_Method: old method of determining percent bleaching from Reef_Check.
Severity_Code: coded range of bleaching severity from Donner et al.10.
Percent_Bleached: percent of coral bleaching.
Number_Bleached_colonies: number of bleached corals from McClanahan et al.20 data source.
Bleaching_intensity: from McClanahan et al.20 data source.
Bleaching_Prevalence_Score: coded range of bleaching prevalence from Safaie et al.21.

6)

Environmental Parameter Information (Environmental_tbl)
Sample_ID: sample ID field from Sample_Event_tbl.
ClimSST: CoRTAD. [Climatological Sea-Surface Temperature (SST)] based on weekly SSTs for the study time frame, created using a harmonics approach.
Temperature_ Kelvin: CoRTAD. SST in Kelvin.
Temperature_Mean: CoRTAD. Mean SST in degrees Celsius.
Temperature_Minimum: CoRTAD. Minimum SST in degrees Celsius.
Temperature_Maximum: CoRTAD. Maximum SST in degrees Celsius.
Temperature_Kelvin_Standard_Deviation: CoRTAD. Standard deviation of SST in Kelvin.
Windspeed: CoRTAD. meters per hour.
SSTA: CoRTAD. (Sea-Surface Temperature Anomaly) weekly SST minus weekly climatological SST.
SSTA_Standard_Deviation: CoRTAD. The Standard Deviation of weekly SSTA in degrees Celsius over the entire period.
SSTA_Mean: CoRTAD. The mean SSTA in degrees Celsius over the entire period.
SSTA_Minimum: CoRTAD. The minimum SSTA in degrees Celsius over the entire period.
SSTA_Maximum: CoRTAD. The maximum SSTA in degrees Celsius over the entire period.
SSTA_Frequency: CoRTAD. (Sea Surface Temperature Anomaly Frequency) number of times over the previous 52 weeks that SSTA  >  = 1 degree Celsius.
SSTA_Frequency_Standard_Deviation: CoRTAD. The standard deviation of SSTA Frequency in degrees Celsius over the entire time period of 40 years.
SSTA_FrequencyMax: CoRTAD. The maximum SSTA Frequency in degrees Celsius over the entire time period.
SSTA_FrequencyMean: CoRTAD. The mean SSTA Frequency in degrees Celsius over the entire time period of 40 years.
SSTA_DHW: CoRTAD. (Sea Surface Temperature Degree Heating Weeks) sum of previous 12 weeks when SSTA  >  = 1 degree Celsius.
SSTA_DHW_Standard_Deviation: CoRTAD. The standard deviation SSTA DHW in degrees Celsius over the entire period.
SSTA_DHWMax: CoRTAD. The maximum SSTA DHW in degrees Celsius over the entire time period of 40 years.
SSTA_DHWMean: CoRTAD. The mean SSTA DHW in degrees Celsius over the entire time period of 40 years.
TSA: CoRTAD. (Thermal Stress Anomaly) weekly SSTs minus the maximum of weekly climatological SSTs in degrees Celsius.
TSA_Standard_Deviation: CoRTAD. The standard deviation of TSA in degrees Celsius over the entire time period of 40 years.
TSA_Minimum: CoRTAD. The minimum TSA in degrees Celsius over the entire time period of 40 years.
TSA_Maximum: CoRTAD. The maximum TSA in degrees Celsius over the entire time period of 40 years.
TSA_Mean: CoRTAD. The mean TSA in degrees Celsius over the entire time period of 40 years.
TSA_Frequency: CoRTAD. The number of times over previous 52 weeks that TSA  >  = 1 degree Celsius.
TSA_Frequency_Standard_Deviation: CoRTAD. The standard deviation of frequency of TSA in degrees Celsius over the entire time period of 40 years.
TSA_FrequencyMax: CoRTAD. The maximum TSA frequency in degrees Celsius over the entire time period of 40 years.
TSA_FrequencyMean: CoRTAD. The mean TSA frequency in degrees Celsius over the entire time period of 40 years.
TSA_DHW: CoRTAD. (Thermal Stress Anomaly Degree Heating Weeks) sum of previous 12 weeks when TSA  >  = 1 degree Celsius.
TSA_DHW_Standard_Deviation: CoRTAD. The standard deviation of TSA DHW in degrees Celsius over the entire time period of 40 years.
TSA_DHWMax: CoRTAD. The maximum TSA DHW in degrees Celsius over the entire time period of 40 years.
TSA_DHWMean: CoRTAD. The mean TSA DHW in degrees Celsius over the entire time period of 40 years.

7)

Author Names (Authors_LUT)
Last_Name: author’s last name.
First_Name: author’s first name.
Middle_Initial: author’s middle initial.

8)

Bleaching Level Information (Bleaching_Level_LUT)
Bleaching_Level: Reef Check data, coral population or coral colony.

9)

City, Town Names (City_Town_Name_LUT)
City_Town_Name: the region, city, or town, where sampling took place.

10)

Country names (Country_Name_LUT)
Country_Name: name of the country where sampling took place.

11)

Data Source Information (Data_Source_LUT)
Data_Source: name of source of original data set.
Sample_Method: Description of the sampling methods used to collect the data. If more than one method was used then we stated that an amalgamation of methods were used to collect the data, and the original papers are found in “Relevant_Papers_tbl”, and can be referenced therein.

12)

Ecoregion Names (Ecoregion_Name_LUT)
Ecoregion_Name: name of Ecoregion from Veron et al.13.

13)

Exposure Type (Exposure_LUT)
Exposure_Type: site exposure to fetch.

14)

Ocean Name Information (Ocean_Name_LUT)
Ocean_Name: name of ocean where sampling took place.

15)

Name of Realm (Realm_Name_LUT)
Realm_Name: name of realm as identified by the Marine Ecoregions of the World (MEOW)12.

16)

State, Island, Province Name (State_Island_Province_Name_LUT)
State_Island_Province_Name, Name of the state, territory (e.g. Guam) or island group (e.g. Hawaiian Islands) where sampling took place.

17)

Substrate Type (Substrate_Type_LUT)
Substrate_Type: type of substrate from Reef Check data.

18)

Relevant Publications (Relevant_Papers_tbl)
Data_Source: source associated with publication.
Author_ID: author ID field from Authors_LUT.
Title: title of published work.
Journal_Name: name of publication journal.
Year_Published: year of publication.
Volume: volume number of journal.
Issue: issue number of journal.
Pages: page range of publication.
URL: hyperlink to publication.
DOI: DOI number of publication.
pdf: pdf attachment of publication.

19)

Severity Index Code (Severity_Code_LUT)
Severity_Code: coded range of bleaching severity from Donner et al.10.

20)

Bleaching Prevalence Code (Bleaching_Prevalence_Score_LUT)

Bleaching_Prevalence_Score: coded range of bleaching prevalence from Safaie et al. 21. More

Here we describe methods for the GAEZ+ 2015 Annual Crop Data, and the GAEZ+ 2015 Monthly Cropland Data. The Annual Crop Data was generated first, then the Monthly Cropland Data was calculated based on the Harvest Area results of the Annual Data (Fig. 1).Fig. 1Schematic overview of annual and monthly data production methods. The GAEZ+ 2015 products described in this paper are in dark blue boxes; publicly available data used are in light blue. Dark blue arrows indicate which data are used in each processing step, and grey arrows from steps to data show which steps result in final GAEZ+ 2015 data products. The processing steps listed here are referred to in the Methods section text.Full size imageGAEZ+ 2015 Annual Crop Data MethodsThe GEAZ+ 2015 Annual Crop Data updates the 2010 GAEZ v4 crop harvest area, yield, and production maps6,7 (identified as Theme 5 in ref. 7) using national-scale data on the change in crop harvested area and livestock numbers from 2010 to 2015, based on statistics for 160 crop groups, and cattle and buffalo, from FAOSTAT5.Three datasets were used to produce GAEZ+ 2015 Annual Crop Data:

1.

FAOSTAT crop production domain: annual, country-level data on crop harvested area (H) and crop production (P) for each crop from the FAOSTAT database (Table 1)Table 1 GAEZ and FAOSTAT crop harmonization.Full size table

2.

GAEZ v46,7 gridded global annual harvested area, yield, and production by crop for the 26 FAOSTAT crops and crop categories at 5-minute resolution

3.

Global Administrative Unit Layer (GAUL 2012)13 data. GAUL 2012 reports the fraction of each global 5-minute grid cell that falls within a given country or disputed territory. There are 275 unique global administrative units.

Step 1. Calculate crop changes from 2010 to 2015 by country:
For each country, we extracted the harvested area (H) and crop production (P) for each of the 160 FAOSTAT crop categories, c, from the FAOSTAT database. We averaged three years (2009–2011) of annual national crop harvested area data to represent 2010 national crop harvest area, H2010, and three years (2014–2016) of annual crop harvested area data to represent 2015 national crop harvest area, H2015, then calculated a ratio, rHc, of 2015 to 2010 harvested areas for each crop c in each country, and equivalently, for crop production:$$r{H}_{c}={H}_{2015}/{H}_{2010}$$
(1)
$$r{P}_{c}={P}_{2015}/{P}_{2010}$$
(2)
This results in 160 rH and rP values per country. If harvest area and production values for a particular crop are zero or unreported in the FAOSTAT data, then rHc and rPc are both set to 1.0 (i.e., no change from 2010 to 2015). Three years of data are averaged (2009 – 2011 and 2014 – 2016) to account for missing data for some country/year combinations and to avoid emphasizing reported outliers.
Step 2. Aggregate FAOSTAT-based ratios to the GAEZ crop categories:
We followed the crop aggregation methods of the GAEZ model to aggregate the FAOSTAT crop list (160 unique crops as of 2019) to 26 crops (see Table 1). For each of the 26 GAEZ crop categories, if there is more than one matching FAOSTAT crop (see Table 1) then we applied an area-weighted average (based on FAOSTAT year 2015 harvested area) of the FAOSTAT crops within each country to the rH and rP values for that crop and country. This results in 26 rH and rP values per country. There was one exception to this: the GAEZ_2010 crop category ‘fodder crops’ was an aggregate of 17 FAOSTAT crops (see Table 1) for which harvest area data are no longer reported on FAOSTAT; i.e., GAEZ_2010 had obtained FAOSTAT data on fodder crops circa 2010, but FAOSTAT no longer provides any data on fodder crops for any year. We assumed that the 2010 to 2015 fractional change in fodder crop harvest area in each country was proportional to the change in the FAOSTAT reported national herd sizes for cattle and buffalo livestock data5 for that country, following the same methodology as for crop harvested area change (see Step 2 below). This method assumes a negligible international trade of fodder crops as indicated by bilateral trade matrices available from FAOSTAT.
Step 3. Apply country-level ratios to grid cells:
Calculated country-level ratios were then applied to each grid cell k, using the GAUL_201213 definitions for which grid cells fall within which countries. Some grid cells are split between two or more countries. In this case, all model output variables for the grid cell are divided between the countries based on the fraction of grid cell area falling within the country i:$${H}_{c,2015}^{k}={H}_{c,2010}^{k}{sum }_{i},{f}_{i}^{k}r{H}_{c,i}$$
(3)
$${P}_{c,2015}^{k}={P}_{c,2010}^{k}{sum }_{i},{f}_{i}^{k}r{P}_{c,i}$$
(4)
where ({H}_{c,2015}^{k}) is the year 2015 harvested area (or production) for crop c in grid cell k; ({f}_{i}^{k}) is the fraction of country i in grid cell k, and rHc,i and rPc,i are the ratios for crop c in country i as calculated in Eqs. 1 and 2. This results in 26 H and P values per grid cell. If the sum of all crop harvest areas exceeds 99% of the grid cell area, all crop harvest areas are reduced equally to fit within 99% of the area.
Special Case: Sudan
FAOSTAT data for years before 2011 report data for Sudan, and for South Sudan and Sudan after 2011. To compute the ratios for these grid cells, we split the 2010 data for Sudan into a virtual ‘North’ Sudan and ‘South_Sudan’, using the data for the year 2012, which was reported for both countries. We then used these generated 2010 data and applied the same methodology as described above to calculate changes in harvested areas and production in all grid cells in both countries.
Special Case: Small regions and islands
Forty-nine countries – generally small regions or islands – had no data reported for crop harvested area by FAOSTAT. We assumed that there was no change in crop harvested area for the grid cells in these countries. Note that many may have had zero ha as previously-reported crop area in GAEZ v4. These countries are (the number following each region is the region’s number in ADM0_CODE in the GAUL_2012 data13):Anguilla (9), Aruba (14), Ashmore_and_Cartier_Islands (16), Azores_Islands (74578), Baker_Island (22), Bassas_da_India (25), Bird_Island (32), Bouvet_Island (36), British_Indian_Ocean_Territory (38), Christmas_Island (54), Clipperton_Island (55), Cocos (Keeling)_Islands (56), Europa_Island (80), French_Southern_and_Antarctic_Territories (88), Glorioso_Island (96), Greenland (98), Guernsey (104), Heard_Island_and_McDonald_Islands (109), Howland_Island (112), Isle_of_Man (120), Jarvis_Island (127), Jersey (128), Johnston_Atoll (129), Juan_de_Nova_Island (131), Kingman_Reef (134), Kuril_islands (136), Madeira_Islands (151), Mayotte (161), Midway_Island (164), Navassa_Island (174), Netherlands_Antilles (176), Norfolk_Island (184), Northern_Mariana_Islands (185), Palmyra_Atoll (190), Paracel_Islands (193), Pitcairn (197), Saint_Helena (207), Scarborough_Reef (216), Senkaku_Islands (218), South_Georgia_and_the_South_Sandwich_Islands (228), Spratly_Islands (230), Svalbard_and_Jan_Mayen_Islands (234), Tromelin_Island (247), Turks_and_Caicos_Islands (251), United_States_Virgin_Islands (258), Wake_Island (265), Gibraltar (95), Holy_See (110), Liechtenstein (146).
Special Case: Disputed Areas
Some grid cells in the GAUL_201213 cell-table database are assigned to nine disputed areas, rather than to specific countries. We assumed that there was no change in crop harvested area or production from 2010 to 2015 for grid cells these disputed areas. These areas are (the number following each region is the region’s number of the ADM0_CODE in the GAUL_201213 data):Abyei (102), Aksai_Chin (2), Arunachal_Pradesh (15), China/India (52), Hala’ib_Triangle (40760), Ilemi_Triangle (61013), Jammu_and_Kashmir (40781), Ma’tan_al-Sarra (40762), Falkland_Islands_(Malvinas) (81).
Step 4. Compute 2015 crop yields:
Crop yields were computed for each crop, c, and grid cell, k, as the ratio of crop production to crop harvest area (if harvest area, Hc,k,2015, is zero, then yield, Yc,k,2015, is set to zero):$${Y}_{c,k,2015}={P}_{c,k,2015}/{H}_{c,k,2015}$$
(5)
The resulting gridded global data are:

A.

GAEZ+ 2015 Crop Harvest Area14

B.

GAEZ+ 2015 Crop Yield15

C.

GAEZ+ 2015 Crop Production16

This new data product consists of 156 data files in geotiff format, one rainfed harvested area file and one irrigated harvested area file for each crop harvest area (1000 ha (107 m2) per 5-minute grid cell), crop production (1000 tonnes (106 kg) per 5-minute grid cell), and crop yield (tonnes per ha (10−1 kg m−2) per 5-minute grid cell), for each of the 26 GAEZ crops or crop categories in Table 1.GAEZ+ 2015 monthly cropland area methodsTwo datasets were used to produce monthly cropland area by crop and by irrigated vs rainfed management. These are:

1.

GAEZ+ 2015 Annual Harvested Area14 (as developed above)

2.

MIRCA2000 cropland area4

Step 5. Harmonize the GAEZ+ 2015 and MIRCA2000 crop lists
The MIRCA20004 cropland product provides monthly growing area grids (gridded physical cropland area) for 26 irrigated and rainfed crops and crop categories, as well as cropping calendars that identify the planting month and harvesting month for each crop (via ‘subcrops’ – see below). However, the MIRCA2000 crop list is not the same as the GAEZ+ 2015 crop list; we matched each crop type in the GAEZ+ 2015 crop list to a crop type in the MIRCA2000 crop list to enable the application of MIRCA2000 crop calendars to GAEZ+ 2015 crops (Table 2). Out of the 26 GAEZ+ 2015 crops, 18 had clear 1:1 matching crop categories within MIRCA2000. The remaining 8 crops were matched based on general crop characteristics, i.e., annual vs. perennial, or to unmatched MIRCA2000 cereals.Table 2 List of GAEZ crop categories used in all GAEZ+ 2015 products, as well as the matching between GAEZ+ 2015 crops and MIRCA20004 crop categories for the purposes of producing GAEZ+ 2015 monthly cropland data.Full size tableAn essential component of the MIRCA2000 cropland dataset is the identification of subcrop categories within each crop category to split crops into areas grown in different seasons, or crops with different planting and harvesting dates within the same season. Up to 5 subcrops can be defined to represent such multi-cropping practices. Below, we use the following notation:HG = annual harvested area from the GAEZ+ 2015 product for a given cropHM = annual harvested area calculated from the MIRCA2000 data for a given cropAM,n = cropland area of MIRCA2000 crop, subcrop n, by monthAG,n = cropland area of GAEZ+ 2015 crop, subcrop n, by monthAG = cropland area of GAEZ+ 2015 crop, by month
Step 6. Apply MIRCA2000 monthly crop calendars to GAEZ+ 2015 annual data
To generate the monthly cropland physical area of GAEZ+ 2015 crops, we followed these steps for each GAEZ crop in each grid cell:

1.

For a given GAEZ crop in a given grid cell, is the area reported >0 for the matching MIRCA2000 crop?

a.

If YES, then use the MIRCA2000 data for the grid cell and crop considered.

b.

If NO, then find the closest grid cell with the matching MIRCA2000 crop category, and apply the MIRCA2000 crop rotation from that grid cell to the given crop/grid cell combination for the following steps.

2.

Does the matching MIRCA2000 crop category (Table 1) have more than 1 subcrop?

a.

If NO, then AG = HG for all months of the cropping season, as defined by the MIRCA2000 crop calendar.

b.

If YES, then for each subcrop category n, apply the ratio of AM,n/HM to HG, then sum the subcrop areas within each month such that:

$${A}_{G}=sum _{n}frac{{A}_{M,n}}{{H}_{M}}{H}_{G}$$

3.

For each month and each grid cell, check if the sum of all crops (irrigated and rainfed) is greater than the 99% of area of the grid cell. We assume that at least 1% of land must be retained as non-cropland for agricultural infrastructure such as roads, buildings, irrigation infrastructure, and other landcovers (e.g. rivers, wetlands).

a.

If NO, then no further processing is done.

b.

If YES, then reduce crop area by the excess value based on a removal order (Table 2). Rainfed crops have higher removal order numbers for the excess truncation (starting with 1) before removing irrigated crops, until the cell area is not exceeded. A large removal number (e.g., 20) indicates that the crop’s land is unlikely to be removed. Large priority numbers are given to the staple crops to ensure these important food producing lands are consistent with FAOSTAT country data.

The maximum monthly amount of physical cropland that was removed by step 3 is 711,543 ha, which is 0.05% of total global cropland physical area.The resulting global gridded data from Step 6 are monthly time series of cropland physical area by crop, subcrop, and production system, called GAEZ+_2015 Monthly Cropland Data17. Combining the MIRCA2000 crop calendar and subcrop rotation information with the GAEZ+ 2015 annual data allows for the representation of crop seasonality; e.g., Fig. 2 shows the aggregate monthly cropland physical area for Rice 1 and Rice 2 (two sub-crops of rice) over the northern hemisphere, clearly illustrating the two main rice-growing seasons.Fig. 2Aggregate monthly cropland physical area for Rice 1 and Rice 2 subcrops from monthly GAEZ+ 2015 over the northern hemisphere shows the two main rice-growing seasons. This seasonality is the result of combining GAEZ+ 2015 annual data with the MIRCA20004 crop calendars and subcrop divisions.Full size image More

1.Ottesen, E. A. et al. Pattern and synchrony of gene expression among sympatric marine microbial populations. Proc. Natl Acad. Sci. USA 110, E488–E497 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
2.Muñoz-Marín, M. D. C. et al. The transcriptional cycle is suited to daytime N2 fixation in the unicellular cyanobacterium “Candidatus Atelocyanobacterium thalassa” (UCYN-A). mBio 10, e02495-18 (2019).PubMed 
PubMed Central 

Google Scholar 
3.Vislova, A., Sosa, O. A., Eppley, J. M., Romano, A. E. & DeLong, E. F. Diel oscillation of microbial gene transcripts declines with depth in oligotrophic ocean waters. Front. Microbiol. 10, 2191 (2019).PubMed 
PubMed Central 

Google Scholar 
4.Harke, M. J. et al. Periodic and coordinated gene expression between a diazotroph and its diatom host. ISME J. 13, 118–131 (2019).CAS 
PubMed 

Google Scholar 
5.Hernández Limón, M. D. et al. Transcriptional patterns of Emiliania huxleyi in the North Pacific Subtropical Gyre reveal the daily rhythms of its metabolic potential.Environ. Microbiol. 22, 381–396 (2020).PubMed 

Google Scholar 
6.Becker, K. W. et al. Daily changes in phytoplankton lipidomes reveal mechanisms of energy storage in the open ocean. Nat. Commun. 9, 5179 (2018).PubMed 
PubMed Central 

Google Scholar 
7.Frischkorn, K. R., Haley, S. T. & Dyhrman, S. T. Coordinated gene expression between Trichodesmium and its microbiome over day–night cycles in the North Pacific Subtropical Gyre. ISME J. 12, 997–1007 (2018).PubMed 
PubMed Central 

Google Scholar 
8.Ottesen, E. A. et al. Ocean microbes. Multispecies diel transcriptional oscillations in open ocean heterotrophic bacterial assemblages. Science 345, 207–212 (2014).CAS 
PubMed 

Google Scholar 
9.Wilson, S. T. et al. Coordinated regulation of growth, activity and transcription in natural populations of the unicellular nitrogen-fixing cyanobacterium Crocosphaera. Nat. Microbiol. 2, 17118 (2017).CAS 
PubMed 

Google Scholar 
10.Saito, M. A. et al. Iron conservation by reduction of metalloenzyme inventories in the marine diazotroph Crocosphaera watsonii. Proc. Natl Acad. Sci. USA 108, 2184–2189 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
11.Strenkert, D. et al. Multiomics resolution of molecular events during a day in the life of Chlamydomonas. Proc. Natl Acad. Sci. USA 116, 2374–2383 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
12.Boysen, A. K. et al. Particulate metabolites and transcripts reflect diel oscillations of microbial activity in the surface ocean. mSystems 6, e00896-20 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
13.White, A. E., Barone, B., Letelier, R. M. & Karl, D. M. Productivity diagnosed from the diel cycle of particulate carbon in the North Pacific Subtropical Gyre: optically derived productivity. Geophys. Res. Lett. 44, 3752–3760 (2017).CAS 

Google Scholar 
14.DeLong, E. F. et al. Community genomics among stratified microbial assemblages in the ocean’s interior. Science 311, 496–503 (2006).CAS 
PubMed 

Google Scholar 
15.Sunagawa, S. et al. Ocean plankton. Structure and function of the global ocean microbiome. Science 348, 1261359 (2015).PubMed 

Google Scholar 
16.Coles, V. J. et al. Ocean biogeochemistry modeled with emergent trait-based genomics. Science 358, 1149–1154 (2017).CAS 
PubMed 

Google Scholar 
17.Walbauer, J. R., Rodrigue, S., Coleman, M. L. & Chisholm, S. W. Transcriptome and proteome dynamics of a light–dark synchronized bacterial cell cycle.PLoS ONE 7, e43432 (2012).
Google Scholar 
18.Steiner, P. A. et al. Highly variable mRNA half-life time within marine bacterial taxa and functional genes. Environ. Microbiol. 21, 3873–3884 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
19.Moran, M. A. et al. Sizing up metatranscriptomics. ISME J. 7, 237–243 (2013).CAS 
PubMed 

Google Scholar 
20.Tamames, J., Cobo-Simón, M. & Puente-Sánchez, F. Assessing the performance of different approaches for functional and taxonomic annotation of metagenomes. BMC Genomics 20, 960 (2019).21.DiTullio, G. R. & Laws, E. A. Diel periodicity of nitrogen and carbon assimilation in five species of marine phytoplankton: accuracy of methodology for predicting N-assimilation rates and N/C composition ratios. Mar. Ecol. Prog. Ser. 32, 123–132 (1986).CAS 

Google Scholar 
22.Granum, E., Kirkvold, S. & Myklestad, S. M. Cellular and extracellular production of carbohydrates and amino acids by the marine diatom Skeletonema costatum: diel variations and effects of N depletion. Mar. Ecol. Prog. Ser. 242, 83–94 (2002).CAS 

Google Scholar 
23.Lacour, T., Sciandra, A., Talec, A., Mayzaud, P. & Bernard, O. Diel variations of carbohydrates and neutral lipids in nitrogen-sufficient and nitrogen-starved cyclostat cultures of Isochrysis sp. J. Phycol. 48, 966–975 (2012).PubMed 

Google Scholar 
24.Follett, C. L., Dutkiewicz, S., Karl, D. M., Inomura, K. & Follows, M. J. Seasonal resource conditions favor a summertime increase in North Pacific diatom–diazotroph associations. ISME J. 12, 1543–1557 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
25.Chen, W.-N. U. et al. Diel rhythmicity of lipid-body formation in a coral-Symbiodinium endosymbiosis. Coral Reefs 31, 521–534 (2012).
Google Scholar 
26.Zhou, X. & Mopper, K. Photochemical production of low-molecular-weight carbonyl compounds in seawater and surface microlayer and their air-sea exchange. Mar. Chem. 56, 201–213 (1997).CAS 

Google Scholar 
27.Durham, B. P. et al. Sulfonate-based networks between eukaryotic phytoplankton and heterotrophic bacteria in the surface ocean.Nat. Microbiol. 4, 1706–1715 (2019).CAS 
PubMed 

Google Scholar 
28.Lambert, S. et al. Rhythmicity of coastal marine picoeukaryotes, bacteria and archaea despite irregular environmental perturbations. ISME J. 13, 388–401 (2019).PubMed 

Google Scholar 
29.Kolody, B. C. et al. Diel transcriptional response of a California Current plankton microbiome to light, low iron, and enduring viral infection. ISME J. 13, 2817–2833 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
30.Aylward, F. O. et al. Microbial community transcriptional networks are conserved in three domains at ocean basin scales. Proc. Natl Acad. Sci. USA 112, 5443–5448 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
31.Rusch, D. B. et al. The Sorcerer II Global Ocean Sampling expedition: northwest Atlantic through eastern tropical Pacific. PLoS Biol. 5, e77 (2007).PubMed 
PubMed Central 

Google Scholar 
32.Bork, P. et al. Tara Oceans studies plankton at planetary scale. Science 348, 873 (2015).CAS 
PubMed 

Google Scholar 
33.Delmont, T. O. et al. Nitrogen-fixing populations of Planctomycetes and Proteobacteria are abundant in surface ocean metagenomes. Nat. Microbiol. 3, 804–813 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
34.Fuhrman, J. A. et al. Annually reoccurring bacterial communities are predictable from ocean conditions. Proc. Natl Acad. Sci. USA 103, 13104–13109 (2006).CAS 
PubMed 
PubMed Central 

Google Scholar 
35.Morris, R. M. et al. Temporal and spatial response of bacterioplankton lineages to annual convective overturn at the Bermuda Atlantic Time‐series Study site. Limnol. Oceanogr. 50, 1687–1696 (2005).CAS 

Google Scholar 
36.Mende, D. R. et al. Environmental drivers of a microbial genomic transition zone in the ocean’s interior. Nat. Microbiol. 2, 1367–1373 (2017).CAS 
PubMed 

Google Scholar 
37.Keeling, P. J. et al. The Marine Microbial Eukaryote Transcriptome Sequencing Project (MMETSP): illuminating the functional diversity of eukaryotic life in the oceans through transcriptome sequencing. PLoS Biol. 12, e1001889 (2014).PubMed 
PubMed Central 

Google Scholar 
38.Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M. & Tanabe, M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 44, D457–D462 (2016).CAS 

Google Scholar 
39.Thaben, P. F. & Westermark, P. O. Detecting rhythms in time series with RAIN. J. Biol. Rhythms 29, 391–400 (2014).PubMed 
PubMed Central 

Google Scholar 
40.Cuhel, R. L., Ortner, P. B. & Lean, D. R. S. Night synthesis of protein by algae. Limnol. Oceanogr. 29, 731–744 (1984).CAS 

Google Scholar 
41.Coesel, S. N. et al. Diel transcriptional oscillations of light-sensitive regulatory elements in open-ocean eukaryotic plankton communities. Proc. Natl Acad. Sci. USA 118, e2011038118 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
42.Bolay, P., Muro-Pastor, M. I., Florencio, F. J. & Klähn, S. The distinctive regulation of cyanobacterial glutamine synthetase. Life (Basel) 8, 52 (2018).CAS 

Google Scholar 
43.Karl, D. M., Church, M. J., Dore, J. E., Letelier, R. M. & Mahaffey, C. Predictable and efficient carbon sequestration in the North Pacific Ocean supported by symbiotic nitrogen fixation. Proc. Natl Acad. Sci. USA 109, 1842–1849 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
44.Berman, T. & Bronk, D. A. Dissolved organic nitrogen: a dynamic participant in aquatic ecosystems. Aquat. Microb. Ecol. 31, 279–305 (2003).
Google Scholar 
45.Lee, C. & Bada, J. L. Amino acids in equatorial Pacific Ocean water. Earth Planet. Sci. Lett. 26, 61–68 (1975).CAS 

Google Scholar 
46.Bada, J. L. & Lee, C. Decomposition and alteration of organic compounds dissolved in seawater. Mar. Chem. 5, 523–534 (1977).CAS 

Google Scholar 
47.Poretsky, R. S., Sun, S., Mou, X. & Moran, M. A. Transporter genes expressed by coastal bacterioplankton in response to dissolved organic carbon. Environ. Microbiol. 12, 616–627 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
48.Berthelot, H. et al. NanoSIMS single cell analyses reveal the contrasting nitrogen sources for small phytoplankton. ISME J. 13, 651–662 (2019).CAS 
PubMed 

Google Scholar 
49.Moore, L. R., Post, A. F., Rocap, G. & Chisholm, S. W. Utilization of different nitrogen sources by the marine cyanobacteria Prochlorococcus and Synechococcus. Limnol. Oceanogr. 47, 989–996 (2002).CAS 

Google Scholar 
50.Hu, S. K., Connell, P. E., Mesrop, L. Y. & Caron, D. A. A hard day’s night: diel shifts in microbial eukaryotic activity in the North Pacific Subtropical Gyre. Front. Mar. Sci. https://doi.org/10.3389/fmars.2018.00351 (2018).51.Hannides, C. C. S., Popp, B. N., Choy, C. A. & Drazen, J. C. Midwater zooplankton and suspended particle dynamics in the North Pacific Subtropical Gyre: a stable isotope perspective. Limnol. Oceanogr. 58, 1931–1946 (2013).CAS 

Google Scholar 
52.Becker, K. W. et al. Combined pigment and metatranscriptomic analysis reveals highly synchronized diel patterns of phenotypic light response across domains in the open oligotrophic ocean.ISME J. 15, 520–533 (2021).CAS 
PubMed 

Google Scholar 
53.Mruwat, N. et al. A single-cell polony method reveals low levels of infected Prochlorococcus in oligotrophic waters despite high cyanophage abundances. ISME J. 15, 41–54 (2021).CAS 
PubMed 

Google Scholar 
54.Chesson, P. L. & Warner, R. R. Environmental variability promotes coexistence in lottery competitive systems. Am. Nat. 117, 923–943 (1981).
Google Scholar 
55.Shmida, A. & Ellner, S. Coexistence of plant species with similar niches. Vegetatio 58, 29–55 (1984).
Google Scholar 
56.Ellner, S. P., Snyder, R. E. & Adler, P. B. How to quantify the temporal storage effect using simulations instead of math. Ecol. Lett. 19, 1333–1342 (2016).PubMed 

Google Scholar 
57.Adler, P. B., Fajardo, A., Kleinhesselink, A. R. & Kraft, N. J. B. Trait-based tests of coexistence mechanisms. Ecol. Lett. 16, 1294–1306 (2013).PubMed 

Google Scholar 
58.Adler, P. B., HilleRisLambers, J., Kyriakidis, P. C., Guan, Q. & Levine, J. M. Climate variability has a stabilizing effect on the coexistence of prairie grasses. Proc. Natl Acad. Sci. USA 103, 12793–12798 (2006).CAS 
PubMed 
PubMed Central 

Google Scholar 
59.Cáceres, C. E. Temporal variation, dormancy, and coexistence: a field test of the storage effect. Proc. Natl Acad. Sci. USA 94, 9171–9175 (1997).PubMed 
PubMed Central 

Google Scholar 
60.Padisák, J. Identification of relevant time-scales in non-equilibrium community dynamics: conclusions from phytoplankton surveys. N. Z. J. Ecol. 18, 169–176 (1994).
Google Scholar 
61.Anderies, J. M. & Beisner, B. E. Fluctuating environments and phytoplankton community structure: a stochastic model. Am. Nat.155, 556–569 (2000).PubMed 

Google Scholar 
62.Wagg, C. et al. Functional trait dissimilarity drives both species complementarity and competitive disparity. Funct. Ecol. 31, 2320–2329 (2017).
Google Scholar 
63.Bligh, E.G. & Dyer, W. J. A rapid method of total lipid extraction and purification. Can. J. Biochem. Physiol. 37, 911–917 (1959).CAS 
PubMed 

Google Scholar 
64.Boysen, A. K., Heal, K. R., Carlson, L. T. & Ingalls, A. E. Best-matched internal standard normalization in liquid chromatography–mass spectrometry metabolomics applied to environmental samples. Anal. Chem. 90, 1363–1369 (2018).CAS 
PubMed 

Google Scholar 
65.MacLean, B. et al. Skyline: an open source document editor for creating and analyzing targeted proteomics experiments. Bioinformatics 26, 966–968 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
66.Fountoulakis, M. & Lahm, H. W. Hydrolysis and amino acid composition analysis of proteins. J. Chromatogr. A 826, 109–134 (1998).CAS 
PubMed 

Google Scholar 
67.Popendorf, K. J., Fredricks, H. F. & Van Mooy, B. A. S. Molecular ion-independent quantification of polar glycerolipid classes in marine plankton using triple quadrupole MS. Lipids 48, 185–195 (2013).CAS 
PubMed 

Google Scholar 
68.Collins, J. R., Edwards, B. R., Fredricks, H. F. & Van Mooy, B. A. S. LOBSTAHS: an adduct-based lipidomics strategy for discovery and identification of oxidative stress biomarkers. Anal. Chem. 88, 7154–7162 (2016).CAS 
PubMed 

Google Scholar 
69.Hummel, J. et al. Ultra performance liquid chromatography and high resolution mass spectrometry for the analysis of plant lipids. Front. Plant Sci. 2, 54 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
70.Smith, C. A., Want, E. J., O’Maille, G., Abagyan, R. & Siuzdak, G. XCMS: processing mass spectrometry data for metabolite profiling using nonlinear peak alignment, matching, and identification. Anal. Chem. 78, 779–787 (2006).CAS 
PubMed 

Google Scholar 
71.Kuhl, C., Tautenhahn, R., Böttcher, C., Larson, T. R. & Neumann, S. CAMERA: an integrated strategy for compound spectra extraction and annotation of liquid chromatography/mass spectrometry data sets. Anal. Chem. 84, 283–289 (2012).CAS 
PubMed 

Google Scholar 
72.Biller, S. J. et al. Prochlorococcus extracellular vesicles: molecular composition and adsorption to diverse microbes.Environ. Microbiol. https://doi.org/10.1111/1462-2920.15834 (2021).Article 
PubMed 

Google Scholar 
73.Aylward, F. O. et al. Diel cycling and long-term persistence of viruses in the ocean’s euphotic zone. Proc. Natl Acad. Sci. USA 114, 11446–11451 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
74.Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
75.Masella, A. P., Bartram, A. K., Truszkowski, J. M., Brown, D. G. & Neufeld, J. D. PANDAseq: paired-end assembler for illumina sequences. BMC Bioinformatics 13, 31 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
76.Joshi, N. & Fass, J. Sickle: A sliding-window, adaptive, quality-based trimming tool for FastQ files. Version 1.33. GitHub https://github.com/najoshi/sickle (2015).77.Kopylova, E., Noé, L. & Touzet, H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics 28, 3211–3217 (2012).CAS 

Google Scholar 
78.Kiełbasa, S. M., Wan, R., Sato, K., Horton, P. & Frith, M. C. Adaptive seeds tame genomic sequence comparison. Genome Res. 21, 487–493 (2011).PubMed 
PubMed Central 

Google Scholar 
79.Li, H. & Durbin, R. Fast and accurate long-read alignment with Burrows–Wheeler transform. Bioinformatics 26, 589–595 (2010).PubMed 
PubMed Central 

Google Scholar 
80.Alexander, H. et al. Functional group-specific traits drive phytoplankton dynamics in the oligotrophic ocean. Proc. Natl Acad. Sci. USA 112, E5972–E5979 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
81.Anders, S., Pyl, P. T. & Huber, W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169 (2015).CAS 

Google Scholar 
82.Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).PubMed 
PubMed Central 

Google Scholar 
83.Meinicke, P. UProC: tools for ultra-fast protein domain classification. Bioinformatics 31, 1382–1388 (2015).CAS 
PubMed 

Google Scholar 
84.Mende, D. R., Boeuf, D. & DeLong, E. F. Persistent core populations shape the microbiome throughout the water column in the North Pacific Subtropical Gyre. Front. Microbiol. 10, 2273 (2019).PubMed 
PubMed Central 

Google Scholar 
85.White, A. E. et al. Phenology of particle size distributions and primary productivity in the North Pacific subtropical gyre (Station ALOHA). J. Geophys. Res. Oceans 120, 7381–7399 (2015).PubMed 
PubMed Central 

Google Scholar 
86.Borchers, H. W. pracma: Practical numerical math functions. R package version 2 https://cran.r-project.org/web/packages/pracma/index.html (2019).87.Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M. & Hornik, K. cluster: Cluster analysis basics and extensions. R package version 1.56 (2012).88.Wehrens, R. & Buydens, L. M. C. Self- and super-organizing maps in R: the Kohonen package. J. Stat. Softw. 21, 1–19 (2007).
Google Scholar 
89.Hennig, C. fpc: Flexible procedures for clustering. R package version 2.2-9 (2010).90.Muratore, D. Code for complex marine microbial communities partition metabolism of scarce resources over the diel cycle. Zenodo https://doi.org/10.5281/zenodo.3817416 (2020). More

How many parameters best describe data in muon spectroscopy?Here we find that the Bayes factor demands the inclusion of more physically-meaningful parameters than the BIC or significance tests. Figure 1a presents some data that might reasonably be fitted with as few as three or as many as 22 physically-meaningful parameters. We find that the Bayes factor encourages the inclusion of all these parameters until the onset of over-fitting. Even though many of them have fitted values that fail significance tests (i.e. are consistent with zero), their omission distorts the fitting results severely.Figure 1Full size imageFigure 1a shows an anti-level-crossing spectrum observed in photo-excited muon-spin spectroscopy26 from an organic molecule27. The data are presented in Fig. 2a of Ref.27 and are given in the SI. These spectra are expected to be Lorentzian peaks. Theory permits optical excitation to affect the peak position, the width and the strength (photosensitivity). In the field region over which the measurements are carried out, there is a background from detection of positrons, which has been subtracted from the data presented27. Wang et al.27 did not attempt to fit the data rigorously; they did report a model-independent integration of the data, which demonstrated a change in area and position.The model that we fit hypothesises one or more Lorentzian peaks, with optional photosensitivity on each fitting parameter and with optional linear backgrounds y = a + bx underlying the peaks, described by the full equation given in the SI, equation (S3). To do a single LS fit to all the data, we extend the data to three dimensions, (x gauss, y asymmetry, z) where z = 0 for data in the dark and z = 1 for photoexcited data. Including all the data in a single LS fit in this way, rather than fitting the dark and photoexcited data separately, simplifies both setting up the fit and doing the subsequent analysis.Figure 1b shows the evolution of the SBIC and the lnBF as the number of fitting parameters in the model is increased. Starting with a single Lorentzian peak, three parameters are required, peak position P, width W and intensity A. Three photosensitivity parameters ΔLP, ΔLW and ΔLA are then introduced successively to the fit, (open and small data points for n = 3–6). The SBIC decreases and the lnMLI scarcely increases. It is only with the inclusion of one background term (n = 7) that any figure of merit shows any substantial increase. There is no evidence here for photosensitivity. The weak peak around 7050 G does not seem worth including in a fit, as it is evidenced by only two or three data points and is scarcely outside the error bars. However, a good fit with two peaks (P1 ~ 7210 G, P2 ~ 7150 G, the subscripts 1 and 2 in accordance with the site labelling of Fig. 2a of Ref.27) can be obtained with just five parameters (P1, P2, A1, A2, W). This gives substantial increases in the SBIC and lnMLI, further increased when W1 and W2 are distinguished and then when the single background term and the three photosensitivity parameters ΔLP2, ΔLW2 and ΔLA2 are successively included (solid or large data points for n = 5–10 in Fig. 1b). The SBIC reaches its maximum here, at n = 10, and then decreases substantially when the other three photosensitivity parameters and the other three background terms are included. These additional parameters fail significance tests as well as decreasing the SBIC (Fig. 1b). Conventionally, the n = 10 fit would be accepted as best. The outcome would be reported as two peaks, with significant photo-sensitivities ΔLP2, ΔLW2 and ΔLA2 for all three of the 7150 G peak parameters, but no photosensitivity for the 7210 G peak (Table 1).Table 1 Photosensitivity results of fitting the data of Fig. 1a with 10, 16 and 19 parameters. Parameter units as implied by Fig. 1a.Full size tableThe Bayes factor gives a very different outcome. From 10 to 16 parameters, the Bayes factor between any two of these seven models is close to unity (Fig. 1b). That is, they have approximately equal probability. The Bayes factor shows that what the conventional n = 10 analysis would report is false. Specifically, it is not the case that ΔLP2, reported as − 14 ± 4 G, has a roughly ({raise0.5exhbox{$scriptstyle 2$} kern-0.1em/kern-0.15em lower0.25exhbox{$scriptstyle 3$}}) probability of lying between − 10 and − 18 G. That is not consistent with the roughly equal probability that it lies in the n = 16 range (− 24 ± 8 G). Table 1 shows that at n = 16, ΔLP2 is the only photosensitivity parameter to pass significance tests. ΔLA2, which had the highest significance level at n = 10, is now the parameter most consistent with zero. The other four are suggestively (about 1({raise0.5exhbox{$scriptstyle 1$} kern-0.1em/kern-0.15em lower0.25exhbox{$scriptstyle 2$}})σ) different from zero.Since the Bayes factor has already radically changed the outcome by encouraging more physically-meaningful parameters, it is appropriate to try the 7050 G peak parameters in the fit. With only 28 data-points, we should be alert to over-fitting. We can include P3 and A3 (n = 18), and ΔLP3 (n = 19), but W3 and ΔLA3 do cause overfitting. Figure 1b shows substantial increases of both the SBIC and the lnMLI for n = 18 to n = 20, where the twentieth parameter is in fact ΔLA3. The symptom of over-fitting that we observe here is an increase in the logarithm of the Occam Factor (lnMLI − lnL), the values of which decrease, − 26.9, − 33.5, − 34.8, and then increase, − 33.4, for n = 16, 18, 19 and 20 respectively. Just as lnL must increase with every additional parameter, so should the Occam factor decrease, as the prior parameter volume should increase more with a new parameter than the posterior parameter volume. So we stop at n = 19. The outcome, Table 1, is that the uncertainties on the n = 16 parameters have decreased markedly. This is due to the better fit, with a substantial increase in lnL corresponding to reduced residuals on all the data. The 7210 G peak 2 now has photosensitivities on all its parameters, significant to at least the 2σ or p value ~ 0.05 level. And the photosensitivities ΔLW2 and ΔLA2, both so significant at n = 10, and already dwindling in significance at n = 16, are both now taking values quite consistent with zero. In the light of Table 1, we see that stopping the fit at n = 10 results in completely incorrect results—misleading fitted values, with certainly false uncertainties.Discriminating between models for the pressure dependence of the GaAs bandgapThe main purpose of this example is to show how the Bayes factor can be used to decide between two models which have equal goodness of fit to the data (equal values of lnL and BIC, as well as p values, etc.). This illustrates the distinction it makes between physically-meaningful and physically meaningless parameters. This example also shows how ML fitting can be used together with the Bayes factor to obtain better results. For details, see SI §7.Figure 2 shows two datasets for the pressure dependence of the bandgap of GaAs (data given in the SI). The original authors published quadratic fits, ({E}_{g}(P)={E}_{0}+bP+c{P}^{2}), with b = 10.8 ± 0.3 meV kbar−1 (Goñi et al.28) and 11.6 ± 0.2 meV kbar−1 (Perlin et al.29). Other reported experimental and calculated values for b ranged from 10.02 to 12.3 meV kbar−130. These discrepancies of about ± 10% were attributed to experimental errors in high-pressure experimentation. However, from a comparison of six such datasets, Frogley et al.30 were able to show that the discrepancies arose from fitting the data with the quadratic formula. The different datasets were reconciled by using the Murnaghan equation of state and supposing the band-gap to vary linearly with the density (see SI, §7, equations (S4) and (S5)30. The curvature c of the quadratic is constant, while the curvature of the density, due to the pressure dependence Bʹ of the bulk modulus B0, decreases with pressure—and the six datasets were recorded over very different pressure ranges, as in Fig. 2. So the fitted values of c, c0, were very different, and the correlation between b and c resulted in the variations in b0.Here, using the Bayes factor, we obtain the same result from a single dataset, that of Goñi et al.28 The two fits are shown in Fig. 2. They are equally good, with values of lnL and SBIC the same to 0.01. The key curvature parameters, c and ({text{B}}^{prime }), are both returned as non-zero by 13.5σ (SI, §7, Table S1), consequently both with p-values less than 10−18. However, c is a physically-meaningless parameter. The tightest constraint we have for setting its range is the values previously reported, ranging from 0 to 60 μeV kbar−2, so we use Δc = 100 μeV kbar−2. In contrast, ({text{B}}^{prime }) is known for GaAs to be 4.4931. For many other materials and from theory the range 4–5 is expected, so we use (Delta {text{B}}^{prime } = 1). The other ranges are same for both models (see SI §7). This difference gives a lnBF of 3.8 in favour of the Murnaghan model against the quadratic, which is strong evidence for it. Moreover, the value of ({text{B}}^{prime }) returned is 4.47 ± 0.33, in excellent agreement with the literature value. Had it been far out of range, the model would have to be rejected. The quadratic model is under no such constraint; indeed, a poor fit might be handled by adding cubic and higher terms ad lib. This justifies adding about 5 to lnBF (see “Background in fitting a carbon nanotube Raman spectrum” section), giving a decisive preference to the Murnaghan model, and the value of b it returns, 11.6 ± 0.3. Note the good agreement with the value from Perlin et al.29 If additionally we fix ({mathrm{B}}^{prime}) at its literature value of 4.4931, lnBF is scarcely improved, because the Occam factor against this parameter is small, but the uncertainty on the pressure coefficient, Ξ/B0, is much improved.When we fit the Perlin data, the Murnaghan fit returns ({text{B}}^{prime }) = 6.6 ± 2.4. This is outside range, and indicates that this data cannot give a reliable value—attempting it is over-fitting. However, it is good to fit this data together with the Goñi data. The Perlin data, very precise but at low pressures only, complement the Goñi data with their lower precision but large pressure range. We notice also that the Perlin data has a proportion of outlier data points. Weighted or rescaled LS fitting can handle the different precisions, but it cannot handle the outliers satisfactorily. Maximum Likelihood fitting handles both issues. We construct lnL using different pdfs P(r) for the two datasets, and with a double-Gaussian pdf for the Perlin data (see equation (S6) in the SI §7). Fixing ({text{B}}^{prime }) at 4.49, fitting with the same Ξ/B0 returns 11.42 ± 0.04 meV kbar−1. Separate Ξ/B0 parameters for the two datasets give an increase of lnL of 4.6, with values 11.28 ± 0.06 and 11.60 ± 0.04 meV kbar−1—a difference in b of 0.32 ± 0.07 meV kbar−1, which is significant at 4½σ. This difference could be due to systematic error, e.g. in pressure calibration. Or it could be real. Goñi et al.28 used absorption spectroscopy to measure the band-gap; Perlin et al.29 used photoluminescence. The increase of the electron effective mass with pressure might give rise to the difference. In any case, it is clear that high-pressure experimentation is much more accurate than previously thought, and that ML fitting exploits the information in the data much better than LS fitting.Figure 2GaAs band-gap. Data for Eg(P) in GaAs from Goñi et al.28 (
) and from Perlin et al.29 (
) are shown after subtraction of the straight line E0 + 8.5P to make the curvature more visible. The Perlin data is expanded × 10 on both axes for clarity. Two least-squares fits to the Goñi data are shown, polynomial (dashed red line) and Murnaghan (solid blue line). (Figure prepared using Mathematica 12.0, www.wolfram.com/mathematica/).Full size imageBackground in fitting a carbon nanotube Raman spectrumThis example demonstrates how the Bayes Factor provides a quantitative answer to the problem, whether we should accept a lower quality of fit to the data if the parameter set is intuitively preferable. It also provides a simple example of a case where the MLI calculated by Eq. (1) is in error and can readily be corrected (see SI §8 Fig. S3).The dataset is a Raman spectrum of the radial breathing modes of a sample of carbon nanotubes under pressure32. The whole spectrum at several pressures is shown with fits in Fig. 1 of Ref.32. The traditional fitting procedure used there was to include Lorentzian peaks for the clear peaks in the spectra, and then to add broad peaks as required to get a good fit, but without quantitative figures of merit and without any attempt to explain the origin of the broad peaks, and therefore with no constraints on their position, widths or intensities. The key issue in the fitting was to get the intensities of the peaks as accurately as possible, to help understand their evolution with pressure. Here, we take a part of the spectrum recorded at 0.23 GPa (the data is given in the SI.) and we monitor the quality of fit and the Bayes factor while parameters are added in four models. This part of the spectrum has seven sharp pseudo-Voigt peaks (Fig. 3a; the two strong peaks are clearly doublets). With seven peak positions Pi, peak widths Wi and peak intensities Ai, and a factor describing the Gaussian content in the pseudo-Voigt peak shape, there are already 22 parameters (for details, see SI §8). This gives a visibly very poor fit, with lnL = − 440, SBIC = − 510 and lnMLI = − 546. The ranges chosen for these parameters for calculating the MLI (see SI §8) are not important because they are used in all the subsequent models, and so they cancel out in the Bayes factors between the models.Figure 3Carbon nanotube Raman spectrum. In (a), the carbon nanotube Raman spectrum is plotted (black datapoints) with a fit (cyan solid line) using the Fourier model. The residuals for four good fits are shown, × 10 and displaced successively downwards (Fourier, Polynomial, Peaks and Tails; all at lnL about − 60, see text). The backgrounds are shown, × 8 (long dashed, chain-dotted, short dashed and solid, respectively. In (b), the evolution of the MLIs is shown against the number of parameters for these four models. (Figure prepared using Mathematica 12.0, www.wolfram.com/mathematica/).Full size imageTo improve the fit, in the Fourier model we add a Fourier background (y=sum {c}_{i}mathrm{cos}ix+{s}_{i}mathrm{sin}ix) (i = 0,..) and in the Polynomial model, we add (y=sum {a}_{i}{x}^{i}) (i = 0,..) for the background. In both, the variable x is centred (x = 0) at the centre of the fitted spectrum and scaled to be ± π or ± 1 at the ends. In the Peaks model we add extra broad peaks as background, invoking extra parameter triplets (Pi, Wi, Ai). These three models all gave good fits; at the stage shown in Fig. 3a they gave lnL values of − 65, − 54 and − 51 and BIC values of − 156, − 153 and − 148 respectively. Thus there is not much to choose between the three models, but it is noteworthy that they give quite different values for the intensities of the weaker peaks, with the peak at 265 cm−1 at 20.5 ± 1.1, 25.5 ± 1.3 and 27 ± 1.7 respectively (this is related to the curvature of the background function under the peak). So it is important to choose wisely.A fourth model was motivated by the observation that the three backgrounds look as if they are related to the sharp peaks, rather like heavily broadened replicas (see Fig. 3a). Accordingly, in the fourth model, we use no background apart from the zeroth term c0 or a0 to account for dark current). Instead, the peak shape is modified, giving it stronger, fatter tails than the pseudo-Voigt peaks (Tails model). This was done by adding to the Lorentzian peak function a smooth function approximating to exponential tails on both sides of the peak position (for details, see SI §8) with widths and amplitudes as fitting parameters. What is added may be considered as background and is shown in Fig. 3a. This model, at the stage of Fig. 3a, returned lnL = − 62, BIC = − 146, and yet another, much smaller value of 15.5 ± 1.0 for the intensity of the 265 cm−1 peak.The Tails model is intuitively preferable to the other three because it does not span the data space—e.g. if there was really were broad peaks at the positions identified by the Peaks model, or elsewhere, the Tails model could not fit them well. That it does fit the data is intuitively strong evidence for its correctness. The Bayes factor confirms this intuition quantitatively. At the stage of Fig. 3a, the lnMLI values are − 251, − 237 and − 223 for the Fourier, Poly and Peaks models, and − 211 for the Tails model. This gives a lnBF value of 12 for the Tails model over the Peaks model—decisive—and still larger lnBF values for these models over the Fourier and Poly models.All models can be taken further, with more fitting parameters. More Fourier or polynomial terms or more peaks can be added, and for the Tails model more parameters distinguishing the tails attached to each of the seven Lorentizian peaks. In this way, the three background models can improve to a lnL ~ − 20; the Tails model does not improve above lnL ~ − 50. However, as seen in Fig. 3b, the MLIs get worse with too many parameters, except when over-fitting occurs, as seen for the Poly model at 35 parameters. The Tails model retains its positive lnBF  > 10 over the other models.The other models can have an indefinite number of additional parameters—more coefficients or more peaks, to fit any data set. It is in this sense that they span the data space. The actual number used is therefore itself a fitting parameter, with an uncertainty perhaps of the order of ± 1, and a range from 0 to perhaps a quarter or a half of the number of data points m. We may therefore penalise their lnMLIs by ~ ln 4 m−1 or about − 5 for a few hundred data points. This takes Tails to a lnBF  > 15 over the other models—overwhelmingly decisive. This quantifies the intuition that a model that is not guaranteed to fit the data, but which does, is preferable to a model that certainly can fit the data because it spans the data space. It quantifies the question, how much worse a quality of fit should we accept for a model that is intuitively more satisfying. Here we accept a loss of − 30 on lnL for a greater gain of + 45 in the Occam factor. It quantifies the argument that the Tails model is the most worthy of further investigation because the fat tails probably have a physical interpretation worth seeking. In this context, it is interesting that in Fig. 3a fat tails have been added only to the 250, 265 and 299 cm−1 peaks; adding fat tails to the others did not improve the fit; however, a full analysis and interpretation is outside the scope of this paper. In the Peaks model it is not probable (though possible) that the extra peaks would have physical meaning. In the other two models it is certainly not the case that their Fourier or polynomial coefficients will have physical meaning. More

Study area and datasetThe study area covers all Estonia located between 57.5(^circ ) N, 21.5(^circ ) E and 59.8(^circ ) N, 28.2(^circ ) E. The study area is relatively flat with no steep slopes and altitudes ranging between 0 and 200m above the sea level. Data about events were collected directly from field books that contained information about the mowing activity’s start and end date and the covered area. Considering the main agricultural areas of the country, we consider 2000 fields in which events are geographically evenly distributed across all Estonia, as shown in Fig. 1. In total, data about 1800 mowing and 200 non-mown events were collected in 2018, based on manual labelling. During manual labelling, the specific mowing days were labelled based on the following: a) information recorded by farmers in field books regarding mowing days, b) domain experts knowledge about the most probable days for mowing based on the climate, weather, and field conditions, c) rapid decrease in the Normalized Difference Vegetation Index (NDVI) and rapid increase in the coherence compared to past measurements. The average field size is 6.0ha, and around 95% of the fields were mown during the year. 90% of the fields are in the range of (0.5-10)ha. The greatest density of the fields is located in Lääne-Viru, Tartu and Jõgeva countries. Grassland parcels vector layer is provided by Estonian Agricultural Registers and Information Board (ARIB)50. The satellite imagery used in the study is from Copernicus program that provides free open Earth observation data to help service providers, public authorities, and international organizations improve European citizens’ quality of life.Figure 1Geographic distribution of events used in this study (This map was created by QGIS version 3.16, which can be accessed on https://qgis.org/en/site/).Full size imageSentinel-1 and Sentinel-2 dataFor Sentinel-1 data, in total, 400 S1A/BSLCIW products acquired between 1st of May 20017 and 30th of October 2018, were processed. 87 products were from relative orbit number (RON)160, 62 from RON131, 84 from RON87, 93 from RON58, and 60 from RON29. These were organised into S1A/S1B 6-day pairs. Sentinel-2 provides high spatial resolution optical imagery to perform terrestrial observations with global coverage of the Earth’s land surface. Sentinel-2 data is provided by the European Space Agency (ESA) together with a cloud mask, which can filter clouds on the image with moderately good accuracy. 400 Sentinel-2A and -2B L2A products acquired between 1 May 2017 and 30 October 2018 were processed. Each Sentinel-2 image is a maximum of three days off from the closest Sentinel-1 image. Only the NDVI was derived from Sentinel-2. NDVI has been widely used in the classification of grassland24,51 and that is mainly due to its ability in limiting spectral noise. The spatial resolution of the derived Sentinel-2 NDVI feature is 10 m.MethodsThe goal of the analysis is to detect mowing events from Sentinel-1 (S-1) and Sentinel-2 (S-2) data. For this, coherence time series were calculated about every field in the database about the event. Average coherence of a field, imaging geometry parameters, imaging time and average NDVI were stored in a database. The database formation process involved preprocessing many satellite images where average coherence and NDVI value was calculated for every parcel for every available date (constrained by image availability and cloud cover). The overall scheme of the proposed methodology is illustrated in Fig. 2. First, the time-series data from S-1 and S-2 images are preprocessed. Then, the most important features are used in a deep neural network to predict mowing events. The model has a reject region option that enables the model to abstain from the prediction in case of uncertainty, which increases trust in the model.We used the Sentinel Application Platform (SNAP) toolbox for processing S-1 data. More specifically, we followed the same following pre-processing steps in16: apply orbit file, back-geocoding (using Shuttle Radar Topography Mission (SRTM) data), coherence calculation, deburst, terrain correction, and reprojection to the local projection (EPSG:3301). Lastly, we resampled the data to 4m resolution to preserve the maximum spatial resolution and square-shaped pixels. Because the study areas’ terrain is relatively flat, there are few topographic distortions in the SAR data. Each swath’s coherence was calculated independently. Only pixels totally inside the parcel boundaries (including the average window used for coherence computation) were utilized to calculate results, and any interference beyond the parcel limits was discarded. Pair-wise coherence was calculated with 6-day time step. The data was stored into a database using a forward-looking convention: coherence regarding date X refers to the coherence between S-1 images over the period between date X and X + 6 days. For preprocessing S-2 data, L1C and L2A Sentinel-2 products were obtained through Copernicus Open Access Hub6. Next, a rule-based cloud mask solution was applied52. Finally, the fourth and eighth bands were extracted to compute NDVI values.Figure 2Flowchart of the proposed approach to detect mowing events.Full size imageFeature extraction from Sentinel-1 dataCoherence is a normalized measure of similarity between two consecutive (same relative orbit) S-1 images. Interferometric 6 day repeat pass coherence in VV polarization (cohvv), and coherence in VH polarization (cohvh) are chosen features as they are shown to be sensitive to changes in vegetation and agricultural events25. The shorter the time interval after the mowing event and the first interferometric acquisition, the higher the coherence value. Generally, up to 24 to 36 days after a mowing event, coherence stays relatively high. Precipitation caused the coherence to drop, which disturbs the detection of a mowing event. The spatial resolution of the S-1 6-day repeat pass interferometric coherence is 70 m. Given two S-1 images (s_{1}) and (s_{2}), coherence is calculated as follows:$$begin{aligned} wp =frac{|langle s_{1}s*_{2}rangle |}{sqrt{langle s_{1}s*_{1}rangle | langle s_{2}s*_{2}rangle |}} end{aligned}$$
(1)
where (|langle s_{1}s*_{2}rangle |) is the absolute value of the spatial average of the complex conjugate product.Coherence between two S-1 images (s_1) and (s_2) reaches its maximum value of 1 when both images have the same position and physical characteristics of the scatters. In contrast, the coherence value declines when the position or properties of the scatters change.Feature extraction from Sentinel-2 dataNDVI is related to the amount of live green vegetation. Generally, NDVI increases and decreases over the season, indicating the natural growth decay of vegetation, while the significant drops in the NDVI indicate an agricultural event such as mowing. NDVI is derived from S2 images and is calculated as follows:$$begin{aligned} NDVI=frac{band8 – band4}{band8 + band4} end{aligned}$$
(2)
Figure 3Typical signature of NDVI and coherence in VV and VH polarisation for non mown field during the year.Full size imageFigure 4Field with single mowing event during the year.Full size imageFigure 5NDVI measurement for a field example with a single mowing event during the season.Full size imageFigures 3, 4 and 5 show different samples of mown and non mown fields. NDVI measurements are green, cohvh and cohvv are blue and black, respectively. For non mown field, the typical signature of NDVI during the year is shown in Fig. 3. For non mown field, the typical signature of NDVI during the season is a half-oval curve; coherence is not stable but remains at almost the same level without apparent trend changes, as shown in Fig. 3. An example of a field with a single mowing event during the season is shown in Fig. 4. A mowing event is characterized by a rapid increase in both cohvh and cohvv and a sharp decrease in NDVI, as observed at day 150 (See Fig. 4). Forty days later, a similar signature is probably not due to a mowing event but likely caused by drought during summer.Notably, NDVI measurements are irregular and relatively sparse. Around 75% of total NDVI measurements are invalid in Estonia, and the percentage is slightly lower in Southern Sweden and Denmark due to cloud cover. The Cloud mask indicates the percentage of cloud coverage and allows the cloudy and cloud-free pixels to be identified. Using the standard cloud mask technique by the European Space Agency (ESA) leads to outliers noticed in the sudden decrease in the NDVI. Figure 5 shows an extreme value of NDVI that is supposed to be an outlier due to high differences to the precedent and subsequent values. The outlier is marked with a yellow dot (NDVI=0.38), nearest previous (NDVI=0.75), and next (NDVI=0.78) measurements are marked with a blue colour.Sentinel-1 and Sentinel-2 data preprocessingTo detect NDVI outliers effectively, a good understanding of the data is needed. NDVI outliers due to cloud mask errors rarely co-occur together, and hence, they can be treated as independent events53. NDVI outliers are usually identified with a sudden drop to almost zero and do not form a sequence. It is enough to look at neighbouring measurements (one before and one after) to detect individual outliers. If the difference between the adjacent measurements is high, this is an outlier signature. Hence, outliers can be handled by iterating through every three consecutive NDVI measurements for a given field and checking the difference between the first and second values and between third and second values. Figure 6 shows the scatter plot of all three consecutive NDVI measurements. The Y-axis shows the difference between third and second NDVI values in a triplet, while X-axis represents the difference between second and first NDVI values in a triplet. Triplets with up to 7 days difference are shown in blue, and triplets from 7 to 14 days are shown in green. The points structure forms a rhombus shape with a small cloud of possible outliers in the upper left corner. To filter outliers from the list of actual mowing events, we only consider triplets within up to 10 days interval (as the mowing event signature can recover in 10 days). Knowing rhombus equation (the centre is approximately in (0, 0), and the side length is around 0.6), the filtering rule can be easily applied as follows:$$begin{aligned} ndvi_3 – 2 cdot ndvi_2 + ndvi_1 ge 0.6 end{aligned}$$
(3)
where ndvi_1, ndvi_2, and ndvi_3 are consecutive NDVI measurements within 10 days interval.All outliers are removed, which represent around 0.1% of NDVI measurements.Figure 6Scatter plot of NDVI triplets.Full size imageSmoothing is an essential pre-processing step for noisy features. In this work, cohvh and cohvv features are smoothed using different techniques, including exponential moving average (EMA), moving average54, and Kalman filter55. Smoothing using moving average is done by taking the averages of raw data sequences. The length of the sequence over which we take the average is called the filter width. Table 1 shows the performance of moving average smoothing technique using different values for the filter width. The results show that the best AUC-ROC of 0.9671 is achieved at a filter size of 7. The Kalman filter produces estimates of the current state variables and their uncertainties. Once the outcome of the subsequent measurement is observed, these estimates are updated using a weighted average, giving more weight to estimates with higher certainty. The AUC-ROC achieved using Kalman filter is 0.962. The EMA is done by taking averages of sequences of data, in addition to assigning weights to every data point. More specifically, as values get older, they are given exponentially decreasing weights. The smoothed cohvh and cohvv EMA for cohvh and cohvv are calculated using a recursive definition (i.e., from its previous value) as follows:$$begin{aligned}&cohvh_sm(cohvh_n, alpha ) = alpha cdot (cohvh_n) + (1 – alpha ) cdot cohvh_sm(cohvh_{n-1}, alpha ) end{aligned}$$
(4)
$$begin{aligned}&cohvv_sm(cohvv_n, alpha ) = alpha cdot (cohvv_n) + (1 – alpha ) cdot cohvv_sm(cohvv_{n-1}, alpha ) end{aligned}$$
(5)
where (cohvh_sm(cohvh_{n-1}, alpha )): exponential moving average for end of (cohvh_{n-1}). (cohvv_sm(cohvv_{n-1}, alpha )): exponential moving average for end of (cohvv_{n-1}). (alpha ): a smoothing parameter.The higher the smoothing parameter, the more it reacts to fluctuations in the original signal. The lower the smoothing parameter, the more the signal is smoothed. Experimentally, we found that the best value for (alpha ) to achieve the best AUC-ROC of 0.968 is (frac{1}{3}) as shown in Table 2. The different smoothing techniques achieve comparable performance. EMA technique was selected as it achieves slightly higher performance.Table 1 Performance of moving average smoothing using different filter width.Full size tableTable 2 Performance of EMA smoothing using different values of (alpha ).Full size tableDerived featuresNew derived features from S-1 and S-2 are extracted to improve the performance of the machine learning model. The features were derived based on the following knowledge about mowing events: coherence tends to increase. In contrast, ndvi tends to decrease after mowing events and, many farmers perform mowing during the same time of the year due to the good weather conditions. Such knowledge was elaborated with the derived features. In the following, we will go through the list of derived features considered in this study. Mixed coherence is derived from S-1 features to capture the overall coherence trend. Mixed coherence is a non-linear combination of cohvh and cohvv and is calculated as follows:$$begin{aligned} Mixed_coh = sqrt{cohvh cdot cohvv} end{aligned}$$
(6)
The date is an important feature for the model to adapt, as it is more likely to have mowing events in the summer rather than in early spring, especially in Estonia. The normalized day of the year is calculated as normalization improves the training process of the neural network. Some methods normalize features during the training process, such as Batch Normalization used in this study56. However, neighbouring batches could have entirely different normalization variables (batch mean and variance). At the same time, DOY is a feature susceptible to small changes, e.g., mowing prediction on day 108 or 109 could have drastically different meaning (weekend or working day, day with sunny weather or day with heavy rain). It implies that unified normalization of the DOY feature before training could help avoid the unwanted impact of Batch normalization and possible gradient computation issues. The normalized day of the year is calculated as follows:$$begin{aligned} t = frac{day_of_year}{365} end{aligned}$$
(7)
where (day_of_year) is the year’s day, which is a number between 1 and 365, January 1st is day 1.In addition, we use another time feature dt to capture the gaps in time series. dt is defined to be the normalized difference in days between the current measurement and the previous one. Normalization was performed with min-max scaling. dt is calculated as follows:$$begin{aligned} dt = frac{diff – min_diff}{max_diff – min_diff} end{aligned}$$
(8)
where (min_diff): the minimum difference in days between two previous consecutive measurements obtained from training data. (max_diff): the maximum difference in days between two previous consecutive measurements obtained from training data.Since mowing is characterized by an increase in the coherence and decline in the NDVI, it is important to capture the difference in the values of features and/or slopes of the features’ curves. In the following, we summarize the list of original and derived features extracted from Sentinel-1 and Sentinel-2 included in this study.

ndvi Normalized difference vegetation index, obtained from Sentinel-2.

cohvv Coherence in VV polarization, Sentinel-1 feature.

cohvh Coherence in VH polarization, Sentinel-1 feature.

t Normalized day of the year when the measurement is obtained.

dt Normalized difference in days between current and previous measurement. The data was interpolated with a daily grid, this feature differentiated between interpolated data and real data by capturing the difference between valid (not interpolated) measurements.

cohvv_sm Smoothed cohvv with exponential mowing average (with parameter (frac{1}{3})).

cohvh_sm Smoothed cohvh with exponential moving average (with parameter (frac{1}{3})).

mixed_coh Harmonic mean of cohvv and cohvh. The harmonic mean is chosen as one of the simplest options of non-linear combination.

ndvi_diff Difference between current and previous NDVI measurements. This feature captures the decrease in the ndvi, which is highly related to mowing detection.

cohvv_sm_diff difference between current and previous (cohvv_sm) measurements. This feature captures the increase in the (cohvv_sm), which is highly related to mowing detection.

cohvh_sm_diff difference between current and previous (cohvh_sm) measurements. This feature captures the increase in the (cohvh_sm), which is highly related to mowing detection.

ndvi_der The slope of the line between previous and current NDVI values.

cohvh_sm_der The slope of the line between previous and current (cohvh_sm) values. This feature captures the change in the smoothed cohvh.

cohvv_sm_der The slope of the line between previous and current (cohvv_sm) values. This feature captures the change in the smoothed cohvv.

Feature selectionThe permutation feature importance measurement was introduced by Breiman57. The importance of a particular feature is measured by the increase in the model’s prediction error after we permuted the values of this feature, which breaks the relationship between the feature and the outcome. A feature is important if shuffling its values increases the model error and is less important otherwise. The importance of features considered in this study is ranked in Table 3. It is notable from Table 3 that the ordinal features are significantly more important than the derived ones. We used backwards elimination to select the optimal subset of features to be used by the machine learning model. More specifically, we start with all the features and then remove the least significant feature at each iteration, which improves the model’s overall performance. We repeat this until no improvement is observed on the removal of features. Figures 7 and 8 show that the end of season accuracy(EOS) and event accuracy, respectively, for training using a different subsets of the most important features. We refer to (F_{x}-F_{y}) to be the set of important features from feature x to feature y in Table 3. Figure 7 shows that using only ndvi and (mixed_{coh}) achieves EOS of 93%. Increasing the number of the most important features to 3 achieves a comparable performance to the best one, as shown in Fig. 7. The results show that using the ndvi and (mixed_{coh}) achieve around 73% event accuracy while increasing the number of features, the performance declines as shown in Fig. 8. As an outcome of the feature selection process, the developed machine learning model used all the 14 features, shown in Table 3, that achieve the highest combined performance.Table 3 Ranking features based on their performance.Full size tableFigure 7End of season accuracy for different number of features.Full size image
Figure 8Event-based accuracy for different number of features.Full size image
Machine learning modelEach record in our dataset represents specific features about a field during one season at a particular time, in addition to the target variable (mown or non mown). In this work, we use a neural network to predict mowing events. We are interested only in observations during the vegetative season, so winter measurements are not included. More specifically, we only include the data in the vegetative season, which is almost the same across all Estonia from April till October (215 days). The dataset is partitioned into 64% for training, 20% for testing and 16% for validation. All training and testing were performed using TensorFlow58 deep learning framework with default parameters. The architecture of the neural network used is shown in Fig. 9. To guarantee a fixed time interval of 1-day, all the missing values in S-1 and S-2 features are interpolated, as shown in Fig. 10. The data is processed in batches of size (64 times 215) (times )14, where 64 is the number of fields considered per patch, 215 is the number of days in the vegetation season in Estonia, 14 is the number of selected features.Figure 9Architecture of the proposed model.Full size imageThe network’s output is a vector of size 215, representing the probability of a mowing event on each day in the vegetation season. The network consists of three one dimension convolution layers. The first and second convolution layers are followed by the Softmax activation function and batch normalization layer, while the third convolution is followed by Sigmoid activation function. The NN hyperparameters required to achieve the model learning process can significantly affect model performance. These hyperparameters include the following56:

Number of epochs represents how many times you want your algorithm to train on your whole dataset.

Loss function represents the prediction error of Neural Network.

Optimizer represents algorithm or method used to change the attributes of the neural network such as weights and learning rate to reduce the loss.

Activation function is the function through which we pass our weighted suown to have a significant output, namely as a vector of probability or a 0–1 output.

Learning rate refers to the step of backpropagation, when parameters are updated according to an optimization function.

Figure 10Time series mowing events before and after linear interpolation.Full size imageA good model uses the optimal combination of these hyperparameters and achieves good generalization capability. The training was performed with the conjugate gradient descent method and the binary cross-entropy loss function. The neural network was trained during 300 epochs; an early stopping was used59. The optimizer used in our model is Nadam optimizer60 with the following parameters: (beta_1=0.9), (beta_2=0.999), (epsilon=None), (schedule_{decay}=0.004), and learning (rate=0.001). Different activation functions such as ReLU, Sigmoid, Linear, and Tanh have been experimentally evaluated on the testing dataset as shown in Fig. 11. The results show that the Softmax activation function achieves the highest combined performance (event accuracy of 72.6% and EOS of 94.5%), as shown in Fig. 11.Figure 11Performance of different activation functions.Full size imageUsing 1D convolution layer acts as a filter that slides on the time dimension allowing the model to predict future mowing events from past events. However, this approach is not suitable for real-time detection of mowing events, but we use it to predict mowing events within a fixed time frame (window). Such a time frame should be greater than half the (1-D) convolution window length.Model evaluationTo evaluate our model, we used two metrics, EOS accuracy and Event-based accuracy. EOS is the accuracy of detecting a mowing event at least once during the season. If the probability of detecting a mowing event at least once during the season is more than 50%, then the field is considered mown, otherwise not mown. Event-based accuracy is used to evaluate how well our model correctly predicts mowing events. The formula for quantifying the binary accuracy is defined as follows:$$begin{aligned} acc = frac{TP + TN}{TP + TN + FP + FN} end{aligned}$$
(9)
where TP is the number of times that the model correctly predicted mowing events, given that the start day of the predicted mowing event is not more than 3 days earlier and not more than 6 days later than the actual start day of the mowing event. Within these 9 days, several mowing events may be predicted. To handle this case, only the first predicted mowing event is considered TP, and every next one is considered an FP. TN is the number of times that the model correctly predicted the absence of mowing events. FP is the number of times that the model incorrectly predicted mowing events. It also includes the number of times that the model correctly predicted mowing events, but the start of the event does not fit into a 9-days time frame with the actual start of some mowing event. FN is the number of times where the model missed actual mowing events.Reject region
Figure 12Calibration plot for proposed model.Full size image
Sometimes the model is not confident enough to give a reliable decision about the state of the field. We cannot expect reliable and confident predictions from inaccurate, incomplete or uncertain data. So, it is better in the cases of uncertainty about the prediction to allow the model to abstain from prediction. In this way, the obtained predictions are more accurate, while human experts could check rejected fields. Given the true positive rate and the true negative rate on the validation set, the reject region technique outputs a probability interval ((t_{low}), (t_{upper})) in which the model abstain prediction, where (t_{low}) and (t_{upper}) are the minimum and maximum probabilities that the model is uncertain about its prediction. Out of this interval, the model is confident about its prediction and predicts afield as mowed if the probability is higher than (t_{upper}) and not mown if the probability is less than (t_{low}). We select (t_{upper}), such that the desired true positive rate is reached. To find (t_{upper}), we sort all positives descending by their predicted probabilities and select the top percentage equal to the true positive rate. We choose (t_{low}) such that the desired true negative rate on validation data is reached. To find (t_{low}), we sort all negatives ascending by their predicted probabilities and select the top percentage equal to the true negative rate.Figure 12 shows the calibration plot for our proposed model. Notably, the predicted probabilities are close to the diagonal, which implies that the model is well-calibrated. More

1.Mordecai, E. A., Ryan, S. J., Caldwell, J. M., Shah, M. M. & LaBeaud, A. D. Climate change could shift disease burden from malaria to arboviruses in Africa. Lancet Planet. Health 4, e416–e423 (2020).PubMed 
PubMed Central 

Google Scholar 
2.W. H. O. Multisectoral approach to the prevention and control of vector-borne diseases (2020).3.Ryan, S. J. et al. Warming temperatures could expose more than 1.3 billion new people to Zika virus risk by 2050. Glob. Change Biol. 27, 84–93 (2021).
Google Scholar 
4.Iwamura, T., Guzman-Holst, A. & Murray, K. A. Accelerating invasion potential of disease vector Aedes aegypti under climate change. Nat. Commun. 11, 2130 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
5.Savage, V. M., Gillooly, J. F., Brown, J. H., West, G. B. & Charnov, E. L. Effects of body size and temperature on population growth. Am. Nat. 163, 429–441 (2004).PubMed 

Google Scholar 
6.Shocket, M. S. et al. Transmission of West Nile and five other temperate mosquito-borne viruses peaks at temperatures between 23 °C and 26 °C. eLife 9, 1–67 (2020).
Google Scholar 
7.Couret, J., Dotson, E. & Benedict, M. Q. Temperature, larval diet, and density effects on development rate and survival of Aedes aegypti (Diptera: Culicidae). PLoS ONE 9, 1–9 (2014).
Google Scholar 
8.Barreaux, A. M. G., Stone, C. M., Barreaux, P. & Koella, J. C. The relationship between size and longevity of the malaria vector Anopheles gambiae (s.s.) depends on the larval environment. Parasites Vectors 11, 485 (2018).PubMed 
PubMed Central 

Google Scholar 
9.Huxley, P. J., Murray, K. A., Pawar, S. & Cator, L. J. The effect of resource limitation on the temperature dependence of mosquito population fitness. Proc. R. Soc. B: Biol. Sci. 288, rspb.2020.3217 (2021).10.Ostfeld, R. S. & Keesing, F. Pulsed resources and community dynamics of consumers in terrestrial ecosystems. Trends Ecol. Evol. 15, 232–237 (2000).CAS 
PubMed 

Google Scholar 
11.Beltran, R. S. et al. Seasonal resource pulses and the foraging depth of a Southern Ocean top predator. Proc. R. Soc. B: Biol. Sci. 288, rspb.2020.2817 (2021).12.Yang, L. H., Bastow, J. L., Spence, K. O. & Wright, A. N. What can we learn from resource pulses? Ecology 89, 621–634 (2008).PubMed 

Google Scholar 
13.Dye, C. Models for the population dynamics of the yellow fever mosquito, Aedes aegypti. J. Animal Ecol. 53, 247 (1984).
Google Scholar 
14.Southwood, T. R., Murdie, G., Yasuno, M., Tonn, R. J. & Reader, P. M. Studies on the life budget of Aedes aegypti in Wat Samphaya, Bangkok, Thailand. Bull. World Health Organ. 46, 211–226 (1972).CAS 
PubMed 
PubMed Central 

Google Scholar 
15.Arrivillaga, J. & Barrera, R. Food as a limiting factor for Aedes aegypti in water-storage containers. J. Vector Ecol. 29, 11–20 (2004).PubMed 

Google Scholar 
16.Barrera, R., Amador, M. & Clark, G. G. Ecological factors influencing Aedes aegypti (Diptera: Culicidae) productivity in artificial containers in Salinas, Puerto Rico. J. Med. Entomol. 43, 484–492 (2006).PubMed 

Google Scholar 
17.Yee, D. A. & Juliano, S. A. Concurrent effects of resource pulse amount, type, and frequency on community and population properties of consumers in detritus-based systems. Oecologia 169, 511–522 (2012).PubMed 

Google Scholar 
18.Subra, R. & Mouchet, J. The regulation of preimaginal populations of Aedes aegypti (L.) (Diptera: Culicidae) on the Kenya coast. Ann. Trop. Med. Parasitol. 78, 63–70 (1984).CAS 
PubMed 

Google Scholar 
19.Amarasekare, P. & Savage, V. A framework for elucidating the temperature dependence of fitness. Am. Nat. 179, 178–191 (2012).PubMed 

Google Scholar 
20.Huey, R. B. & Kingsolver, J. G. Climate warming, resource availability, and the metabolic meltdown of ectotherms. Am. Nat. 194, 6 (2019).21.García-Carreras, B. et al. Role of carbon allocation efficiency in the temperature dependence of autotroph growth rates. Proc. Natl Acad. Sci. USA 115, E7361–E7368 (2018).PubMed 
PubMed Central 

Google Scholar 
22.Smith, T. P., Clegg, T., Bell, T. & Pawar, S. Systematic variation in the temperature dependence of bacterial carbon use efficiency. Ecol. Lett. 24, 2123–2133 (2021).PubMed 

Google Scholar 
23.Lehmann, P. et al. Complex responses of global insect pests to climate warming. Front. Ecol. Environ. 18, 141–150 (2020).
Google Scholar 
24.Amarasekare, P. Effects of climate warming on consumer-resource interactions: a latitudinal perspective. Front. Ecol. Evol. 7, 1–15 (2019).25.Amarasekare, P. & Simon, M. W. Latitudinal directionality in ectotherm invasion success. Proc. R. Soc. B: Biol. Sci. 287, 20191411 (2020).
Google Scholar 
26.Diagne, C. et al. High and rising economic costs of biological invasions worldwide. Nature 592, 571–576 (2021).CAS 
PubMed 

Google Scholar 
27.Cross, W. F., Hood, J. M., Benstead, J. P., Huryn, A. D. & Nelson, D. Interactions between temperature and nutrients across levels of ecological organization. Glob. Change Biol. 21, 1025–1040 (2015).
Google Scholar 
28.Mordecai, E. A. et al. Thermal biology of mosquito‐borne disease. Ecol. Lett. 22, 1690–1708 (2019).PubMed 
PubMed Central 

Google Scholar 
29.Thomas, M. K. et al. Temperature-nutrient interactions exacerbate sensitivity to warming in phytoplankton. Glob. Change Biol. 23, 3269–3280 (2017).
Google Scholar 
30.Siegel, P., Baker, K. G., Low‐Décarie, E. & Geider, R. J. High predictability of direct competition between marine diatoms under different temperatures and nutrient states. Ecol. Evol. 10, 7276–7290 (2020).PubMed 
PubMed Central 

Google Scholar 
31.Bestion, E., García-Carreras, B., Schaum, C.-E., Pawar, S. & Yvon-Durocher, G. Metabolic traits predict the effects of warming on phytoplankton competition. Ecol. Lett. 21, 655–664 (2018).PubMed 
PubMed Central 

Google Scholar 
32.Jackson, C. flexsurv: A Platform for Parametric Survival Modeling in R. J. Stat. Softw. 70, 1–33 (2016).
Google Scholar 
33.Bellows, T. S. The descriptive properties of some models for density dependence. J. Animal Ecol. 50, 139–156 (1981).
Google Scholar 
34.Orcutt, J. D. & Porter, K. G. The synergistic effects of temperature and food concentration of life history parameters of Daphnia. Oecologia 63, 300–306 (1984).PubMed 

Google Scholar 
35.Huey, R. B. & Berrigan, D. Temperature, demography, and ectotherm fitness. Am. Nat. 158, 204–210 (2001).CAS 
PubMed 

Google Scholar 
36.Caswell, H. A general formula for the sensitivity of population growth rate to changes in life history parameters. Theor. Popul. Biol. 14, 215–230 (1978).CAS 
PubMed 

Google Scholar 
37.Kammenga, J. E., Busschers, M., Straalen, N. M., Van, Jepson, P. C. & Bakker, J. Stress induced fitness reduction is not determined by the most sensitive life-cycle trait. Funct. Ecol. 10, 106 (1996).
Google Scholar 
38.Cator, L. J. et al. The role of vector trait variation in vector-borne disease dynamics. Front. Ecol. Evol. 8, 1–25 (2020).
Google Scholar 
39.Juliano, S. A. Species introduction and replacement among mosquitoes: interspecific resource competition or apparent competition? Ecology 79, 255 (1998).
Google Scholar 
40.Shapiro, L. L. M., Murdock, C. C., Jacobs, G. R., Thomas, R. J. & Thomas, M. B. Larval food quantity affects the capacity of adult mosquitoes to transmit human malaria. Proc. R. Soc. B: Biol. Sci. 283, 20160298 (2016).
Google Scholar 
41.Carvajal-Lago, L., Ruiz-López, M. J., Figuerola, J. & Martínez-de la Puente, J. Implications of diet on mosquito life history traits and pathogen transmission. Environ. Res. 195, 110893 (2021).CAS 
PubMed 

Google Scholar 
42.Reiner, R. C. et al. A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970-2010. J. R. Soc. Interface 10, 20120921–20120921 (2013).PubMed 
PubMed Central 

Google Scholar 
43.Farjana, T., Tuno, N. & Higa, Y. Effects of temperature and diet on development and interspecies competition in Aedes aegypti and Aedes albopictus. Med.Vet. Entomol. 26, 210–217 (2012).CAS 
PubMed 

Google Scholar 
44.Kooijman, S. A. L. M. Dynamic energy and mass budgets in biological systems. (Cambridge University Press, 2000).45.Merritt, R. W., Dadd, R. H. & Walker, E. D. Feeding behaviour, natural food, and nutritional relationships and larval mosquitoes. Annu. Rev. Entomol. 37, 349–376 (1992).46.Craine, J. M., Fierer, N. & McLauchlan, K. K. Widespread coupling between the rate and temperature sensitivity of organic matter decay. Nat. Geosci. 3, 854–857 (2010).CAS 

Google Scholar 
47.Smith, T. P. et al. Community-level respiration of prokaryotic microbes may rise with global warming. Nat. Commun. 10, 5124 (2019).PubMed 
PubMed Central 

Google Scholar 
48.Yee, D. A., Kaufman, M. G. & Juliano, S. A. The significance of ratios of detritus types and micro-organism productivity to competitive interactions between aquatic insect detritivores. J. Animal Ecol. 76, 1105–1115 (2007).
Google Scholar 
49.Chouaia, B. et al. Delayed larval development in Anopheles mosquitoes deprived of Asaia bacterial symbionts. BMC Microbiol. 12, S2 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
50.Souza, R. S. et al. Microorganism-based larval diets affect mosquito development, size and nutritional reserves in the yellow fever mosquito Aedes aegypti (Diptera: Culicidae). Front. Physiol. 10, 1–24 (2019).
Google Scholar 
51.Dickson, L. B. et al. Carryover effects of larval exposure to different environmental bacteria drive adult trait variation in a mosquito vector. Sci. Adv. 3, e1700585 (2017).PubMed 
PubMed Central 

Google Scholar 
52.Hery, L. et al. Natural variation in physicochemical profiles and bacterial communities associated with Aedes aegypti breeding sites and larvae on Guadeloupe and French Guiana. Microbial Ecol. 81, 93–109 (2021).CAS 

Google Scholar 
53.Liikanen, A., Murtoniemi, T., Tanskanen, H., Väisänen, T. & Martikainen, P. J. Effects of temperature and oxygen availability on greenhouse gas and nutrient dynamics in sediment of a eutrophic mid-boreal lake. Biogeochemistry 59, 269–286 (2002).CAS 

Google Scholar 
54.Lister, B. C. & Garcia, A. Climate-driven declines in arthropod abundance restructure a rainforest food web. Proc. Natl Acad. Sci. USA 115, E10397–E10406 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
55.Du, E. et al. Global patterns of terrestrial nitrogen and phosphorus limitation. Nat. Geosci. 13, 221–226 (2020).CAS 

Google Scholar 
56.Briegel, H. Metabolic relationship between female body size, reserves, and fecundity of Aedes aegypti. J. Insect Physiol. 36, 165–172 (1990).
Google Scholar 
57.Steinwascher, K. Relationship between pupal mass and adult survivorship and fecundity for Aedes aegypti. Environ. Entomol. 11, 150–153 (1982).
Google Scholar 
58.Trisos, C. H., Merow, C. & Pigot, A. L. The projected timing of abrupt ecological disruption from climate change. Nature 580, 496–501 (2020).CAS 
PubMed 

Google Scholar 
59.Parmesan, C. Ecological and evolutionary responses to recent climate change. Annu. Rev. Ecol., Evol. Syst. 37, 637–669 (2006).
Google Scholar 
60.Taheri, S., Naimi, B., Rahbek, C. & Araújo, M. B. Improvements in reports of species redistribution under climate change are required. Sci. Adv. 7, eabe1110 (2021).PubMed 
PubMed Central 

Google Scholar 
61.Bargielowski, I. E., Lounibos, L. P. & Carrasquilla, M. C. Evolution of resistance to satyrization through reproductive character displacement in populations of invasive dengue vectors. Proc. Natl Acad. Sci. USA 110, 2888–2892 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
62.Arguez, A. et al. NOAA’s 1981–2010 U.S. climate normals: an overview. Bull. Am. Meteorol. Soc. 93, 1687–1697 (2012).
Google Scholar 
63.Caswell, H. Matrix population models construction, analysis, and interpretation. Nat. Resource Model. (Sinauer Associates, 1989).64.Birch, L. C. The intrinsic rate of natural increase of an insect population. J. Animal Ecol. 17, 15 (1948).
Google Scholar 
65.Cole, L. C. The population consequences of life history phenomena. Q. Rev. Biol. 29, 103–137 (1954).CAS 
PubMed 

Google Scholar 
66.R. Core Team. R: A language and environment for statistical computing. (2018).67.Stubben, C. & Milligan, B. Estimating and analyzing demographic models using the popbio Package in R. J. Stat. Softw. 22, 1–23 (2007).
Google Scholar 
68.Therneau, T. A Package for Survival Analysis in R. (2021).69.Agnew, P., Hide, M., Sidobre, C. & Michalakis, Y. A minimalist approach to the effects of density-dependent competition on insect life-history traits. Ecol. Entomol. 27, 396–402 (2002).
Google Scholar 
70.Honěk, A. Intraspecific variation in body size and fecundity in insects: a general relationship. Oikos 66, 483 (1993).
Google Scholar 
71.Livdahl, T. P. & Sugihara, G. Non-linear interactions of populations and the importance of estimating per capita rates of change. J. Animal Ecol. 53, 573 (1984).
Google Scholar 
72.Juliano, S. A. & Lounibos, L. P. Ecology of invasive mosquitoes: effects on resident species and on human health. Ecol. Lett. 8, 558–574 (2005).PubMed 
PubMed Central 

Google Scholar 
73.van den Heuvel, M. J. The effect of rearing temperature on the wing length, thorax length, leg length and ovariole number of the adult mosquito, Aedes aegypti (L.). Trans. R. Entomol. Soc. Lond. 115, 197–216 (1963).
Google Scholar 
74.Farjana, T. & Tuno, N. Effect of body size on multiple blood feeding and egg retention of Aedes aegypti (L.) and Aedes albopictus (Skuse) (Diptera: Culicidae). Med. Entomol. Zool. 63, 123–131 (2012).
Google Scholar 
75.Skalski, J. R., Millspaugh, J. J., Dillingham, P. & Buchanan, R. A. Calculating the variance of the finite rate of population change from a matrix model in Mathematica. Environ. Model. Softw. 22, 359–364 (2007).
Google Scholar 
76.Hope, R. M. Rmisc: Rmisc: Ryan Miscellaneous. (2013).77.Caswell, H., Naiman, R. J. & Morin, R. Evaluating the consequences of reproduction in complex salmonid life cycles. Aquaculture 43, 123–134 (1984).
Google Scholar 
78.de Kroon, H., Plaisier, A., van Groenendael, J. & Caswell, H. Elasticity: the relative contribution of demographic parameters to population growth rate. Ecology 67, 1427–1431 (1986).
Google Scholar 
79.Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Statist. Softw. 67, (2015).80.Padfield, D., O’Sullivan, H. & Pawar, S. rTPC and nls.multstart: a new pipeline to fit thermal performance curves in R. Methods Ecol. Evol. 12, 1138–1143 (2021).
Google Scholar 
81.Lactin, D. J., Holliday, N. J., Johnson, D. L. & Craigen, R. Improved rate model of temperature-dependent development by arthropods. Environ. Entomol. 24, 68–75 (1995).
Google Scholar 
82.Kamykowski, D. & McCollum, S. A. The temperature acclimatized swimming speed of selected marine dinoflagellates. J. Plankton Res. 8, 275–287 (1986).
Google Scholar  More