Ecology

1.Li, J. et al. Development, progress and future prospects in cryobiotechnology of Lilium spp. Plant Methods 15, 125 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
2.Jin, L., Zhang, Y., Yan, L., Guo, Y. & Niu, L. Phenolic compounds and antioxidant activity of bulb extracts of six Lilium species native to China. Molecules 17, 9361–9378 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
3.Zhao, B., Zhang, J., Guo, X. & Wang, J. Microwave-assisted extraction, chemical characterization of polysaccharides from Lilium davidii var. unicolor Salisb and its antioxidant activities evaluation. Food Hydrocoll. 31, 346–356 (2013).CAS 
Article 

Google Scholar 
4.Wu, Z. et al. In vitro study of the growth, development and pathogenicity responses of Fusarium oxysporum to phthalic acid, an autotoxin from Lanzhou lily. World J. Microb. Biot. 31, 1227–1234 (2015).CAS 
Article 

Google Scholar 
5.Wu, Z. et al. Identification of autotoxins from root exudates of Lanzhou lily (Lilium davidii var. unicolor). Allelopathy J. 35, 35–48 (2015).
Google Scholar 
6.Vaitauskienė, K., Šarauskis, E., Naujokienė, V. & Liakas, V. The influence of free-living nitrogen-fixing bacteria on the mechanical characteristics of different plant residues under no-till and strip-till conditions. Soil Till. Res. 154, 91–102 (2015).Article 

Google Scholar 
7.Doornbos, R. F., van Loon, L. C. & Bakker, P. A. H. M. Impact of root exudates and plant defense signaling on bacterial communities in the rhizosphere. A review. Agron. Sustain. Dev. 32, 227–243 (2012).Article 

Google Scholar 
8.Liu, X. et al. Long-term greenhouse cucumber production alters soil bacterial community structure. J. Soil Sci. Plant Nut. 20, 306–321 (2020).Article 
CAS 

Google Scholar 
9.Liu, X. et al. Microbial community diversities and taxa abundances in soils along a seven-year gradient of potato monoculture using high throughput pyrosequencing approach. PLoS ONE 9, e86610 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
10.Li, X., Ding, C., Zhang, T. & Wang, X. Fungal pathogen accumulation at the expense of plant-beneficial fungi as a consequence of consecutive peanut monoculturing. Soil Biol. Biochem. 72, 11–18 (2014).Article 
CAS 

Google Scholar 
11.Zhao, Q. et al. Long-term coffee monoculture alters soil chemical properties and microbial communities. Sci. Rep. 8, 6116 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
12.Shang, Q. et al. Illumina-based analysis of the rhizosphere microbial communities associated with healthy and wilted Lanzhou lily (Lilium davidii var. unicolor) plants grown in the field. World J. Microb. Biot. 32, 95 (2016).Article 
CAS 

Google Scholar 
13.Wu, K. et al. Effects of bio-organic fertilizer plus soil amendment on the control of tobacco bacterial wilt and composition of soil bacterial communities. Biol. Fert. Soils 50, 961–971 (2014).Article 

Google Scholar 
14.Baćmaga, M., Wyszkowska, J. & Kucharski, J. The influence of chlorothalonil on the activity of soil microorganisms and enzymes. Ecotoxicology 27, 1188–1202 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
15.Ding, H. et al. Influence of chlorothalonil and carbendazim fungicides on the transformation processes of urea nitrogen and related microbial populations in soil. Environ. Sci. Pollut. R. 26, 31133–31141 (2019).CAS 
Article 

Google Scholar 
16.Passari, A. K. et al. Biocontrol of Fusarium wilt of Capsicum annuum by rhizospheric bacteria isolated from turmeric endowed with plant growth promotion and disease suppression potential. Eur. J. Plant Pathol. 150, 831–846 (2018).CAS 
Article 

Google Scholar 
17.Wang, L. et al. Application of bioorganic fertilizer significantly increased apple yields and shaped bacterial community structure in orchard soil. Microb. Ecol. 73, 404–416 (2017).CAS 
PubMed 
Article 

Google Scholar 
18.Bi, Y. et al. Differential effects of two earthworm species on Fusarium wilt of strawberry. Appl. Soil Ecol. 126, 174–181 (2018).Article 

Google Scholar 
19.Zhao, F. et al. Vermicompost improves microbial functions of soil with continuous tomato cropping in a greenhouse. J. Soil Sediment 20, 380–391 (2020).CAS 
Article 

Google Scholar 
20.Mendes, R. et al. Deciphering the rhizosphere microbiome for disease-suppressive bacteria. Science 332, 1097–1100 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
21.Hermans, S. M. et al. Bacteria as emerging indicators of soil condition. Appl. Environ. Microb. 83, 13 (2017).Article 

Google Scholar 
22.Brown, M. E. & Chang, M. C. Y. Exploring bacterial lignin degradation. Curr. Opin. Chem. Biol. 19, 1–7 (2014).CAS 
PubMed 
Article 

Google Scholar 
23.Tian, J., Pourcher, A., Bouchez, T., Gelhaye, E. & Peu, P. Occurrence of lignin degradation genotypes and phenotypes among prokaryotes. Appl. Microbiol. Biot. 98, 9527–9544 (2014).CAS 
Article 

Google Scholar 
24.She, S. et al. Significant relationship between soil bacterial community structure and incidence of bacterial wilt disease under continuous cropping system. Arch. Microbiol. 199, 267–275 (2017).CAS 
PubMed 
Article 

Google Scholar 
25.Shen, Z. et al. Soils naturally suppressive to banana Fusarium wilt disease harbor unique bacterial communities. Plant Soil 393, 21–33 (2015).CAS 
Article 

Google Scholar 
26.Feng, X. Y. The taxonomic characteristics of various genera of terrestrial earthworms in China. Chin. J. Zool. 1, 44–47 (1985).
Google Scholar 
27.Wu, W. R. Studies on the Germplasm Resources and Quality Estimation of Dilong (Eartheworm) (Guangzhou University of Chinese Medicine, 2008).
Google Scholar 
28.Zhou, L. et al. Effects of lily/maize intercropping on rhizosphere microbial community and yield of Lilium davidii var. unicolor. J. Basic Microb. 58, 892–901 (2018).Article 

Google Scholar 
29.Shi, G. Y. et al. Soil fungal diversity loss and appearance of specific fungal pathogenic communities associated with the consecutive replant problem (CRP) in lily. Front. Microbiol. 11, 1649 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
30.Blouin, M. et al. A review of earthworm impact on soil function and ecosystem services. Eur. J. Soil Sci. 64, 161–182 (2013).Article 

Google Scholar 
31.Doran, J. W. & Zeiss, M. R. Soil health and sustainability: Managing the biotic component of soil quality. Appl. Soil Ecol. 15, 3–11 (2000).Article 

Google Scholar 
32.Kaneda, S., Ohkubo, S., Wagai, R. & Yagasaki, Y. Soil temperature and moisture-based estimation of rates of soil aggregate formation by the endogeic earthworm Eisenia japonica (Michaelsen, 1892). Biol. Fert. Soils 52, 789–797 (2016).CAS 
Article 

Google Scholar 
33.Bottinelli, N. et al. Earthworms accelerate soil porosity turnover under watering conditions. Geoderma 156, 43–47 (2010).ADS 
Article 

Google Scholar 
34.Eriksen-Hamel, N. S. & Whalen, J. K. Impacts of earthworms on soil nutrients and plant growth in soybean and maize agroecosystems. Agric. Ecosyst. Environ. 120, 442–448 (2007).CAS 
Article 

Google Scholar 
35.Rousk, J., Brookes, P. C. & Baath, E. The microbial PLFA composition as affected by pH in an arable soil. Soil Biol. Biochem. 42, 516–520 (2010).CAS 
Article 

Google Scholar 
36.Stiles, C. M. & Murray, T. D. Infection of field-grown winter wheat by cephalosporium gramineum and the effect of soil pH. Phytopathology 86, 177–183 (1996).Article 

Google Scholar 
37.Weyman-Kaczmarkowa, W. & Pedziwilk, Z. The development of fungi as affected by pH and type of soil, in relation to the occurrence of bacteria and soil fungistatic activity. Microbiol. Res. 155, 107–112 (2000).CAS 
PubMed 
Article 

Google Scholar 
38.Salmon, S. Earthworm excreta (mucus and urine) affect the distribution of springtails in forest soils. Biol. Fert. Soils 34, 304–310 (2001).CAS 
Article 

Google Scholar 
39.García-Montero, L. G., Valverde-Asenjo, I., Grande-Ortíz, M. A., Menta, C. & Hernando, S. Impact of earthworm casts on soil pH and calcium carbonate in black truffle burns. Agroforest. Syst. 87, 815–826 (2013).Article 

Google Scholar 
40.Bending, G. D., Turner, M. K. & Jones, J. E. Interactions between crop residue and soil organic matter quality and the functional diversity of soil microbial communities. Soil Biol. Biochem. 34, 1073–1082 (2002).CAS 
Article 

Google Scholar 
41.Wang, H. et al. Effects of long-term application of organic fertilizer on improving organic matter content and retarding acidity in red soil from China. Soil Till. Res. 195, 104382 (2019).Article 

Google Scholar 
42.Blouin, M., Sery, N., Cluzeau, D., Brun, J. J. & Bédécarrats, A. Balkanized research in ecological engineering revealed by a bibliometric analysis of earthworms and ecosystem services. Environ. Manage. 52, 309–320 (2013).ADS 
PubMed 
Article 

Google Scholar 
43.Wu, Y., Zhang, N., Wang, J. & Sun, Z. An integrated crop-vermiculture system for treating organic waste on fields. Eur. J. Soil Biol. 51, 8–14 (2012).Article 

Google Scholar 
44.Basker, A., Macgregor, A. N. & Kirkman, J. H. Influence of soil ingestion by earthworms on the availability of potassium in soil: An incubation experiment. Biol. Fert. Soils 14, 300–303 (1992).CAS 
Article 

Google Scholar 
45.Canellas, L. P., Olivares, F. L., Okorokova-Facanha, A. L. & Facanha, A. R. Humic acids isolated from earthworm compost enhance root elongation, lateral root emergence, and plasma membrane H+-ATPase activity in maize roots. Plant Physiol. 130, 1951–1957 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
46.Dvorak, J. et al. Sensing microorganisms in the gut triggers the immune response in Eisenia andrei earthworms. Dev. Comp. Immunol. 57, 67–74 (2016).CAS 
PubMed 
Article 

Google Scholar 
47.Ma, X., Xing, M., Wang, Y., Xu, Z. & Yang, J. Microbial enzyme and biomass responses: Deciphering the effects of earthworms and seasonal variation on treating excess sludge. J. Environ. Manage. 170, 207–214 (2016).CAS 
PubMed 
Article 

Google Scholar 
48.Groffman, P. M. et al. Earthworms increase soil microbial biomass carrying capacity and nitrogen retention in northern hardwood forests. Soil Biol. Biochem. 87, 51–58 (2015).CAS 
Article 

Google Scholar 
49.Gomez-Brandon, M., Lazcano, C., Lores, M. & Dominguez, J. Detritivorous earthworms modify microbial community structure and accelerate plant residue decomposition. Appl. Soil Ecol. 44, 237–244 (2010).Article 

Google Scholar 
50.Dempsey, M. A., Fisk, M. C. & Fahey, T. J. Earthworms increase the ratio of bacteria to fungi in northern hardwood forest soils, primarily by eliminating the organic horizon. Soil Biol. Biochem. 43, 2135–2141 (2011).CAS 
Article 

Google Scholar 
51.Wolfarth, F., Schrader, S., Oldenburg, E. & Weinert, J. Contribution of the endogeic earthworm species Aporrectodea caliginosa to the degradation of deoxynivalenol and Fusarium biomass in wheat straw. Mycotoxin Res. 27, 215–220 (2011).CAS 
PubMed 
Article 

Google Scholar 
52.Paul, B. K., Lubbers, I. M. & van Groenigen, J. W. Residue incorporation depth is a controlling factor of earthworm-induced nitrous oxide emissions. Glob. Change Biol. 18, 1141–1151 (2012).ADS 
Article 

