Philippa Kaur

1.
Cryan, J. F. & Dinan, T. G. Mind-altering microorganisms: the impact of the gut microbiota on brain and behaviour. Nat. Rev. Neurosci. 13, 701–712 (2012).
CAS  PubMed  Article  Google Scholar 
2.
Sherwin, E., Bordenstein, S. R., Quinn, J. L., Dinan, T. G. & Cryan, J. F. Microbiota and the social brain. Science (80-) 366, eaar2016 (2019).
CAS  Article  Google Scholar 

3.
Heijtz, R. D. et al. Normal gut microbiota modulates brain development and behavior. Proc. Natl. Acad. Sci. 108, 3047–3052 (2011).
ADS  CAS  Article  Google Scholar 

4.
Foster, J. A. & McVey Neufeld, K.-A. Gut–brain axis: how the microbiome influences anxiety and depression. Trends Neurosci. 36, 305–312 (2013).
CAS  PubMed  Article  Google Scholar 

5.
Clarke, G. et al. The microbiome-gut–brain axis during early life regulates the hippocampal serotonergic system in a sex-dependent manner. Mol. Psychiatry 18, 666–673 (2013).
CAS  PubMed  Article  Google Scholar 

6.
Desbonnet, L., Clarke, G., Shanahan, F., Dinan, T. G. & Cryan, J. F. Microbiota is essential for social development in the mouse. Mol. Psychiatry 19, 146–148 (2014).
CAS  PubMed  Article  Google Scholar 

7.
Hoban, A. E. et al. The microbiome regulates amygdala-dependent fear recall. Mol. Psychiatry 23, 1134–1144 (2018).
CAS  PubMed  Article  Google Scholar 

8.
Magnusson, K. R. et al. Relationships between diet-related changes in the gut microbiome and cognitive flexibility. Neuroscience 300, 128–140 (2015).
CAS  PubMed  Article  Google Scholar 

9.
Ogbonnaya, E. S. et al. Adult Hippocampal Neurogenesis Is Regulated by the Microbiome. Biol. Psychiat. 78, e7–e9 (2015).
PubMed  Article  Google Scholar 

10.
Gareau, M. G. et al. Bacterial infection causes stress-induced memory dysfunction in mice. Gut 60, 307–317 (2011).
PubMed  Article  Google Scholar 

11.
Stilling, R. M. et al. The neuropharmacology of butyrate: the bread and butter of the microbiota-gut–brain axis?. Neurochem. Int. 99, 110–132 (2016).
CAS  PubMed  Article  Google Scholar 

12.
Davidson, G. L., Raulo, A. & Knowles, S. C. L. Identifying microbiome-mediated behaviour in wild vertebrates. Trends Ecol. Evol. https://doi.org/10.1016/j.tree.2020.06.014 (2020).
Article  PubMed  Google Scholar 

13.
Davidson, G. L., Cooke, A. C., Johnson, C. N. & Quinn, J. L. The gut microbiome as a driver of individual variation in cognition and functional behaviour. Philos. Trans. R. Soc. B Biol. https://doi.org/10.1098/rstb.2017.0286 (2018).
Article  Google Scholar 

14.
Morand-Ferron, J., Cole, E. F. & Quinn, J. L. Studying the evolutionary ecology of cognition in the wild: a review of practical and conceptual challenges. Biol. Rev. 91, 367–389 (2016).
PubMed  Article  Google Scholar 

15.
Stephens, D. W. & Krebs, J. R. Foraging Theory (Princeton University Press, Princeton, 2019).https://doi.org/10.2307/j.ctvs32s6b
Google Scholar 

16.
De Filippo, C. et al. Impact of diet in shaping gut microbiota revealed by a comparative study in children from Europe and rural Africa. Proc. Natl. Acad. Sci. U. S. A. 107, 14691–14696 (2010).
ADS  PubMed  PubMed Central  Article  Google Scholar 

17.
Gillingham, M. A. F. et al. Offspring microbiomes differ across breeding sites in a panmictic species. Front. Microbiol. https://doi.org/10.3389/fmicb.2019.00035 (2019).
Article  PubMed  PubMed Central  Google Scholar 

18.
Costa, S., Lopes, I., Proença, D. N., Ribeiro, R. & Morais, P. V. Diversity of cutaneous microbiome of Pelophylax perezi populations inhabiting different environments. Sci. Total Environ. 572, 995–1004 (2016).
ADS  CAS  PubMed  Article  Google Scholar 

19.
Knutie, S. A., Chaves, J. A. & Gotanda, K. M. Human activity can influence the gut microbiota of Darwin’s finches in the Galapagos Islands. Mol. Ecol. 28, 2441–2450 (2019).
PubMed  Article  Google Scholar 

20.
Amato, K. R. et al. The gut microbiota appears to compensate for seasonal diet variation in the Wild Black Howler Monkey (Alouatta pigra). Microb. Ecol. 69, 434–443 (2014).
PubMed  Article  CAS  PubMed Central  Google Scholar 

21.
Hicks, A. L. et al. Gut microbiomes of wild great apes fluctuate seasonally in response to diet. Nat. Commun. 9, 1786 (2018).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

22.
Maurice, C. F. et al. Marked seasonal variation in the wild mouse gut microbiota. ISME J. 9, 2423–2434 (2015).
CAS  PubMed  PubMed Central  Article  Google Scholar 

23.
Lozupone, C. A., Stombaugh, J. I., Gordon, J. I., Jansson, J. K. & Knight, R. Diversity, stability and resilience of the human gut microbiota. Nature 489, 220–230 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

24.
Pan, D. & Yu, Z. Intestinal microbiome of poultry and its interaction with host and diet. Gut Microbes 5, 108–119 (2013).
PubMed  PubMed Central  Article  Google Scholar 

25.
Teyssier, A. et al. Diet contributes to urban-induced alterations in gut microbiota: experimental evidence from a wild passerine. Proc. R. Soc. B Biol. Sci. https://doi.org/10.1098/rspb.2019.2182 (2020).
Article  Google Scholar 

26.
David, L. A. et al. Diet rapidly and reproducibly alters the human gut microbiome. Nature 505, 559–563 (2014).
ADS  CAS  PubMed  Article  Google Scholar 

27.
Clarke, S. F. et al. Exercise and associated dietary extremes impact on gut microbial diversity. Gut 63, 1913–1920 (2014).
CAS  PubMed  Article  Google Scholar 

28.
Fava, F. et al. The type and quantity of dietary fat and carbohydrate alter faecal microbiome and short-chain fatty acid excretion in a metabolic syndrome ‘at-risk’ population. Int. J. Obes. 37, 216–223 (2013).
CAS  Article  Google Scholar 

29.
Wu, G. D. et al. Comparative metabolomics in vegans and omnivores reveal constraints on diet-dependent gut microbiota metabolite production. Gut 65, 63–72 (2016).
CAS  PubMed  Article  Google Scholar 

30.
Zimmer, J. et al. A vegan or vegetarian diet substantially alters the human colonic faecal microbiota. Eur. J. Clin. Nutr. 66, 53–60 (2012).
CAS  PubMed  Article  Google Scholar 

31.
Youngblut, N. D. et al. Host diet and evolutionary history explain different aspects of gut microbiome diversity among vertebrate clades. Nat. Commun. 10, 2200 (2019).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

32.
Hird, S. M., Sánchez, C., Carstens, B. C. & Brumfield, R. T. Comparative gut microbiota of 59 neotropical bird species. Front. Microbiol. 6, 1403 (2015).
PubMed  PubMed Central  Article  Google Scholar 

33.
Kartzinel, T. R., Hsing, J. C., Musili, P. M., Brown, B. R. P. & Pringle, R. M. Covariation of diet and gut microbiome in African megafauna. Proc. Natl. Acad. Sci. U. S. A. 116, 23588–23593 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

34.
Bolnick, et al. The ecology of individuals: incidence and implications of individual specialization. Am. Nat. 161, 1–28 (2003).
MathSciNet  PubMed  Article  Google Scholar 

35.
Li, W., Dowd, S. E., Scurlock, B., Acosta-Martinez, V. & Lyte, M. Memory and learning behavior in mice is temporally associated with diet-induced alterations in gut bacteria. Physiol. Behav. 96, 557–567 (2009).
CAS  PubMed  Article  PubMed Central  Google Scholar 

36.
Ezra-Nevo, G., Henriques, S. F. & Ribeiro, C. The diet-microbiome tango: how nutrients lead the gut brain axis. Curr. Opin. Neurobiol. 62, 122–132 (2020).
CAS  PubMed  Article  PubMed Central  Google Scholar 

37.
Psaltopoulou, T. et al. Mediterranean diet, stroke, cognitive impairment, and depression: a meta-analysis. Ann. Neurol. 74, 580–591 (2013).
PubMed  Article  PubMed Central  Google Scholar 

38.
Carlson, A. L. et al. Infant gut microbiome associated with cognitive development. Biol. Psychiatry 83, 148–159 (2018).
PubMed  Article  PubMed Central  Google Scholar 

39.
Dunn, J. C., Cole, E. F. & Quinn, J. L. Personality and parasites: Sex-dependent associations between avian malaria infection and multiple behavioural traits. Behav. Ecol. Sociobiol. 65, 1459–1471 (2011).
Article  Google Scholar 

40.
Cole, E. F., Morand-Ferron, J., Hinks, A. E. & Quinn, J. L. Cognitive ability influences reproductive life history variation in the wild. Curr. Biol. 22, 1808–1812 (2012).
CAS  PubMed  Article  Google Scholar 

41.
Seed, A. & Mayer, C. Problem Solving. in APA handbook of comparative psychology: Perception, learning, and cognition, Vol. 2 601–625 (American Psychological Association, 2017).

42.
Cole, E. F., Cram, D. L. & Quinn, J. L. Individual variation in spontaneous problem-solving performance among wild great tits. Anim. Behav. 81, 491–498 (2011).
Article  Google Scholar 

43.
Morand-Ferron, J., Cole, E. F., Rawles, J. E. C. & Quinn, J. L. Who are the innovators? A field experiment with 2 passerine species. Behav. Ecol. 22, 1241–1248 (2011).
Article  Google Scholar 

44.
Quinn, J. L., Cole, E. F., Reed, T. E. & Morand-Ferron, J. Environmental and genetic determinants of innovativeness in a natural population of birds. Philos. Trans. R. Soc. Biol. B Sci. 371, 20150184 (2016).
Article  CAS  Google Scholar 

45.
Ducatez, S., Clavel, J. & Lefebvre, L. Ecological generalism and behavioural innovation in birds: technical intelligence or the simple incorporation of new foods?. J. Anim. Ecol. 84, 79–89 (2015).
PubMed  Article  Google Scholar 

46.
Reader, S. M. & MacDonald, K. Environmental variability and primate behavioural flexibility. Anim. Innov. https://doi.org/10.1093/acprof:oso/9780198526223.003.0004 (2012).
Article  Google Scholar 

47.
Biard, C. et al. Growing in cities: an urban penalty for wild birds? A study of phenotypic differences between urban and rural great tit chicks (Parus major). Front. Ecol. Evol. https://doi.org/10.3389/fevo.2017.00079 (2017).
Article  Google Scholar 

48.
Teyssier, A. et al. Inside the guts of the city: urban-induced alterations of the gut microbiota in a wild passerine. Sci. Total Environ. 612, 1276–1286 (2018).
ADS  CAS  PubMed  Article  Google Scholar 

49.
Escallón, C., Belden, L. K. & Moore, I. T. The cloacal microbiome changes with the breeding season in a wild bird. Integr. Org. Biol. https://doi.org/10.1093/iob/oby009 (2019).
Article  Google Scholar 

50.
Waite, D. W. & Taylor, M. W. Characterizing the avian gut microbiota: membership, driving influences, and potential function. Front. Microbiol. https://doi.org/10.3389/fmicb.2014.00223 (2014).
Article  PubMed  PubMed Central  Google Scholar 

51.
Singh, R. K. et al. Influence of diet on the gut microbiome and implications for human health. J. Transl. Med. 15, 73 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

52.
Knutie, S. A. Food supplementation affects gut microbiota and immunological resistance to parasites in a wild bird species. J. Appl. Ecol. 57, 536–547 (2020).
CAS  Article  Google Scholar 

53.
Caporaso, J. G. et al. QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

54.
Veľký, M., Kaňuch, P. & Krištín, A. Food composition of wintering great tits (Parus major): habitat and seasonal aspects. Folia Zool. 60, 228–236 (2011).
Article  Google Scholar 

55.
Phillips, J. N., Berlow, M. & Derryberry, E. P. The effects of landscape urbanization on the gut microbiome: an exploration into the gut of urban and rural white-crowned Sparrows. Front. Ecol. Evol. 6, 148 (2018).
Article  Google Scholar 

