Ecology

1.
Minnis, A. M. & Lindner, D. L. Phylogenetic evaluation of geomyces and allies reveals no close relatives of Pseudogymnoascus destructans, comb. Nov., in bat hibernacula of eastern North America. Fungal Biol. 117, 638–649 (2013).
PubMed  Article  PubMed Central  Google Scholar 
2.
Blehert, D. S. et al. Bat white-nose syndrome: An emerging fungal pathogen?. Science 323, 227. https://doi.org/10.1126/science.1163874 (2009).
CAS  Article  PubMed  PubMed Central  Google Scholar 

3.
Chaturvedi, V. & Chaturvedi, S. Editorial: What is in a Name? A proposal to use geomycosis instead of white nose syndrome (WNS) to describe bat infection caused by geomyces destructans. Mycopathologia 171, 231–233. https://doi.org/10.1007/s11046-010-9385-3 (2011).
Article  PubMed  PubMed Central  Google Scholar 

4.
Frick, W. F., Puechmaille, S. J. & Willis, C. K. Bats in the Anthropocene: Conservation of Bats in a Changing World 245–262 (Springer, Berlin, 2016).
Google Scholar 

5.
Bandouchova, H. et al. Alterations in the health of hibernating bats under pathogen pressure. Sci. Rep. 8, 6067 (2018).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

6.
Zukal, J. et al. White-nose syndrome without borders: Pseudogymnoascus destructans infection tolerated in Europe and Palearctic Asia but not in North America. Sci. Rep. 6, 19829 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

7.
Drees, K. P. et al. Phylogenetics of a fungal invasion: Origins and widespread dispersal of white-nose syndrome. mBio 8, 11941–11917 (2017).
Article  Google Scholar 

8.
Leopardi, S., Blake, D. & Puechmaille, S. J. White-nose syndrome fungus introduced from Europe to North America. Curr. Biol. 25, 217–219 (2015).
Article  CAS  Google Scholar 

9.
Palmer, J. M. et al. Molecular characterization of a heterothallic mating system in Pseudogymnoascus destructans, the fungus causing white-nose syndrome of bats. G3 Genes Genomes Genet. 4, 1755–1763 (2014).
Google Scholar 

10.
Trivedi, J. N. Population genomics and mutational history of the invasive, epidemic clone of Pseudogymnoascus destructans, causal agent of White-nose Syndrome in bats (University of Toronto, Toronto, 2017).
Google Scholar 

11.
Rajkumar, S. S. et al. Clonal genotype of Geomyces destructans among bats with white nose syndrome, New York, USA. Emerg. Infect. Dis 17, 1273–1276 (2011).
PubMed  PubMed Central  Article  Google Scholar 

12.
Lorch, J. M. et al. First detection of bat white-nose syndrome in western North America. mSphere 1, e00148-e1116 (2016).
CAS  PubMed  PubMed Central  Article  Google Scholar 

13.
Forsythe, A., Giglio, V., Asa, J. & Xu, J. Phenotypic divergence along geographic gradients reveals potential for rapid adaptation of the White-nose syndrome pathogen, Pseudogymnoascus destructans, North America. Appl. Environ. Microbiol. 84, e00863-e1818. https://doi.org/10.1128/aem.00863-18 (2018).
CAS  Article  PubMed  PubMed Central  Google Scholar 

14.
Khankhet, J. et al. Clonal expansion of the Pseudogymnoascus destructans genotype in North America is accompanied by significant variation in phenotypic expression. PLoS ONE 9, e104684 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

15.
Cryan, P. M., Meteyer, C. U., Boyles, J. G. & Blehert, D. S. Wing pathology of white-nose syndrome in bats suggests life-threatening disruption of physiology. BMC Biol. 8, 1–8. https://doi.org/10.1186/1741-7007-8-135 (2010).
Article  Google Scholar 

16.
Meteyer, C. U. et al. Histopathologic criteria to confirm white-nose syndrome in bats. J. Vet. Diagn. Invest. 21, 411–414. https://doi.org/10.1177/104063870902100401 (2009).
Article  PubMed  PubMed Central  Google Scholar 

17.
Pikula, J. et al. White-nose syndrome pathology grading in nearctic and palearctic bats. PLoS ONE 12, e0180435 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

18.
Warnecke, L. et al. Inoculation of bats with European Geomyces destructans supports the novel pathogen hypothesis for the origin of white-nose syndrome. Proc. Natl. Acad. Sci. 109, 6999–7003 (2012).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

19.
Flieger, M. et al. Vitamin B2 as a virulence factor in Pseudogymnoascus destructans skin infection. Sci. Rep. 6, 33200 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

20.
Hayman, D. T. S., Pulliam, J. R. C., Marshall, J. C., Cryan, P. M. & Webb, C. T. Environment, host, and fungal traits predict continental-scale white-nose syndrome in bats. Sci. Adv. 2, e1500831. https://doi.org/10.1126/sciadv.1500831 (2016).
ADS  Article  PubMed  PubMed Central  Google Scholar 

21.
Warnecke, L. et al. Pathophysiology of white-nose syndrome in bats: A mechanistic model linking wing damage to mortality. Vol. 9 (2013).

22.
Wibbelt, G. in Emerging and Epizootic Fungal Infections in Animals 289–307 (Springer, Berlin, 2018).
Google Scholar 

23.
Achterman, R. R. & White, T. C. Dermatophyte virulence factors: Identifying and analyzing genes that may contribute to chronic or acute skin infections. Int. J. Microbiol. 20, 12 (2011).
Google Scholar 

24.
Chinnapun, D. Virulence factors involved in pathogenicity of dermatophytes. Walailak J. Sci. Technol. (WJST) 12, 573–580 (2015).
Google Scholar 

25.
Pannkuk, E. L., Risch, T. S. & Savary, B. J. Isolation and identification of an extracellular subtilisin-like serine protease secreted by the bat pathogen Pseudogymnoascus destructans. PLoS ONE 10, e0120508. https://doi.org/10.1371/journal.pone.0120508 (2015).
CAS  Article  PubMed  PubMed Central  Google Scholar 

26.
Raudabaugh, D. B. & Miller, A. N. Nutritional capability of and substrate suitability for Pseudogymnoascus destructans, the causal agent of bat white-nose syndrome. PLoS ONE 8, e78300. https://doi.org/10.1371/journal.pone.0078300 (2013).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

27.
Mascuch, S. J. et al. Direct detection of fungal siderophores on bats with white-nose syndrome via fluorescence microscopy-guided ambient ionization mass spectrometry. PLoS ONE 10, e0119668. https://doi.org/10.1371/journal.pone.0119668 (2015).
CAS  Article  PubMed  PubMed Central  Google Scholar 

28.
van Burik, J. A. H. & Magee, P. T. Aspects of fungal pathogenesis in humans. Annu. Rev. Microbiol. 55, 743–772 (2001).
PubMed  Article  Google Scholar 

29.
Donaldson, M. E. et al. Growth medium and incubation temperature alter the Pseudogymnoascus destructans transcriptome: Implications in identifying virulence factors. Mycologia 110, 300–315 (2018).
CAS  PubMed  Article  Google Scholar 

30.
Field, K. A. et al. The white-nose syndrome transcriptome: activation of anti-fungal host responses in wing tissue of hibernating little brown Myotis. PLoS Pathog. 11, e1005168. https://doi.org/10.1371/journal.ppat.1005168 (2015).
CAS  Article  PubMed  PubMed Central  Google Scholar 

31.
Reeder, S. M. et al. Pseudogymnoascus destructans transcriptome changes during white-nose syndrome infections. Virulence 8, 1695–1707 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

32.
Lorch, J. M. et al. Experimental infection of bats with Geomyces destructans causes white-nose syndrome. Nature 480, 376 (2011).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

33.
Lorch, J. M. et al. A culture-based survey of fungi in soil from bat hibernacula in the eastern United States and its implications for detection of Geomyces destructans, the causal agent of bat white-nose syndrome. Mycologia 105, 237–252. https://doi.org/10.3852/12-207 (2012).
CAS  Article  PubMed  PubMed Central  Google Scholar 

34.
Meyer, A. D., Stevens, D. F. & Blackwood, J. C. Predicting bat colony survival under controls targeting multiple transmission routes of white-nose syndrome. J. Theor. Biol. 409, 60–69 (2016).
MathSciNet  CAS  PubMed  MATH  Article  PubMed Central  Google Scholar 

35.
Gargas, A., Trest, M., Christensen, M., Volk, T. J. & Blehert, D. Geomyces destructans sp. nov. associated with bat white-nose syndrome. Mycotaxon 108, 147–154 (2009).
Article  Google Scholar 

36.
Chaturvedi, V. et al. Morphological and molecular characterizations of psychrophilic fungus Geomyces destructans from New York bats with white nose syndrome (WNS). PLoS ONE 5, e10783 (2010).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

37.
Verant, M. L., Boyles, J. G., Waldrep, W. Jr., Wibbelt, G. & Blehert, D. S. Temperature-dependent growth of Geomyces destructans, the fungus that causes bat white-nose syndrome. PLoS ONE 7, e46280 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

38.
Palmer, J. M., Drees, K. P., Foster, J. T. & Lindner, D. L. Extreme sensitivity to ultraviolet light in the fungal pathogen causing white-nose syndrome of bats. Nat. Commun. 9, 35. https://doi.org/10.1038/s41467-017-02441-z (2018).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

39.
Campbell, L. J., Walsh, D. P., Blehert, D. S. & Lorch, J. M. Long-term survival of Pseudogymnoascus destructans at elevated temperatures. J. Wildlife Dis. 56, 278–287 (2020).
Article  Google Scholar 

40.
Reynolds, H. T. & Barton, H. A. Comparison of the white-nose syndrome agent Pseudogymnoascus destructans to cave-dwelling relatives suggests reduced saprotrophic enzyme activity. PLoS ONE 9, e86437. https://doi.org/10.1371/journal.pone.0086437 (2014).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

41.
Smyth, C., Schlesinger, S., Overton, B. & Butchkoski, C. The alternative host hypothesis and potential virulence genes in Geomyces destructans. Bat Res. News 54, 17–24 (2013).
Google Scholar 