Google Scholar 
53.Chen, Y., Chang, S. K. C., Chen, J., Zhang, Q. & Yu, H. Characterization of microbial community succession during vermicomposting of medicinal herbal residues. Bioresource Technol. 249, 542–549 (2018).CAS 
Article 

Google Scholar 
54.Tao, J. et al. Effects of earthworms on soil enzyme activity in an organic residue amended rice-wheat rotation agro-ecosystem. Appl. Soil Ecol. 42, 221–226 (2009).Article 

Google Scholar 
55.Bertrand, M. et al. Earthworm services for cropping systems. A review. Agron. Sustain. Dev. 35, 553–567 (2015).CAS 
Article 

Google Scholar 
56.Jayasinghe, B. A. T. D. & Parkinson, D. Earthworms as the vectors of actinomycetes antagonistic to litter decomposer fungi. Appl. Soil Ecol. 43, 1–10 (2009).Article 

Google Scholar 
57.Fierer, N., Bradford, M. A. & Jackson, R. B. Toward an ecological classification of soil bacteria. Ecology 88, 1354–1364 (2007).PubMed 
Article 

Google Scholar 
58.Kopecky, J. et al. Actinobacterial community dominated by a distinct clade in acidic soil of a waterlogged deciduous forest. Fems Microbiol. Ecol. 78, 386–394 (2011).CAS 
PubMed 
Article 

Google Scholar 
59.Sun, J., Zhang, Q., Zhou, J. & Wei, Q. Pyrosequencing technology reveals the impact of different manure doses on the bacterial community in apple rhizosphere soil. Appl. Soil Ecol. 78, 28–36 (2014).Article 

Google Scholar 
60.Bull, A. T., Stach, J. E. M., Ward, A. C. & Goodfellow, M. Marine actinobacteria: perspectives, challenges, future directions. Anton. Leeuw. Int. J. G. 87, 65–79 (2005).Article 

Google Scholar 
61.Wang, Q. et al. Long-term fertilization changes bacterial diversity and bacterial communities in the maize rhizosphere of Chinese Mollisols. Appl. Soil Ecol. 125, 88–96 (2018).ADS 
Article 

Google Scholar 
62.Jones, R. T. et al. A comprehensive survey of soil acidobacterial diversity using pyrosequencing and clone library analyses. ISME J. 3, 442–453 (2009).CAS 
PubMed 
Article 

Google Scholar 
63.Liu, L. et al. Neorhizobium lilium sp. nov., an endophytic bacterium isolated from Lilium pumilum bulbs in Hebei province. Arch. Microbiol. 202, 609–616 (2020).CAS 
PubMed 
Article 

Google Scholar 
64.Vessey, J. K. Plant growth promoting rhizobacteria as biofertilizers. Plant Soil 255, 571–586 (2003).CAS 
Article 

Google Scholar 
65.Taylor, W. J. & Draughon, F. A. Nannocystis exedens: A potential biocompetitive agent against Aspergillus flavus and Aspergillus parasiticus. J. Food Protect. 64, 1030–1034 (2001).CAS 
Article 

Google Scholar 
66.Viswanathan, R. & Samiyappan, R. Induced systemic resistance by fluorescent pseudomonads against red rot disease of sugarcane caused by Colletotrichum falcatum. Crop Prot. 21, 1–10 (2002).Article 

Google Scholar 
67.Wolters, V. Invertebrate control of soil organic matter stability. Biol. Fert. Soils 31, 1–19 (2000).MathSciNet 
CAS 
Article 

Google Scholar 
68.Manna, M. C., Jha, S., Ghosh, P. K. & Achaya, C. L. Comparative efficacy of three epigeic earthworms under different deciduous forest litters decomposition. Bioresource Technol. 88, 197–206 (2003).CAS 
Article 

Google Scholar 
69.Felten, D. & Emmerling, C. Earthworm burrowing behaviour in 2D terraria with single-and multi-species assemblages. Biol. Fert. Soils 45, 789–797 (2009).Article 

Google Scholar 
70.Wang, Z. et al. Soil protist communities in burrowing and casting hotspots of different earthworm species. Soil Biol. Biochem. 144, 107774 (2020).CAS 
Article 

Google Scholar 
71.Danielle, J., Yvan, C. & Daniel, C. Interactions between earthworm species in artificial soil cores assessed through the 3D reconstruction of the burrow systems. Geoderma 102, 123–137 (2001).Article 

Google Scholar 
72.Gomes, N. C. M. et al. Dynamics of fungal communities in bulk and maize rhizosphere soil in the tropics. Appl. Environ. Microb. 69, 3758–3766 (2003).ADS 
CAS 
Article 

Google Scholar 
73.Bao, S. D. Analytical Methods of Soil Agrochemistry 3rd edn. (China Agricultural Press, 2000).
Google Scholar 
74.Chen, S., Zhou, Y., Chen, Y. & Gu, J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884–i890 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
75.Magoč, T. & Salzberg, S. L. FLASH: Fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27, 2957–2963 (2011).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
76.Edgar, R. C. UPARSE: Highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996–998 (2013).CAS 
PubMed 
Article 

Google Scholar 
77.Stackebrandt, E. & Goebel, B. M. Taxonomic note: A place for DNA-DNA reassociation and 16S rRNA sequence analysis in the present species definition in bacteriology. Int. J. Syst. Evol. Micr. 44, 846–849 (1994).CAS 
Article 

Google Scholar 
78.Wang, Q., Garrity, G. M., Tiedje, J. M. & Cole, J. R. Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261–5267 (2007).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar  More

1.IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (eds Pörtner, H.-O. et al.) (IPCC, 2019).2.Doney, S. C. et al. Climate change impacts on marine ecosystems. Ann. Rev. Mar. Sci. 4, 11–37 (2012).Article 

Google Scholar 
3.Bindoff, N. L. et al. in Special Report on the Ocean and Cryosphere in a Changing Climate (eds Pörtner, H.-O. et al.) Ch. 5 (IPCC, 2019).4.Griffith, G. P., Fulton, E. A., Gorton, R. & Richardson, A. J. Predicting interactions among fishing, ocean warming, and ocean acidification in a marine system with whole-ecosystem models. Conserv. Biol. 26, 1145–1152 (2012).Article 

Google Scholar 
5.Fu, C. et al. Risky business: the combined effects of fishing and changes in primary productivity on fish communities. Ecol. Modell. 368, 265–276 (2018).Article 

Google Scholar 
6.Tittensor, D. P. et al. Integrating climate adaptation and biodiversity conservation in the global ocean. Sci. Adv. https://doi.org/10.1126/sciadv.aay9969 (2019).7.IPBES: Summary for Policymakers. In Global Assessment Report on Biodiversity and Ecosystem Services (eds Díaz, S. et al.) (IPBES Secretariat, 2019).8.Boyce, D. G., Lotze, H. K., Tittensor, D. P., Carozza, D. A. & Worm, B. Future ocean biomass losses may widen socioeconomic equity gaps. Nat. Commun. 11, 2235 (2020).CAS 
Article 

Google Scholar 
9.Payne, M. R. et al. Uncertainties in projecting climate-change impacts in marine ecosystems. ICES J. Mar. Sci. 73, 1272–1282 (2016).Article 

Google Scholar 
10.Eyring, V. et al. Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. Geosci. Model Dev. 9, 1937–1958 (2016).Article 

Google Scholar 
11.Tittensor, D. P. et al. A protocol for the intercomparison of marine fishery and ecosystem models: Fish-MIP v1.0. Geosci. Model Dev. 11, 1421–1442 (2018).Article 

Google Scholar 
12.Lotze, H. K. et al. Global ensemble projections reveal trophic amplification of ocean biomass declines with climate change. Proc. Natl Acad. Sci. USA 116, 12907–12912 (2019).CAS 
Article 

Google Scholar 
13.Bryndum-Buchholz, A. et al. Twenty-first-century climate change impacts on marine animal biomass and ecosystem structure across ocean basins. Glob. Change Biol. 25, 459–472 (2019).Article 

Google Scholar 
14.Bryndum-Buchholz, A. et al. Differing marine animal biomass shifts under 21st century climate change between Canada’s three oceans. Facets 5, 105–122 (2020).Article 

Google Scholar 
15.Bryndum-Buchholz, A. et al. Climate-change impacts and fisheries management challenges in the North Atlantic Ocean. Mar. Ecol. Prog. Ser. 648, 1–17 (2020).Article 

Google Scholar 
16.Ruane, A. C. et al. The vulnerability, impacts, adaptation and climate services advisory board (VIACS AB v1.0) contribution to CMIP6. Geosci. Model Dev. 9, 3493–3515 (2016).Article 

Google Scholar 
17.Kwiatkowski, L. et al. Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections. Biogeosciences 17, 3439–3470 (2020).CAS 
Article 

Google Scholar 
18.Séférian, R. et al. Tracking improvement in simulated marine biogeochemistry between CMIP5 and CMIP6. Curr. Clim. Change Rep. 6, 95–119 (2020).Article 

Google Scholar 
19.Meehl, G. A. et al. Context for interpreting equilibrium climate sensitivity and transient climate response from the CMIP6 Earth system models. Sci. Adv. 6, eaba1981 (2020).Article 

Google Scholar 
20.Tebaldi, C. et al. Climate model projections from the Scenario Model Intercomparison Project (ScenarioMIP) of CMIP6. Earth Syst. Dyn. 12, 253–293 (2021).Article 

Google Scholar 
21.Heneghan, R. F. et al. Disentangling diverse responses to climate change among global marine ecosystem models. Prog. Oceanogr. 198, 102659 (2021).Article 

Google Scholar 
22.Zelinka, M. D. et al. Causes of higher climate sensitivity in CMIP6 models. Geophys. Res. Lett. 47, e2019GL085782 (2020).Article 

Google Scholar 
23.Kwiatkowski, L. et al. Emergent constraints on projections of declining primary production in the tropical oceans. Nat. Clim. Change 7, 355–358 (2017).CAS 
Article 

Google Scholar 
24.Cabré, A., Marinov, I. & Leung, S. Consistent global responses of marine ecosystems to future climate change across the IPCC AR5 Earth system models. Clim. Dyn. 45, 1253–1280 (2015).Article 

Google Scholar 
25.Laufkötter, C. et al. Drivers and uncertainties of future global marine primary production in marine ecosystem models. Biogeosciences 12, 6955–6984 (2015).Article 

Google Scholar 
26.Doney, S. C. Plankton in a warmer world. Nature 444, 695–696 (2006).CAS 
Article 

Google Scholar 
27.Rykaczewski, R. R. & Dunne, J. P. Enhanced nutrient supply to the California Current Ecosystem with global warming and increased stratification in an Earth system model. Geophys. Res. Lett. 37, L21606 (2010).Article 

Google Scholar 
28.Laufkötter, C., John, J. G., Stock, C. A. & Dunne, J. P. Temperature and oxygen dependence of the remineralization of organic matter. Glob. Biogeochem. Cycles 31, 1038–1050 (2017).Article 
CAS 

Google Scholar 
29.Coll, M. et al. Advancing global ecological modeling capabilities to simulate future trajectories of change in marine ecosystems. Front. Mar. Sci. 7, 741 (2020).Article 

Google Scholar 
30.Hawkins, E. & Sutton, R. The potential to narrow uncertainty in regional climate predictions. Bull. Am. Meteorol. Soc. 90, 1095–1107 (2009).Article 

Google Scholar 
31.Frölicher, T. L., Rodgers, K. B., Stock, C. A. & Cheung, W. W. L. Sources of uncertainties in 21st century projections of potential ocean ecosystem stressors. Glob. Biogeochem. Cycles 30, 1224–1243 (2016).Article 
CAS 

Google Scholar 
32.Gaines, S. D. et al. Improved fisheries management could offset many negative effects of climate change. Sci. Adv. 4, eaao1378 (2018).Article 

Google Scholar 
33.The State of World Fisheries and Aquaculture 2020 (FAO, 2020).34.Dahlke, F. T., Wohlrab, S., Butzin, M. & Pörtner, H.-O. Thermal bottlenecks in the life cycle define climate vulnerability of fish. Science 369, 65–70 (2020).CAS 
Article 