56.
Rosshart, S. P. et al. Wild mouse gut microbiota promotes host fitness and improves disease resistance. Cell 171, 1015–1028 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

57.
Griffin, A. S. & Guez, D. Innovation and problem solving: a review of common mechanisms. Behav. Process. 109, 121–134 (2014).
Article  Google Scholar 

58.
Alcock, J., Maley, C. C. & Aktipis, C. A. Is eating behavior manipulated by the gastrointestinal microbiota? Evolutionary pressures and potential mechanisms. BioEssays 36, 940–949 (2014).
PubMed  PubMed Central  Article  Google Scholar 

59.
Maniscalco, J. W. & Rinaman, L. Vagal interoceptive modulation of motivated behavior. Physiology 33, 151–167 (2018).
CAS  PubMed  PubMed Central  Article  Google Scholar 

60.
Bruce-Keller, A. J. et al. Obese-type gut microbiota induce neurobehavioral changes in the absence of obesity. Biol. Psychiatry 77, 607–615 (2015).
PubMed  Article  PubMed Central  Google Scholar 

61.
Greyson-Gaito, C. J. et al. Into the wild: microbiome transplant studies need broader ecological reality. Proc. R. Soc. B Biol. Sci. https://doi.org/10.1098/rspb.2019.2834 (2020).
Article  Google Scholar 

62.
Roager, H. M. & Dragsted, L. O. Diet-derived microbial metabolites in health and disease. Nutr. Bull. 44, 216–227 (2019).
Article  Google Scholar 

63.
Möhle, L. et al. Ly6Chi monocytes provide a link between antibiotic-induced changes in gut microbiota and adult hippocampal neurogenesis. Cell Rep. 15, 1945–1956 (2016).
PubMed  Article  CAS  PubMed Central  Google Scholar 

64.
Cryan, J. F. et al. The microbiota-gut-dbrain axis. Physiol. Rev. 99, 1877–2013 (2019).
CAS  PubMed  Article  PubMed Central  Google Scholar 

65.
Heintz-Buschart, A. & Wilmes, P. Human gut microbiome: function matters. Trends Microbiol. 26, 563–574 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

66.
O’Connor, R. J. Identification guide to European Passerines L. Svensson. Auk 102, (1985).

67.
Khan, G., Kangro, H. O., Coates, P. J. & Heath, R. B. Inhibitory effects of urine on the polymerase chain reaction for cytomegalovirus DNA. J. Clin. Pathol. 44, 360–365 (1991).
CAS  PubMed  PubMed Central  Article  Google Scholar 

68.
Eisenhofer, R. et al. Contamination in low microbial biomass microbiome studies: issues and recommendations. Trends Microbiol. 27, 105–117 (2019).
CAS  PubMed  Article  Google Scholar 

69.
Perrins, C. M. Tits and their caterpillar food supply. Ibis (Lond. 1859) 133, 49–54 (1991).
Article  Google Scholar 

70.
Serrano-Davies, E., O’Shea, W. & Quinn, J. L. Individual foraging preferences are linked to innovativeness and personality in the great tit. Behav. Ecol. Sociobiol. 71, 161 (2017).
Article  Google Scholar 

71.
Aplin, L. M., Sheldon, B. C. & McElreath, R. Conformity does not perpetuate suboptimal traditions in a wild population of songbirds. Proc. Natl. Acad. Sci. U. S. A. 114, 7830–7837 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

72.
O’Shea, W., Serrano-Davies, E. & Quinn, J. L. Do personality and innovativeness influence competitive ability? An experimental test in the great tit. Behav. Ecol. 28, 1435–1444 (2017).
Article  Google Scholar 

73.
Shutt, J. D. et al. Gradients in richness and turnover of a forest passerine’s diet prior to breeding: a mixed model approach applied to faecal metabarcoding data. Mol. Ecol. 29, 1199–1213 (2020).
PubMed  Article  Google Scholar 

74.
Crouch, N. M. A., Lynch, V. M. & Clarke, J. A. A re-evaluation of the chemical composition of avian urinary excreta. J. Ornithol. 161, 17–24 (2020).
Article  Google Scholar 

75.
Fouhy, F. et al. Perinatal factors affect the gut microbiota up to four years after birth. Nat. Commun. 10, 1517 (2019).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

76.
Konstantinidis, K. T. & Tiedje, J. M. Genomic insights that advance the species definition for prokaryotes. Proc. Natl. Acad. Sci. U. S. A. 102, 2567–2572 (2005).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

77.
McMurdie, P. J. & Holmes, S. Phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8, e61217 (2013).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

78.
R Core Team. R: A language and environment for statistical computing. R Found. Stat. Comput. Vienna, Austria. https://www.R-project.org/ (2014).

79.
Bokulich, N. A. et al. Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing. Nat. Methods 10, 57–59 (2013).
CAS  PubMed  Article  Google Scholar 

80.
Di Rienzi, S. C. et al. The human gut and groundwater harbor non-photosynthetic bacteria belonging to a new candidate phylum sibling to Cyanobacteria. Elife https://doi.org/10.7554/eLife.01102 (2013).
Article  PubMed  PubMed Central  Google Scholar 

81.
Bates, D., Maechler, M., Bolker, B. & Walker, S. lme4: linear mixed-effects models using Eigen and S4. R package version 1.1–7, https://CRAN.R-project.org/package=lme4. R Packag. version (2014).

82.
Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest package: tests in linear mixed effects models. J. Stat. Softw. https://doi.org/10.18637/jss.v082.i13 (2017).
Article  Google Scholar 

83.
Zakrzewski, M. et al. Calypso: a user-friendly web-server for mining and visualizing microbiome-environment interactions. Bioinformatics https://doi.org/10.1093/bioinformatics/btw725 (2017).
Article  PubMed  Google Scholar  More

1.
Ripple, W. J. et al. Status and ecological effects of the world’s largest carnivores. Science 343, 1241484–1241484 (2014).
PubMed  Article  CAS  Google Scholar 
2.
Loucks, C. J. et al. Giant Pandas in a Changing Landscape (American Association for the Advancement of Science, Washington, 2001).
Google Scholar 

3.
Derocher, A. E. et al. Rapid ecosystem change and polar bear conservation. Conserv. Lett. 6, 368–375 (2013).
Google Scholar 

4.
Liu, F. et al. Human–wildlife conflicts influence attitudes but not necessarily behaviors: factors driving the poaching of bears in China. Biol. Conserv. 144, 538–547 (2011).
Article  Google Scholar 

5.
McKenney, E. A., Koelle, K., Dunn, R. R. & Yoder, A. D. The ecosystem services of animal microbiomes. Mol. Ecol. 27, 2164–2172 (2018).
CAS  PubMed  Article  Google Scholar 

6.
Nicholson, J. K. et al. Host-gut microbiota metabolic interactions. Science 336, 1262–1267 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

7.
Hill, M. J. Intestinal flora and endogenous vitamin synthesis. Eur. J. Cancer Prev. Off. J. Eur. Cancer Prev. Organ. ECP 6, S43–S45 (1997).
Article  Google Scholar 

8.
Hooper, L. V., Littman, D. R. & Macpherson, A. J. Interactions between the microbiota and the immune system. Science 336(6086), 1268–1273 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

9.
Turnbaugh, P. J. et al. A core gut microbiome in obese and lean twins. Nature 457, 480–484 (2009).
ADS  CAS  PubMed  Article  Google Scholar 

10.
Hauffe, H. C. & Barelli, C. Conserve the germs: the gut microbiota and adaptive potential. Conserv. Genet. 20, 19–27 (2019).
Article  Google Scholar 

11.
Dominianni, C. et al. Sex, body mass index, and dietary fiber intake influence the human gut microbiome. PLoS ONE 10, e0124599 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

12.
McKenney, E. A., Rodrigo, A. & Yoder, A. D. Patterns of gut bacterial colonization in three primate species. PLoS ONE 10, e0124618 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

13.
Barelli, C. et al. Habitat fragmentation is associated to gut microbiota diversity of an endangered primate: implications for conservation. Sci. Rep. 5, 14862 (2015).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

14.
Phillips, C. D. et al. Microbiome analysis among bats describes influences of host phylogeny, life history, physiology and geography: microbiome analysis among bats. Mol. Ecol. 21, 2617–2627 (2012).
PubMed  Article  Google Scholar 

15.
Hooper, L. V., Midtvedt, T. & Gordon, J. I. How host-microbial interactions shape the nutrient environment of the mammalian intestine. Annu. Rev. Nutr. 22, 283–307 (2002).
CAS  PubMed  Article  Google Scholar 

16.
Turnbaugh, P. J. et al. An obesity-associated gut microbiome with increased capacity for energy harvest. Nat. Lond. 444, 1027–1031 (2006).
ADS  Article  Google Scholar 

17.
Clayton, J. B. et al. Captivity humanizes the primate microbiome. Proc. Natl. Acad. Sci. 113, 10376–10381 (2016).
CAS  PubMed  Article  Google Scholar 

18.
Cheng, Y. et al. The Tasmanian devil microbiome: implications for conservation and management. Microbiome 3, 76 (2015).
PubMed  PubMed Central  Article  Google Scholar 

19.
McKenzie, V. J. et al. The effects of captivity on the mammalian gut microbiome. Integr. Comp. Biol. 57, 690–704 (2017).
PubMed  PubMed Central  Article  Google Scholar 

20.
Borgström, B., Dahlqvist, A., Lundh, G. & Sjövall, J. Studies of intestinal digestion and absorption in the human1. J. Clin. Invest. 36, 1521–1536 (1957).
PubMed  PubMed Central  Article  Google Scholar 

21.
Thomson, A. B. R. et al. Normal physiology, part 1. Dig. Dis. Sci. 48, 19 (2003).
Google Scholar 

22.
Amato, K. R. Co-evolution in context: the importance of studying gut microbiomes in wild animals. Microbiome Sci. Med. 1, 10–29 (2013).
Article  Google Scholar 

23.
Stevens, C. E. & Hume, I. D. Comparative Physiology of the Vertebrate Digestive System (Cambridge University Press, Cambridge, 1995).
Google Scholar 

24.
Lafferty, D. J. R., Belant, J. L. & Phillips, D. L. Testing the niche variation hypothesis with a measure of body condition. Oikos 124, 732–740 (2015).
Article  Google Scholar 

25.
Baruch-Mordo, S. et al. Stochasticity in natural forage production affects use of urban areas by black bears: implications to management of human-bear conflicts. PLoS ONE 9, e85122 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

26.
Ayres, L. A., Chow, L. S. & Graber, D. M. Black bear activity patterns and human induced modifications in sequoia national park. Bears Biol. Manag. 6, 151–154 (1986).
Google Scholar 

27.
Enders, M. S. & Vander Wall, S. B. Black bears Ursus americanus are effective seed dispersers, with a little help from their friends. Oikos 121, 589–596 (2012).
Article  Google Scholar 

28.
Pritchard, G. T. & Robbins, C. T. Digestive and metabolic efficiencies of grizzly and black bears. Can. J. Zool. 68, 1645–1651 (1990).
Article  Google Scholar 

29.
Nelson, R. A. et al. Behavior, biochemistry, and hibernation in black, grizzly, and polar bears. Bears Biol. Manag. 5, 284–290 (1983).
Google Scholar 

30.
Brody, A. J. & Pelton, M. R. Seasonal changes in digestion in black bears. Can. J. Zool. 66, 1482–1484 (1988).
Article  Google Scholar 

31.
Hellgren, E. C. Ecology and Physiology of a Black Bear (Ursus americanus) Population in the Great Dismal Swamp and Reproduction Physiology in the Captive Female Black Bear (Virginia Polytechnic Institute and State University, Blacksburg, 1988).
Google Scholar 

32.
Fowler, N. L., Belant, J. L., Wang, G. & Leopold, B. D. Ecological plasticity of denning chronology by American black bears and brown bears. Glob. Ecol. Conserv. 20, e00750 (2019).
Article  Google Scholar 

33.
Samson, C. & Huot, J. Reproductive biology of female black bears in relation to body mass in early winter. J. Mammal. 76, 68–77 (1995).
Article  Google Scholar 

34.
Garshelis, D. L., Scheick, B. K., Doan-Crider, D. L., Beecham, J. J. & Obbard, M. E. Ursus americanus. The IUCN Red List of Threatened Species 2016: e.T41687A114251609 (2016).

35.
Sundin, O. H. et al. The human jejunum has an endogenous microbiota that differs from those in the oral cavity and colon. BMC Microbiol. 17, 160 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

36.
Hayashi, H. Molecular analysis of jejunal, ileal, caecal and recto-sigmoidal human colonic microbiota using 16S rRNA gene libraries and terminal restriction fragment length polymorphism. J. Med. Microbiol. 54, 1093–1101 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

37.
Xiao, Y. et al. Comparative biogeography of the gut microbiome between Jinhua and Landrace pigs. Sci. Rep. 8, 1–10 (2018).
Article  CAS  Google Scholar 