42.
Chaturvedi, V., DeFiglio, H. & Chaturvedi, S. Phenotype profiling of white-nose syndrome pathogen Pseudogymnoascus destructans and closely-related Pseudogymnoascus pannorum reveals metabolic differences underlying fungal lifestyles. F1000Research 7, 2 (2018).
Article  CAS  Google Scholar 

43.
Vanderwolf, K. J., Malloch, D., McAlpine, D. F. & Forbes, G. J. A world review of fungi, yeasts, and slime molds in caves. Int. J. Speleol. 42, 9 (2013).
Article  Google Scholar 

44.
Wilson, M. B., Held, B. W., Freiborg, A. H., Blanchette, R. A. & Salomon, C. E. Resource capture and competitive ability of non-pathogenic Pseudogymnoascus spp. and P. destructans, the cause of white-nose syndrome in bats. PLoS ONE 12, e0178968. https://doi.org/10.1371/journal.pone.0178968 (2017).
CAS  Article  PubMed  PubMed Central  Google Scholar 

45.
Gabriel, K. T., Neville, J. J., Pierce, G. E. & Cornelison, C. T. Lipolytic activity and the utilization of fatty acids associated with bat sebum by Pseudogymnoascus destructans. Mycopathologia 184, 625–636. https://doi.org/10.1007/s11046-019-00381-4 (2019).
CAS  Article  PubMed  PubMed Central  Google Scholar 

46.
Park, M., Do, E. & Jung, W. H. Lipolytic enzymes involved in the virulence of human pathogenic fungi. Mycobiology 41, 67–72. https://doi.org/10.5941/myco.2013.41.2.67 (2013).
CAS  Article  PubMed  PubMed Central  Google Scholar 

47.
Carlini, C. R. & Ligabue-Braun, R. Ureases as multifunctional toxic proteins: A review. Toxicon 110, 90–109 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

48.
Cox, G. M., Mukherjee, J., Cole, G. T., Casadevall, A. & Perfect, J. R. Urease as a virulence factor in experimental cryptococcosis. Infect. Immun. 68, 443–448 (2000).
CAS  PubMed  PubMed Central  Article  Google Scholar 

49.
Vylkova, S. & Lorenz, M. C. Modulation of phagosomal pH by Candida albicans promotes hyphal morphogenesis and requires Stp2p, a regulator of amino acid transport. PLoS Pathog. 10, e1003995. https://doi.org/10.1371/journal.ppat.1003995 (2014).
CAS  Article  PubMed  PubMed Central  Google Scholar 

50.
Vylkova, S. Environmental pH modulation by pathogenic fungi as a strategy to conquer the host. PLoS Pathog. 13, e1006149. https://doi.org/10.1371/journal.ppat.1006149 (2017).
CAS  Article  PubMed  PubMed Central  Google Scholar 

51.
Shawcross, D. L. et al. Ammonia impairs neutrophil phagocytic function in liver disease. Hepatology 48, 1202–1212 (2008).
CAS  PubMed  Article  Google Scholar 

52.
O’Donoghue, A. J. et al. Destructin-1 is a collagen-degrading endopeptidase secreted by Pseudogymnoascus destructans, the causative agent of white-nose syndrome. PNAS 112, 7478–7483. https://doi.org/10.1073/pnas.1507082112 (2015).
ADS  CAS  Article  PubMed  Google Scholar 

53.
Marroquin, C. M., Lavine, J. O. & Windstam, S. T. Effect of humidity on development of Pseudogymnoascus destructans, the causal agent of bat white-nose syndrome. Northeastern Nat. 24, 54–64 (2017).
Article  Google Scholar 

54.
Kolařík, M. et al. Geosmithia associated with bark beetles and woodborers in the western USA: Taxonomic diversity and vector specificity. Mycologia 109, 185–199. https://doi.org/10.1080/00275514.2017.1303861 (2017).
CAS  Article  PubMed  Google Scholar 

55.
Garland, J. L. Analytical approaches to the characterization of samples of microbial communities using patterns of potential C source utilization. Soil Biol. Biochem. 28, 213–221. https://doi.org/10.1016/0038-0717(95)00112-3 (1996).
CAS  Article  Google Scholar 

56.
Dobranic, J. K. & Zak, J. C. A microtiter plate procedure for evaluating fungal functional diversity. Mycologia 91, 756–765 (1999).
Article  Google Scholar 

57.
Harch, B. D., Correll, R. L., Meech, W., Kirkby, C. A. & Pankhurst, C. E. Using the Gini coefficient with BIOLOG substrate utilisation data to provide an alternative quantitative measure for comparing bacterial soil communities. J. Microbiol. Methods 30, 91–101. https://doi.org/10.1016/s0167-7012(97)00048-1 (1997).
CAS  Article  Google Scholar 

58.
Sobek, E. A. & Zak, J. C. The Soil FungiLog procedure: Method and analytical approaches toward understanding fungal functional diversity. Mycologia 95, 590–602 (2003).
CAS  PubMed  Article  PubMed Central  Google Scholar 

59.
Kouker, G. & Jaeger, K.-E. Specific and sensitive plate assay for bacterial lipases. Appl. Environ. Microbiol. 53, 211–213 (1987).
CAS  PubMed  PubMed Central  Article  Google Scholar 

60.
Lupan, D. M. & Nziramasanga, P. Collagenolytic activity of Coccidioides immitis. Infect. Immun. 51, 360–361 (1986).
CAS  PubMed  PubMed Central  Article  Google Scholar 

61.
Saleh-Rastin, N., Petersen, M., Coleman, S. & Hubbell, D. The rhizosphere and plant growth 188–188 (Springer, Berlin, 1991).
Google Scholar 

62.
NziramasangaM, P. & Lupan, D. Elastase activity of Coccidioides immitis. J. Med. Microbiol. 19, 109–114 (1985).
Article  Google Scholar 

63.
Dietz, M. & Kalko, E. K. Seasonal changes in daily torpor patterns of free-ranging female and male Daubenton’s bats (Myotis daubentonii). J. Comp. Physiol. B. 176, 223–231 (2006).
PubMed  Article  Google Scholar 

64.
Sephton-Clark, P. C. S. & Voelz, K. In Advances in applied microbiology (eds Sima, S. & Geoffrey, M. G.) 117–157 (Academic Press, New York, 2018).
Google Scholar 

65.
Hammer, O., Harper, D. A. T. & Ryan, P. D. PAST: Paleontological statistics software package for education and data analysis. Palaeontol. Electron. 4, 1–9 (2001).
Google Scholar 

66.
Martínková, N. et al. Increasing incidence of Geomyces destructans fungus in bats from the Czech Republic and Slovakia. PLoS ONE 5, e13853 (2010).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

67.
Větrovský, T., Kolařík, M., Žifčáková, L., Zelenka, T. & Baldrian, P. The rpb2 gene represents a viable alternative molecular marker for the analysis of environmental fungal communities. Mol. Ecol. Resour. 16, 388–401 (2016).
PubMed  Article  CAS  Google Scholar 

68.
Crous, P. et al. Fungal planet description sheets: 558–624. Persoonia 38, 240 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

69.
Hubka, V. et al. A reappraisal of Aspergillus section Nidulantes with descriptions of two new sterigmatocystin-producing species. Plant Syst. Evol. 302, 1267–1299 (2016).
CAS  Article  Google Scholar 

70.
Kubátová, A., Hujslová, M., Frisvad, J. C., Chudíčková, M. & Kolařík, M. Taxonomic revision of the biotechnologically important species Penicillium oxalicum with the description of two new species from acidic and saline soils. Mycol. Progr. 18, 215–228 (2019).
Article  Google Scholar 

71.
Gabrielová, A. et al. The oomycete Pythium oligandrum can suppress and kill the causative agents of dermatophytoses. Mycopathologia 183, 751–764 (2018).
PubMed  PubMed Central  Article  CAS  Google Scholar  More

1.
FAO. In World Food Summit — Food for All Ch. 2 (FAO, 1996).
2.
Challinor, A. J. et al. Nat. Clim. Change 4, 287–291 (2014).
Article  Google Scholar 

3.
Myers, S. S. et al. Annu. Rev. Publ. Health 38, 259–277 (2017).
Article  Google Scholar 

4.
Woolfe, J. A. Sweet Potato: An Untapped Food Resource (Cambridge Univ. Press, 1992).

5.
Heider, B. et al. Nat. Clim. Change https://doi.org/10.1038/s41558-020-00924-4 (2020).

6.
Meyer, R. S. & Purugganan, M. D. Nat. Rev. Genet. 14, 840–852 (2013).
CAS  Article  Google Scholar 

7.
Pironon, S. et al. Nat. Clim. Change 9, 758–763 (2019).
Article  Google Scholar 

8.
Borrell, J. S. et al. Environ. Exp. Bot. 170, 103872 (2020).
CAS  Article  Google Scholar 

9.
Wu, S. et al. Nat. Commun. 9, 4580 (2018).
Article  Google Scholar 

10.
Khoury, C. K. et al. Front. Plant Sci. 6, 251 (2015).
Article  Google Scholar 

11.
Shiotani, I., Huang, Z. Z., Sakamoto, S. & Miyazaki, T. Acta Hortic. 380, 388–398 (1994).
Article  Google Scholar 

12.
Dempewolf, H. et al. Crop Sci. 57, 1070–1082 (2017).
Article  Google Scholar 

13.
Castañeda-Álvarez, N. P. et al. Nat. Plants 2, 16022 (2016).
Article  Google Scholar 

14.
Cámara-Leret, R., Fortuna, M. A. & Bascompte, J. Proc. Natl Acad. Sci. USA 116, 9913–9918 (2019).
Article  Google Scholar 

15.
McCouch, S. et al. Nature 499, 23–24 (2013).
CAS  Article  Google Scholar  More