Google Scholar 
35.Stuart-Smith, R. D., Edgar, G. J. & Bates, A. E. Thermal limits to the geographic distributions of shallow-water marine species. Nat. Ecol. Evol. 1, 1846–1852 (2017).Article 

Google Scholar 
36.Carozza, D. A., Bianchi, D. & Galbraith, E. D. Metabolic impacts of climate change on marine ecosystems: implications for fish communities and fisheries. Glob. Ecol. Biogeogr. 28, 158–169 (2019).Article 

Google Scholar 
37.du Pontavice, H., Gascuel, D., Reygondeau, G., Stock, C. & Cheung, W. W. L. Climate-induced decrease in biomass flow in marine food webs may severely affect predators and ecosystem production. Glob. Change Biol. 27, 2608–2622 (2021).Article 

Google Scholar 
38.Piroddi, C. et al. Effects of nutrient management scenarios on marine food webs: a pan-European assessment in support of the marine strategy framework directive. Front. Mar. Sci. 8, 179 (2021).Article 

Google Scholar 
39.Maury, O. An overview of APECOSM, a spatialized mass balanced ‘Apex Predators ECOSystem Model’ to study physiologically structured tuna population dynamics in their ecosystem. Prog. Oceanogr. 84, 113–117 (2010).Article 

Google Scholar 
40.Maury, O. & Poggiale, J. C. From individuals to populations to communities: a dynamic energy budget model of marine ecosystem size-spectrum including life history diversity. J. Theor. Biol. 324, 52–71 (2013).Article 

Google Scholar 
41.Carozza, D. A., Bianchi, D. & Galbraith, E. D. The ecological module of BOATS-1.0: a bioenergetically-constrained model of marine upper trophic levels suitable for studies of fisheries and ocean biogeochemistry. Geosci. Model Dev. 9, 1545–1565 (2016).Article 

Google Scholar 
42.Carozza, D. A. et al. Formulation, general features and global calibration of a bioenergetically-constrained fishery model. PLoS ONE 12, e0169763 (2017).Article 
CAS 

Google Scholar 
43.Cheung, W. W. L. et al. Building confidence in projections of the responses of living marine resources to climate change. ICES J. Mar. Sci. https://doi.org/10.1093/icesjms/fsv250 (2016).Article 

Google Scholar 
44.Cheung, W. W. L., Dunne, J., Sarmiento, J. L. & Pauly, D. Integrating ecophysiology and plankton dynamics into projected maximum fisheries catch potential under climate change in the Northeast Atlantic. ICES J. Mar. Sci. 68, 1008–1018 (2011).Article 

Google Scholar 
45.Blanchard, J. L. et al. Potential consequences of climate change for primary production and fish production in large marine ecosystems. Phil. Trans. R. Soc. B 367, 2979–2989 (2012).Article 

Google Scholar 
46.Christensen, V. et al. The global ocean is an ecosystem: simulating marine life and fisheries. Glob. Ecol. Biogeogr. 24, 507–517 (2015).Article 

Google Scholar 
47.Gascuel, D., Guénette, S. & Pauly, D. The trophic-level-based ecosystem modelling approach: theoretical overview and practical uses. ICES J. Mar. Sci. 68, 1403–1416 (2011).Article 

Google Scholar 
48.Petrik, C. M., Stock, C. A., Andersen, K. H., van Denderen, P. D. & Watson, J. R. Bottom-up drivers of global patterns of demersal, forage, and pelagic fishes. Prog. Oceanogr. 176, 102124 (2019).Article 

Google Scholar 
49.Jennings, S. & Collingridge, K. Predicting consumer biomass, size-structure, production, catch potential, responses to fishing and associated uncertainties in the world’s marine ecosystems. PLoS ONE 10, e0133794 (2015).Article 
CAS 

Google Scholar 
50.Heneghan, R. F. et al. A functional size-spectrum model of the global marine ecosystem that resolves zooplankton composition. Ecol. Modell. 435, 109265 (2020).CAS 
Article 

Google Scholar 
51.Dunne, J. P. et al. GFDL’s ESM2 global coupled climate–carbon Earth system models. Part I: physical formulation and baseline simulation characteristics. J. Clim. 25, 6646–6665 (2012).Article 

Google Scholar 
52.Dunne, J. P. et al. Carbon Earth system models. Part II: carbon system formulation and baseline simulation characteristics. J. Clim. 26, 2247–2267 (2013).Article 

Google Scholar 
53.Dufresne, J.-L. et al. Climate change projections using the IPSL-CM5 Earth system model: from CMIP3 to CMIP5. Clim. Dyn. 40, 2123–2165 (2013).Article 

Google Scholar 
54.Dunne, J. P. et al. The GFDL Earth System Model Version 4.1 (GFDL-ESM 4.1): overall coupled model description and simulation characteristics. J. Adv. Model. Earth Syst. 12, e2019MS002015 (2020).
Google Scholar 
55.Krasting, J. P. et al. NOAA-GFDL GFDL-ESM4 Model Output Prepared for MIP6 CMIP Historical Version 20190726 (Earth System Grid Federation, 2018); https://doi.org/10.22033/ESGF/CMIP6.859756.John, J. G. et al. NOAA-GFDL GFDL-ESM4 Model Output Prepared for CMIP6 ScenarioMIP ssp585 Version 20180701 (Earth System Grid Federation, 2018); https://doi.org/10.22033/ESGF/CMIP6.870657.Boucher, O. et al. Presentation and evaluation of the IPSL-CM6A-LR climate model. J. Adv. Model. Earth Syst. 12, e2019MS002010 (2020).Article 

Google Scholar 
58.Boucher, O. et al. IPSL IPSL-CM6A-LR Model Output Prepared for CMIP6 CMIP Version 20180727 (Earth System Grid Federation, 2018); https://doi.org/10.22033/ESGF/CMIP6.153459.Boucher, O. et al. IPSL IPSL-CM6A-LR Model Output Prepared for CMIP6 CMIP Historical Version 20180103 (Earth System Grid Federation, 2018); https://doi.org/10.22033/ESGF/CMIP6.5195 More

1.Fahad, S. et al. Crop production under drought and heat stress: Plant responses and management options. Front. Plant Sci. 8, 1147 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
2.Engelbrecht, B. M., Kursar, T. A. & Tyree, M. T. Drought effects on seedling survival in a tropical moist forest. Trees 19, 312–321 (2005).Article 

Google Scholar 
3.Michaelian, M., Hogg, E. H., Hall, R. J. & Arsenault, E. Massive mortality of aspen following severe drought along the southern edge of the Canadian boreal forest. Glob. Change Biol. 17, 2084–2094 (2011).ADS 
Article 

Google Scholar 
4.Lehner, B., Döll, P., Alcamo, J., Henrichs, T. & Kaspar, F. Estimating the impact of global change on flood and drought risks in Europe: A continental, integrated analysis. Clim. Change 75, 273–299 (2006).ADS 
Article 

Google Scholar 
5.Allen, C. D. et al. A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. For. Ecol. Manag. 259, 660–684 (2010).Article 

Google Scholar 
6.IPCC. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. (Cambridge University Press, 2013).7.Nayyar, H., Kaur, S., Singh, S. & Upadhyaya, H. D. Differential sensitivity of Desi (small-seeded) and Kabuli (large-seeded) chickpea genotypes to water stress during seed filling: Effects on accumulation of seed reserves and yield. J. Sci. Food Agric. 86, 2076–2082 (2006).CAS 
Article 

Google Scholar 
8.Bodner, G., Nakhforoosh, A. & Kaul, H.-P. Management of crop water under drought: A review. Agron. Sustain. Dev. 35, 401–442 (2015).Article 

Google Scholar 
9.Osakabe, Y., Osakabe, K., Shinozaki, K. & Tran, L.-S.P. Response of plants to water stress. Front. Plant Sci. 5, 86 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
10.Anjum, S. A. et al. Morphological, physiological and biochemical responses of plants to drought stress. Afr. J. Agric. Res. 6, 2026–2032 (2011).
Google Scholar 
11.Farooq, M., Wahid, A., Kobayashi, N., Fujita, D. & Basra, S. Sustainable Agriculture 153–188 (Springer, Berlin, 2009).Book 

Google Scholar 
12.Schneider, F. & Don, A. Root-restricting layers in German agricultural soils. Part I: extent and cause. Plant Soil 442, 433–451 (2019).CAS 
Article 

Google Scholar 
13.Gupta, A., Rico-Medina, A. & Caño-Delgado, A. I. The physiology of plant responses to drought. Science 368, 266–269 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
14.Basu, S., Ramegowda, V., Kumar, A. & Pereira, A. Plant adaptation to drought stress. F1000Research 5 (2016).15.Coskun, D., Britto, D. T., Huynh, W. Q. & Kronzucker, H. J. The role of silicon in higher plants under salinity and drought stress. Front. Plant Sci. 7, 1072 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
16.Schaller, J., Puppe, D., Kaczorek, D., Ellerbrock, R. & Sommer, M. Silicon cycling in soils revisited. Plants 10, 295 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Epstein, E. The anomaly of silicon in plant biology. Proc. Natl. Acad. Sci. U. S. A. 91, 11–17 (1994).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Ma, J. F. & Yamaji, N. Silicon uptake and accumulation in higher plants. Trends Plant Sci. 11, 392–397 (2006).CAS 
PubMed 
Article 

Google Scholar 
19.Gong, H. J., Chen, K. M., Chen, G. C., Wang, S. M. & Zhang, C. L. Effects of silicon on growth of wheat under drought. J. Plant Nutr. 26, 1055–1063 (2003).CAS 
Article 

Google Scholar 
20.Hattori, T. et al. Application of silicon enhanced drought tolerance in Sorghum bicolor. Physiol. Plant. 123, 459–466 (2005).CAS 
Article 

Google Scholar 
21.Chen, W., Yao, X., Cai, K. & Chen, J. Silicon alleviates drought stress of rice plants by improving plant water status, photosynthesis and mineral nutrient absorption. Biol. Trace Elem. Res. 142, 67–76 (2011).CAS 
PubMed 
Article 

Google Scholar 
22.Ibrahim, M. A., Merwad, A.-R.M. & Elnaka, E. A. Rice (Oryza Sativa L.) tolerance to drought can be improved by silicon application. Commun. Soil Sci. Plant Anal. 49, 945–957 (2018).CAS 
Article 

Google Scholar 
23.Alzahrani, Y., Kuşvuran, A., Alharby, H. F., Kuşvuran, S. & Rady, M. M. The defensive role of silicon in wheat against stress conditions induced by drought, salinity or cadmium. Ecotox. Environ. Safe. 154, 187–196 (2018).CAS 
Article 

Google Scholar 
24.Meunier, J. D. et al. Effect of phytoliths for mitigating water stress in durum wheat. New Phytol. 215, 229–239 (2017).CAS 
PubMed 
Article 

Google Scholar 
25.Schaller, J., Cramer, A., Carminati, A. & Zarebanadkouki, M. Biogenic amorphous silica as main driver for plant available water in soils. Sci. Rep. 10, 2424 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
26.Schaller, J., Frei, S., Rohn, L. & Gilfedder, B. S. Amorphous silica controls water storage capacity and phosphorus mobility in soils. Front. Environ. Sci. 8, 94 (2020).Article 

Google Scholar 
27.Neu, S., Schaller, J. & Dudel, E. G. Silicon availability modifies nutrient use efficiency and content, C:N: P stoichiometry, and productivity of winter wheat (Triticum aestivum L.). Sci. Rep. 7, 40829 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
28.Schaller, J. et al. Silicon increases the phosphorus availability of Arctic soils. Sci. Rep. 9, 449 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
29.Munné-Bosch, S. & Alegre, L. Die and let live: Leaf senescence contributes to plant survival under drought stress. Funct. Plant Biol. 31, 203–216 (2004).PubMed 
Article 