38.
Xue, Z. et al. The bamboo-eating giant panda harbors a carnivore-like gut microbiota, with excessive seasonal variations. MBio 6, e00022-e115 (2015).
CAS  PubMed  PubMed Central  Google Scholar 

39.
Rojas, C. A., Holekamp, K. E., Winters, A. D. & Theis, K. R. Body-site specific microbiota reflect sex and age-class among wild spotted hyenas. FEMS Microbiol. Ecol. 96(2), fiaa007 (2020).
PubMed  Article  Google Scholar 

40.
Sommer, F. et al. The gut microbiota modulates energy metabolism in the hibernating brown bear Ursus arctos. Cell Rep. 14, 1655–1661 (2016).
CAS  PubMed  Article  Google Scholar 

41.
Schwab, C. & Gänzle, M. Comparative analysis of fecal microbiota and intestinal microbial metabolic activity in captive polar bears. Can. J. Microbiol. 57, 177–185 (2011).
CAS  PubMed  Article  Google Scholar 

42.
Zhu, L., Wu, Q., Dai, J., Zhang, S. & Wei, F. Evidence of cellulose metabolism by the giant panda gut microbiome. Proc. Natl. Acad. Sci. 108, 17714–17719 (2011).
ADS  CAS  PubMed  Article  Google Scholar 

43.
Borbón-García, A., Reyes, A., Vives-Flórez, M. & Caballero, S. Captivity shapes the gut microbiota of Andean bears: insights into health surveillance. Front. Microbiol. 8, 13–16 (2017).
Article  Google Scholar 

44.
Song, C. et al. Comparative analysis of the gut microbiota of black bears in China using high-throughput sequencing. Mol. Genet. Genomics 292, 407–414 (2017).
CAS  PubMed  Article  Google Scholar 

45.
McKenney, E. A., Maslanka, M., Rodrigo, A. & Yoder, A. D. Bamboo specialists from two mammalian orders (primates, carnivora) share a high number of low-abundance gut microbes. Microb. Ecol. 76, 272–284 (2018).
PubMed  Article  Google Scholar 

46.
Bollinger, R. R., Barbas, A. S., Bush, E. L., Lin, S. S. & Parker, W. Biofilms in the large bowel suggest an apparent function of the human vermiform appendix. J. Theor. Biol. 249, 826–831 (2007).
CAS  Article  Google Scholar 

47.
Smith, H. F. et al. Comparative anatomy and phylogenetic distribution of the mammalian cecal appendix. J. Evol. Biol. 22, 1984–1999 (2009).
CAS  PubMed  Article  Google Scholar 

48.
Sanders, N. L., Bollinger, R. R., Lee, R., Thomas, S. & Parker, W. Appendectomy and clostridium difficile colitis: relationships revealed by clinical observations and immunology. World J. Gastroenterol. WJG 19, 5607–5614 (2013).
PubMed  Article  Google Scholar 

49.
Merchant, R. et al. Association between appendectomy and clostridium difficile infection. J. Clin. Med. Res. 4, 17–19 (2012).
PubMed  PubMed Central  Google Scholar 

50.
Greene, L. K. & McKenney, E. A. The inside tract: the appendicular, cecal, and colonic microbiome of captive aye-ayes. Am. J. Phys. Anthropol. 166, 960–967 (2018).
PubMed  Article  Google Scholar 

51.
Tilg, H. & Kaser, A. Gut microbiome, obesity, and metabolic dysfunction. J. Clin. Invest. 121, 2126–2132 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

52.
Ley, R. E. Obesity and the human microbiome. Curr. Opin. Gastroenterol. 26, 5–11 (2010).
PubMed  Article  Google Scholar 

53.
Cani, P. D. et al. Changes in gut microbiota control metabolic endotoxemia-induced inflammation in high-fat diet-induced obesity and diabetes in mice. Diabetes 57, 1470–1481 (2008).
CAS  PubMed  Article  Google Scholar 

54.
Scher, J. U. et al. Decreased bacterial diversity characterizes the altered gut microbiota in patients with psoriatic arthritis, resembling dysbiosis in inflammatory bowel disease. Arthritis Rheumatol. 67, 128–139 (2015).
CAS  PubMed  PubMed Central  Article  Google Scholar 

55.
Dietz, R. et al. Ursidibacter maritimus gen. nov., sp. nov. and Ursidibacter arcticus sp. nov., two new members of the family Pasteurellaceae isolated from the oral cavity of bears. Int. J. Syst. Evol. Microbiol. 65, 3683–3689 (2015).
PubMed  Article  CAS  Google Scholar 

56.
Christensen, H. & Bisgaard, M. Taxonomy and biodiversity of members of Pasteurellaceae. In Pasteurellaceae: Biology, Genomics and Molecular Aspects (eds Kuhnert, P. & Christensen, H.) 1–26 (Caister Academic Press, Norfolk, 2008).
Google Scholar 

57.
Ma, J. et al. High-fat maternal diet during pregnancy persistently alters the offspring microbiome in a primate model. Nat. Commun. 5, 3889 (2014).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

58.
Yasuda, K. et al. Biogeography of the intestinal mucosal and lumenal microbiome in the rhesus macaque. Cell Host Microbe 17, 385–391 (2015).
CAS  PubMed  PubMed Central  Article  Google Scholar 

59.
Carthey, A. J. R., Blumstein, D. T., Gallagher, R. V., Tetu, S. G. & Gillings, M. R. Conserving the holobiont. Funct. Ecol. https://doi.org/10.1111/1365-2435.13504 (2020).
Article  Google Scholar 

60.
Cappa, F., Laut, J., Nov, O., Giustiniano, L. & Porfiri, M. Activating social strategies: face-to-face interaction in technology-mediated citizen science. J. Environ. Manag. 182, 374–384 (2016).
Article  Google Scholar 

61.
Budde, M. et al. Participatory sensing or participatory nonsense? Mitigating the effect of human error on data quality in citizen science. Proc. ACM Interact. Mob. Wearable Ubiquitous Technol. 1, 1–23 (2017).
Article  Google Scholar 

62.
McKenney, E. A., Greene, L. K., Drea, C. M. & Yoder, A. D. Down for the count: cryptosporidium infection depletes the gut microbiome in Coquerel’s sifakas. Microb. Ecol. Health Dis. 28, 1335165 (2017).
PubMed  PubMed Central  Google Scholar 

63.
Caporaso, J. G. et al. Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J. 6, 1621–1624 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

64.
Bolyen, E. et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. https://doi.org/10.1038/s41587-019-0209-9 (2019).
Article  PubMed  PubMed Central  Google Scholar 

65.
Callahan, B. J. et al. DADA2: high-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581 (2016).
CAS  PubMed  PubMed Central  Article  Google Scholar 

66.
Quast, C. et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596 (2012).
PubMed  PubMed Central  Article  CAS  Google Scholar 

67.
Bokulich, N. A. et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome 6, 90 (2018).
PubMed  PubMed Central  Article  Google Scholar 

68.
Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772–780 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

69.
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, 2018).
Google Scholar 

70.
Allaire, J. RStudio: Integrated Development Environment for R 770 (RStudio, Boston, 2012).
Google Scholar 

71.
Anderson, M. J. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46 (2001).
Google Scholar 

72.
Anderson, M. J. Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62, 245–253 (2006).
MathSciNet  PubMed  MATH  Article  Google Scholar 

73.
Segata, N. et al. Metagenomic biomarker discovery and explanation. Genome Biol. 12, R60 (2011).
PubMed  PubMed Central  Article  Google Scholar  More

We consider a two-layer network. One layer is the slave layer, which corresponds to the original network over which one wants to impose the desired dynamics (i.e. a given evolution compatible with the equations of motion). The other layer is the master layer, which is identical to the slave layer, but starts from a different initial condition (i.e. the one generating the specific desired dynamics towards which the state of the slave layer is to be steered), and evolves autonomously. In applications, the master layer may just be an experimental recording or a simulation of the original system—as long as it can be coupled to the slave network its physical nature is irrelevant. Our control method consists then in establishing directed inter-layer links from nodes in the master layer to their counterparts in the slave layer. Once they are established, these links remain in place as more nodes are connected in sequential control steps. At each step the selected node is the one whose pinning causes the most rapid approach towards inter-layer synchronization (i.e. the imposition of the evolution followed by the master layer on the slave layer). While the two layers have to be identical, the nodes (i.e. the dynamical units) and links (the coupling structure connecting the dynamical systems) on each layer can be completely different, as we will see below. This is thus a generalization of the method proposed in Ref.6.
We illustrate our method by applying it to networks of identical chaotic oscillators, and leave the applicability to more challenging real-world systems to the next section. Specifically, we consider networks of (N=50) nodes whose topology is that of a mixed random graph, i.e. containing both bidirectional and unidirectional links. These graphs are realizations of the configuration model22 with the in-degree (k_text {in}) (i.e. the number of links pointing to a given node) and the out-degree (k_text {out}) (i.e. the number of links emanating from a given node) uniformly distributed in ({5,6,ldots ,45}). Each node evolves autonomously in time as a chaotic Rössler oscillator, which we simply denote as ({dot{mathbf{r}}} = mathbf{f}(mathbf{r})), where (mathbf{r} = (x,y,z)^text {T}) and ({dot{x}} = -y – z, {dot{y}} = x + a y, {dot{z}} = b + z (x – c)), with parameters (a=0.2), (b=0.2) and (c=7). Nodes are coupled quadratically via their z variables, a nonlinear coupling form that was previously considered in Ref.23.
Before the first control step is applied (prior to the creation of the first inter-layer connection) both master and slave layers evolve spontaneously as follows

$$begin{aligned} {dot{mathbf{r}}}_i = mathbf{f}(mathbf{r}_i) + sigma _1 sum _{j=1}^N D_{ji} (z_j^2-z_i^2) = mathbf{f}(mathbf{r}_i) + sigma _1 sum _{j=1}^N {mathcal {L}}_{ji} z_j^2. end{aligned}$$
(1)

where (D_{ji} = 1) if there is a directed link from node j to node i, and is zero otherwise (for bidirectional links (D_{ij} = D_{ji})). As we do not consider self-links, the diagonal terms vanish, i.e. (D_{ii} = 0 forall , i), and the in-degree of node i is (k_{text {in},i} = sum _{j} D_{ji}). The graph can thus be alternatively represented by the Laplacian matrix ({mathcal {L}}_{ji} =D_{ji} – k_{text {in},i} delta _{ji}). The vector field (mathbf{f}(mathbf{r}_i)) governs the dynamics of node i, which would evolve autonomously (if uncoupled from its neighbors) simply as ({dot{mathbf{r}}}_i = mathbf{f}(mathbf{r}_i)), and the parameter (sigma _1) is the intra-layer coupling strength
When the control procedure starts, each node i in the master network keeps evolving according to the dynamics in Eq. (1), ({dot{mathbf{r}}}^M_i = mathbf{f}(mathbf{r}^M_i) + sigma _1 sum _{j} {mathcal {L}}_{ji} (z^M_j)^2). In the slave layer dynamics, however, one has to consider an additional term which accounts for the inter-layer coupling from the master layer (without loss of generality, we here take a linear coupling through the y variable). One has

$$begin{aligned} {dot{mathbf{r}}}^S_i = mathbf{f}(mathbf{r}^S_i) + sigma _1 sum _{j} {mathcal {L}}_{ji} (z^S_j)^2 + sigma _2 chi _i (y^M_i – y^S_i). end{aligned}$$
(2)

Here (chi _i) is a binary variable that is one if there is a link coupling node i in the master layer to node i in the slave layer (i.e. if the targeting procedure includes a pinning action from master to slave at node i) and is zero otherwise. The parameter (sigma _2) is the inter-layer coupling strength. We emphasize that the coupling that is linear (in fact, diffusive) is the externally-imposed inter-layer coupling, which does not restrict in any way the form of the (intra-layer) couplings between the nodes of the system under study. Such diffusive inter-layer coupling is chosen as it is the simplest form that makes the inter-layer synchronization manifold into an invariant set of the dynamics (for a detailed mathematical treatment of invariant sets and related concepts, see e.g. Ref.24).
Figure 1

Controlling the dynamics of a mixed network with uniform (k_text {in}) and (k_text {out}) distributions comprising (N=50) nonlinearly-coupled Rössler oscillators with intra-layer coupling (sigma _1 = 0.01) and inter-layer coupling (sigma _2 = 1). (Top) Maximum Lyapunov exponent (lambda _text {max}) (main panel) and synchronization error (inset) as functions of the targeting step. (Bottom) Influence index (k_text {out}/k_text {in}) of the node that is pinned at each targeting step. The curves are averages of 20 different network realizations. A 4th-order Runge-Kutta method with a step of 0.01 time units has been employed for the numerical integration of the systems of (3 N=150) ordinary differential equations corresponding to each layer.