Epidemiological data
No officially reported line list was available for cases in China. We used a standardized protocol43 to extract individual-level data from December 1, 2019, to March 30, 2020. Sources were mainly official reports from provincial, municipal or national health governments. Data included basic demographics (age and sex), travel histories and key dates (dates of onset of symptoms, hospitalization and confirmation). Data were entered by a team of data curators on a rolling basis, and technical validation and geo-positioning protocols were applied continuously to ensure validity. A detailed description of the methodology is available22. Lastly, total numbers were matched with officially reported data from China and other government reports. Daily case counts from Italian provinces (n = 107) were extracted from the Presidenza del Consiglio dei Ministri Dipartimento della Protezione Civile (https://github.com/pcm-dpc/COVID-19).
Estimating epidemic peakedness
Epidemic peakedness was estimated for each prefecture by calculating the inverse Shannon entropy of the distribution of COVID-19 cases. Inverse Shannon entropy was used to fit time series of other respiratory infections (influenza)13. The inverse Shannon entropy of incidence for a given prefecture in 2020 is then given by (v_j = left( { – mathop {sum}nolimits_i {p_{ij}log p_{ij}} } right)^{ – 1}). Because vj is a function of incidence distribution in each location rather than raw incidence, it is invariant under differences in overall reporting rates between cities or attack rates. We then assessed how peakedness ({it{v}} propto mathop {sum}nolimits_j {v_j}) varied across geographic areas in China. As an alternative measure of temporal clustering of cases, we estimated the proportion of cases at the peak ± 1 d (Extended Data Fig. 2).
Proxies for COVID-19 interventions using within-city human mobility data from China
Estimates of within-city reductions of human mobility between the period before and after the lockdown was implemented on January 23, 2020, were extracted from Lai et al.36. Daily measures of human mobility were extracted from the Baidu Qianxi web platform to estimate the proportion of daily movement within prefectures in China. Relative mobility volume was available from January 2, 2020, to January 25, 2020. For each city, change in relative mobility was defined by mi=mil(lockdown)/mib(baseline), where mi is defined as mobility in prefecture i. Baidu’s mapping service is estimated to have a 30% market share in China, and more data can be found5,6.
Data on drivers of transmission of COVID-19
Prefecture-specific population counts and densities were derived from the 2020 Gridded Population of The World, a modeled continuous surface of population estimated from national census data and the United Nations World Population Prospectus44. Population counts are defined at a 30-arc-second resolution (approximately 1 km × 1 km at the equator) and extracted within administrative 2 level cartographic boundaries defined by the National Bureau of Statistics of China. Lloyd’s mean crowding, (frac{{left[ {mathop {sum }nolimits_i left( {q_i – 1} right)q_i} right]}}{{mathop {sum }nolimits_i q_i}}), was estimated for each prefecture, where qi represents the population count of each non-zero pixel within a prefecture’s boundary and the resulting value estimates an individual’s mean number of expected neighbors13,45. When fitting the models, we consider the numerator (left[ {mathop {sum}nolimits_i {left( {q_i – 1} right)q_i} } right]), which we refer to as ‘contacts’, and the denominator (mathop {sum}nolimits_i {q_i}) (that is, population size) as separate predictors. We note that a negative slope for ‘contacts’ and a positive slope for ‘population’ support a negative coefficient for Lloyd’s mean crowding.
Daily temperature (°F), relative humidity (%) and atmospheric pressure (Pa) at the centroid of each prefecture was provided by The Dark Sky Company via the Dark Sky API and aggregated across a variety of data sources. Specific humidity (kg/kg) was then calculated using the R package humidity16. Meteorological variables for each prefecture were then averaged across the entirety of the study period.
Statistical analysis
We normalized the values of epidemic peakedness between 0 and 1 and, for all non-zero values, fit a generalized linear model of the form

$${mathrm{log}}left( {Y_j} right)sim beta _0 + beta _1{mathrm{log}}left( {C_j} right) + beta _2{mathrm{log}}left({h_j}right) + beta _3{mathrm{log}}left( {P_j} right) + beta _4{mathrm{log}}left( {f_j} right) + beta _5{mathrm{log}}left( {t_j} right)$$

where, for each prefecture j, Y is the scaled inverse Shannon entropy measure of epidemic peakedness derived from the COVID-19 time series; C is the mean number of contacts26,46; h is the mean specific humidity over the reporting period in kg/kg; P is the estimated population density; f is the relative change in population flows within each prefecture; and t is daily mean temperature.
Projecting epidemic peakedness in cities around the world
We selected 310 urban centers from the European Commission Global Human Settlement Urban Centre Database and their included cartographic boundaries47. To ensure global coverage, up to the five most populous cities in each country were selected from the 1,000 most populous urban centers recorded in the database. Population count, crowding and meteorological variables were then estimated following identical procedures used to calculate these variables in the Chinese prefectures. Weather measurements were averaged over the 2-month period starting on February 1, 2020.
The parameters from the model of epidemic peakedness predicted by humidity, crowding and population size (Supplementary Table 1, model 6) were used to estimate relative peakedness in the 310 urban centers. A full list of predicted epidemic peakedness values can be found in Supplementary Table 3.
Global human mobility data
We used the Google COVID-19 Aggregated Mobility Research Dataset, which contains anonymized relative mobility flows aggregated over users who have turned on the Location History setting, which is off by default. This is similar to the data used to show how busy certain types of places are in Google Maps, helping identify when a local business tends to be the most crowded. The mobility flux is aggregated per week, between pairs of approximately 5-km2 cells worldwide, and for the purpose of this study aggregated for 310 cities worldwide. We calculated both mobility within each city’s shapefile and mobility coming into each city. For each city, change in relative mobility was defined by mi = mil(April)/mib(December), where mi is defined as mobility in city i.
To produce this data set, machine learning was applied to log data to automatically segment it into semantic trips48. To provide strong privacy guarantees, all trips were anonymized and aggregated using a differentially private mechanism49 to aggregate flows over time (https://policies.google.com/technologies/anonymization). This research is done on the resulting heavily aggregated and differentially private data. No individual user data were ever manually inspected; only heavily aggregated flows of large populations were handled.
All anonymized trips were processed in aggregate to extract their origin and destination location and time. For example, if users traveled from location a to location b within time interval t, the corresponding cell (a, b, t) in the tensor would be n ± err, where err is Laplacian noise. The automated Laplace mechanism adds random noise drawn from a zero-mean Laplace distribution and yields (𝜖, δ)-differential privacy guarantee of 𝜖 = 0.66 and δ = 2.1 × 10−29. The parameter 𝜖 controls the noise intensity in terms of its variance, whereas δ represents the deviation from pure 𝜖-privacy. The closer they are to zero, the stronger the privacy guarantees. Each user contributes, at most, one increment to each partition. If they go from a region a to another region b multiple times in the same week, they contribute only once to the aggregation count.
These results should be interpreted in light of several important limitations. First, the Google mobility data are limited to smartphone users who have opted in to Google’s Location History feature, which is off by default. These data might not be representative of the population as whole, and, furthermore, their representativeness might vary by location. Importantly, these limited data are viewed only through the lens of differential privacy algorithms, specifically designed to protect user anonymity and obscure fine detail. Moreover, comparisons across, rather than within, locations are descriptive only because these regions can differ in substantial ways.
Simulating epidemic dynamics
We simulated a simple stochastic SIR model of infection spread on weighted networks created to represent hierarchically structured populations. Individuals were first assigned to households using the distribution of household sizes in China (data from the United Nations Population Division; mean, 3.4 individuals). Households were then assigned to ‘neighborhoods’ of ~100 individuals, and all neighborhood members were connected with a lower weight. A randomly chosen 10% of individuals were given ‘external’ connections to individuals outside the neighborhood. The total population size was n = 1,000. Simulations were run for 300 d, and averages were taken over 20 iterations. The SIR model used a per-contact transmission rate of 𝛽 = 0.15 per day and recovery rate 𝛾 = 0.1 per day. For the simulations without interventions, the weights were wHH = 1, wNH = 0.01 and wEX = 0.001 for the crowded prefecture and wEX = 0.0001 for the less crowded prefecture. For the simulations with interventions, the household and neighborhood weights were the same, but we used wEX = 0.01 for the crowded prefecture and wEX = 0.001 for the ‘sparse’ prefecture. The intervention reduced the weight of all connections outside the household by 75%.
Reporting Summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.
Pandolfi, J. M. et al. Global trajectories of the long-term decline of coral reef ecosystems. Science 301, 955–958 (2003).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 
2.
Death, G., Fabricius, K. E., Sweatman, H. & Puotinen, M. The 27-year decline of coral cover on the Great Barrier Reef and its causes. Proc. Natl. Acad. Sci. USA 109, 1–5 (2012).
Article  Google Scholar 

3.
Burke, L., Reytar, K., Spalding, M. & Perry, A. Reefs at Risk Revisited (World Resources Institute, Washington, D.C., 2011).
Google Scholar 

4.
Mumby, P. J. et al. Fishing, trophic cascades, and the process of grazing on coral reefs. Science 311, 98–101 (2006).
ADS  CAS  PubMed  Article  Google Scholar 

5.
McLeod, E., Salm, R., Green, A. & Almany, J. Designing marine protected area networks to address the impacts of climate change. Front. Ecol. Environ. 7, 362–370 (2009).
Article  Google Scholar 

6.
McClanahan, T. R., Graham, N. A. J., Calnan, J. M. & MacNeil, M. A. Toward pristine biomass: reef fish recovery in coral reef marine protected areas in Kenya. Ecol. Appl. 17, 1055–1067 (2007).
PubMed  Article  Google Scholar 

7.
Stockwell, B., Jadloc, C. R. L., Abesamis, R. A., Alcala, A. C. & Russ, G. R. Trophic and benthic responses to no-take marine reserve protection in the Philippines. Mar. Ecol. Prog. Ser. 389, 1–15 (2009).
ADS  Article  Google Scholar 

8.
Rotjan, R. D. & Lewis, S. M. Impact of coral predators on tropical reefs. Mar. Ecol. Prog. Ser. 367, 73–91 (2008).
ADS  Article  Google Scholar 

9.
Shaver, E. C., Burkepile, D. E. & Silliman, B. R. Local management actions can increase coral resilience to thermally-induced bleaching. Nat Ecol Evol 2, 1075–1079 (2018).
PubMed  Article  Google Scholar 

10.
Hoeksema, B. W., Scott, C. & True, J. D. Dietary shift in corallivorous Drupella snails following a major bleaching event at Koh Tao Gulf of Thailand. Coral Reefs 32, 423–428 (2013).
ADS  Article  Google Scholar 

11.
Moerland, M. S., Scott, C. M. & Hoeksema, B. W. Prey selection of corallivorous muricids at Koh Tao (Gulf of Thailand) four years after a major coral bleaching event. Contrib. Zool. 85, 291–309 (2019).
Article  Google Scholar 

12.
Kayal, M. et al. Predator crown-of-thorns starfish (Acanthaster planci) outbreak, mass mortality of corals, and cascading effects on reef fish and benthic communities. PLoS ONE 7, e47363 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

13.
Boucher, L. M. Coral predation by muricid gastropods of the genus Drupella at Enewetak Marshall Islands. Bull. Mar. Sci. 38, 9–11 (1986).
Google Scholar 

14.
Moyer, J. T., Emerson, W. K. & Ross, M. Massive destruction of scleractinian corals by the muricid gastropod, Drupella Japan and the Philippines. Nautilus 96, 69–82 (1982).
Google Scholar 

15.
McClanahan, T. R. Coral-eating snail Drupella cornus population increases in Kenyan coral reef lagoons. Mar. Ecol. Prog. Ser. 115, 131–138 (1994).
ADS  Article  Google Scholar 

16.
Knowlton, N., Lang, J. C. & Keller, B. D. Case study of natural population collapse: post-hurricane predation on Jamaican staghorn corals. Smithson. Contrib. Mar. Sci. 1–3, 1–25 (1990).
Article  Google Scholar 

17.
Grober-Dunsmore, R., Bonito, V. & Frazer, T. K. Potential inhibitors to recovery of Acropora palmata populations in St. John, US Virgin Islands. Mar. Ecol. Prog. Ser. 321, 123–132 (2006).
ADS  Article  Google Scholar 

18.
Brodie, J., Fabricius, K., Death, G. & Okaji, K. Are increased nutrient inputs responsible for more outbreaks of crown-of-thorns starfish? An appraisal of the evidence. Mar. Pollut. Bull. 51, 266–278 (2005).
CAS  PubMed  Article  Google Scholar 

19.
McClanahan, T. R. Kenyan coral reef-associated gastropod fauna: a comparison between protected and unprotected reefs. Mar. Ecol. Prog. Ser. 53, 11–20 (1989).
ADS  Article  Google Scholar 

20.
Dulvy, N. K., Freckleton, R. P. & Polunin, N. V. C. Coral reef cascades and the indirect effects of predator removal by exploitation. Ecol. Lett. 7, 410–416 (2004).
Article  Google Scholar 

21.
Sweatman, H. No-take reserves protect coral reefs from predatory starfish. Curr. Biol. 18, 598–599 (2008).
Article  CAS  Google Scholar 

22.
Gardner, T. A., Côté, I. M., Gill, J. A., Grant, A. & Watkinson, A. R. Long-term region-wide declines in Caribbean corals. Science 301, 958–960 (2003).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

23.
Aronson, R. B. & Precht, W. F. White-band disease and the changing face of Caribbean coral reefs. Hydrobiologia 460, 25–38 (2001).
Article  Google Scholar 

24.
Jackson, J. B. C., Donovan, M. K., Cramer, K. L. & Lam, W. Status and Trends of Caribbean Coral Reefs: 1970-2012. Global Coral Reef Monitoring Network, IUCN, Gland. Switz. (2014).

25.
Williams, D. E., Miller, M. W., Bright, A. J. & Cameron, C. M. Removal of corallivorous snails as a proactive tool for the conservation of acroporid corals. PeerJ 2, e680 (2014).
PubMed  PubMed Central  Article  Google Scholar 

26.
Sutherland, K. P., Shaban, S., Joyner, J. L., Porter, J. W. & Lipp, E. K. Human pathogen shown to cause disease in the threatened eklhorn coral Acropora palmata. PLoS ONE 6, e23468 (2011).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

27.
Gignoux-Wolfsohn, S. A., Marks, C. J. & Vollmer, S. V. White Band Disease transmission in the threatened coral Acropora cervicornis. Sci. Rep 2, 804 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

28.
Shaver, E. C. et al. Effects of predation and nutrient enrichment on the success and microbiome of a foundational coral. Ecology 98, 830–938 (2017).
PubMed  Article  Google Scholar 

29.
Williams, D. E. & Miller, M. W. Coral disease outbreak: pattern, prevalence and transmission in Acropora cervicornis. Mar. Ecol. Prog. Ser. 301, 119–128 (2005).
ADS  Article  Google Scholar 

30.
Baums, I. B., Miller, M. W. & Szmant, A. M. Ecology of a corallivorous gastropod, Coralliophila abbreviata, on two scleractinian hosts, I: population structure of snails and corals. Mar. Biol. 142, 1083–1091 (2003).
Article  Google Scholar 

31.
Baums, I. B., Miller, M. W. & Szmant, A. M. Ecology of a corallivorous gastropod, Coralliophila abbreviata, on two seleractinian hosts. II. Feeding, respiration and growth. Mar. Biol. 142, 1093–1101 (2003).
Article  Google Scholar 

32.
National Marine Fisheries Service. Recovery Plan: Elkhorn coral (Acropora palmata) and staghorn coral (A. cervicornis). https://www.fisheries.noaa.gov/resource/document/recovery-plan-elkhorn-coral-acropora-palmata-and-staghorn-coral-cervicornis (2015).

33.
Clements, C. S. & Hay, M. E. Overlooked coral predators suppress foundation species as reefs degrade. Ecol. Appl. 28, 1673–1682 (2018).
PubMed  PubMed Central  Article  Google Scholar 

34.
Sharp, W. C. & Delgado, G. A. Predator-prey interactions between the corallivorous snail Coralliophila abbreviata and the carnivorous deltoid rock snail Thais deltoidea. Biol. Bull. 229, 129–133 (2015).
CAS  PubMed  Article  Google Scholar 

35.
Humann, P. & Deloach, N. Reef coral identification: Florida Caribbean Bahamas (New World Publications Inc, Jacksonville, 2014).
Google Scholar 

36.
National Marine Sanctuaries. Florida Keys National Marine Sanctuary Regulations. National Oceanic and Atmospheric Administration. https://floridakeys.noaa.gov/regs/ (2015).

37.
Bartholomew, A. et al. Influence of marine reserve size and boundary length on the initial response of exploited reef fishes in the Florida Keys National Marine Sanctuary USA. Landsc. Ecol. 23, 55–65 (2008).
Article  Google Scholar 

38.
Kramer, K. L. & Heck, K. L. Top-down trophic shifts in Florida Keys patch reef marine protected areas. Mar. Ecol. Prog. Ser. 349, 111–123 (2007).
ADS  Article  Google Scholar 

39.
Miller, M. W. Corallivorous snail removal: evaluation of impact on Acropora palmata. Coral Reefs 19, 293–295 (2001).
Article  Google Scholar 

40.
Johnston, L. & Miller, M. W. Variation in life-history traits of the corallivorous gastropod Coralliophila abbreviata on three coral hosts. Mar. Biol. 150, 1215–1225 (2006).
Article  Google Scholar 

41.
Lutes, D. L. & Szmant, A. M. Evidence for an indirect link between the Florida Keys spiny lobster fishery and the decline of Acropora palmata. Proceeding in Benthic Ecology Meeting (2005).

42.
Randall, J. E. Food habits of reef fishes of the West Indies. Stud. Trop. Oceanogr. 5, 665–847 (1967).
Google Scholar 

43.
Goldberg, W. M. A note on the feeding behavior of the snapping shrimp Synalpheus Fritzmuelleri Coutière (Decapoda, Alpheidae). Crustaceana 21, 318–320 (1971).
Article  Google Scholar 

44.
Ayotte, P., Mccoy, K., Williams, I. & Zamzow, J. Coral Reef Ecosystem Division Standard Operating Procedures: Data Collection for Rapid Ecological Assessment Fish Surveys. https://repository.library.noaa.gov/view/noaa/702 (2011).

45.
Venables, W. N. & Ripley, B. D. Modern Applied Statistics with S 4th edn. (Springer, Berlin, 2002).
Google Scholar 

46.
R Core Team. foreign: Read Data Stored by ‘Minitab’, ‘S’, ‘SAS’, ‘SPSS’, ‘Stata’, ‘Systat’, ‘Weka’, ‘dBase’, … R package version 0.8–78. https://CRAN.R-project.org/package=foreign (2020).

47.
R Core Team. R: A Language and Environment for Statistical Computing. (2020).

48.
Wickham, H. et al. Welcome to the Tidyverse. J. Open Source Softw. 4, 1686 (2019).
ADS  Article  Google Scholar 

49.
Office of National Marine Sanctuaries. Florida Keys National Marine Sanctuary Condition Report 2011. https://sanctuaries.noaa.gov/science/condition/fknms/ (2011).

50.
Roberts, C. M. & Hawkins, J. E. How small can a marine reserve be and still be effective?. Coral Reefs 16, 150 (1997).
ADS  Article  Google Scholar 

51.
Cox, C. Monitoring Caribbean Spiny Lobsters in the Florida Keys National Marine Sanctuary, 1997–2002. Florida Keys National Marine Sanctuary. (2006).

52.
Cox, C. & Hunt, J. H. Change in size and abundance of Caribbean spiny lobsters Panulirus argus in a marine reserve in the Florida Keys National Marine Sanctuary, USA. Mar. Ecol. Prog. Ser. 294, 227–239 (2005).
ADS  Article  Google Scholar 

53.
Mumby, P. J. et al. Trophic cascade facilitates coral recruitment in a marine reserve. Proc. Natl. Acad. Sci. USA 104, 8362–8367 (2007).
ADS  CAS  PubMed  Article  Google Scholar 

54.
Ott, B. & Lewis, J. B. The importance of the gastropod Coralliophila abbreviata (Lamarck) and the polychaete Hermodice carunculata (Pallas) as coral reef predators. Can. J. Zool. 50, 1651–1656 (1972).
Article  Google Scholar 

55.
Cumming, R. L. Predation on reef-building corals: multiscale variation in the density of three corallivorous gastropods Drupella spp.. Coral Reefs 18, 147–157 (1999).
Article  Google Scholar 

56.
Ling, S. D., Johnson, C. R., Frusher, S. D. & Ridgway, K. R. Overfishing reduces resilience of kelp beds to climate-driven catastrophic phase shift. Proc. Natl. Acad. Sci. U.S.A. 106, 22341–22345 (2009).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar  More

Study area and animals
This study was conducted on the Shiretoko Peninsula (43°50′–44°20′ N, 144°45′–145°20′ E; Fig. 1), eastern Hokkaido, Japan. We conducted a monitoring survey of adult female bears inhabiting a region extending from the Horobetsu-Iwaobetsu area (44°5′ N, 145°1′ E; Fig. 1) to the Rusha area (44°11′ N, 145°11′ E), and divided the bears into three groups according to their home ranges and levels of human habituation. The home ranges were determined by the core areas revealed by radiotracking data for captured bears, and by the locations of DNA sampling (collection of hair, feces, and skin samples) for the remaining bears. The level of human habituation was determined based on the bears’ behavioral responses to humans, as follows: very high, no aversive reaction when approached by humans within 20 m; high, no aversive reaction to humans  > 20 m away who move to within 20 m; moderate, movement away when approached within 50 m; and low, rare appearance before humans and rapid movement away upon awareness of humans. For bears inhabiting the area patrolled daily by park managers (i.e., the area excluding the Rusha area, as described below), this classification was done by two park managers with long experience with those bears from 1997 to 2018, based on patrol records that described the response of the bears when they approached to chase the bears away. We reviewed multiple randomly selected records for each bear and determined the habituation level based on the response most frequently observed. Alternatively, for bears inhabiting the Rusha area, outside the area patrolled for management purposes, two persons (one in common with the former case) determined the level, based on the response when they approached the bears multiple times during the long-term survey conducted between 2006 and 2018 (detailed in our previous reports31). All bears whose behavioral responses were not available, e.g., bears identified by genetic sampling that did not appear in front of humans, were categorized as low habituation level.
The first group comprised bears living around the front portion of the national park, including the Horobetsu-Iwaobetsu area, Shiretoko Five Lakes (44°7′ N, 145°5′ E), and the Utoro Plateau (44°3′ N, 144°59′ E), which is adjacent to the town of Utoro (44°4′ N, 144°59′ E; Fig. 1), a gateway community. Group 1 comprised 23 females, 14 of which were captured and released with ear tags for radiotracking in the Horobetsu-Iwaobetsu area between 1998 and 2013 and 9 of which were determined by genetic sampling to have been alive and present in this area from 1999 to 2018 [2 with distinctive characteristics (e.g., chest marking) were periodically identified in this area between 2013 and 2018]. This area is the main place where park visitors experience activities in nature. In the last two decades, bear sightings have been increasingly common in this area. When a bear appears near a road or development, park managers rush to the site and chase the bear away by shooting it with non-lethal rubber slugs and cracker shells. Thus, some bears are moderately habituated, but are not allowed to act freely during daylight hours. The second group comprised bears living between the Idashubetsu River (44°9′ N, 145°6′ E; Fig. 1) and the Rusha area, about 12.4–21.8 km from Utoro. Group 2 comprised 26 females, 4 of which were captured and released with ear tags for radiotracking in this area between 1993 and 2015 and 22 of which were determined by genetic sampling to have been alive and present in this area from 2000 to 2015 (13 were periodically identified in the area between 2013 and 2018). Kamuiwakka Falls (44°9′ N, 145°8′ E; Fig. 1) is the final tourist destination in this region; visitors are not allowed to pass through from Kamuiwakka Falls to the Rusha area. The third group comprised human-habituated bears living around the Rusha area, a special wildlife protection area with no public access without permission, and no human residence except for one fishermen’s settlement. Because the fishermen have not excluded bears from the settlement area in the last few decades, the bears have become habituated to the existence of humans and ignore people, which enables direct observation at close range (Fig. S1). Group 3 comprised 13 very highly habituated bears that had been visually identified and frequently observed in this area from the late 1990s to 201828,31,32,33. Among them, six bears were captured and released with ear tags for radiotracking between 2013 and 201831. Non- and less-habituated bears (low–high habituation levels) observed and sampled in the Rusha area were allocated to Group 2. All collared bears were captured live in accordance with the Guidelines for Animal Care and Use of Hokkaido University and all procedures were approved by the Animal Care and Use Committee of the Graduate School of Veterinary Medicine, Hokkaido University (Permit Number: 1152 and 15009, 17005, 18-0083). The protocols for capture received annual approval from the Ministry of the Environment, Japan, and the Hokkaido Government through research permit applications.
Sample collection
During 1998–2018, we collected genetic samples from bears throughout the peninsula, including the towns of Shari, Rausu, and Shibetsu (Fig. 1), using multiple methods. Most samples were from bears that were dead due to nuisance control or hunting, or that were captured for research purposes. Blood, hair, liver, and/or muscle samples were obtained. Hunting is limited from October to January and is not allowed inside the national park. When a bear is killed, all hunters are required to report the reason why it was killed (i.e., nuisance control or hunting), date, location, sex, estimated age (based on body size), body mass, and size (if available), and to provide the genetic samples described above. Due to the strong relationship between park managers and hunters in the peninsula, poaching or hunting without a report are very unlikely to have occurred over the past two decades. In the current analysis, nuisance killing and hunting were included as human-caused mortality because, in most situations, hunting occurs near human settlements, especially around agricultural land, to mitigate the human–bear conflict in autumn, which makes it reasonable to consider that the likelihood of mortality by hunting is also influenced by the bears’ behavioral traits, including habituation level. We also obtained: (1) hair samples collected from hair traps, including fence traps and tree-rub traps placed in several locations (including the Horobetsu-Iwaobetsu and Rusha areas) during 2010–2018; (2) skin tissues collected by biopsy dart sampling during 2011–2018; and (3) fecal samples collected during 2009–2018. The sampling methods have been detailed in our previous reports30,31,32,34. For bears captured or killed between 1998 and 2018, age was estimated by counting the dental cementum annuli35.
Human-caused mortality of the offspring of adult females
To determine how often offspring of adult females were killed due to human activities after separation from their mothers at 1.5–2.5 years of age31, the number of offspring (aged ≥ 1 year) produced by those females during the follow-up period must be known. However, we were unable to determine the exact number of offspring because the females, except for habituated bears in the Rusha area, were not monitored annually by direct observation. Therefore, the number of offspring produced by each female was estimated using the following formula: bear-years (follow-up years at ages ≥ 5 years) × reproductive rate (young born/year/reproductive adult female aged ≥ 5 years) × cub survival rate (from approximately 0.5 to 1.5 years of age). We assumed an equal sex ratio at birth. The latter two parameters were investigated in a long-term, individual-based monitoring survey in the Rusha area; values obtained were 0.74 and 0.63, respectively31. The bear-years were calculated from the year of first identification at ≥ 5 years of age (the earliest age of first reproduction in the Rusha area31) to the year of last identification. The initiation year was determined based on: (1) direct observation of a bear whose birth year was known, (2) dental age estimation at the time of capture or death, or (3) reproductive history (i.e., the age at first identification with cubs was assumed to be ≥ 5 years). The latter parameter was determined by not only direct observation of a female with cubs, but also indirect evidence of birth experience based on parentage analysis (described below). For example, if a female was assigned to the mother of 3-year-old bears (estimated by dental examination) killed in 2012, the mother was assumed to be ≥ 5 years old in 2009. The completion year was determined based on: (1) the year of last observation or death or (2) the reproductive history, revealed by parentage analysis. In the latter case, for example, if a female was assigned to the mother of 3-year-old bears in 2012, we assumed that the female survived at least until separation of the offspring (i.e., 2010), considering that a cub is not likely to survive its first year without the mother. Because genetic analysis of killed bears was initiated in 1998, the earliest year for the start of bear-year calculation was set at 1997, considering that offspring born in 1997 were potentially killed at 1 year of age in 1998. The latest year for bear-year calculation was set at 2015, as most human-caused mortality of offspring occurs at 1–3 years of age; thus, dead and alive bears born in 2015 can be evaluated in 2016–2018. The number of offspring (aged ≥ 1 year) killed by humans after separation from the mother was calculated by parentage analysis. Bears killed with their mothers were not included in the calculation. The age of offspring was estimated by counting the dental cementum annuli, for bears without birth records based on direct observation. Offspring sex was confirmed by genetic methods (described below) and external characteristics.
DNA extraction and genotyping
DNA was extracted using the DNeasy Blood & Tissue Mini Kit (Qiagen Inc., Tokyo, Japan) for blood and tissue samples, the DNA Extractor FM Kit (Wako, Osaka, Japan) or Isohair Easy (Nippon Gene, Inc., Tokyo, Japan) for hair samples, and the QIAamp DNA Stool Mini Kit (Qiagen Inc.) for feces samples, according to the manufacturers’ protocols. Twenty-one microsatellite markers and one sex marker, amelogenin, were analyzed by multiplex PCR assay under conditions described previously (primers are listed in31). Allele size was determined using an ABI PRISM 310 genetic analyzer (Life Technologies Japan Ltd.). In addition, a mitochondrial DNA haplotype analysis targeting the control region was conducted to select candidate mothers of each individual in the parentage analysis (for details, see30).
Parentage analysis
The parentage analysis was performed using a likelihood-based approach with CERVUS version 3.0.736. The simulation parameters were 10,000 cycles, 150 candidate mothers and fathers per offspring, 40% of candidate parents sampled, and 1% of loci mistyped32. In the first step of the CERVUS analysis, we assigned a parent pair. The confidence level was set at 80%, and no mismatching was allowed in a parents–offspring combination (i.e., mother–father–offspring trio). One mismatch was allowed in a parents–offspring combination obtained at a ≥ 95% confidence level when the same mother and father were selected as the most likely parents (≤ 1 mismatch per pair) in maternity and paternity assignment analyses, respectively. If a parent pair could not be assigned due to a low ( More

1.
MEA, M. E. A. Millennium ecosystem assessment. Ecosystems and Human Well-Being: Biodiversity Synthesis, Published by World Resources Institute, Washington, DC (2005).
2.
Cherlet, M. et al. World Atlas of Desertification: Rethinking Land Degradation and Sustainable Land Management. (Publications Office of the European Union, 2018).

3.
Huang, J., Yu, H., Guan, X., Wang, G. & Guo, R. Accelerated dryland expansion under climate change. Nature Climate Change 6, 166 (2016).
ADS  Article  Google Scholar 

4.
Gaitán, J. J. et al. Biotic and abiotic drivers of topsoil organic carbon concentration in drylands have similar effects at regional and global scales. Ecosystems 22, 1445–1456 (2019).
Article  Google Scholar 

5.
Middleton, N., Stringer, L., Goudie, A. & Thomas, D. The forgotten billion: MDG achievement in the drylands. (UNCCD Secretariat, 2011).

6.
Reynolds, J. F. et al. Global Desertification: Building a. Science for Dryland Development. Science 316, 847–851 (2007).
CAS  PubMed  Google Scholar 

7.
Bestelmeyer, B. T. et al. Land management in the American Southwest: A State-and-Transition approach to Ecosystem Complexity. Environmental Management 34, 38–51 (2004).
Article  Google Scholar 

8.
Maestre, F. T. et al. Structure and functioning of dryland ecosystems in a changing world. Annual review of ecology, evolution, and systematics 47, 215–237 (2016).
Article  Google Scholar 

9.
Maestre, F., Salguero-Gómez, R. & Quero, J. It’s getting hotter in here: determining and projecting the impacts of global change on dryland ecosystems and on the people living in them. Philosophical Transactions of the Royal Society B: Biological Sciences 367, 3062–3075 (2012).
Article  Google Scholar 

10.
Tongway, D. & Hindley, N. L. Landscape Function Analysis: Procedures for monitoring and assessing landscapes. With special reference to Minesites and Rangelands. Vol. 1 (CSIRO, 2004).

11.
Sankey, T. T., Leonard, J. M. & Moore, M. M. Unmanned aerial vehicle− Based rangeland monitoring: examining a century of vegetation changes. Rangeland Ecology & Management 72, 858–863 (2019).
Article  Google Scholar 

12.
Watson, I. W., Novelly, P. E. & Thomas, P. W. E. Monitoring changes in pastoral rangelands – the Western Australian Rangeland Monitoring System (WARMS). The Rangeland Journal 29, 191–205 (2007).
Article  Google Scholar 

13.
White, A. et al. AUSPLOTS rangelands survey protocols manual. (University of Adelaide Press, 2012).

14.
Guerin, G. R. et al. Opportunities for integrated ecological analysis across inland Australia with standardised data from Ausplots Rangelands. PloS one 12 (2017).

15.
Maestre, F. T. et al. Plant species richness and ecosystem multifunctionality in global drylands. Science 335, 214–218 (2012).
ADS  CAS  Article  Google Scholar 

16.
Herrick, J. E., Van Zee, J. W., Havstad, K. M., Burkett, L. M. & Whitford, W. G. Monitoring Manual for Grassland, Shrubland and Savanna Ecosystems. Vol. I (USDA-ARS Jornada Experimental Range, 2005).

17.
Oliva, G. et al. Monitoring drylands: The MARAS system. Journal of Arid Environments 161, 55–63 (2019).
ADS  Article  Google Scholar 

18.
Oliva, G., Escobar, J., Siffredi, G., Salomone, J. & Buono, G. In Monitoring Patagonian Rangelands: The MARAS System. Monitoring Science and Technology Symposium. Denver CO. (eds C. Aguirre-Bravo, P. Pellicane, D. Burns, & S. Draggan) 188–193 (U.S. Dept. Agriculture, Forest Service, 2006).

19.
Oliva, G. et al. Manual para la instalación y lectura de monitores MARAS. Vol. 1 (PNUD, 2011).

20.
Bran, D. et al. Regiones Ecológicas Homogéneas de la Patagonia Argentina., (INTA, 2005).

21.
Oliva, G. et al. Installation Manual for MARAS Monitors: Environmental monitoring for arid and semiarid lands (PNUD, 2011).

22.
Hernández, F., Ríos, C. & Perotto-Baldivieso, H. L. Evolutionary history of herbivory in the Patagonian steppe: The role of climate, ancient megafauna, and guanaco. Quaternary Science Reviews 220, 279–290 (2019).
ADS  Article  Google Scholar 

23.
Borrelli, P. In Ganadería ovina sustentable en la Patagonia Austral Vol. Cap 5 (eds P Borrelli & G Oliva) 131-162 (INTA, 2001).

24.
Cornforth, I. S. & Sinclair, A. G. Fertiliser recommendations for pastures and crops in New Zealand. (New Zealand Ministry of Agriculture, 1984).

25.
Elissalde, N., Escobar, J. & Nakamatsu, V. Inventario y evaluación de pastizales naturales de la zona arida y semiarida de la Patagonia. (INTA Trelew, 2002).

26.
McLaren, C. A. Dry Sheep Equivalents for comparing different classes of livestock. 4 (Department of Primary Industries, State of Victoria, Victoria, 1997).

27.
INIA. Revisión y análisis de las bases históricas y científicas del uso de la equivalencia ovino:bovino “Hacia una nueva equivalencia para ser utilizada en Uruguay”. (INIA, 2012).

28.
SRM, G. R. S. C. A Glossary of terms used in range management: a definition of terms commonly used in range management. (Society for Range Management, 1989).

29.
Oliva, G. et al. The MARAS dataset, vegetation and soil characteristics of dryland rangelands across Patagonia. figshare https://doi.org/10.6084/m9.figshare.c.4789113 (2020).

30.
Tongway, D. Rangeland soil condition assessment manual. (CSIRO. Division of Wildlife and Ecology, 1994).

31.
Daget, P. & Poissonet, J. Une methode d’analyse phytologique des prairies. Ann Argr. France 22, 5–41 (1971).
Google Scholar 

32.
Halloy, S. & Barratt, B. I. P. Patterns of abundance and morphology as indicators of ecosystem status: A meta-analysis. Ecological Complexity 4, 128–147 (2007).
Article  Google Scholar 

33.
Zuloaga, F., Morrone, O. & Belgrano, M. Catálogo de las Plantas Vasculares del Cono Sur. Versión base de datos en sitio web. (Instituto Darwinion, 2009).

34.
Ludwig, J. A., Wilcox, B. P., Breshears, D. D., Tongway, D. J. & Imeson, A. C. Vegetation patches and runoff–erosion as interacting ecohydrological processes in semiarid landscapes. Ecology 86, 288–297 (2005).
Article  Google Scholar 

35.
Walkley, A. & Black, I. A. An examination of the Degtjareff method for determining soil organic matter, and a proposed modification of the chromic acid titration method. Soil science 37, 29–38 (1934).
ADS  CAS  Article  Google Scholar 

36.
Schulte, E. & Hoskins, B. Recommended soil organic matter tests. Recommended Soil Testing Procedures for the North Eastern USA. Northeastern Regional Publication, 52-60 (1995).

37.
López, C., Rial, P., Elissalde, N., Llanos, E. & Behr, S. Grandes paisajes de la Patagonia Argentina. (INTA, 2005).

38.
Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. International journal of climatology 37, 4302–4315 (2017).
ADS  Article  Google Scholar 

39.
Elzinga, C. L., Salzer, D. W. & Willoughby, J. W. Measuring & Monitoring Plant Populations. (BLM National Business Center, 1998).

40.
Magurran, A. E. Measuring Biological Diversity. (Blackwell Publishing, 2004).

41.
Oliva, G. et al. Estado de los Recursos Naturales de la Patagonia Sur 66 (INTA CRPATSU, Trelew, 2017).

42.
Gaitan, J. J. et al. Vegetation structure is as important as climate for explaining ecosystem function across Patagonian rangelands. Journal of Ecology 102, 1419–1428 (2014).
Article  Google Scholar 

43.
Gaitán, J. J. et al. Evaluating the performance of multiple remote sensing indices to predict the spatial variability of ecosystem structure and functioning in Patagonian steppes. Ecological Indicators 34, 181–191 (2013).
Article  Google Scholar 

44.
Gaitán, J. J. et al. Plant species richness and shrub cover attenuate drought effects on ecosystem functioning across Patagonian rangelands. Biology letters 10, 20140673 (2014).
Article  Google Scholar 

45.
Gaitán, J. J. et al. Aridity and Overgrazing Have Convergent Effects on Ecosystem Structure and Functioning in Patagonian Rangelands. Land Degradation & Development 29, 210–218 (2017).
Article  Google Scholar 

46.
Oliva, G., et al.) 1115-1117 (IRC 2016).

47.
Domínguez Díaz, E., Oliva, G. E., Báez Madariaga, J., Suárez Navarro, Á. & Pérez Castillo, C. Efectos del pastoreo holístico sobre la estructura y composición vegetal en praderas naturalizadas de uso ganadero, provincia de Última Esperanza, región de Magallanes, Chile. Anales del Instituto de la Patagonia 46, 17–28 (2018).
Article  Google Scholar 

48.
Borrelli, P. et al. Estándar para la regeneración y la sustentabilidad de los pastizales (GRASS). The Nature Conservancy, OVIS 21 (2013).

49.
T. exchange. RWS, Responsible Wool Standard. 73 (London, 2016).

50.
Cabrera, A. Fitogeografia de la República Argentina. Boletín de la Sociedad Argentina de Botánica 14, 1–42 (1971).
Google Scholar 

51.
León, R., Bran, D., Collantes, M., Paruelo, J. & Soriano, A. Grandes unidades de vegetación de la Patagonia extra andina. Ecologia Austral 8, 125–144 (1998).
Google Scholar 

52.
Luebert, F. & Pliscoff, P. Sinopsis bioclimática y vegetacional de Chile. (Santiago de Chile: Editorial Universitaria, 2006). More

Study area and period
Our study area is the tropics between 15°N–35°S6. We divided the study area into three continents and studied them separately: South America, Africa, and Australasia. Australasia includes Australia and southeast Asia, but excludes southern India. Our results are generated on 0.25° spatial resolution. We classify a cell as forest if it contained at least 50% tree cover (‘forest cover’ in this manuscript) in 1999 according to the dataset from ref. 41. The moisture recycling simulations were carried out for 2003–2014 (‘recent climate’), for which a consistent set of input data was available (see also ref. 11). ‘Late 21st century’ refers to 2071–2100.
Local-scale forest hysteresis
Previous research has shown that tropical forests may have local-scale tipping points at certain mean annual rainfall levels, but are also affected by the seasonality of that rainfall5,6,8. Local-scale tipping points for forest were determined using tree cover data following a method from ref. 6. Using potential analysis42, an empirical stability landscape (as in Fig. 1a) is constructed based on spatial distributions of tree cover against environmental variables such as mean annual rainfall for each continent separately. For each value of the environmental variable, the probability density of tree cover was determined using the MATLAB function ksdensity with a bandwidth of 5%. We applied Gaussian weights to the environmental variable with a standard deviation of 0.05 times the length of the axis of the environmental variable. Local maxima of the resulting probability density function are interpreted as stable states, where we ignored local maxima below a threshold value of 0.004. We used Landsat tree cover data for 2000 on 30 m resolution downloaded for every 0.01°43. We masked out human-used areas, water bodies, and bare ground using the ESA GlobCover land cover dataset for 2009 on 300 m resolution (values 11–30 and ≥190). From the resulting dataset we randomly sampled one million locations for each continent and used them to construct the stability landscapes6 against mean annual rainfall and average MCWD. MCWD is the cumulative difference between evapotranspiration and rainfall using monthly averages of those fluxes calculated for each calendar year44. It is set to zero when monthly rainfall exceeds monthly evapotranspiration and becomes more negative with an increasing water deficit. Following ref. 11, for both mean annual rainfall and MCWD, we took monthly data from the GLDAS 2.0 dataset45 for 1970–1999 so the 30-year period leading up to the land-cover sample (for the year 2000) was used.
Forest evapotranspiration
To estimate the fraction of evapotranspiration attributable to forest cover we used the large-scale hydrological model PCR-GLOBWB, run at 0.5° resolution46. Per grid cell, the model simulates evapotranspiration for a range of land-cover types. Here, we are specifically interested in the evapotranspiration of forests, or ‘tall natural vegetation’ in PCR-GLOBWB. Note that we here account for both forest transpiration and canopy interception evaporation instead of, as in ref. 11, only transpiration.
PCR-GLOBWB computes the water balance in two soil layers and a groundwater layer. Soil type, fractional area of saturated soil, and the spatiotemporal distribution of groundwater depth are accounted for (see refs. 46,47). It includes six land-cover types, with spatially varying parameters46: tall and short natural vegetation, pasture, rainfed crops, and paddy and non-paddy irrigated crops. The model was forced with WATCH Forcing Data ERA-Interim precipitation, temperature, and reference potential evapotranspiration for 1979–201448. We used monthly evapotranspiration output of PCR-GLOBWB, implying that we assume that forest component of evapotranspiration remains equal within each month. For detailed model descriptions and validation, we refer to earlier studies11,46,49,50.
Atmospheric moisture tracking
As an essential step in estimating the forest–rainfall feedback, we determined where the moisture from enhanced evapotranspiration precipitates again by using atmospheric moisture tracking. The method for atmospheric moisture tracking resembles that in ref. 11. Apart from the expansion of the study area to the entire tropics, a notable difference is that we here used ERA5 reanalysis data rather than ERA-Interim, meaning that the simulations were based on finer resolution input data (i.e. on 0.25° instead of 0.75° for spatial resolution, and 1 h instead of 3 h for temporal resolution). ERA5 has better performance than ERA-Interim regarding wind fields and rainfall, especially in the tropics51,52,53. Below we summarize the method (see also ref. 11).
We used a Lagrangian method of moisture tracking that is based on previous studies11,54,55,56 that track parcels of evaporated moisture forward through the atmosphere to their subsequent precipitation location. Moisture particles that enter the atmosphere are assigned a random location within the 0.25° grid cell and random starting height in the atmosphere scaled with the humidity profile, and their trajectories are then tracked through the atmosphere. The trajectories are forced with the three-dimensional ERA5 reanalysis estimates of wind speed and direction, which were linearly interpolated at every time step of 0.25 h. Water particles in the atmosphere have an equal probability of raining out, regardless of vertical position. Rainfall A (mm per time step) at location x,y and time t that has evaporated from any location of release in any cell is ref. 56

$$A_{x,y,t} = P_{x,y,t}frac{{W_{{mathrm{parcel}},t}E_{{mathrm{source}},t}}}{{{mathrm{TPW}}_{x,y,t}}},$$
(1)

where P is rainfall in mm per time step, Wparcel is the water in the tracked parcel in mm, Esource is its fraction of water that evapotranspired from the source, and TPW is the precipitable water in the atmospheric water column in mm. Every time step, the amount of water in the parcel is updated based on evapotranspiration ET into the parcel and rainfall P from it:

$$W_{{mathrm{parcel}},t} = W_{{mathrm{parcel}},t – 1} + ({mathrm{ET}}_{x,y,t} – P_{x,y,t})frac{{W_{{mathrm{parcel}},t – 1}}}{{{mathrm{TPW}}_{x,y,t}}}.$$
(2)

The fraction of water in the parcel that has evapotranspired from the source then becomes

$$E_{{mathrm{source}},t} = frac{{E_{{mathrm{source}},t – 1}W_{{mathrm{parcel}},t – 1} – A_{x,y,t}}}{{W_{{mathrm{parcel}},t}}}.$$
(3)

Thus, the amount of water that was tracked from the source location decreases with precipitation along its trajectory. Parcels were followed until either less than 5% of its original amount was left in the atmosphere, or the tracking time was 30 days. Any moisture remaining in the parcel when the trajectories end is assumed to rain out over non-land areas, thus not contributing to our analysis. We analysed each continent separately for reasons of computability. However, small moisture flows between forests in different continents can be expected, as has been simulated for flows from Africa to the Amazon57. Over all land points, ERA5 hourly evapotranspiration is linearly interpolated to every 0.25-h time step. The moisture flow mij in mm per month linking evapotranspiration in cell i to rainfall in cell j where (left[ {x,y} right] , {it{epsilon }} , j) over the course of a given month becomes

$$m_{ij} = mathop {sum}limits_{t = 0}^{{mathrm{month}}} {A_{j,t}} cdot frac{{{mathrm{ET}}_{i,t}}}{{W_{i,t}}},$$
(4)

where ETi,t is the evapotranspiration in mm per time step, and Wi,t is the tracked amount of water from source cell i at time step t.
At continental scales, evaporated moisture can rain down and re-evaporate multiple times. This also means that forest evapotranspiration can enhance rainfall multiple times. We accounted for this ‘cascading moisture recycling’ following refs. 11,14, in which the rainfall attributed to an upwind forest source is subsequently tracked upon re-evaporation. After six re-evapotranspiration cycles, cascading moisture recycling has decreased to practically zero11. Therefore, following ref. 11, seven iterations of cascading moisture recycling were performed.
Hysteresis experiments
We determined the hysteresis of tropical forests through a series of iterative runs; each one started either from a fully forested continent or from a fully deforested continent. We simulated the hypothetical mean annual rainfall levels across the (tropical part of the) continent given this initial condition, that is, rainfall without any forest evapotranspiration or rainfall including evapotranspiration from an entirely forested continent. Next, based on the empirical bifurcation diagrams (i.e. nonforest, bistable forests, and stable forests in each continent are determined based on the bistability ranges shown in Supplementary Figs. 1−3), we determined either the minimal distribution of tropical forest (in case of a no-forest initial condition, based on the higher end of the bistability range) or the maximal distribution (in case of a fully forested initial condition, based on the lower end of the bistability range) at these rainfall levels. Thus, in the simulations with an empty initial condition, only stable forests (green in Fig. 1) would regrow; in those with a full initial condition, both stable and bistable forests (green and yellow in Fig. 1) would disappear. Because the resulting new distribution of forest would generate different levels of rainfall, the procedure was repeated with the respective forest distribution as initial condition. This occurred until rainfall levels had (practically) stabilized between iterative runs (up to three iterations).
We assumed that no other vegetation type replaces the forest in order to show the theoretically possible distributions of tropical forests. This may lead to an overestimation of the actual effects of forest on rainfall, especially if forests would be replaced by highly transpiring crops58. Furthermore, land-cover changes will alter wind patterns and therefore the expected coupling between forests through evapotranspiration and rainfall59. Fossil fuel emissions not only alter the climate, but the emitted CO2 also fertilizes plants. This increases trees’ water-use efficiency, reducing their water demand, but it also increases biomass production60. The effects of this CO2 fertilization on the water cycle might be small61, but its net effects on tropical forest hysteresis remains uncertain.
For display of Fig. 2, the resolution of rainfall values was increased by a factor of 2, to 0.125° and smoothed using a two-dimensional Gaussian kernel with a standard deviation of 0.5°.
Climate-change scenario
As the estimate of late 21st-century rainfall conditions, we used the rainfall output from the SSP5-8.5 scenario simulations by seven available CMIP6 models62. These models are BCC-CSM2-MR, CanESM5, CNRM-CM6-1, CNRM-ESM2, IPSL-CM6A-LR, MRI-ESM2.0, and UKESM1.0-LL. We took the mean across the model outputs for the annual rainfall values for 2071–2100. The scenario is considered a severe climate-change scenario. Because the moisture tracking model is forced with atmospheric reanalysis data, we assumed that (forest-induced) moisture flows in the scenario are the same as in the period of our atmospheric simulations (2003–2014).
We assumed that a tipping point from a forested to a nonforested state occurs when mean annual rainfall in a forested cell (forest cover ≥ 50%) is currently (2003–2014) above the lower tipping point (Supplementary Figs. 1–3), but is reduced to below the lower tipping point in the climate-change scenario. Similarly, a tipping point from a nonforested to a forested state occurs when mean annual rainfall in a nonforested cell (forest cover  More

Scopes and databases
The CEAP dataset comprises all the thermal power plants operating in China, totalling 2,714 plants (or 6,267 units), from 2014 to 2017 in 26 provinces and 4 municipalities (except Hong Kong, Macao, Taiwan and Tibet; Table 1). The thermal power plants produce electricity by combusting a variety of fossil energies, which fall into 4 categories: coal, gas plus oil, biomass and others (detailed in Table 2).
Table 1 China’s thermal power plants in CEAP.
Full size table

Table 2 Fuel type descriptions.
Full size table

The CEAP dataset integrates two databases, i.e., the CEMS data and unit-specific information. The CEMS data—the direct, real-time measurements of stack gas concentrations of PM, SO2 and NOX from China’s power plant stacks—are monitored by China’s CEMS network and reported to the China Ministry of Ecology and Environment (MEE; http://www.envsc.cn/). The CEMS data are recorded on a source and hourly basis. In total, the CEMS dataset covers 4,622 emission sources (i.e., power plant stacks) associated with 5,606 units (accounting for 98% of China’s thermal power capacity), 35,064 hours from 2014 to 2017, and 3 air pollutants (i.e., PM, SO2 and NOX) for each source-hour sample (Table 3). The MEE has also provided stack-specific information (regarding latitude and longitude, heights, temperature, diameter, etc.; http://permit.mee.gov.cn/).
Table 3 CEMS coverage of China’s thermal power plant units or stacks in CEAP.
Full size table

Unit-specific information is also derived from the MEE, involving activity levels (energy consumption and power generation), operating capacities, geographic allocations and pollution control equipment (particularly the types and removal efficiencies) at a yearly frequency. Due to data availability, the unit information is available only until 2016, and the activity levels for 2017 are projected following the overall trends in provincial thermal power generation between 2016 and 2017 (which are available in the China Energy Statistical Yearbooks26), under the assumption that new units constructed in 2017 have the same structures of installed capacities, energy uses and regions as those of the existing units in 2016.
With a combination of the two datasets, the CEAP dataset provides nationwide, plant-level, dynamic PM, SO2 and NOX emissions from China’s thermal power plants from 2014 to 2017. Relative to existing inventories, the CEAP dataset is innovative in that it incorporates comprehensive real CEMS-measured emission data, avoiding the use of average emission factors and the associated operational assumptions and uncertain parameters.
Pre-processing of CEMS data
We have been exclusively granted access to the data from China’s CEMS network. Generally, the CEMS consists of a sampling system (for filtering and sampling flue gas), an online analytical component (for monitoring flue gas parameters, particularly emission concentrations) and a data processing system (for collecting, processing and reporting monitoring data)27,28. According to the GB13223-2003 regulation29, the CEMS network should cover all power plant furnaces that burn coal (except stoker and spreader stoker) and oil and generate >65 tons of steam each hour, as well as those that burn pulverized coal and gas. Thus, some power plants have not yet been incorporated into the CEMS network (accounting for 3–4% of the total thermal power capacity from 2014 to 2017) because their furnaces did not meet the requirements necessary to install a CEMS. For the power plants outside the CEMS network, we assume their stack concentrations are similar to the averages of the units with similar fuel types and similar regions within the CEMS network.
To guarantee the reliability of CEMS data, China’s government has made great efforts in developing specific regulations and technical guidelines for power plants and local entities to follow and supervise, respectively24,28,30,31,32. These official documents elaborate on all the processes required to regulate the CEMS network, including not only CEMS installation, operation, inspection, maintenance and repair but also CEMS data collection, processing, reporting, analysis and storage28,32,33. Since 2014, all state-monitored companies have been mandated to report their CEMS data to the local governments through a series of online platforms for different provinces (listed in Supplementary Table 1). Local entities have random onsite inspections to check the truthfulness of the reported results on at least a quarterly basis23,24,28,32,34; this system enables a comparison of CEMS data across different firms to explore potential outliers and abnormalities and prevent data manipulation28,35. Then, the governments release the inspection results to the public through the same online platforms (listed in Supplementary Table 1)24,36,37. Severe financial penalties and criminal punishments can be imposed on firms that adopt data manipulation (in terms of deleting, distorting and forging CEMS data, for example)38,39.
The malfunction of CEMS monitors may also introduce large uncertainty to CEMS data during the processes of operation (indication errors, span drift, zero drift, etc.), maintenance (particularly the failure to perform calibration and reference tests) and data reporting (invalid data communication, data missing, etc.)24,28. Accordingly, each power plant is required to make at least one A-, B- and C-grade overhaul for 32–80, 14–50 and 9–30 days per 4–6, 2–3 and 1 year(s), respectively, as well as one D-grade overhaul (if needed) for 5–15 days per year, to check, maintain and upgrade its technologies, thereby reducing measurement uncertainty40. During these overhauls, CEMS operators conduct CEMS calibration (i.e., zero and span calibration), maintenance procedures (e.g., examining and cleaning major CEMS components and replacing or upgrading parts, if necessary, such as optical lens, filter and sampling meter) and a reference test (i.e., relative accuracy test audit). Furthermore, third-party operators examine CEMS operation and maintenance routines, to guarantee standardized CEMS operation and facilitate improvement in CEMS data accuracy27,28,31. All the related activities should be documented according to standardized requirement contents27,28. Even with the aforementioned efforts, there is still a small proportion of nulls and outliers in the CEMS database, which represent 1% and 0.1% of the total operating hours, respectively, from 2014 to 2017. We treat these samples seriously by following the relevant official documents, which have been released by China’s government. Table 4 provides the treatment methods for nulls or zeros, which can be divided into 3 types based on duration. On the one hand, we consider nulls and/or zeros that span at least 5 successive days as a downtime or overhaul and omit them in the estimation, according to the regulation27. On the other hand, missing data lasting  24 hours are assumed around the valid values near the time and set to the monthly averages27:

$${hat{C}}_{s,i,y,m,h}={bar{C}}_{s,i,y,m,bullet }$$
(1)

where ({C}_{s,i,y,m,h}) denotes the stack gas concentrations of pollutant s emitted by unit i for year y, month m and hour h (i.e., the actual measurements monitored by the CEMS network), defined as the amount of pollutants per unit of emitted stack gas (g m−3)41,42; ({widehat{C}}_{s,i,y,m,h}) is the imputation for the missing data ({C}_{s,i,y,m,h}); ({bar{C}}_{s,i,y,m,.}) is the mean of the hourly valid values for the same pollutant, unit, year and month as ({C}_{s,i,y,m,h}). In contrast, the missing data for 1–24 hour(s) are interpolated with the arithmetic averages of the two nearest valid points before and after them27,43:

$${widehat{C}}_{s,i,y,m,h}=frac{{C}_{s,i,y,m,h-l}+{C}_{s,i,y,m,h+q}}{2}$$
(2)

where ({C}_{s,i,y,m,h-l}) and ({C}_{s,i,y,m,h{rm{+}}q}) represent the nearest last known measurements (l hour(s) before) and next known measurements (q hour(s) after), respectively, for the missing data ({C}_{s,i,y,m,h}), namely, the series data ({C}_{s,i,y,m,h-l+1}),…, ({C}_{s,i,y,m,h}),…, ({C}_{s,i,y,m,h+q-1}) are all missing values. Furthermore, we treat the measurements that are out of the measurement ranges of CEMS instruments (outside of which the data are unreliable30,44; detailed in Supplementary Table 2) as abnormal data and process them in a similar way to nulls according to the official regulation27.
Table 4 Treatment methods for nulls and the relevant official documents.
Full size table

CEMS-based estimation of emission factors and absolute emissions
The introduction of real CEMS-monitored measurements provides a direct estimation for emission factors on a source and hourly basis, avoiding the use of average emission factors with many assumptions and uncertain parameters17,42,44.

$$E{F}_{s,i,y,m,h}={C}_{s,i,y,m,h}{V}_{i,y}$$
(3)

In Eq. (3), (E{F}_{s,i,y,m,h}) indicates the emission factor, defined as the amount of emissions per unit of fuel use (in g kg−1 for solid or liquid fuel and in g m−3 for gas fuel), and ({V}_{i,y}) is the theoretical flue gas rate, defined as the expected volume of flue gas per unit of fuel use under standard production conditions (m3 kg−1 for solid or liquid fuel and m3 m−3 for gas fuel)42, which was estimated by the China Pollution Source Census (2011)45 based on sufficient field measurements (detailed in Table 5). Based on Eq. (3), abated emission factors can be directly obtained even without the use of removal efficiencies and the relevant parameters, because CEMS monitors the gas concentrations at stacks after the effect of control equipment (if any).
Table 5 Theoretical flue gas rate.
Full size table

Notably, recent clean air policies (particularly different emissions standards) target stack concentrations, such that a large proportion of missing data exist regarding other measurements (particularly flue gas rates, with missing data accounting for 34.62%, 31.91%, 29.97% and 42.96% of the total samples in 2014, 2015, 2016 and 2017, respectively). Accordingly, we introduce theoretical flue gas rates into the estimation to avoid significant underestimation of the actual volume when there are too many missing data values46. In addition, the adoption of theoretical flue gas rates can address flue gas leakage, a common problem in power plants that greatly distorts the real flue gas volume46. The theoretical flue gas rates are derived from the China Pollution Source Census, with values varying across operating capacities, fuel types and boiler types42,45. Thus, the actual volume of flue gas is computed in terms of the theoretical flue gas rate times actual fuel consumption.
The absolute emissions of PM, SO2 and NOX from individual power plants can be estimated in terms of the emission factors times the activity levels21:

$${E}_{s,i,y,m}=E{F}_{s,i,y,m}{A}_{i,y,m}$$
(4)

where ({E}_{s,i,y,m}) represents the air pollution emissions (g); and ({A}_{i,y,m}) is the activity data, i.e., the amount of fuel use (kg for solid or liquid fuel and m3 for gas fuel). In the CEAP dataset, power plant emissions are estimated on a monthly basis (the smallest scale for activity data), in which the yearly unit-level activity data are allocated at the monthly scale using the monthly province-level thermal power generation as weights16:

$${A}_{i,y,m}=frac{{F}_{{p}_{i},y,m}}{{sum }_{m=1}^{12}{F}_{{p}_{i},y,m}}{A}_{i,y}$$
(5)

where ({F}_{{p}_{i},y,m}) denotes the thermal power generation by province Pi, which is obtained from the Chinese Energy Statistics Yearbook26, and ({p}_{i}) indicates the province where unit i is located. More