Google Scholar 
30.Joshi, S. et al. Improved wheat growth and yield by delayed leaf senescence using developmentally regulated expression of a cytokinin biosynthesis gene. Front. Plant Sci. 10, 1285 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
31.Rivero, R. M. et al. Delayed leaf senescence induces extreme drought tolerance in a flowering plant. Proc. Natl. Acad. Sci. 104, 19631–19636 (2007).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Chapman, E. A., Orford, S., Lage, J. & Griffiths, S. Capturing and selecting senescence variation in wheat. Front. Plant Sci. 12, 638738 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
33.Camp, P. J., Huber, S. C., Burke, J. J. & Moreland, D. E. Biochemical changes that occur during senescence of wheat leaves: I. Basis for the reduction of photosynthesis. Plant Physiol. 70, 1641–1646 (1982).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Lopes, M. S. & Reynolds, M. P. Stay-green in spring wheat can be determined by spectral reflectance measurements (normalized difference vegetation index) independently from phenology. J. Exp. Bot. 63, 3789–3798 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Rodriguez, D. et al. Modelling the response of wheat canopy assimilation to atmospheric CO2 concentrations. New Phytol. 150, 337–346 (2001).Article 

Google Scholar 
36.Moureaux, C. et al. Carbon balance assessment of a Belgian winter wheat crop (Triticum aestivum L.). Glob. Change Biol. 14, 1353–1366 (2008).ADS 
Article 

Google Scholar 
37.Gao, X. P., Zou, C. Q., Wang, L. J. & Zhang, F. S. Silicon improves water use efficiency in maize plants. J. Plant Nutr. 27, 1457–1470 (2004).CAS 
Article 

Google Scholar 
38.Agarie, S., Uchida, H., Agata, W., Kubota, F. & Kaufman, P. B. Effects of silicon on transpiration and leaf conductance in rice plants (Oryza sativa L.). Plant Prod. Sci. 1, 89–95 (1998).Article 

Google Scholar 
39.Dakora, F. D. & Nelwamondo, A. Silicon nutrition promotes root growth and tissue mechanical strength in symbiotic cowpea. Funct. Plant Biol. 30, 947–953 (2003).CAS 
PubMed 
Article 

Google Scholar 
40.Steudle, E. & Peterson, C. A. How does water get through roots?. J. Exp. Bot. 49, 775–788 (1998).CAS 

Google Scholar 
41.Gao, X., Zou, C., Wang, L. & Zhang, F. Silicon decreases transpiration rate and conductance from stomata of maize plants. J. Plant Nutr. 29, 1637–1647 (2006).CAS 
Article 

Google Scholar 
42.Vandegeer, R. K. et al. Silicon deposition on guard cells increases stomatal sensitivity as mediated by K+ efflux and consequently reduces stomatal conductance. Physiol. Plant. 171, 358–370 (2020).PubMed 
Article 
CAS 

Google Scholar 
43.Flexas, J. & Medrano, H. Drought-inhibition of photosynthesis in C3 plants: Stomatal and non-stomatal limitations revisited. Ann. Bot. 89, 183–189 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
44.Kuhla, J., Pausch, J. & Schaller, J. Effect on soil water availability, rather than silicon uptake by plants, explains the beneficial effect of silicon on rice during drought. Plant Cell Environ. (2020).45.Schaller, J. et al. Silicon accumulation in rice plant aboveground biomass affects leaf carbon quality. Plant Soil 444, 399–407 (2019).CAS 
Article 

Google Scholar 
46.IUSS_Working_Group_Wrb. (Food and Agriculture Organization of the United Nations Rome, 2015).47.Schüller, H. Die CAL-Methode, eine neue Methode zur Bestimmung des pflanzenverfügbaren Phosphates in Böden. Z. Pflanzenernähr. Bodenkd. 123, 48–63 (1969).Article 

Google Scholar 
48.Huth, V. et al. Divergent NEE balances from manual-chamber CO2 fluxes linked to different measurement and gap-filling strategies: A source for uncertainty of estimated terrestrial C sources and sinks?. J. Plant Nutr. Soil Sci. 180, 302–315 (2017).CAS 
Article 

Google Scholar 
49.Livingston, G. & Hutchinson, G. Enclosure-based measurement of trace gas exchange: applications and sources of error. Biog. Trace Gases Meas. Emiss. Soil Water 51, 14–51 (1995).
Google Scholar 
50.Hoffmann, M. et al. Detecting small-scale spatial heterogeneity and temporal dynamics of soil organic carbon (SOC) stocks: A comparison between automatic chamber-derived C budgets and repeated soil inventories. Biogeosciences 14, 1003–1019 (2017).ADS 
CAS 
Article 

Google Scholar 
51.Davidson, E., Savage, K., Verchot, L. & Navarro, R. Minimizing artifacts and biases in chamber-based measurements of soil respiration. Agric. For. Meteorol. 113, 21–37 (2002).ADS 
Article 

Google Scholar 
52.Kutzbach, L. et al. CO2 flux determination by closed-chamber methods can be seriously biased by inappropriate application of linear regression. Biogeosciences 4, 1005–1025 (2007).ADS 
CAS 
Article 

Google Scholar 
53.Langensiepen, M., Kupisch, M., van Wijk, M. T. & Ewert, F. Analyzing transient closed chamber effects on canopy gas exchange for optimizing flux calculation timing. Agric. For. Meteorol. 164, 61–70 (2012).ADS 
Article 

Google Scholar 
54.Webb, E. K., Pearman, G. I. & Leuning, R. Correction of flux measurements for density effects due to heat and water vapour transfer. Q. J. R. Meteorol. Soc. 106, 85–100 (1980).ADS 
Article 

Google Scholar 
55.Leiber-Sauheitl, K., Fuß, R., Voigt, C. & Freibauer, A. High greenhouse gas fluxes from grassland on histic gleysol along soil carbon and drainage gradients. Biogeosci. Discuss. 10 (2013).56.Adamsen, F. et al. Measuring wheat senescence with a digital camera. Crop Sci. 39, 719–724 (1999).Article 

Google Scholar 
57.Idso, S., Pinter, P. Jr., Jackson, R. & Reginato, R. Estimation of grain yields by remote sensing of crop senescence rates. Remote Sens. Environ. 9, 87–91 (1980).ADS 
Article 

Google Scholar 
58.Kandel, T. P., Elsgaard, L. & Lærke, P. E. Measurement and modelling of CO2 flux from a drained fen peatland cultivated with reed canary grass and spring barley. Gcb Bioenergy 5, 548–561 (2013).CAS 
Article 

Google Scholar 
59.Görres, C.-M., Kutzbach, L. & Elsgaard, L. Comparative modeling of annual CO2 flux of temperate peat soils under permanent grassland management. Agr. Ecosyst. Environ. 186, 64–76 (2014).Article 
CAS 

Google Scholar  More

Global statistical analysisOur first attempt to identify plausible effects of meteorological covariates on COVID-19 spread applied a comparative regression analysis. To this end, we focused on the exponential onset of the disease, as it is the epidemic phase that allows for a better comparison between countries or regions, without the confounding effect of intervention policies. We first determined, for each of the spatial units (either countries or NUTS (nomenclature of territorial units for statistics) 2 regions), the day in which 20 or more cumulative cases were officially reported. We then fitted the first-order polynomial function f(t) = x0 + rt for the next 20 days of log-transformed data, where t represents time (in days) and ({{x}_0}) is the value at initial condition t = 0. The r parameter can be understood as the exponential growth rate, and is then used to estimate the basic reproduction number (R0) using the estimated serial interval T for COVID-19 of 4.7 days53, such that R0 = 1 + rT (ref. 54). (We note that we are interested here in the relationship between the reproductive number and not in the actual inference of R0.) Once R0 was obtained for all our spatial units, we filtered our meteorological data to match the same fitting period (with a 10-day negative delay to account for an incubation and reporting lapse) for every spatial unit. To compute a single average of the meteorological variables per regional unit, we computed a weighted average on the basis of the population contribution of each grid cell to the total population of the region. We did so to have an aggregated value that would better represent the impact of these factors on the population transmission of COVID-19, as the same variation in weather in a high-density urban area is more likely to contribute to a change in population-level transmission than that of an unpopulated rural area. We then averaged the daily values of temperature and AH for each country and computed univariate linear models for each of these variables as predictors of R0. Given the somewhat arbitrary criteria to select the dates to estimate the R0 in each country, a sensitivity analysis was run to test the robustness of the regressions to changes in the related parameters. We tested 70 different combinations of two parameters: the total number of days used for the fit (18–27) and the threshold of cumulative COVID-19 cases used to select the initial day of the fit (15–45). We also calculated the weather averages by shifting the selected dates accordingly. Then, a linear model for each of the estimates was fitted for both T and AH. A summary of the distribution of parameter estimates (the regression slope coefficients and the R2 of the models) is shown in Extended Data Fig. 3.Bivariate time-series analysis with scale-dependent correlationsTo examine associations between cases and climate factors in more detail, SDC was performed on the daily time series of both COVID-19 incidence and a given meteorological variable. SDC is an optimal method for identifying dynamical couplings in short and noisy time series20,21. In general, Spearman correlations between incidence and a meteorological time series assess whether there is a monotonic relation between the variables. SDC analysis was specifically developed to study transitory associations that are local in time at a specified temporal scale corresponding to the size of the time intervals considered (s). The two-way implementation (TW-SDC) is a bivariate method that computes non-parametric Spearman rank correlations between two time series, for different pairs of time intervals along these series. Different window sizes (s) can be used to examine increasingly finer temporal resolution. The results are sensitive to the value of this window size, s, with expected significant and highest correlation values at the scale of the transient coupling between variables. Correlation values decrease in magnitude as window size increases, and averages are computed over too long a time interval. Values can also decrease and become non-significant for small windows when correlations are spurious. Here, the method was applied for windows of different length (from s = 75 to 14 days) and, despite a weekly cycle showing up in some cases for small s, results removing this cycle were robust. We therefore did not remove this cycle.The results are typically displayed in a figure with the following subplots: (1) the two time series, to the left and top of the matrix of correlation values, respectively; (2) the matrix or grid of correlation values itself in the center, with significant correlations colored in blue when positive and in red when negative, with rows and columns corresponding to the temporal localization of the moving window along the time series on the left and top, respectively; (3) a time series at the bottom, below this grid, with the highest significant correlations for a given time (vertically, and therefore for the variable that acts as the driver, here the meteorological time series). To read the results, one starts at the diagonal and moves vertically down from it to identify a given lag for which significant correlations are found (the closest to the main diagonal). In some of the SDC figures, the time intervals with high local correlations are highlighted with boxes. These intervals alternate with other ones (left blank) for which no significant correlation is found. All colored areas correspond to significance levels of at least P  fs/fr, where fs is the sampling rate and fr the minimum frequency. Another strategy is that M be large enough that the M-lagged vector incorporates the temporal scale of the time series that is of interest. The larger the M, the more detailed the resulting decomposition of the signal. In particular, the most detailed decomposition is achieved when the embedding dimension is approximately equal to half of the total signal length. A compromise must be reached, however, as a large M implies increased computation, and too large a value may produce mixing of components. SSA is especially well suited for separating components corresponding to different frequencies in nonlinear systems. Here, we applied it to remove the weekly cycle.MSDC analysisMSDC provides a scan of the SDC analyses over a range of different scales (here, S from 5 to 100 days at 5-day intervals), by selecting the maximum correlation values (positive or negative) closer to the diagonal. The goal is to consider the evolution of transient correlations at all scales pooled together in a single analysis. The MSDC plot displays time on the x axis and scale (S) on the y axis, and positive and negative correlations either jointly or separately. The rationale behind MSDC is that correlations at very small scales can occur by chance because of coincident similar patterns, but that as one moves up to larger scales (by increasing S), the correlation patterns that are spurious tend to vanish, whereas those reflecting mechanistic links increase in strength. This increase in correlation values should occur up to the real scale of interaction, decreasing afterwards. By ‘real’, we mean here the temporal scale covering the extent of the interaction between the driver and the response process (in this case, the response of disease transmission to a given climate factor). Thus, continuity of the same sign correlations together with transitions to larger values are indicative of causal effects, whereas the rapid vanishing of small-scale significant correlations signals spurious ones.Process-based modelDescriptionThe dynamical model is a discrete stochastic model that incorporates seven different compartments: S, E, I, C, Q, R and D. The model structure is illustrated in Fig. 4. The transition probabilities of the stochastic model are based on the corresponding rates of the transitions between classes in the deterministic (mean-field) model (specified in Fig. 4b). These probabilities are defined as follows. P(e) = (1.0 − exp(−β dt)) is the probability of infection exposure of the susceptible class, where β = (1/N)(βII + βQQ) is the infection rate (of the deterministic model). P(i) = (1.0 − exp(−γ dt)) is the probability that an new exposed individual becomes infectious, where γ denotes the incubation rate. P(r) = (1.0 − exp(−Λ dt)) is the recovery probability, where λ0(1 − exp(λ1t)) is the (deterministic) recovery rate. P(p) = (1.0 − exp(−α dt)) is the protection probability, where α = α0exp(α1t). P(d) = (1.0 − exp(−K dt)) is the mortality probability, with K = k0exp(k1t). P(re) = (1.0 − exp(−τ dt)) is the release probability from confinement, where τ = τ0exp(τ1t). Finally, P(q) = (1.0 − exp(−δ  dt)) is the detection probability, where δ is the quarantine rate (for example, at which infected individuals are isolated from the rest of the population).In the model, both infected non-detected and infected detected individuals can infect susceptible ones. In the model incorporating temperature in the transmission rate, the respective values of βI and βQ are calculated as follows:$${beta }_{I}(t)={beta }_{I},T_{mathrm{inv}}(t);quad {beta }_{Q}(t)={beta }_{Q},T_{mathrm{inv}}(t)$$where (T_{mathrm{inv}}=fleft(frac{1-T(t)}{bar{T}}right)), with (bar{T}) corresponding to the overall mean of the temperature time series and f(·) to a Savitzky–Golay filter, used to smooth the temperature series with a window size of 50 data points and a polynomial order of 3. When the infection rate is constant, we simply omit the temperature term. For further comparison, in a third model, β is specified with a sinusoidal function of period equal to 12 months and an estimated phase.The number of individuals transitioning from compartment i to j at time t are determined by means of binomial distributions P(Xi,P(y)), where Xi corresponds to one of the compartments S, E, I, Q, R, D, C, and P(y) to the respective transition probability defined above. Thus,