Full size image

The results of applying our method to the network of Rössler chaotic oscillators are shown in Fig. 1. Two observables are employed to characterize the inter-layer synchronization between master and slave as more and more inter-layer links are established in successive targeting steps. One is the maximum Lyapunov exponent, (lambda _text {max}), computed from the dynamics of the slave network linearized around that of the master network as in Ref.6. For a review of the theory and numerical computation of Lyapunov spectra, see e.g. Refs.25,26. The other observable is the synchronization error, which is the time average of the Euclidean distance in phase space ({mathbb {R}}^{N m}) (N is the number of nodes, m is the phase space dimensionality of the node dynamics—in our case (N=50), (m = 3)) between the full state of the master layer and that of the slave layer, (lim _{Trightarrow infty } frac{1}{T} int _0^T sqrt{ sum _{i=1}^{N} (x^M_i(t) – x^S_i(t))^2 } dt). In practice, T is finite, but orders of magnitude larger than the characteristic timescales of oscillation (thus the numerical convergence to an asymptotic value is guaranteed). In the top panel of Fig. 1, we show the maximum Lyapunov exponent (lambda _text {max}) as a function of the targeting step, which is seen to progressively decrease as more and more nodes are pinned. Analogous results in terms of the synchronization error are reported in the inset, which shows how the synchronization error becomes zero when (lambda _text {max}) becomes negative.
The maximum Lyapunov exponent is also used to identify the node to be targeted at each step: of all the nodes that remain unconnected to their counterparts in the other layer, the one that, when a master-slave connection is established, leads to the largest decrease in (lambda _text {max}) is targeted next. An exploration of possible correlations between the resulting targeting sequence (i.e. the ordered list of nodes that are targeted at successive steps) and local topological properties yields a remarkable correspondence between the targeting sequence position and the ranking of nodes in terms of their influence index (k_text {out}/k_text {in}), as shown in the lower panel of Fig. 1. This index is large when a node has a privileged position for influencing other nodes, while receiving very little influence from the rest of the network. No such correlations are observed for connectivity indices that are insensitive to the directionality of connections, such as ((k_text {in}+k_text {out})/2), while correlations only based on (k_text {out}) or (k_text {in}) give considerably poorer results that those shown in the figure. Other measures of connection directionality that we have inspected, such as ((k_text {out} – k_text {in})/(k_text {out} + k_text {in})), show weaker correlations with the targeting sequence than the influence index does. While these results are based on networks with uniform distributions of (k_text {in}) and (k_text {out}), which have been chosen precisely because a large variety of possible degree values is desirable, a strong correlation between the influence-index ranking and the targeting ranking is also observed for Barabási-Albert scale-free networks27 and Erdös-Rényi random graph28 topologies, as shown in Sect. A of the Supplementary Information.
This correlation is most clearly seen for small values of the intra-layer coupling strength, such as the value (sigma _1 = 0.01) considered in Fig. 1. For larger values of (sigma _1), which make inter-layer synchronization possible with a very small number of steps, the correlation is less strong, while no obvious correlation between the targeting sequence and local topological properties are found for very large (sigma _1), see again Sect. A of the Supplementary Information. This might be related to the enhanced contribution of next-nearest neighbors and other relative distant nodes as the coupling strength is increased. Despite its limited range of validity, this correlation is nontheless remarkable, as it is very robust, and quite different from the situation observed in undirected networks, where the topological observable correlating with the targeting sequence is the degree6. On the other hand, there is an intriguing parallel between the correlation reported in the lower panel of Fig. 1 and the fact that, in undirected networks, nodes with a higher dynamic vulnerability are those with less influence from the rest of the network, followed by those that have the strongest ability to influence the rest29. In fact, both aspects of a node position are combined in the influence index in the case of directed or mixed networks. More

A waterborne, abiotic, and other indirectly transmitted (WAIT) model for the dynamics of emergent viral outbreaks
Several models have been engineered to explore aspects of COVID-19 dynamics. For example, models have been used to investigate the role of social distancing20,42, social mixing43, the importance of undocumented infections44, the role of mobility in the early spread of disease in China45, and the potential for contact tracing as a solution46. Only a few notable models of SARS-CoV-2 transmission incorporate features of indirect or environmental transmission40,46 and none consider the dynamical properties of viral free-living survival in the environment. Such a model structure would provide an avenue towards exploring how variation in free-living survival influences disease outbreaks. Indirect transmission includes those routes where pathogen is spread through means other than from person to person, and includes transmission through environmental reservoirs. Environmental transmission models are aplenty in the literature and serve as a theoretical foundation for exploring similar concepts in newer, emerging viruses1,2,3,4,5,6,7,8,9,10.
Here, we parameterize and validate an SEIR-W model: Susceptible (S), Exposed (E), Infectious (I), Recovered (R), and WAIT (W) model. Here W represents the environmental component of the transmission cycle during the early stage of the SARS CoV-2 pandemic. This component introduces more opportunities for infection, and complex dynamics resulting from viral persistence in the environment. In this framework, both indirect and direct transmission occur via mass-action, “random” encounters.
This model is derived from a previously developed framework called “WAIT”—which stands for Waterborne, Abiotic, and other Indirectly Transmitted—that incorporates an environmental reservoir where a pathogen remains in the environment and “waits” for hosts to interact with it11,12. The supplementary information contains a much more rigorous discussion of the modeling details. In the main text, we provide select details.
Building the SEIR-W model framework for SARS-CoV-2
Here W represents the environmental component of the early stage of the SARS CoV-2 pandemic (Fig. S1). This environmental compartment refers to reservoirs that people may have contact with on a daily basis, such as doorknobs, appliances, and non-circulating air indoors. The W compartment of our model represents the fraction of these environmental reservoirs that house some sufficiently transmissible amount of infectious virus. We emphasize that the W compartment is meant to only represent reservoirs that are common sites for interaction with people. Thus, inclusion of the W compartment allows us to investigate the degree to which the environment is infectious at any given point, and its impact on the transmission dynamics of SARS CoV-2.
Model parameters are described in detail in Table 1. The system of equations in the proposed mathematical model corresponding to these dynamics are defined in Eqs. (1)–(6):
Table 1 Model population definitions and initial values denoted with subscript 0 for each state variable.
Full size table

$$frac{dS}{dt}=mu (N-S)-left(frac{{beta }_{A}{I}_{A}+{beta }_{S}{I}_{S}}{N}+{beta }_{W}Wright)S$$
(1)

$$frac{dE}{dt}=left(frac{{beta }_{A}{I}_{A}+{beta }_{S}{I}_{S}}{N}+{beta }_{W}Wright)S-(epsilon +mu )E$$
(2)

$$frac{d{I}_{A}}{dt}=epsilon E-(omega +mu ){I}_{A}$$
(3)

$$frac{d{I}_{S}}{dt}=(1-p)omega {I}_{A}-(nu +{mu }_{S}){I}_{S}$$
(4)

$$frac{dR}{dt}=pomega {I}_{A}+nu {I}_{S}-mu R$$
(5)

$$frac{dW}{dt}=left(frac{{sigma }_{A}{I}_{A}+{sigma }_{S}{I}_{S}}{N}right)left(1-Wright)-kW$$
(6)

Infection trajectories
In addition to including a compartment for the environment (W), our model also deviates from traditional SEIR form by splitting the infectious compartment into an IA-compartment (A for asymptomatic), and an IS-compartment (S for symptomatic). As we discuss below, including asymptomatic (or sub-clinical) transmission is both essential for understanding how environmental—as opposed to simply unobserved or hidden—transmission affects the ecological dynamics of pathogens and also for analyzing SARS-CoV-2. The former represents an initial infectious stage (following the non-infectious, exposed stage), from which individuals will either move on to recovery directly (representing those individuals who experienced mild to no symptoms) or move on to the IS-compartment (representing those with a more severe response). Finally, individuals in the IS-compartment will either move on to recovery or death due to the infection. This splitting of the traditional infectious compartment is motivated by mounting evidence of asymptomatic transmission of SARS CoV-244,47,48,49,50. Thus, we consider two trajectories for the course of the disease, similar to those employed in prior studies42: (1) E → IA → R and (2) E → IA → IS → R (or death). More precisely, once in the E state, an individual will transition to the infectious state IA, at a per-person rate of ε. A proportion p will move from IA to the recovered state R (at a rate of p ⍵). A proportion (1—p) of individuals in the IA state will develop more severe systems and transition to Is (at a rate of (1—p) ⍵). Individuals in the Is state recover at a per-person rate of ν or die at a per-person rate μS. In each state, normal mortality of the individual occurs at the per-person rate μ and newly susceptible (S) individuals enter the population at a rate μN. The important differences between these two trajectories are in how likely an individual is to move down one path or another, how infectious individuals are (both for people and for the environment), how long individuals spend in each trajectory, and how likely death is along each trajectory.
Clarification on the interactions between hosts and reservoirs
The model couples the environment and people in two ways: (1) people can deposit the infectious virus onto environmental reservoirs (e.g. physical surfaces, and in the case of aerosols, the ambient air) and (2) people can become infected by interacting with these reservoirs. While most of our study is focused on physical surfaces, we also include data and analysis of SARS-CoV-2 survival in aerosols. While aerosols likely play a more significant role in person-to-person transmission, they also facilitate an indirect means of transmitting. For example, because SARS-CoV-2 can remain suspended in the air, other individuals can become infected without ever having to be in especially close physical proximity to the aerosol emitter (only requires that they interact with the same stagnant air, containing infectious aerosol particles)51. That is, a hypothetical infectious person A may produce aerosols, leave a setting, and those aerosols may infect a susceptible individual B who was never in close proximity to person A. In the transmission event between person A and person B, aerosol transmission functions in a similar fashion to surface transmission, where aerosols may be exchanged in the same room where infected individuals were, rather than exchanging infectious particles on a surface.
In our model, indirect infection via aerosols is encoded into the terms associated with the W component, just as the different physical surfaces are. Alternatively, aerosol transmission that leads to direct infection between hosts is encoded in the terms associated with direct infection between susceptible individuals and those infected (see section entitled Infection Trajectories).
Environmental reservoirs infect people through the βW term (Eqs. 1 and 2), a proxy for a standard transmission coefficient, corresponding specifically to the probability of successful infectious transmission from the environment reservoir to a susceptible individual (the full rate term being βWW·S). Hence, the βW factor is defined as the fraction of people who interact with the environment daily, per fraction of the environment, times the probability of transmitting infection from environmental reservoir to people. The factor βWW (where W is the fraction of environmental reservoirs infected) represents the daily fraction of people that will interact with the infected portion of the environment and become infected themselves. The full term βWW·S is thus the total number of infections caused by the environment per day.
In an analogous manner, we model the spread of infection to the environment with the two terms σAIA·(1—W) / N and σSIS·(1—W) / N representing deposition of infection to the environment by asymptomatic individuals, in the former, and symptomatic individuals, in the latter. In this case, σA (and analogously for σS) gives the fraction of surfaces/reservoirs that interact with people at least once per day, times the probability that a person (depending on whether they are in the IA or the IS compartment) will deposit an infectious viral load to the reservoir. Thus, σAIA/ N and σSIS/ N (where N is the total population of people) represent the daily fraction of the environment that interacts with asymptomatic and symptomatic individuals, respectively. Lastly, the additional factor of (1—W) gives the fraction of reservoirs in the environment that have the potential for becoming infected, and so σAIA·(1—W) / N (and analogously for IS) gives the fraction of the environment that becomes infected by people each day. We use W to represent a fraction of the environment, although one could also have multiplied the W equation by a value representing the total number of reservoirs in the environment (expected to remain constant throughout the course of the epidemic, assuming no intervention strategies).
Parameter values estimation
Table 1 displays information on the population definitions and initial values in the model. Tables 2 and 3 contain the fixed and estimated values and their sources (respectively). The model’s estimated parameters are based on model fits to 17 countries with the highest cumulative COVID-19 cases (of the 181 total countries affected) as of 03/30/2020, who have endured outbreaks that had developed for at least 30 days following the first day with ≥ 10 cumulative infected cases within each country14 (See supplementary information Tables S1–S3). In addition, we compare country fits of the SEIR-W model to fits with a standard SEIR model. Lastly, we compare how various iterations of these mathematical models compare to one another with regards to the general model dynamics. For additional details, see the supplementary information.
Table 2 Fixed parameter values estimated based on available published literature.
Full size table

Table 3 Estimated parameter values, averaged across countries.
Full size table

Basic reproductive ratios (({mathcal{R}}_{0}))
We can express the ({mathcal{R}}_{0}) (Eq. 7) in a form that makes explicit the contributions from the environment and from person-to-person interactions. In this way, the full ({mathcal{R}}_{0}) is observed to comprise two ({mathcal{R}}_{0}) sub-components: one the number of secondary infections caused by a single infected person through person-to-person contact alone (Rp) and the other is the number of secondary infections caused by exchanging infection with the environment (Re).