e(t) = P(S(t), P(e)), new exposed individuals at time t

p(t) = P(S(t), P(p)), protected individuals at time t

i(t) = P(E(t), P(i)), new infected not detected individuals at time t

q(t) = P(I(t), P(q)), new infected and detected individuals at time t

r(t) = P(Q(t), P(r)), total recovered individuals at time t

d(t) = P(Q(t), P(d)), total dead individuals at time t

re(t) = P(C(t), P(re)), individuals released from confinement at time t

Then, the final dynamics are given by the following equations:$$S(t)=S(t-{rm{d}}t)-e(t)-p(t)+re(t)$$$$E(t)=E(t-{rm{d}}t)+e(t)-i(t)$$$$I(t)=I(t-{rm{d}}t)+i(t)-q(t)$$$$Q(t)=Q(t-{rm{d}}t)+q(t)-r(t)-d(t)$$$$R(t)=R(t-{rm{d}}t)+r(t)$$$$D(t)=D(t-{rm{d}}t)+d(t)$$$$C(t)=C(t-{rm{d}}t)+p(t)-re(t)$$CalibrationThe model was implemented using Python and calibrated by means of the least squares algorithm of the scipy library. The error function minimized with this algorithm was obtained from the normalized residuals on the basis of total cases (Q + R + D) and deaths (D).To search parameter space, we ran 100 calibrations starting from different initial choices of parameter combinations. The tolerance for termination in the change of the cost function was set to 1 × 10−10. Tolerance for termination by the norm of the gradient was also set to 1 × 10−10, and the tolerance for termination by the change of the independent variables was set to 1 × 10−10. The solver was the lsmr method (which is suitable for problems with sparse and large Jacobian matrices) with a differential step of 1 × 10−5. With this configuration, each fitting run usually converged after ~500 iterations.ValidationTo compare the model including an effect of T in the transmission rate to those without it, we calculated the chi-square, Akaike information criterion (AIC) and Bayesian information criterion (BIC) indices for the residuals obtained from the optimization process. The resulting values are shown in Supplementary Table 1.Our choice of T to modulate the infection rate (β) instead of AH underlies the fact that the temporal dynamics of both factors roughly follow the same shape, with the advantage that T shows less oscillatory behavior than AH. This fact adds stability to the model when the inverse relationship is used in the calculation of β (Supplementary Information). This selection is further reinforced by the results from the SDC analyses, which yielded larger correlations for temperature, even when penalizing for the larger autocorrelation structure.Our choice to modulate β using T instead of AH follows from the fact that the temporal dynamics of both climate variables present roughly the same shape, with the advantage that T exhibits weaker oscillations. This less fluctuating pattern provides stability to the model fitting when the inverse relationship is used in the calculation of β (Supplementary Information). Additionally, the transient correlations obtained with SDC yielded higher values for T than for AH (even when accounting for concurrent levels of autoregression in the two variables). More

Study sites and samplingThe study was conducted in Garhwal region (Western Himalaya) from 2016 to 2017 at eight Rhododendron arboreum rich areas in four hill districts (Chamoli, Tehri, Pauri and Rudraprayag). Voucher specimen of Rhododendron arboreum collected and have been deposited in the Herbarium, Botany department, HNB Garhwal University (specimen no. GUH 8510)6. Identification of R. arboreum has been done through A Field Guide book authored by Rai et al.7. Since it is a wild species and flowers have been collected for our research and field study under the permission from competent authority of State Forest Department, Govt. of Uttarakhand. According to IUCN’s Red List Categories and Criteria, globally Rhododendron arboreum comes under Least Concern (LC) category8. These sites are situated between 30°08′47″ to 30°24′06″ N latitude and 78°25′05″ to 79°12′39″ E longitude with altitudes from 1820 m asl in Nandasain and 2270 m asl in Jadipani (Table 1; Fig. 1). All sites were well stocked (mean stand density ≥ 500 tree/ha) with Rhododendron arboreum trees mixed with Quercus leucotrichophora. We referred these resource rich sites as R. arboreum habitats (Table 1). Stratified random sampling method (i.e. stand density and CBH class’s strata) were carried out these eight sites. Total sampled area 0.2 ha in each site; two sample plots (size of each plot is 0.1 ha or 31.62 × 31.62 m) nested within 0.2 ha in each site were laid out for trees enumeration. Sample size (number of R. arboreum tree) for a total population in each site were 166 in Phadkhal; 110 in Khirsu; 104 in Khadpatiya; 166 in Ghimtoli; 80 in Jadipani; 74 in Ranichauri; 74 in Nandasain and 96 in Nauti. Out of the standing trees in sample plots, flower bearing trees were 96 in Phadkhal; 90 in Khirsu; 102 in Khadpatiya; 126 in Ghimtoli; 64 in Jadipani; 58 in Ranichauri; 68 in Nandasain and 82 in Nauti, and without flower or smaller trees were 70 in Phadkhal; 20 in Khirsu; 02 in Khadpatiya; 40 in Ghimtoli; 16 in Jadipani; 16 in Ranichauri; 06 in Nandasain and 14 in Nauti. The individuals of all tree species in each plot were recorded along with their CBH (circumference at breast height, 1.3 m above from the ground). Individuals were categorized as mature trees (≥ 31 cm CBH), saplings (11–30 cm CBH) and seedlings (≤ 10 cm CBH)9. Further all the tree individuals have been grouped into 8 CBH classes: (A) 5–15 cm, (B) 16–25 cm, (C) 26–35 cm, (D) 36–45 cm, (E) 46–55 cm, (F) 56–65 cm, (G) 66–75 cm, (H) 76–85 cm. Recorded data were used for the analysis of density10.Table 1 Physical characteristics of study sites in four districts of Garhwal region.Full size tableFigure 1Locations of Rhododendron arboreum study sites in Garhwal region (ARC GIS software 10.5 version was used for map preparation. The map was created by Mr. Raman Patel, Research scholar, Dept. of Geology, HNB Garhwal University, Srinagar, Uttarakhand, India).Full size imageFlower yield estimationFlower yield (kg/tree) was estimated during full bloom (flowering season/harvest season February–April 2017). In each sample plot, numbers of flower bearing trees varied from 29–63 trees/0.1 ha. At each site of 0.2 ha sample plot, total 40 trees, 05 flower bearing trees in each of the 08 CBH classes were marked for estimation of flower yield. The number of main branches, the number of sub- branches/offshoots per main branches (i.e. average per five randomly selected main branches per tree), and the amount of flower per sub-branches/offshoot (i.e. the average per five offshoots from the low, middle and upper canopy of each tree) were counted form marked individuals. This way flower yield/tree was calculated9,11,12.The flowers from all CBH classes in each site were mixed and weighted in 5 lots of 1 kg each. The number of flowers in each lot was then counted and the mean value (400.0 ± 9.56) was considered as a standard for conversion into kilograms. Based on this conversion flower yield kg/tree was obtained. Flower yield data were pooled and mean yield (kg/tree) for each CBH class (A–H) calculated. For each site, flower yield in kg/0.2 ha was obtained by multiplying flower yield/tree by the density of flower bearing trees/0.2 ha. The total yield kg/ha for each site was calculated as total yield = (yield/ha) × density of flower bear trees/ha9,11,12.Extraction/harvesting and marketing trendsFlower extraction and collection were totally dependent on market availability and accessibility of site; one of the selected sites (Ranichauri) was easily accessible, while Phadkhal, Khirsu and Jadipani were moderately accessible. Khadpatiya, Ghimtoli, Nandasain and Nauti sites were far-flung from market (Table 1). The highest extraction was recorded between second week of February and first week of April. During this period, data was obtained for three consecutive days at each site.Questionnaire based survey was carried out in selective forest fringe villages. Across the sites, total sixteen villages were selected for questionnaire survey, three villages each in Jadipani, Ranichauri, Nandasain and Nauti sites, while one village each in Phadkhal, Khirsu, Khadpatiya and Ghimtoli. In each village 15 families were randomly chosen for semi-structured questionnaire survey.Considering the market availability for trading of the R. arboreum flower products, Nandasain and Nauti sites are located nearest to local market whereas Khadpatiya and Ghimtoli sites are farthest from local market. As far as the access to resources is concerned, four sites represent open and easy access to resource and four sites represent open and moderate access of resource (Table 1). During questionnaire survey, villagers were asked about the number of persons involved in resource collection (hereafter referred as collectors), age of collectors, timing of collection (early morning and late evening) etc. Ten individuals in each group (adults and children) were randomly interviewed on their harvest load to generate data on the average collection per individual, the number of days spent in flower collection, and the total income generated through this activity.Squash/juice making factories are generally located nearby urban centers; local NGOs and small entrepreneurs are engaged in this work. These peoples purchase flower from the collectors or middleman for preparation of value product (squash). Collectors of each families (varied from n = 15 in Nanadasain to n = 31 in Jadipani) and buyers (n = 5 each site) were contacted to obtain information on the benefits accrued. The income values are given in Indian rupees (USD 1 = Rs. 68.00, 2017 exchange rates). Projections of potential (probable/-could generate) income (with flower processed into juice or squash) were made. The involvement of rural inhabitants as flowers collectors and the income that subsequently accrued (within a 10 km radius of fringe area) was also analyzed for sixteen villages across the sites. One adult member from each household was contacted in a village to collect information on involvement of flower collection/extraction.Juice/squash preparation methods and value-added productsThe collected flowers are graded for their size and healthiness and the stamens are separated from petals by laborers in the juice processing unit. Petals are cleaned washed with tap water and grinded into small pieces. The petal mass is retained in the water and then boiled for one hour. The slurry (aqueous solution) obtained in this process is left at room temperature for cooling and when it get cold, filtered through linen cloth. The filtrate solution is the pure juice of the flower. For the preparation of squash from the pure juice, about 2 kg of sugar is boiled in one liter of water. Further one liter of pure juice and a small quantity of citric acid (10 g/2 kg sugar) are added to this solution. The mixture is boiled again for 30 min and then left to cool at room temperature13. The obtained solution known as squash is then filtered through linen cloth and stored into containers and bottles for marketing. For long term storage and good test and aroma small amount of sodium benzoate and vanilla or kawra is also mixed in the squash.Cost–benefit analysis of value- added productsThe cost–benefit analysis of value added products prepared from the R. arboreum flowers was calculated in Rs./day which includes labour charges of workers involved in flower collection and materials/items required for preparation of different value added products viz: sugar, preservatives, essence, plastic containers/bottles, packaging materials etc. Labour charge was calculated on the basis of existing daily wages as per market rates. The monetary output was calculated as per the current market rates of the products (Table 2). The cost- benefit analysis of the squash product prepared from the flowers was calculated as Rs./day which includes: (i) Man days incumbent for the flowers extraction from the forest and for the preparation of squash product, (ii) Essential items such as sugar, preservatives etc. and their monetary equivalents, (iii) The total quantity of squash product and their monetary equivalents.Table 2 Market cost in rupees (Rs.) of essential commodity in the preparation of R. arboreum juice/squash in Garhwal region.Full size tableStatistical analysisData failed to meet the assumption of normality (Shapiro–Wilk test) as well as homogeneity (Levene statistic); therefore, a non-parametric test (i.e. Independent–Samples Kruskal–Wallis test) was applied for one-way ANOVA. However, to find the interaction of site and cbh on flower production (yield), the same data set was subjected to two-way analysis using univariate analysis. To find if (?) flower yield depends on tree diameter or not, data of actual cbh and flower yield per tree were used to determine a correlation (Pearson Correlation Coefficient) between them. In case of correlation found significant then regression equation was developed to predict flower production based on tree diameter. All analysis were performed using IBM-SPSS 16.0 version14.Ethics approval and consent to participateAll necessary approval, free prior informed consent, permit, and certification were secured. This was done to adhere to the ethical standards of human participation in scientific research. This study was approved by Research and Consultancy Cell (Ethics Committee) of HNB Garhwal University, Srinagar Garhwal, Uttarakhand, India. All the methods were performed in accordance with the relevant guidelines and regulations. More

1.Dairy Production and Products: Milk and Milk Products (FAO, 2013); http://www.fao.org/dairy-production-products/production/dairy-animals/cattle/en/2.Background: Corn and Other Feedgrains (USDA ERS, 2018); https://www.ers.usda.gov/topics/animal-products/dairy/background/3.National Agricultural Statistics Service (US Department of Agriculture); https://www.nass.usda.gov/index.php4.Capper, J. L., Cady, R. A. & Bauman, D. E. The environmental impact of dairy production: 1944 compared with 2007. J. Anim. Sci. 87, 2160–2167 (2009).CAS 
Article 

Google Scholar 
5.Niles, M. T. & Wiltshire, S. Tradeoffs in US dairy manure greenhouse gas emissions, productivity, climate, and manure management strategies. Environ. Res. Commun 1, 075003 (2019).Article 

Google Scholar 
6.Field, T. G. & Taylor, R. E. Scientific Farm Animal Production: An Introduction, Eleventh Edition (Pearson, 2018).7.Fuquay, J. W. Heat stress as it affects animal production. J. Anim. Sci. 52, 164–174 (1981).CAS 
Article 

Google Scholar 
8.St-Pierre, N. R., Cobanov, B. & Schnitkey, G. Economic losses from heat stress by US livestock industries. J. Dairy Sci. 86, E52–E77 (2003).Article 

Google Scholar 
9.Kadzere, C. T., Murphy, M. R., Silanikove, N. & Maltz, E. Heat stress in lactating dairy cows: a review. Livest. Prod. Sci. 77, 59–91 (2002).Article 

Google Scholar 
10.Bouraoui, R., Lahmar, M., Majdoub, A., Djemali, M. & Belyea, R. The relationship of temperature–humidity index with milk production of dairy cows in a Mediterranean climate. Anim. Res. 51, 479–491 (2002).Article 

Google Scholar 
11.West, J. W. Effects of heat-stress on production in dairy cattle. J. Dairy Sci. 86, 2131–2144 (2003).CAS 
Article 

Google Scholar 
12.Vitali, A. et al. Seasonal pattern of mortality and relationships between mortality and temperature–humidity index in dairy cows. J. Dairy Sci. 92, 3781–3790 (2009).CAS 
Article 

Google Scholar 
13.Pragna, P. et al. Heat stress and dairy cow: impact on both milk yield and composition. Int. J. Dairy Sci. 12, 1–11 (2017).CAS 
Article 

Google Scholar 
14.Hoffmann, I. Climate change and the characterization, breeding and conservation of animal genetic resources. Anim. Genet. 41, 32–46 (2010).Article 

Google Scholar 
15.Qi, L., Bravo-Ureta, B. E. & Cabrera, V. E. From cold to hot: a preliminary analysis of climatic effects on the productivity of Wisconsin dairy farms. AgEconSearch https://doi.org/10.22004/ag.econ.172411 (2014).16.Bohmanova, J., Misztal, I. & Cole, J. B. Temperature–humidity indices as indicators of milk production losses due to heat stress. J. Dairy Sci. 90, 1947–1956 (2007).CAS 
Article 

Google Scholar 
17.Field, C. B. et al. Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation (IPCC, 2021); https://www.ipcc.ch/report/managing-the-risks-of-extreme-events-and-disasters-to-advance-climate-change-adaptation/18.Mueller, N. D. et al. Cooling of US Midwest summer temperature extremes from cropland intensification. Nat. Clim. Chang. 6, 317–322 (2016).ADS 
MathSciNet 
Article 

Google Scholar 
19.Seneviratne, S. I., Donat, M. G., Mueller, B. & Alexander, L. V. No pause in the increase of hot temperature extremes. Nat. Clim. Chang. 4, 161–163 (2014).ADS 
Article 

Google Scholar 
20.Dairy 2014: Dairy Cattle Management Practices in the United States, 2014 (USDA, APHIS, NAHMS, 2016); https://www.aphis.usda.gov/animal_health/nahms/dairy/downloads/dairy14/Dairy14_dr_PartI_1.pdf21.Mondaca, M. R. & Cook, N. B. Modeled construction and operating costs of different ventilation systems for lactating dairy cows. J. Dairy Sci. 102, 896–908 (2019).CAS 
Article 

Google Scholar 
22.Ferreira, F. C., Gennari, R. S., Dahl, G. E. & De Vries, A. Economic feasibility of cooling dry cows across the United States. J. Dairy Sci. 99, 9931–9941 (2016).CAS 
Article 

Google Scholar 
23.Hayhoe, K. et al. Emissions pathways, climate change, and impacts on California. Proc. Natl Acad. Sci. USA 101, 12422–12427 (2004).ADS 
CAS 
Article 

Google Scholar 
24.Klinedinst, P. L., Wilhite, D. A., Hahn, L. G. & Hubbard, K. G. The potential effects of climate change on summer seasonal dairy cattle milk production and reproduction. Clim. Chang. 23, 21–36 (1993).ADS 
Article 

Google Scholar 
25.Mauger, G., Bauman, Y., Nennich, T. & Salathé, E. Impacts of climate change on milk production in the United States. Prof. Geogr. 67, 121–131 (2015).Article 

Google Scholar 
26.Key, N. & Sneeringer, S. Potential effects of climate change on the productivity of U.S. dairies. Am. J. Agric. Econ. 96, 1136–1156 (2014).Article 

Google Scholar 
27.Ortiz-Bobea, A., Knippenberg, E. & Chambers, R. G. Growing climatic sensitivity of U.S. agriculture linked to technological change and regional specialization. Sci. Adv. 4, eaat4343 (2018).ADS 
Article 

Google Scholar 
28.Butler, E. E., Mueller, N. D. & Huybers, P. Peculiarly pleasant weather for US maize. Proc. Natl Acad. Sci. USA 115, 11935–11940 (2018).CAS 
Article 

Google Scholar 
29.Lobell, D. B., Schlenker, W. & Costa-Roberts, J. Climate trends and global crop production since 1980. Science 333, 616–620 (2011).ADS 
CAS 
Article 

Google Scholar 
30.Tigchelaar, M., Battisti, D. S., Naylor, R. L. & Ray, D. K. Future warming increases probability of globally synchronized maize production shocks. Proc. Natl Acad. Sci. U. S. A. 115, 6644–6649 (2018).ADS 
Article 

Google Scholar 
31.PRISM Climate Data (Oregon State Univ., 2019); http://www.prism.oregonstate.edu/32.Daly, C. et al. Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int. J. Climatol. https://doi.org/10.1002/joc.1688 (2008).33.National Research Council. Nutrient Requirements of Dairy Cattle, Seventh Revised Edition (National Academies Press, 2001).34.Auldist, M. J., Walsh, B. J. & Thomson, N. A. Seasonal and lactational influences on bovine milk composition in New Zealand. J. Dairy Res. 65, 401–411 (1998).CAS 
Article 

Google Scholar 
35.Lobell, D. B. Climate change adaptation in crop production: beware of illusions. Glob. Food Sec. 3, 72–76 (2014).Article 

Google Scholar 
36.Mukherjee, D., Bravo-Ureta, B. E. & De Vries, A. Dairy productivity and climatic conditions: econometric evidence from South-eastern United States. Aust. J. Agric. Resour. Econ. 57, 123–140 (2013).Article 

Google Scholar 
37.Milk Cost of Production Estimates: Cost-of-Production Estimates-2016 Base (USDA ERS, 2021); https://www.ers.usda.gov/data-products/milk-cost-of-production-estimates/milk-cost-of-production-estimates/#Milk38.Liang, X. Z. et al. Determining climate effects on US total agricultural productivity. Proc. Natl Acad. Sci. USA 114, E2285–E2292 (2017).CAS 
Article 

Google Scholar 
39.Malikov, E., Miao, R. & Zhang, J. Distributional and temporal heterogeneity in the climate change effects on U.S. agriculture. J. Environ. Econ. Manage. 104, 102386 (2020).Article 

Google Scholar 
40.MacDonald, J. M., Law, J. & Mosheim, R. Consolidation in U.S. Dairy Farming Economic Research Report No. 274 (ERS, USDA, 2020); https://www.ers.usda.gov/publications/pub-details/?pubid=9890041.Hemme, T. & Otte, J. Pro-Poor Livestock Policy Initiative Status and Prospects for Smallholder Milk Production a Global Perspective (Food and Agriculture Organization of the United Nations, 2010).42.Osei-Amponsah, R. et al. Heat stress impacts on lactating cows grazing Australian summer pastures on an automatic robotic dairy. Animals 10, 869 (2020).Article 

Google Scholar 
43.Chang-Fung-Martel, J., Harrison, M. T., Rawnsley, R., Smith, A. P. & Meinke, H. The impact of extreme climatic events on pasture-based dairy systems: a review. Crop Pasture Sci 68, 1158 (2017).Article 

Google Scholar 
44.Livestock Hot Weather Stress. Operations Manual (NOAA, 1976); https://scirp.org/reference/referencespapers.aspx?referenceid=191321645.Pinheiro J., Bates D., Debroy S. S. D. Linear and nonlinear mixed effects models, R package nlme version 3.1-152 (2021).46.Conley, T. G. GMM estimation with cross sectional dependence. J. Econom. 92, 1–45 (1999).MathSciNet 
Article 

Google Scholar 
47.Borchers, H. W. pracma: practical numerical math functions, version 2.2.9.1–393 (2019).48.Colin Cameron, A., Gelbach, J. B. & Miller, D. L. Robust inference with multiway clustering. J. Bus. Econ. Stat. 29, 238–249 (2011).MathSciNet 
Article 

Google Scholar 
49.Zeileis, A., Köll, S. & Graham, N. Various versatile variances: an object-oriented implementation of clustered covariances in R. J. Stat. Softw. https://doi.org/10.18637/jss.v095.i01 (2020). More