$${R}_{0}=frac{{R}_{p} + sqrt{{R}_{p}^{2} + 4 {R}_{e}^{2}}}{2}$$
(7)

where Rp and Re are defined in Eqs. (8a) and (8b)

$${R}_{p}=frac{varepsilon ({beta }_{A} ({mu }_{S }+ nu ) + {beta }_{S}(1 – p) omega )}{(mu + varepsilon )(mu + omega )({mu }_{S} + nu )}$$
(8a)

$${R}_{e}^{2}=frac{varepsilon { beta }_{W} ({sigma }_{A} ({mu }_{S }+ nu ) + {sigma }_{S}(1 – p) omega )}{k (mu + varepsilon )(mu + omega )({mu }_{S} + nu )}$$
(8b)

Note that when Rp = 0, the ({mathcal{R}}_{0}) reduces to Re and when Re = 0, the ({mathcal{R}}_{0}) reduces to Rp. Thus, when person-to-person transmission is set to zero, the ({mathcal{R}}_{0}) consists only of terms associated with transmission from the environment, and when transmission from the environment is set to zero, the ({mathcal{R}}_{0}) consists only of infection directly between people. When both routes of transmission are turned on, the two ({mathcal{R}}_{0})-components combine in the manner in Eq. (7).
While Re represents the component of the ({mathcal{R}}_{0}) formula associated with infection from the environment, the square of this quantity Re2 represents the expected number of people who become infected in the two-step infection process: people → environment → people, representing the flow of infection from people to the environment, and then from the environment to people. Thus, while Rp gives the expected number of people infected by a single infected person when the environmental transmission is turned off, Re2 gives the expected number of people infected by a single infected person by way of the environmental route exclusively, with no direct person-to-person transmission. Also note that Re2/(Re2 + Rp) can be used to measure the extent of transmission that is mediated by the environment exclusively. This proportion can be used as a proxy for how important environmental transmission is in a given setting. Elaboration on formulas 8a–b—and associated derivation-discussions—appear in the supplementary information. More

Case 1: spectacled eiders, from Wang et al.42
The experiment
This case is based on captive feeding trials conducted on 8 adult spectacled eiders, Somateria fischeri, which were maintained on an initial diet containing 1% Atlantic surf clam, 3% Antarctic krill, 88% Mazuri sea duck formula, 4% blue mussel and 4% Atlantic silverside, for 69 days prior to the start of the feeding experiment. After this, on day 0, a biopsy sample of the synsacral adipose tissue was obtained from each eider. With the FA data of the adipose tissue the authors calculated CCs. Feeding trials started on Day 0, and spectacled eiders were switched to diet A, consisting of 56% krill and 44% Mazuri sea duck formula for 21 days. On Day 21, eiders were biopsied again and switched to diet B consisting of 48% Mazuri formula and 52% silverside. On Day 50, a final biopsy sample was collected (Fig. 2A). FA turnover was considered near complete by 69 days.
Figure 2

Spectacled eiders case. (A) Feeding experiment: spectacled eiders (n = 8) spent 69 days on the initial diet described in the figure; after this, on day 0 eiders were biopsied and switched to diet A. On day 21 they were biopsied again and switched to diet B. After 29 days on diet B, eiders were biopsied on day 50. (B) Plots for diet estimations of spectacled eiders fed different combined diets. The true diet is indicated in each plot by the blue asterisks. CCs were calculated using FA data of day 0. Images by Alicia Guerrero.

Full size image

The model
This analysis is based on FA data from day 0, 21 and 50. We used the CCs calculated in this study after eiders were maintained on the same initial diet for 69 days. All sources were significantly different from each other (PERMANOVA, F4 = 9936.9, P = 0.001). ‘Day’ was set as fixed factor in the model.
Diet predictions
Based on FA data of day 0, MixSIAR estimated a contribution of 0.1% clam, 4% krill, 88% Mazuri, 6% mussel, and 1% silverside (Fig. 2B). After the first shift of diet, MixSIAR estimations changed to 28% krill and 62% Mazuri for FA data obtained on day 21. The final biopsy sample on day 50 produced estimations of 62% Mazuri, and 27% silverside.
Case 2: Steller’s eiders, from Wang et al.42
The experiment
This case corresponds to a feeding trial conducted simultaneously to the previous case by Wang et al.42, although the diet of Steller’s eiders, Polysticta stelleri, differed slightly. For 69 days prior to the start of the feeding trial, 8 adult Steller’s eiders were maintained on an initial diet containing 1% clam, 1% Antarctic krill, 88% Mazuri sea duck formula, 7% mussel, and 3% silverside. CCs were calculated after a biopsy was extracted to each eider on day 0. At the start of the feeding experiment, on day 0, Steller’s eiders were switched to diet A, containing 66% krill and 34% Mazuri formula. Then, on day 21, they were switched to diet B, consisting of 34% Mazuri formula and 66% silverside. Biopsy samples were collected on days 0, 21 and 50 (Fig. 3A). FA turnover was considered near complete by 69 days.
Figure 3

Steller’s eiders case. (A) Feeding experiment: eiders (n = 8) were maintained on the initial diet described in the figure for 69 days after which biopsy samples were collected (day 0) and eiders were switched to diet A. After 21 days, eiders were biopsied again and switched to diet B. On day 50, after 29 days on diet B, eiders were biopsied one last time. (B) Plots of diet estimation for Steller’s eiders fed different combined diets. Diet estimates are based on biopsy samples collected on days 0, 21 and 50. The true diet is indicated in each plot by red asterisks. CCs were calculated using FA data of day 0. Images by Alicia Guerrero.

Full size image

The model
We used the CCs calculated in the same study after Steller’s eiders were maintained on the same diet for 69 days. All sources were significantly different from each other (PERMANOVA, F4 = 10,079, P = 0.001). ‘Day’ was set as fixed factor in the model.
Diet predictions
Based on FA data of Day 0, MixSIAR estimated a contribution of 0.2% clam, 2% krill, 90% Mazuri, 4% mussel, and 3% silverside (Fig. 3B). After the first diet switch, MixSIAR estimations changed to 46% krill and 27% Mazuri. The final biopsy sample on day 50 produced estimations of 9% krill, 28% Mazuri, and 46% silverside.
Case 3: Atlantic salmon, from Budge et al.43
The experiment
For 22 weeks, tank-reared juvenile Atlantic salmon, Salmo salar (n = 132), were fed one of four different formulated feeds based on two marine oils: 100% krill oil, 100% herring oil, a mixture of 70:30 herring to krill oil, or a mixture of 30:70 herring to krill oil. Muscle samples were analysed for FAs after the 22-week experiment, which allowed the calculation of CCs. In this experiment, two sets of CCs were calculated: one derived from salmon fed the diet based on 100% herring oil and another from salmon fed the diet based on 100% krill oil (Fig. 4A). Unlike the previous two cases, where CCs were derived from consumers eating a mixed diet, here CCs were obtained from consumers feeding on a single type of source: either herring or krill oil. This allowed us to evaluate whether the source used to calculate the CCs affected dietary predictions. Additionally, we calculated a combined CC (an average value between CCs derived from krill and herring diets) and run a separate model. FA turnover was considered complete after 22 weeks on the same diet.
Figure 4