Toxicity of spirotetramat to F. occidentalis eggsThe indoor toxicity of spirotetramat to 0-h-, 12-h-, and 24-h-old eggs of F. occidentalis using egg dipping method and leaf dipping method was shown in Table 1 after egg hatching was observed for 144 h. The results suggested that the LC50 value gradually decreased as the egg age increased. The two methods have the same trend. And in the leaf dipping method, according to the confidence limits analysis, 0-h-old eggs are significantly more sensitive to spirotretramat than 24-h-old eggs.Table 1 Toxicity of spirotetramat to F. occidentalis eggs.Full size tableEgg external shape observationsExternal morphology of normally developed isolated 0-h-old eggs of F. occidentalis in the control treatment (Fig. 1a) were compared with those treated with spirotetramat (Fig. 1b, c). After spirotetramat treatment, some eggs appeared darker, yellowish-brown, and with embryonic development abnormalities (Fig. 1b); some of the egg embryo cells treated with spirotetramat appeared atrophied, and there were obvious gaps between the outer and the inner egg embryo cells, compared with the control treatment (Fig. 1c). Some eggs ruptured at the top after the egg shell was treated with spirotetramat, and the internal egg embryo cells flowed out, thus failing to form a complete embryo (Fig. 1d); the full egg embryo cells formed a control treatment.Figure 1The effect of spirotetramat on external morphology of isolated 0-h-old eggs of F. occidentalis. (a) 0-h-old isolated eggs of normal developing thrips; (b) yellowish-brown, developmentally deformed eggs after spirotetramat treatment; (c) eggs with shrunken oocytes after spirotetramat treatment; (d) eggs with apical rupture of the eggshell after spirotetramat treatment.Full size imageThe effect of spirotetramat on external morphology of live F. occidentalis 0-h-old eggs (Fig. 2b, c) was compared with normal development in the control treatment (Fig. 2a). Similar to the effect on external morphology of isolated 0-h-old eggs of F. occidentalis, the eggs treated with spirotetramat also showed abnormal embryonic development (Fig. 2b), egg embryo cell atrophy (Fig. 2c) and the phenomenon of rupture of the egg shell and outflow of embryo cells.Figure 2The effect of spirotetramat on external morphology of living 0-h-old eggs of F. occidentalis. (a) External morphology of live 0-h-old eggs of normal developing thrips; (b) developmentally deformed eggs after spirotetramat treatment; (c) eggs with shrunken oocytes after spirotetramat treatment; (d) eggs with apical rupture of the eggshell after spirotetramat treatment.Full size imageCompared with the control (Fig. 3a), the isolated 24-h-old eggs treated with spirotetramat (Fig. 3b) did not show obvious external morphological differences. After spirotetramat treatment, the eggs were still white and plump and with no embryonic deformities, egg cell atrophy or egg shell rupture, and could still develop normally. There were clear red eye spots on the head end, and embryo movement was clearly seen under the super-depth microscope (Fig. 3b). Similarly, the live 24-h-old eggs in the control treatment (Fig. 4a) and those treated with spirotetramat (Fig. 4b) showed no obvious external morphological differences, and the eggs developed normally.Figure 3The effect of spirotetramat on external morphology of isolated 24-h-old eggs of F. occidentalis. (a) Normally developing 24-h-isolated eggs in the control group; (b) 24-h-old eggs after spirotetramat treatment.Full size imageFigure 4The effect of spirotetramat on external morphology living of 24-h-old eggs of F. occidentalis. (a) Normally developing 24-h-live eggs in the control group; (b) live 24-h-old eggs after spirotetramat treatment.Full size imageEffect of egg hatchingThe 0-h-old eggs of F. occidentalis treated with spirotetramat did not hatch normally, and the mortality rate was 100% (Fig. 5). Among them, 77 eggs eventually showed rupture of the egg shell, the internal egg embryo cells flowed out and they did not hatch; 23 eggs showed no changes in external morphology, but did not hatch after continuous observation for 144 h, and showed no developmental phenomena such as embryo movement under a super-depth microscope, which was regarded as egg death. In the control treatment, 96 eggs hatched normally, and only six eggs did not rupture but did not hatch normally and were considered dead.Figure 5Effect of spirotetramat on hatching rate of F. occidentalis 0-h-old eggs.Full size imageThere was no significant difference between the 24-h-old eggs of F. occidentalis treated with spirotetramat and the control treatment. After spirotetramat treatment, 93 eggs hatched normally, and the shells of seven eggs were not ruptured (Fig. 6). Any eggs not hatched after 144 h of continuous observation were considered dead. In the control treatment, 95 eggs hatched normally and five eggs did not rupture but did not hatch normally, and so were considered dead.Figure 6The effect of spirotetramat on hatching rate of F. occidentalis 24-h-old eggs.Full size imageSEM observationsThe F. occidentalis eggs in the control treatment were kidney-shaped, with regular egg morphology, smooth surfaces and no folds or protrusions (Fig. 7a). At 24 h after spirotetramat treatment, part of the egg shells treated with spirotetramat had fallen off the chorion, and the embryonic material was exposed (Fig. 7b). The surface of the egg shell was uneven and severely wrinkled (Fig. 7c). The pores of some eggs treated with spirotetramat were sunken down and shrunken (Fig. 7d). Spirotetramat treatment of 0-h-old eggs affect clearly egg shells, resulting in shrinkage of egg shells, ovarian depression and egg malformations, and destroyed the egg shell structure. Thus, normal embryonic development was affected, and disrupted normal hatching.Figure 7The effect of spirotetramat on the surface of egg shells of F. occidentalis 0-h-old eggs. (a) 0-h-old eggs in the control treatment; (b) eggs shells were shed 24 h after treatment with spirotetramat; (c) the surface of the egg shell was uneven and severely wrinkled; (d) the pores of some eggs were sunken down and shrunken after treatment with spirotetramat.Full size imageThe shells of eggs treated with spirotetramat (Fig. 8b) showed no significant difference compared with controls (Fig. 8a). The eggs of the two groups of F. occidentalis were regular in shape, with smooth surfaces and without folds or protrusions. Thus, development of 24-h-old eggs showed some resistance to spirotetramat. Spirotetramat did not destroy the egg shell surface structure of 24-h-old eggs, indicating a high resistance to spirotetramat.Figure 8The effect of spirotetramat on the egg shell surface of F. occidentalis 24-h-old eggs. (a) 24-h-old eggs in the control treatment; (b) 24-h-old eggs in the spirotetramat treatment.Full size imageTEM observationsThe TEM observations showed that the egg structure of the control treatment was complete, the protoplasm and yolk were clearly observed inside the egg and the yolk was packed in the void of the protoplasm network (Fig. 9a). The egg shell structure was clear, and the outer and inner egg shell were clearly observed, as was the yolk membrane and the dense layer structure (Fig. 9c). Eggs treated with spirotetramat were flocculent, and no clear internal material was observed. The protoplasm and yolk structure were blurred, and flocculation in the protoplasm appeared to agglomerate and form blocks (Fig. 9b). The egg shell structure was unclear, and no clear outer egg shell, inner egg shell, yolk membrane and lamellar structures were observed. The egg shell was also filled with many flocs (Fig. 9d).Figure 9The effect of spirotetramat on the structure of F. occidentalis 0-h-old eggs. (a) and (b) 0-h-old eggs in the control treatment; (c) and (d) 0-h-old eggs in the spirotetramat treatment.Full size imageEffect on embryonic developmentThe initial eggs of the control group were kidney-shaped, white and full of vitellin (Fig. 10a). After 12 h of development, the eggs were larger and of oval shape (Fig. 10b). After 24 h of development, the egg had increased in volume, a partially transparent region appeared in the embryo and the embryo had transparent top follicles (Fig. 10c). After 36 h of development, some yolk granules disappeared and eggs became smooth and translucent (Fig. 10d). After 48 h of development, the insect outline was visible within the egg, a pair of antennae were visible on the head and a red eye point was clearly observed on the head during the blastokinesis phenomenon (Fig. 10e). After 60 h of development, embryo color deepened, the eye point was clearer and the head, femur, tibia and tarsus were clear (Fig. 10f). After 72 h of development, the egg shell began to break at the head, the tail constantly jittered, internal fluid flowed and the larva hatched from the top of the egg (Fig. 10g).Figure 10The embryonic development process of control 0-h-old eggs of F. occidentalis. (a) Control initial eggs; (b) eggs after 12 h of development; (c) eggs after 24 h of development; (d) eggs after 36 h of development; (e) eggs after 48 h of development; (f) eggs after 60 h of development; (g) eggs hatching as larvae after 72 h of development.Full size imageEggs of F. occidentalis were initially white, kidney-shaped and full of vitellin (Fig. 11a). Following treatment with spirotetramat, after 12 h of development, the eggs became large and oval, and the embryo was a pale brown color (Fig. 11b). After 24 h of development, color of the egg deepened to dark brown. There was a gap between the egg and the egg shell, and a small amount of spillage appeared at the end of the egg (Fig. 11c). After 36 h of development, the egg shell ruptured, material flowed out of the egg and embryo development did not proceed (Fig. 11d).Figure 11Effects of spirotetramat on development of 0-h-old eggs of F. occidentalis. (a) Frankliniella occidentalis initial eggs; (b) eggs developing 12 h after spirotetramat treatment; (c) eggs developing 24 h after spirotetramat treatment; (d) eggs developing 36 h after spirotetramat treatment.Full size imageIn the control treatment, the egg volume increased at 24 h, the embryo had a partially transparent area and there was a transparent follicle on the top of the embryo (Fig. 12a). After 12 h of development, some of the yolk particles disappeared and the egg body was smooth and translucent (Fig. 12b). After 24 h of development, the body outline, a pair of antennae and red eye spots were visible, and there was obvious embryo movement (Fig. 12c). At 36 h of development, the head, leg segments, tibia and tarsus were apparent (Fig. 12d). After 48 h of development, the embryo moved violently, internal body fluid flowed and the larva was ready to hatch (Fig. 12e). After 60 h of development, the larva emerged from its shell (Fig. 12f).Figure 1224-h-old eggs embryo development of F. occidentalis in control treatment. (a) Control 24-h-old eggs; (b) eggs after 12 h of development; (c) eggs after 24 h of development; (d) eggs after 36 h of development; (e) eggs after 48 h of development; (f) eggs hatching as larvae after 60 h of development.Full size imageThe 24 h old eggs of F. occidentalis showed enlarged volume, and there were transparent follicles on the top of the embryo (Fig. 13a). After 24 h eggs were treated with spirotetramat, they developed for 12–36 h, and the developmental status was the same as that of the control. The embryos developed normally, and there was no egg body discoloration or egg shell rupture (Fig. 13b–d). After 48 h of development, hairy scales appeared on the surface of the egg shell, and the egg body turned yellowish-brown in color, but the egg shell was not broken and no internal material overflow was seen (Fig. 13e). After 60 h of development, larvae hatched normally (Fig. 13f).Figure 13The effect of spirotetramat on embryonic development of F. occidentalis 24-h-old eggs. (a) Frankliniella occidentalis 24-h-old eggs; (b) eggs developing 12 h after spirotetramat treatment; (c) eggs developing 24 h after spirotetramat treatment; (d) eggs developing 36 h after spirotetramat treatment; (e) eggs developing 48 h after spirotetramat treatment; (f) eggs hatching as larvae after 60 h of development.Full size image More

1.Gartner, T., Mulligan, J., Schmidt, R. & Gunn, J. Natural Infrastructure (World Resources Institute, 2013).2.McDonald, R. I. et al. Water on an urban planet: urbanization and the reach of urban water infrastructure. Glob. Environ. Change 27, 96–105 (2014).Article 

Google Scholar 
3.Vorosmarty, C. J. et al. Global threats to human water security and river biodiversity. Nature 467, 555–561 (2010).CAS 
Article 

Google Scholar 
4.Grill, G. et al. Mapping the world’s free-flowing rivers. Nature 569, 215–221 (2019).CAS 
Article 

Google Scholar 
5.Tessler, Z. D. et al. Profiling risk and sustainability in coastal deltas of the world. Science 349, 638–643 (2015).CAS 
Article 

Google Scholar 
6.Palmer, M. A. Water resources: beyond infrastructure. Nature 467, 534–535 (2010).CAS 
Article 

Google Scholar 
7.Michalak, A. M. Study role of climate change in extreme threats to water quality. Nature 535, 349–350 (2016).CAS 
Article 

Google Scholar 
8.McDonald, R. I., Weber, K. F., Padowski, J., Boucher, T. & Shemie, D. Estimating watershed degradation over the last century and its impact on water-treatment costs for the world’s large cities. Proc. Natl Acad. Sci. USA 113, 9117–9122 (2016).CAS 
Article 

Google Scholar 
9.Romulo, C. L. et al. Global state and potential scope of investments in watershed services for large cities. Nat. Commun. 9, 4375 (2018).Article 
CAS 