Atlantic salmon case. (A) Feeding experiment: Atlantic salmon fed formulated feeds based on either solely herring (n = 36) or krill oil (n = 28), or in proportions of 70:30 or 30:70 herring to krill oil (n = 34 each), for 22 weeks. (B) Plots of diet proportions estimated using MixSIAR for models using CCs derived from salmon fed a herring-oil diet (HO-CC), a krill-oil diet (KO-CC) or CCs averaged from these two treatments (Combined-CC). The true diet is indicated in each plot by black asterisks. Salmon image designed by Creazilla (https://creazilla.com).

Full size image

The model
We ran three independent models to estimate the diet of salmon fed each of the four diets: one using the set of CCs derived from herring oil, another using the CCs derived from krill oil, and one using the combined CCs. For each model, we used krill oil, herring oil and initial diet (commercial feed given to salmon prior to the experiment) as sources, and set ‘diet group’ as fixed factor. Significantly different FA compositions were found for the three sets of sources: those multiplied by herring oil CCs (PERMANOVA, F2 = 4008.7, P = 0.003), those multiplied by krill oil CCs (PERMANOVA, F2 = 3354.7, P = 0.002), and those multiplied by the combined CCs (PERMANOVA, F2 = 3677.2, P = 0.003).
Diet predictions
When we used CCs derived from salmon fed on a diet supplemented with herring oil only, MixSIAR correctly estimated the contribution of herring oil in the consumers (salmon) diet (98%). However, the contribution of herring was slightly overestimated where the consumers had been fed a mixture of herring and krill oil (Fig. 4B), and where the salmon’s diet was based on krill oil the contribution of krill oil was slightly underestimated (89%).
We found the opposite trend when we used CCs derived from salmon that had been fed a diet where krill oil had been the lipid source. Here, MixSIAR correctly estimated the contribution of krill oil in the salmon’s diet when they had been fed a diet based on krill oil (98%) or a mixture of 70:30 herring to krill oil (33% krill) or 30:70 herring to krill oil (70% krill); however, when herring oil had been the only dietary source this dietary contribution was underestimated (81%) (Fig. 4B).
Our dietary estimates were less biased when we used CC values that had been derived from the average between the herring- and krill-oil treatments (Combined-CCs). For example, we estimated herring contribution to be 90% of the diet when the actual diet was supplemented only with herring oil, and an estimate of 95% krill contribution when the actual diet was supplemented with krill oil only, and when the actual diet was a combination of herring (70%) and krill (30%), the estimations were 71% and 27%, and where the actual contribution of herring was 30% and 70% of krill, the estimated diets were 34% herring and 65% krill (Fig. 4B).
Case 4: tufted puffin nestlings, from Williams et al. 44
The experiment
Tufted puffin, Fratercula cirrhata, nestlings (n = 6) underwent an experimental feeding trial in their own burrows. Chicks were fed by their parents for approximately 10 days since hatching. When they were estimated to be 10-days old, the access to the burrows was blocked, so adults could not continue feeding their chicks. Through another access hole excavated by the researchers, chicks began being fed Pacific herring once a day, for 27 days. To infer the diet of free-living puffin nestlings during the first 10 days after hatching, wire screens were placed at burrow entrances to collect whole fish dropped by the parents. The species identified, in descending order (by mass), were Pacific sandlance, capelin, Pacific sandfish, salmonid and Pacific cod. An adipose tissue sample was collected on days 10 (start of the feeding trial), 19, 28 and 37. On day 37, assuming complete FA turnover after 27 days on a single prey diet (herring), the researchers calculated CCs (Fig. 5A). The data used to run this model included day 10, which represents the unknown diet provided by the parents, days 19, 28, and 37 which represent the herring diet at different extents. FA turnover was considered “close to, but not entirely complete” after 27 days44.
Figure 5

source even though it was not part of the nestlings’ diet. The red asterisks in each plot represent the potential diet fed by the parents (and used as priors in the first model). (D) Plots of diet estimations for tufted puffins fed herring, based on their FA profiles of days 19, 28 and 37. The true diet is indicated in each plot by red asterisks. CCs were calculated from tufted puffins fed herring, using FAs from day 37. Images by Alicia Guerrero.

Tufted puffins case. (A) Nestlings (n = 6) were fed by their parents for approximately 10 days since hatching. After this, they were fed herring for another 27 days as the entrance to their burrows was blocked and parents could not feed their chicks. (B) Non-metric multidimensional scaling plots for FAs obtained from chicks at different stages of the experiment and their sources. When FAs of sources were multiplied by their respective CCs, source (herring) and consumer (chicks, day 37) overlap in the plot. (C) Plots of the three models run to estimate the diet of nestlings on day 10. From left to right: Model using informative priors based on meals brought by the parents after the burrows were blocked; the same model without informative priors; and a third model without informative priors but including herring as

Full size image

The model
We used the CCs derived from these chicks feeding on Pacific herring for 27 days. For biopsy samples collected on day 10, we conducted three separate analyses: the first model excluded herring as potential prey, and incorporated informative priors based on the amount of different fish (% by mass) dropped by the parents at the burrow entrances; the second model was exactly the same but without prior information; and the third model was run without priors, and included herring as potential prey in order to determine whether this prey could be identified as absent.
For days 19, 28 and 37, we ran another analysis and included herring as source and no informative priors. ‘Day’ was set as a fixed factor in this model. All sources had significantly different FA compositions (PERMANOVA, F5 = 430.9, P = 0.001).
Diet predictions
In this example we included a non-metric dimensional scaling plot (Fig. 5B) to evaluate the effect of applying CCs to sources. Day 0 biopsies show greater variation than those of successive days, in both plots. When using original sources and consumer FA values, chicks are segregated from all the sources, although the similarity of the FAs increases toward the FA of herring as days pass, but they do not match. When CCs were applied to sources, the FA values of herring and chicks from day 37 overlap, indicating that they have the same FA compositions.
For day 0 (Fig. 5C), the estimated diet contributions were similar to meals brought by the parents when the model included informative priors, where the main dietary sources were sandlance (72%) and capelin (15%). Whereas estimates from the model without informative priors misrepresented the diet, as capelin was wrongly identified as the main dietary source (61%), cod the second most important prey (26%), and sandlance was incorrectly estimated to be only 6% of the diet. The third model including herring again identified capelin and cod as the main contributors (56% and 23%, respectively) whereas herring was identified as the least important prey, with 0.9% of contribution.
For day 19 (Fig. 5D), herring was identified as the main source, with 60% of contribution, followed by capelin and sandlance although with greater variation. From day 19 to 37, the contribution of herring increases from 60 to 97%, respectively.
Case 5. Harp seals, from Kirsch et al.45
The experiment
This study evaluated the effect of a low-fat diet on blubber FAs of harp seals, Pagophilus groenlandicus. Only for this experiment, the fat content of the different sources was available. Juvenile harp seals (n = 5) had been maintained on a diet of Atlantic herring (≥ 9% fat) for approximately 1 year prior to the feeding trial. On day 0, a full-depth blubber sample was collected from the posterior flank of each animal. For 30 days, seals were kept on a diet consisting solely of Atlantic pollock (1.7% fat). Blubber biopsies were taken again on days 14 and 30 (Fig. 6A). FA turnover was not considered complete after 30 days on the same diet, and authors suggest that a longer period on the diet, or higher intakes of fat, would be needed to accomplish it.
Figure 6

Diet estimates for juvenile harp seals fed a low-fat prey. (A) Feeding experiment: For a year prior to the feeding experiment, harp seals (n = 5) had been eating only Atlantic herring, a prey with a high-fat content. During the feeding experiment, seals were fed Atlantic pollock, a low-fat prey, for 30 days. (B) Plots of estimates derived from MixSIAR models based on FA data of whole blubber cores from days 0, 14 and 30. The black asterisks in the plots indicate the true diet. CCs correspond to harp seals fed herring, calculated using FAs from day 0. Images: fish by Lukas Guerrero Zambra, harp seal by Alicia Guerrero.

Full size image

The model
Since seals had been feeding on the same source for a year, we calculated CCs using FA data from day 0. Thus, consumer values were divided by Atlantic herring values producing CCs that were applied to both herring and pollock FA data. Sources were significantly different (PERMANOVA, F1 = 3680.4, P = 0.001) and ‘day’ was set as fixed factor in the model.
Diet predictions
For all three data sets (day 0, 14 and 30) the main contributor to the diet was Atlantic herring. The predicted proportion of Atlantic herring decreased only slightly from 99% on day 0 to 95% on day 30. Consequently, the contribution of pollock increased from 1% on day 0, to 5% on day 30 (Fig. 6B).
Case 6: Harbour seals, from Nordstrom et al.46
The experiment
Estimations are based on a feeding experiment by Nordstrom et al.46. Prior to the feeding study, juvenile harbour seals, Phoca vitulina, were fed a homogenate of 2:1 Pacific herring to salmon oil for approximately three weeks, and then fed only Pacific herring for four to six days. The feeding experiment consisted of three diets: one group of seals was fed only herring for 42 days (n = 3); the second group was fed only surf smelt for the same period (n = 6); and a third group (n = 7) was fed smelt for 21 days and then only herring for 21 days (Fig, 7A). Whole blubber core samples were collected on days 0, 21 and 42 for each group. Complete FA turnover was estimated to occur after at least 55 days on the same diet, although it could extend well beyond if turnover rate slowed with time.
The model
Since prey FA data was not provided in the same study, we used Pacific herring, surf smelt, and salmon FA values from Huynh and Kitts47, which had significantly different FA composition (PERMANOVA, F2 = 87.7, P = 0.001). CCs were derived from other harbour seals on a Pacific herring diet for over a year, in the same study46. We estimated the diet of the three groups of harbour seals, based on samples collected on day 42, setting ‘diet group’ as fixed factor in our model.
Diet predictions
Overall, diet differences were evident among groups, and the direction of the change was consistent with the shifts in diet. Estimates for harbour seals fed exclusively Pacific herring for 42 days, correctly showed that diet was predominantly based on herring (94.7%). For seals fed surf smelt for 42 days; however, estimates showed that surf smelt only accounted for 26.6% of the diet whereas herring remained to be the main component. For seals fed surf smelt for 21 days and then herring for another 21 days, MixSIAR again identified herring as the main component, with 90.9%, whereas surf smelt was only 3% (Fig. 7B).
Figure 7

Diet estimates for harbour seals. (A) Feeding experiment: For 42 days, seals were fed the following diets: solely herring (n = 3), solely surf smelt (n = 6), or surf smelt for the first 21 days and then herring for the remaining 21 days (n = 7). Prior to the feeding experiments they had been fed a mixture of herring and salmon. (B) Plots for MixSIAR diet estimations, based on blubber FAs obtained on day 42. The red asterisks indicate the true diet. CCs were obtained from harbour seals (other individuals) fed herring for over a year. Image by Alicia Guerrero.

Full size image More

1.
Running, S. W. A measurable planetary boundary for the biosphere. Science 337, 1458–1459 (2012).
ADS  CAS  PubMed  Article  Google Scholar 
2.
Haberl, H. et al. Quantifying and mapping the human appropriation of net primary production in Earth’s terrestrial ecosystems. Proc. Natl Acad. Sci. USA 104, 12942–12945 (2007).
ADS  CAS  PubMed  Article  Google Scholar 

3.
Xia, J. et al. Spatio-temporal patterns and climate variables controlling of biomass carbon stock of global grassland ecosystems from 1982 to 2006. Remote Sens. 6, 1783 (2014).
ADS  Article  Google Scholar 

4.
Grace, J. B. et al. Integrative modelling reveals mechanisms linking productivity and plant species richness. Nature 529, 390–393 (2016).
ADS  CAS  PubMed  Article  Google Scholar 

5.
Del Grosso, S. et al. Global potential net primary production predicted from vegetation class, precipitation, and temperature. Ecology 89, 2117–2126 (2008).
PubMed  Article  Google Scholar 

6.
Galloway, J. N. et al. Nitrogen cycles: past, present, and future. Biogeochemistry 70, 153–226 (2004).
CAS  Article  Google Scholar 

7.
Chadwick, O. A., Derry, L. A., Vitousek, P. M., Huebert, B. J. & Hedin, L. O. Changing sources of nutrients during four million years of ecosystem development. Nature 397, 491–497 (1999).
ADS  CAS  Article  Google Scholar 

8.
McNaughton, S. J., Oesterheld, M., Frank, D. A. & Williams, K. J. Ecosystem-level patterns of primary productivity and herbivory in terrestrial habitats. Nature 341, 142–144 (1989).
ADS  CAS  PubMed  Article  Google Scholar 

9.
Foley, J. A. et al. Solutions for a cultivated planet. Nature 478, 337–342 (2011).
ADS  CAS  PubMed  Article  Google Scholar 

10.
White, R. P., Murray, S. & Rohweder, M. Pilot Analysis of Global Ecosystems (PAGE): Grassland Ecosystems, 70 (World Resources Institute, Washington, DC, 2000).

11.
Ripple, W. J. et al. Collapse of the world’s largest herbivores. Sci. Adv. 1, e1400103 (2015).

12.
Intergovernmental Panel on Climate Change. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, 151 (Geneva, Switzerland, 2014).

13.
Canfield, D. E., Glazer, A. N. & Falkowski, P. G. The evolution and future of earth’s nitrogen cycle. Science 330, 192–196 (2010).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

14.
Stevens, C. J. Nitrogen in the environment. Science 363, 578–580 (2019).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

15.
Stevens, C. J. et al. Anthropogenic nitrogen deposition predicts local grassland primary production worldwide. Ecology 96, 1459–1465 (2015).
Article  Google Scholar 

16.
Reyer, C. P. O. et al. A plant’s perspective of extremes: terrestrial plant responses to changing climatic variability. Glob. Change Biol. 19, 75–89 (2013).
ADS  Article  Google Scholar 

17.
Knapp, A. K. et al. Consequences of more extreme precipitation regimes for terrestrial ecosystems. BioScience 58, 811–821 (2008).
MathSciNet  Article  Google Scholar 

18.
LeBauer, D. S. & Treseder, K. K. Nitrogen limitation of net primary productivity in terrestrial ecosystems is globally distributed. Ecology 89, 371–379 (2008).
PubMed  Article  PubMed Central  Google Scholar 

19.
Fay, P. A. et al. Grassland productivity limited by multiple nutrients. Nat. Plants 1, 15080 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

20.
Borer, E. T. et al. Herbivores and nutrients control grassland plant diversity via light limitation. Nature 508, 517–520 (2014).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

21.
Knapp, A. K. & Seastedt, T. R. Detritus accumulation limits productivity of tallgrass prairie: the effects of its plant litter on ecosystem function make the tallgrass prairie unique among North American biomes. BioScience 36, 662–668 (1986).
Article  Google Scholar 

22.
Volterra, V. Variations and fluctuations of the numbers of individuals in animal species living together. Nature 118, 558–560 (1926).
ADS  Article  Google Scholar 

23.
Crawley, M. J. Herbivory: The Dynamics of Animal-Plant Interactions, Vol. 10, 437 (University of California Press, 1983).

24.
Oksanen, L. & Oksanen, T. The logic and realism of the hypothesis of exploitation ecosystems. Am. Nat. 155, 703–723 (2000).
PubMed  Article  PubMed Central  Google Scholar 

25.
Gruner, D. S. et al. A cross-system synthesis of consumer and nutrient resource control on producer biomass. Ecol. Lett. 11, 740–755 (2008).
PubMed  Article  PubMed Central  Google Scholar 

26.
Hillebrand, H. et al. Consumer versus resource control of producer diversity depends on ecosystem type and producer community structure. Proc. Natl Acad. Sci. USA 104, 10904–10909 (2007).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

27.
Worm, B., Lotze, H. K., Hillebrand, H. & Sommer, U. Consumer versus resource control of species diversity and ecosystem functioning. Nature 417, 848–851 (2002).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

28.
Turkington, R. Top-down and bottom-up forces in mammalian herbivore – vegetation systems: an essay review. Botany 87, 723–739 (2009).
Article  Google Scholar 

29.
DeAngelis, D. L. Dynamics of Nutrient Cycling and Food Webs (Chapman and Hall, London, 1992).

30.
Arditi, R. & Ginzburg, L. R. Coupling in predator-prey dynamics: ratio-dependence. J. Theor. Biol. 139, 311–326 (1989).
Article  Google Scholar 

31.
Chase, J. M., Leibold, M. A., Downing, A. L. & Shurin, J. B. The effects of productivity, herbivory, and plant species turnover in grassland food webs. Ecology 81, 2485–2497 (2000).
Article  Google Scholar 

32.
Leibold, M. A. Resource edibility and the effects of predators and productivity on the outcome of trophic interactions. Am. Nat. 134, 922–949 (1989).
Article  Google Scholar 

33.
Endara, M. J. & Coley, P. D. The resource availability hypothesis revisited: a meta-analysis. Funct. Ecol. 25, 389–398 (2011).
Article  Google Scholar 

34.
Milchunas, D. G. & Lauenroth, W. K. Quantitative effects of grazing on vegetation and soils over a global range of environments. Ecol. Monogr. 63, 327–366 (1993).
Article  Google Scholar 

35.
Polis, G. A. & Strong, D. R. Food web complexity and community dynamics. Am. Nat. 147, 813–846 (1996).
Article  Google Scholar 

36.
Murdoch, W. Community structure, population control, and competition – a critique. Am. Nat. 100, 219–226 (1966).
Article  Google Scholar 

37.
Borer, E. T. et al. Finding generality in ecology: a model for globally distributed experiments. Methods Ecol. Evolution 5, 65–73 (2014).
Article  Google Scholar 

38.
Anderson, T. M. et al. Herbivory and eutrophication mediate grassland plant nutrient responses across a global climatic gradient. Ecology 99, 822–831 (2018).

39.
Zhu, D. et al. The large mean body size of mammalian herbivores explains the productivity paradox during the Last Glacial Maximum. Nat. Ecol. Evol. 2, 640–649 (2018).
PubMed  PubMed Central  Article  Google Scholar 

40.
Hillebrand, H. Top-down versus bottom-up control of autotrophic biomass – a meta-analysis on experiments with periphyton. J. North Am. Benthol. Soc. 21, 349–369 (2002).
Article  Google Scholar 

41.
Frank, R. & Merle, L. F. Effects of annual applications of low N fertilizer rates on a mixed grass prairie. J. Range Manag. 36, 359–362 (1983).
Article  Google Scholar 

42.
Olofsson, J. et al. Long-term experiments reveal strong interactions between lemmings and plants in the Fennoscandian highland tundra. Ecosystems 17, 606–615 (2014).
Article  Google Scholar 

43.
Lemaire, G., Jeuffroy, M.-H. & Gastal, F. Diagnosis tool for plant and crop N status in vegetative stage: theory and practices for crop N management. Eur. J. Agron. 28, 614–624 (2008).
CAS  Article  Google Scholar 

44.
Hillebrand, H. et al. Herbivore metabolism and stoichiometry each constrain herbivory at different organizational scales across ecosystems. Ecol. Lett. 12, 516–527 (2009).
PubMed  Article  PubMed Central  Google Scholar 

45.
Hempson, G. P., Illius, A. W., Hendricks, H. H., Bond, W. J. & Vetter, S. Herbivore population regulation and resource heterogeneity in a stochastic environment. Ecology 96, 2170–2180 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

46.
Sickman, J. O. et al. Quantifying atmospheric N deposition in dryland ecosystems: a test of the Integrated Total Nitrogen Input (ITNI) method. Sci. Total Environ. 646, 1253–1264 (2019).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

47.
Yahdjian, L., Gherardi, L. & Sala, O. E. Nitrogen limitation in arid-subhumid ecosystems: a meta-analysis of fertilization studies. J. Arid Environ. 75, 675–680 (2011).
ADS  Article  Google Scholar 

48.
Koerner, S. E. et al. Plant community response to loss of large herbivores differs between North American and South African savanna grasslands. Ecology 95, 808–816 (2014).
PubMed  Article  Google Scholar 

49.
Frank, D. A., McNaughton, S. J. & Tracy, B. F. The ecology of the Earth’s grazing ecosystems: profound functional similarities exist between the Serengeti and Yellowstone. Bioscience 48, 513–521 (1998).
Article  Google Scholar 

50.
Augustine, D. J. & McNaughton, S. J. Interactive effects of ungulate herbivores, soil fertility, and variable rainfall on ecosystem processes in a semi-arid savanna. Ecosystems 9, 1242–1256 (2006).
CAS  Article  Google Scholar 

51.
Ritchie, M. E., Tilman, D. & Knops, J. M. H. Herbivore effects on plant and nitrogen dynamics in oak savanna. Ecology 79, 165–177 (1998).
Article  Google Scholar 

52.
Pastor, J., Dewey, B., Naiman, R. J., McInnes, P. F. & Cohen, Y. Moose browsing and soil fertility in the boreal forests of Isle Royale National Park. Ecology 74, 467–480 (1993).
Article  Google Scholar 

53.
Grellmann, D. Plant responses to fertilization and exclusion of grazers on an Arctic tundra heath. Oikos 98, 190–204 (2002).
Article  Google Scholar 

54.
Hartley, S. E. & Mitchell, R. J. Manipulation of nutrients and grazing levels on heather moorland: changes in Calluna dominance and consequences for community composition. J. Ecol. 93, 990–1004 (2005).
Article  Google Scholar 

55.
Lind, E. M. et al. Increased grassland arthropod production with mammalian herbivory and eutrophication: a test of mediation pathways. Ecology 98, 3022–3033 (2017).
PubMed  Article  PubMed Central  Google Scholar 

56.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).
Article  Google Scholar 

57.
Trabucco, A., Zomer, R. J., Bossio, D. A., van Straaten, O. & Verchot, L. V. Climate change mitigation through afforestation/reforestation: a global analysis of hydrologic impacts with four case studies. Agric. Ecosyst. Environ. 126, 81–97 (2008).
Article  Google Scholar 

58.
Dentener, F. J. Global Maps of Atmospheric Nitrogen Deposition, 1860, 1993, and 2050. Oak Ridge National Laboratory Distributed Active Archive Center. https://doi.org/10.3334/ORNLDAAC/830 (2006).

59.
Borer, E. T. et al. Environmental Data Initiative https://doi.org/10.6073/pasta/a318fe0fb11eb43c1a2c8233b2e3494f (2020). More

1.
Odum, E. P. & Barrett, G. W. Fundamentals of Ecology (Saunders, Philadelphia, 1953).
Google Scholar 
2.
Elton, C. S. Ecology of Invasions by Animals And Plants (Methuen, London, 1958).
Google Scholar 

3.
MacArthur, R. Fluctuations of animal populations and a measure of community stability. Ecology 36, 533–536 (1955).
Article  Google Scholar 

4.
Paine, R. T. Food web complexity and species diversity. Am. Nat. 100, 65–75 (1966).
Article  Google Scholar 

5.
Landi, P., Minoarivelo, H. O., Å, Brännström, Hui, C. & Dieckmann, U. Complexity and stability of ecological networks: a review of the theory. Popul. Ecol. 60, 319–345 (2018).
Article  Google Scholar 

6.
McCann, K. S. The diversity–stability debate. Nature 405, 228–233 (2000).
CAS  PubMed  Article  Google Scholar 

7.
May, R. M. Will a large complex system be stable? Nature 238, 413–414 (1972).
ADS  CAS  PubMed  Article  Google Scholar 

8.
May, R. M. Stability in multispecies community models. Math. Biosci. 12, 59–79 (1971).
MathSciNet  MATH  Article  Google Scholar 

9.
Namba, T. Multi-faceted approaches toward unravelling complex ecological networks. Popul. Ecol. 57, 3–19 (2015).
Article  Google Scholar 

10.
Justus, J. A case study in concept determination: Ecological diversity. in Philosophy of Ecology, Handbook of the Philosophy of Science, Vol. 11, (eds deLaplante, K., Brown, B. & Peacock, K. A.) (North-Holland, Amsterdam, 2011) pp. 147–168.

11.
Grilli, J., Rogers, T. & Allesina, S. Modularity and stability in ecological communities. Nat. Commun. 7, 1–10 (2016).
Article  CAS  Google Scholar 

12.
Allesina, S. & Tang, S. The stability–complexity relationship at age 40: a random matrix perspective. Popul. Ecol. 57, 63–75 (2015).
Article  Google Scholar 

13.
Allesina, S. & Tang, S. Stability criteria for complex ecosystems. Nature 483, 205–208 (2012).
ADS  CAS  PubMed  Article  Google Scholar 

14.
Hutchinson, M. C. Seeing the forest for the trees: putting multilayer networks to work for community ecology. Funct. Ecol. 33, 206–217 (2019).
Article  Google Scholar 

15.
Pilosof, S., Porter, M. A., Pascual, M. & Kéfi, S. The multilayer nature of ecological networks. Nat. Ecol. Evolut. 1, 1–9 (2017).
Article  Google Scholar 

16.
Stone, L. The feasibility and stability of large complex biological networks: a random matrix approach. Sci. Rep. 8, 1–12 (2018).
CAS  Article  Google Scholar 

17.
Gibbs, T., Grilli, J., Rogers, T. & Allesina, S. Effect of population abundances on the stability of large random ecosystems. Phys. Rev. E 98, 022410 (2018).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

18.
Fyodorov, Y. V. & Khoruzhenko, B. A. Nonlinear analogue of the may-wigner instability transition. Proc. Natl Acad. Sci. USA 113, 6827–6832 (2016).
MathSciNet  CAS  PubMed  MATH  Article  Google Scholar 

19.
Gross, T., Rudolf, L., Levin, S. A. & Dieckmann, U. Generalized models reveal stabilizing factors in food webs. Science 325, 747–750 (2009).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

20.
Berlow, E. L. et al. Interaction strengths in food webs: issues and opportunities. J. Anim. Ecol. 73, 585–598 (2004).
Article  Google Scholar 

21.
Chesson, P. & Huntly, N. The roles of harsh and fluctuating conditions in the dynamics of ecological communities. Am. Nat. 150, 519–553 (1997).
CAS  PubMed  Article  Google Scholar 

22.
McNaughton, S. J. Diversity and stability of ecological communities: a comment on the role of empiricism in ecology. Am. Nat. 111, 515–525 (1977).
Article  Google Scholar 

23.
Thébault, E. & Loreau, M. Trophic interactions and the relationship between species diversity and ecosystem stability. Am. Nat. 166, E95–E114 (2005).
PubMed  Article  Google Scholar 

24.
Tilman, D. The ecological consequences of changes in biodiversity: a search for general principles. Ecology 80, 1455–1474 (1999).
Google Scholar 

25.
Gravel, D., Massol, F. & Leibold, M. A. Stability and complexity in model meta-ecosystems. Nat. Commun. 7, 12457 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

26.
Turing, A. M. The chemical basis of morphogenesis. Bull. Math. Biol. 52, 153–197 (1990).
CAS  PubMed  Article  Google Scholar 

27.
Hanski, I. Metapopulation dynamics. Nature 396, 41–49 (1998).
ADS  CAS  Article  Google Scholar 

28.
Hanski, I. Habitat connectivity, habitat continuity, and metapopulations in dynamic landscapes. Oikos 87, 209–219 (1999).

29.
Hassell, M. P., Comins, H. N. & May, R. M. Spatial structure and chaos in insect population dynamics. Nature 353, 255–258 (1991).
ADS  Article  Google Scholar 

30.
Hassell, M. P., Comins, H. N. & May, R. M. Species coexistence and self-organizing spatial dynamics. Nature 370, 290–292 (1994).
ADS  Article  Google Scholar 

31.
Levin, S. A. & Segel, L. A. Hypothesis for origin of planktonic patchiness. Nature 259, 659 (1976).
ADS  Article  Google Scholar 

32.
Brechtel, A., Gramlich, P., Ritterskamp, D., Drossel, B. & Gross, T. Master stability functions reveal diffusion-driven pattern formation in networks. Phys. Rev. E 97, 032307 (2018).
ADS  CAS  PubMed  Article  Google Scholar 

33.
Cross, M. & Hohenberg, P. Pattern formation outside of equilibrium. Rev. Mod. Phys. 65, 851–1112 (1993).
ADS  CAS  MATH  Article  Google Scholar 

34.
Grilli, J. et al. Feasibility and coexistence of large ecological communities. Nat. Commun. 8, 14389 (2017).
PubMed Central  Article  CAS  PubMed  Google Scholar 

35.
Allesina, S. et al. Predicting the stability of large structured food webs. Nat. Commun. 6, 1–6 (2015).
Article  CAS  Google Scholar 

36.
Neubert, M. G., Kot, M. & Lewis, M. A. Dispersal and pattern formation in a discrete-time predator-prey model. Theor. Popul. Biol. 48, 7–43 (1995).
MATH  Article  Google Scholar 

37.
Rietkerk, M. & Van de Koppel, J. Regular pattern formation in real ecosystems. Trends Ecol. Evolut. 23, 169–175 (2008).
Article  Google Scholar 

38.
HilleRisLambers, R., Rietkerk, M., van den Bosch, F., Prins, H. H. T. & de Kroon, H. Vegetation pattern formation in semi-arid grazing systems. Ecology 82, 50–61 (2001).
Article  Google Scholar 

39.
Levin, S. Dispersion and population interactions. Am. Nat. 108, 207–228 (1974).
Article  Google Scholar 

40.
Murray, J. D. Mathematical Biology II: Spatial Models and Biomedical Applications (Springer: New York, 2001).
Google Scholar 

41.
Baron, J. W. & Galla, T. Stochastic fluctuations and quasipattern formation in reaction-diffusion systems with anomalous transport. Phys. Rev. E 99, 052124 (2019).
ADS  CAS  PubMed  Article  Google Scholar 

42.
Kuramoto, Y. Diffusion-induced chaos in reaction systems. Progr. Theor. Phys. Suppl. 64, 346–367 (1978).
ADS  CAS  Article  Google Scholar 

43.
Pascual, M. Diffusion-induced chaos in a spatial predator–prey system. Proc. R. Soc. Lond. Ser. B 251, 1–7 (1993).
ADS  Article  Google Scholar 

44.
Rietkerk, M., Dekker, S. C., de Ruiter, P. C. & van de Koppel, J. Self-organized patchiness and catastrophic shifts in ecosystems. Science 305, 1926–1929 (2004).
ADS  CAS  PubMed  Article  Google Scholar 

45.
Rietkerk, M. & van de Koppel, J. Regular pattern formation in real ecosystems. Trends Ecol. Evolut. 23, 169–175 (2008).
Article  Google Scholar 

46.
van de Koppel, J. et al. Experimental evidence for spatial self-organization and its emergent effects in mussel bed ecosystems. Science 322, 739–742 (2008).
ADS  PubMed  Article  CAS  Google Scholar 