Google Scholar 
10.Tellman, B. et al. Opportunities for natural infrastructure to improve urban water security in Latin America. PLoS ONE 13, e0209470 (2018).Article 

Google Scholar 
11.United Nations World Water Assessment Programme/UN-Water The United Nations World Water Development Report 2018: Nature-Based Solutions for Water (UNESCO, 2018).12.Palmer, M. A., Liu, J., Matthews, J. H., Mumba, M. & D’Odorico, P. Manage water in a green way. Science 349, 584–585 (2015).CAS 
Article 

Google Scholar 
13.Ziv, G., Baran, E., Nam, S., Rodríguez-Iturbe, I. & Levin, S. A. Trading-off fish biodiversity, food security, and hydropower in the Mekong River Basin. Proc. Natl Acad. Sci. USA 109, 5609–5614 (2012).CAS 
Article 

Google Scholar 
14.Harrison, I. J. et al. Protected areas and freshwater provisioning: a global assessment of freshwater provision, threats and management strategies to support human water security. Aquat. Conserv. Mar. Freshw. Ecosyst. 26, 103–120 (2016).Article 

Google Scholar 
15.The World Database on Protected Areas (IUCN and UNEP-WCMC, 2017); http://www.protectedplanet.net16.Huber-Stearns, H. R., Goldstein, J. H., Cheng, A. S. & Toombs, T. P. Institutional analysis of payments for watershed services in the western United States. Ecosyst. Serv. 16, 83–93 (2015).Article 

Google Scholar 
17.Moran, E. F., Lopez, M. C., Moore, N., Müller, N. & Hyndman, D. W. Sustainable hydropower in the 21st century. Proc. Natl Acad. Sci. USA 115, 11891–11898 (2018).CAS 
Article 

Google Scholar 
18.Zheng, H. et al. Benefits, costs, and livelihood implications of a regional payment for ecosystem service program. Proc. Natl Acad. Sci. USA 110, 16681–16686 (2013).CAS 
Article 

Google Scholar 
19.Adamowicz, W. et al. Assessing ecological infrastructure investments. Proc. Natl Acad. Sci. USA 116, 201802883 (2019).Article 
CAS 

Google Scholar 
20.McDonald R. I. Conservation for Cities: How to Plan & Build Natural Infrastructure (Island Press, 2015).21.Grill, G. et al. An index-based framework for assessing patterns and trends in river fragmentation and flow regulation by global dams at multiple scales. Environ. Res. Lett. 10, 015001 (2015).Article 

Google Scholar 
22.Poff, N. L. & Schmidt, J. C. How dams can go with the flow. Science 353, 1099–1100 (2016).CAS 
Article 

Google Scholar 
23.Liu, J. & Yang, W. Integrated assessments of payments for ecosystem services programs. Proc. Natl Acad. Sci. USA 110, 16297–16298 (2013).CAS 
Article 

Google Scholar 
24.Muller, M., Biswas, A., Martin-Hurtado, R. & Tortajada, C. Built infrastructure is essential. Science 349, 585–586 (2015).CAS 
Article 

Google Scholar 
25.Veldkamp, T. I. E. et al. Water scarcity hotspots travel downstream due to human interventions in the 20th and 21st century. Nat. Commun. 8, 15697 (2017).CAS 
Article 

Google Scholar 
26.Cohen, S., Kettner, A. J. & Syvitski, J. P. M. Global suspended sediment and water discharge dynamics between 1960 and 2010: continental trends and intra-basin sensitivity. Glob. Planet. Change 115, 44–58 (2014).Article 

Google Scholar 
27.Dottori, F. et al. Development and evaluation of a framework for global flood hazard mapping. Adv. Water Resour. 94, 87–102 (2016).Article 

Google Scholar 
28.Byers L. et al. A Global Database of Power Plants (World Resources Institute, 2018); https://www.wri.org/publication/global-power-plant-database29.Liu, J. Integration across a metacoupled world. Ecol. Soc. 22, 29 (2017).Article 

Google Scholar 
30.Vercruysse, K., Grabowski, R. C. & Rickson, R. J. Suspended sediment transport dynamics in rivers: multi-scale drivers of temporal variation. Earth Sci. Rev. 166, 38–52 (2017).Article 

Google Scholar 
31.Wu, X.-X., Gu, Z.-J., Luo, H., Shi, X.-Z. & Yu, D.-S. Analyzing forest effects on runoff and sediment production using leaf area index. J. Mt. Sci. 11, 119–130 (2014).Article 

Google Scholar 
32.Wang, Y. et al. Annual runoff and evapotranspiration of forestlands and non-forestlands in selected basins of the Loess Plateau of China. Ecohydrology 4, 277–287 (2011).CAS 
Article 

Google Scholar 
33.Bilotta, G. S. & Brazier, R. E. Understanding the influence of suspended solids on water quality and aquatic biota. Water Res. 42, 2849–2861 (2008).CAS 
Article 

Google Scholar 
34.Stickler, C. M. et al. Dependence of hydropower energy generation on forests in the Amazon Basin at local and regional scales. Proc. Natl Acad. Sci. USA 110, 9601–9606 (2013).CAS 
Article 

Google Scholar 
35.Maltby, E. & Acreman, M. C. Ecosystem services of wetlands: pathfinder for a new paradigm. Hydrol. Sci. J. 56, 1341–1359 (2011).Article 

Google Scholar 
36.Shuster, W. D., Bonta, J., Thurston, H., Warnemuende, E. & Smith, D. R. Impacts of impervious surface on watershed hydrology: a review. Urban Water J. 2, 263–275 (2005).Article 

Google Scholar 
37.Borrelli, P. et al. Land use and climate change impacts on global soil erosion by water (2015–2070). Proc. Natl Acad. Sci. USA 117, 21994–22001 (2020).CAS 
Article 

Google Scholar 
38.Blöschl, G. et al. Changing climate both increases and decreases European river floods. Nature 573, 108–111 (2019).Article 
CAS 

Google Scholar 
39.Symes, W. S., Rao, M., Mascia, M. B. & Carrasco, L. R. Why do we lose protected areas? Factors influencing protected area downgrading, downsizing and degazettement in the tropics and subtropics. Glob. Change Biol. 22, 656–665 (2016).Article 

Google Scholar 
40.Kremen, C. & Merenlender, A. M. Landscapes that work for biodiversity and people. Science 362, eaau6020 (2018).Article 
CAS 

Google Scholar 
41.Dinerstein, E. et al. A global deal for nature: guiding principles, milestones, and targets. Sci. Adv. 5, eaaw2869 (2019).CAS 
Article 

Google Scholar 
42.Liu, J. et al. China’s environment on a metacoupled planet. Annu. Rev. Environ. Resour. 43, 1–34 (2018).CAS 
Article 

Google Scholar 
43.Viña, A., McConnell, W. J., Yang, H., Xu, Z. & Liu, J. Effects of conservation policy on China’s forest recovery. Sci. Adv. 2, e1500965 (2016).Article 

Google Scholar 
44.Chen, C. et al. China and India lead in greening of the world through land-use management. Nat. Sustain. 2, 122–129 (2019).Article 

Google Scholar 
45.Ouyang, Z. et al. Improvements in ecosystem services from investments in natural capital. Science 352, 1455–1459 (2016).CAS 
Article 

Google Scholar 
46.Vörösmarty, C. J. et al. Ecosystem-based water security and the Sustainable Development Goals (SDGs). Ecohydrol. Hydrobiol. 18, 317–333 (2018).Article 

Google Scholar 
47.Liu, J. et al. Nexus approaches to global sustainable development. Nat. Sustain. 1, 466–476 (2018).Article 

Google Scholar 
48.Flörke, M., Schneider, C. & McDonald, R. I. Water competition between cities and agriculture driven by climate change and urban growth. Nat. Sustain. 1, 51–58 (2018).Article 

Google Scholar 
49.McDonald, R. I. et al. Urban growth, climate change, and freshwater availability. Proc. Natl Acad. Sci. USA 108, 6312–6317 (2011).CAS 
Article 

Google Scholar 
50.Willner, S. N., Otto, C. & Levermann, A. Global economic response to river floods. Nat. Clim. Change 8, 594–598 (2018).Article 

Google Scholar 
51.Cattaneo, A., Nelson, A. & McMenomy, T. Global mapping of urban–rural catchment areas reveals unequal access to services. Proc. Natl Acad. Sci. USA 118, e2011990118 (2021).CAS 
Article 

Google Scholar 
52.Zarfl, C., Lumsdon, A. E., Berlekamp, J., Tydecks, L. & Tockner, K. A global boom in hydropower dam construction. Aquat. Sci. 77, 161–170 (2015).Article 

Google Scholar 
53.Schneider, A., Friedl, M. A. & Potere, D. A new map of global urban extent from MODIS satellite data. Environ. Res. Lett. 4, 044003 (2009).Article 

Google Scholar 
54.Lehner, B., Verdin, K. & Jarvis, A. New global hydrography derived from spaceborne elevation data. EOS 89, 93–94 (2008).55.Yang, H. et al. A global assessment of the impact of individual protected areas on preventing forest loss. Sci. Total Environ. 777, 145995 (2021).CAS 
Article 

Google Scholar 
56.Smith, A. et al. New estimates of flood exposure in developing countries using high-resolution population data. Nat. Commun. 10, 1814 (2019).Article 
CAS 

Google Scholar 
57.Best, J. Anthropogenic stresses on the world’s big rivers. Nat. Geosci. 12, 7–21 (2019).CAS 
Article 

Google Scholar 
58.Hanasaki, N. et al. An integrated model for the assessment of global water resources—Part 1: model description and input meteorological forcing. Hydrol. Earth Syst. Sci. 12, 1007–1025 (2008).Article 

Google Scholar 
59.Bondeau, A. et al. Modelling the role of agriculture for the 20th century global terrestrial carbon balance. Glob. Change Biol. 13, 679–706 (2007).Article 

Google Scholar 
60.Pokhrel, Y. N. et al. Incorporation of groundwater pumping in a global Land Surface Model with the representation of human impacts. Water Resour. Res. 51, 78–96 (2015).Article 

Google Scholar 
61.Wada, Y., Wisser, D. & Bierkens, M. F. P. Global modeling of withdrawal, allocation and consumptive use of surface water and groundwater resources. Earth Syst. Dyn. 5, 15–40 (2014).Article 

Google Scholar 
62.Müller Schmied, H. et al. Variations of global and continental water balance components as impacted by climate forcing uncertainty and human water use. Hydrol. Earth Syst. Sci. 20, 2877–2898 (2016).Article 

Google Scholar 
63.Sheffield, J., Goteti, G. & Wood, E. F. Development of a 50-year high-resolution global dataset of meteorological forcings for land surface modeling. J. Clim. 19, 3088–3111 (2006).Article 

Google Scholar 
64.Dirmeyer, P. A. et al. GSWP-2: multimodel analysis and implications for our perception of the land surface. Bull. Am. Meteorol. Soc. 87, 1381–1398 (2006).Article 

Google Scholar 
65.Weedon, G. P. et al. The WFDEI meteorological forcing data set: WATCH forcing data methodology applied to ERA-Interim reanalysis data. Water Resour. Res. 50, 7505–7514 (2014).Article 

Google Scholar 
66.Bingham, H. C. et al. Sixty years of tracking conservation progress using the World Database on Protected Areas. Nat. Ecol. Evol. 3, 737–743 (2019).Article 

Google Scholar 
67.ArcGIS Desktop: Release 10.3.1 (Environmental Systems Research Institution, 2015).68.Domisch, S., Amatulli, G. & Jetz, W. Near-global freshwater-specific environmental variables for biodiversity analyses in 1 km resolution. Sci. Data 2, 150073 (2015).CAS 
Article 

Google Scholar 
69.Bennett, G. & Ruef, F. Alliances for Green Infrastructure: State of Watershed Investment 2016 (Forest Trends, 2016).70.R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2017).71.Wellman, B. & Frank, K. in Social Capital: Theory and Research (eds Lin, N. et al.) 233–273 (Routledge, 2001).72.Bates, D., Machler, M., Bolker, B. M. & Walker, S. C. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar  More