47.
Meron, E. Pattern-formation approach to modelling spatially extended ecosystems. Ecol. Model. 234, 70–82 (2012).
Article  Google Scholar 

48.
Liu, Q. et al. Pattern formation at multiple spatial scales drives the resilience of mussel bed ecosystems. Nat. Commun. 5, 5234 (2014).
ADS  CAS  PubMed  Article  Google Scholar 

49.
Karig, D. et al. Stochastic turing patterns in a synthetic bacterial population. Proc. Natl Acad. Sci. USA 115, 6572–6577 (2018).
ADS  CAS  PubMed  Article  Google Scholar 

50.
Lengyel, I. & Epstein, I. A chemical approach to designing turing patterns in reaction-diffusion systems. Proc. Natl Acad. Sci. USA 89, 3977–3979 (1992).
ADS  CAS  PubMed  MATH  Article  Google Scholar 

51.
Castets, V., Dulos, E., Boissonade, J. & De Kepper, P. Experimental evidence of a sustained standing turing-type nonequilibrium chemical pattern. Phys. Rev. Lett. 64, 2953–2956 (1990).
ADS  CAS  PubMed  Article  Google Scholar 

52.
Barabás, G., Michalska-Smith, M. J. & Allesina, S. Self-regulation and the stability of large ecological networks. Nat. Ecol. Evolut. 1, 1870–1875 (2017).
Article  Google Scholar 

53.
Tao, T. & Vu, V. Random matrices: universality of local eigenvalue statistics up to the edge. Commun. Math. Phys. 298, 549–572 (2010).
ADS  MathSciNet  MATH  Article  Google Scholar 

54.
Tao, T. et al. Random matrices: universality of ESDs and the circular law. Ann. Probab. 38, 2023–2065 (2010).
MathSciNet  MATH  Article  Google Scholar 

55.
O’Sullivan, J. D., Knell, R. J. & Rossberg, A. G. Metacommunity-scale biodiversity regulation and the self-organised emergence of macroecological patterns. Ecol. Lett. 22, 1428–1438 (2019).
PubMed  Article  Google Scholar 

56.
Gotelli, N. J. et al. Community-level regulation of temporal trends in biodiversity. Sci. Adv. 3, e1700315 (2017).
ADS  PubMed  PubMed Central  Article  Google Scholar 

57.
Magurran, A. E. et al. Divergent biodiversity change within ecosystems. Proc. Natl Acad. Sci. USA 115, 1843–1847 (2018).
CAS  PubMed  Article  Google Scholar 

58.
Brown, J. H., Ernest, S. K. M., Parody, J. M. & Haskell, J. P. Regulation of diversity: maintenance of species richness in changing environments. Oecologia 126, 321–332 (2001).
ADS  PubMed  Article  Google Scholar 

59.
Parody, J. M., Cuthbert, F. J. & Decker, E. H. The effect of 50 years of landscape change on species richness and community composition. Glob. Ecol. Biogeogr. 10, 305–313 (2001).
Article  Google Scholar 

60.
Haake, F., Izrailev, F., Lehmann, N., Saher, D. & Sommers, H.-J. Statistics of complex levels of random matrices for decaying systems. Z. Phys. B Condens. Matter 88, 359–370 (1992).
ADS  Article  Google Scholar 

61.
Sommers, H.-J., Crisanti, A., Sompolinsky, H. & Stein, Y. Spectrum of large random asymmetric matrices. Phys. Rev. Lett. 60, 1895–1898 (1988).
ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

62.
O’Rourke, S. et al. Low rank perturbations of large elliptic random matrices. Electron. J. Probab. 19, 1–65 (2014).

63.
Süli, E. & Mayers, D. F. An Introduction to Numerical Analysis (Cambridge University Press, 2003).

64.
Baron, J. W. & Galla, T. Dispersal-induced instability in complex ecosystems. GitHub https://doi.org/10.5281/zenodo.4068257 (2020). More

Basic characters of the six chloroplast genomes
The cp genomes of P. cicutarrifolia, P. hubeiensis, P. jiugongshanensis, P. merrilliana and P. ranunculoides (GenBank accessions: MT268974, MT268976, MT937162, MT268977, MT268978) were reported for the first time here, and that of P. filchnerae was downloaded from NCBI (MK88869821).
The sequencing coverage of our five newly assembled cp genomes was from 923 to 6237 (Figure S1). The six cp genomes possessed typical quadripartite structure: IRa, IRb, LSC and SSC (Table 1), and they exhibited the same gene order, no gene rearrangement or inversion occurred (Figure S2). The physical map of the cp genome of P. hubeiensis was shown in Fig. 1. The GC content was ~ 37%. The newly sequenced genomes ranged from 150,187 bp to 151,972 bp, harboring 113 genes: four ribosomal RNA genes, 29 tRNA genes and 80 protein-coding genes, and among them 14 genes was duplicated in IRa and IRb (Table 1). Due to presence of multiple stop codons, the gene infA was pseudogenized in the five newly sequenced species. The open reading frame (ORF) in accD of P. filchnerae (MK888698) was truncated to be only 1305 bp compared with 1455 or 1464 bp ORF of other five species. Lee et al.34 identified five conserved amino acid sequence motifs in accD gene. Conserved amino acid sequence motifs IV and V were absent in accD of P. filchnerae. Therefore, accD was nonfunctional in P. filchnerae.
Table 1 Basic characteristics of cp genomes of the six Primula species (Pc: P. cicutarrifolia; Pf: P. filchnerae; Ph: P. hubeiensis; Pj: P. jiugongshanensis; Pm: P. merrilliana; Pr: P. ranunculoides).
Full size table

Figure 1

Physical map of the P. hubeiensis chloroplast genome.

Full size image

SSRs and repeats
Five categories of SSRs were identified for the six species (Table 2). The least number of SSRs was 41 for P. ranunculoides and the most 59 for P. merrilliana. Three types of SSRs were detected for P. filchnerae, and in the rest species four types could be found. While mono-, di- and tetra-nucleotide repeats existed across all the six species, tri- and penta-inucleotide repeats resided in three and two species respectively. Mono- and dinucleotide repeats accounted for the vast majority of SSRs (65.1% for P. cicutariifolia, 87.5% for P. filchnerae, 69.0% for P. hubeiensis, 62.8% for P. jiugongshanensis, 72.9% for P. merrilliana, 73.2% for P. ranunculoides). Most or all mono- repeats were A/T repeats including 10 to 16 nucleotides. The number of repeat units ranged from five to eight for dinucleotide repeats. The tri- and penta-nucleotide SSRs consisted of four motifs, and tetra-nucleotide SSRs of four to five repeat units.
Table 2 Types and numbers of SSRs in the cp genomes of six Primula species, the numbers in the bracket indicating total number of SSRs (Pc: P. cicutarrifolia; Pf: P. filchnerae; Ph: P. hubeiensis; Pj: P. jiugongshanensis; Pm: P. merrilliana; Pr: P. ranunculoides).
Full size table

Except the largest repeat for each genome (i.e. IRs), a total of 183 repeat pairs (three types: forward (F), reverse (R), and palindromic repeats (P)) were detected in the six genomes (Fig. 2), which ranged from 30 to 137 bp in length. Palindromic repeats were the most common, accounting for 55.2% (101 of 183), followed by forward repeats (44.3%, 81 of 183). No complement repeats were identified in all species and one pair of reverse repeats existed specifically in P. ranunculoides. In the six species, 96.7% (177 of 183 repeat pairs) repeats were 30–59 bp in length, consistent with the length reported in other Primula species20. The longest repeat (137 bp) was found in P. cicutariifolia, and this species contained the most repeats (44 pairs), while P. jiugongshanensis had the least (24 pairs).
Figure 2

Types and numbers of repeat pairs in the cp genomes of six Primula species (Pc: P. cicutarrifolia; Pf: P. filchnerae; Ph: P. hubeiensis; Pj: P. jiugongshanensis; Pm: P. merrilliana; Pr: P. ranunculoides).

Full size image

IR/SC boundary
The IR/SC boundary regions of the six Primula cp genomes were compared, and the IR/SC junction regions showed slight differences in the length of organization genes flanking the junctions or the distance between the junctions and the organization genes (Fig. 3). The genes spanning or flanking the junction of LSC/IRb, IRb/SSC, SSC/IRa and IRa/LSC were rps19/rpl2, ndhF, ycf1, rpl2/trnH, respectively. IR expansion and contraction was observed. P. cicutarrifolia had the smallest size of IR but largest size of both LSC and SSC; though largest size of IR was detected in P. filchnerae, the LSC or SSC was not the smallest in this species. The gene trnH was located in LSC, 0–24 bp away from the IRa/LSC border. The largest extensions of ycf1 into both SSC and IRa occurred in P. filchnerae (4566 bp and 1023 bp, respectively) and ycf1 of P. filchnerae were the longest among the six species. The gene ndhF was utterly situated in SSC and 108 bp distant from the IRb/SSC junction in P. cicutarrifolia; in the rest five species the fragment size of ndhF in SSC was largest in P. hubeiensis (2194 bp). In P. cicutarrifolia, P. jiugongshanensis and P. merrilliana, rps19 and rpl2 were located in the upstream and downstream of the LSC/IRb junction, respectively; rps19 ran across the LSC/IRb junction in P. filchnerae, P. hubeiensis, P. ranunculoides with 161, 62, 56 bp extension in IRb, respectively.
Figure 3

LSC/IR, and SSC/IR border regions of the six Primula cp genomes.

Full size image

Divergent hotspots in the Primula chloroplast genome
As indicated by the value of Pi, the nucleotide variability of the 22 Primula species (Table S1) was evaluated by DnaSP 6.1231 using noncoding sequences (intron and intergenic spacer) or protein coding sequences (CDS) at least 200 bp long. The variation level of DNA polymorphorism was 0.00444–0.11369 for noncoding sequences or 0.00094–0.05036 for CDSs. For the CDSs, the highest Pi value were detected for ycf1 (0.05036), followed by matK (0.04878), rpl22 (0.04364), ndhF (0.03975), rps8 (0.03658), ndhD (0.03455), ccsA (0.03292), rpl33 (0.0303), rps15 (0.03022), and rpoC2 (0.02954). These markers had higher Pi than rbcL (0.02149). Obviously, the gene ycf1 exhibited the greatest diversity and harbored the most abundant variation. The ten most divergent regions among noncoding regons included trnH (GUG)-psbA (0.11369), trnW (CCA)-trnP (UGG) (0.09463), rpl32-trnL (UAG) (0.09337), ndhC-trnV (UAC) (0.09148), ccsA-ndhD (0.08745), ndhG-ndhI (0.08363), trnK (UUU)-rps16 (0.08334), trnM (CAU)-atpE (0.08273), trnS (GGA)-rps4 (0.08028), and trnC (GCA)-petN (0.07971). No intron ranked among the top ten variable noncoding regions.
Phylogenetic analysis
The ML tree of 22 Primula species was constructed with RAxML32 (Fig. 4), based on the whole cp genomes. The six pinnate-leaved Primula species did not form a monophyly, but separated into two distant clades. P. filchnerae grouped with P. sinensis, the other five species clustered together and constituted the clade Sect. Ranunculoides with 100% bootstrap. In the ML tree, Sect. Proliferae exhibited monophyly, while species of Sect. Crystallophlomis separated into different clades.
Figure 4

ML phylogenetic tree of Primula species based on cp genomes. Bootstrap support at nodes are all 100%.

Full size image

The topology of the ML tree based on ycf1 (Figure S3) was consistent with that based on whole cp genomes (Fig. 4), except that the clade formed by P. veris and P. knuthiana were sister to the clade consisting of Sects. Auganthus, Obconicolisteri, Carolinella and Monocarpicae instead of being sister to the clade of Sects. Proliferae, Ranunculoides and Crystallophlomis.
We also constructed both ML and NJ tree of 71 Primula species based on the concatenation of three common barcoding markers (ITS, matK and rbcL). Only the results of NJ analysis (Fig. 5) showed consistency with those of Yan et al.12, Liu et al.35, and ML analysis based on whole cp genomes (Fig. 4). The six pinnate-leaved Primula species were separated into two distantly related groups. The clade consisting of P. filchnerae and P. sinensis (Sect. Auganthus) was sister to the clade formed by Sects. Carolinella, Obconicolisteri, Monocarpicae, Cortusoides, Malvacea, Pycnoloba. The five pinnatisect-leaved species P. cicutarrifolia, P. hubeiensis, P. jiugonshanensis, P. merrilliana and P. ranunculoides (Sect. Ranunculoides) comprised a 100% supported clade, which was sister to the group containing Sects. Crystallophlomis, Petiolares, Proliferae, Amethystina. Sect. Carolinella and Sect. Crystallophlomis, and Sect. Malvacea were polyphyletic.
Figure 5

NJ bootstrap consensus tree of Primula based on concatenation of ITS, matK and rbcL.

Full size image More