Ecology

Our results are consistent with a non-industrial gut harboring a more resilient ecology with respect to SCFA production, while the industrial gut ecology would be vulnerable to disruption of such pathways, yet the pattern is complex and nuanced. The increased gene abundance in non-industrial populations and overall ratio of acetate:butyrate:propionate generally agrees with previous studies of SCFAs5,12. Similarly, the higher genus-level diversity of bacteria encoding acetate, compared to the other SCFAs, is expected and matches studies that have documented the taxa that encode different SCFAs13,17,19. The overall high richness, high diversity at Hill numbers 1 and 2, and high Gini-Simpson indices found in non-industrial populations at the genus level indicates a highly diverse and evenly distributed production of SCFAs. From an ecological perspective, uneven production of SCFA dominated by a few bacteria in industrial gut microbiomes means lower functional diversity and less redundancy, which ultimately leads to an expectation of decreased resilience. In other words, this study finds that industrial gut microbiomes are at a higher risk of reduced SCFA production because SCFA synthesis is dominated by only a few genera. Given the lower resilience, factors that disrupt the gut ecology are expected to have a more extreme consequence to those living an industrial lifestyle.
While there is an overall trend of increased genus-level functional diversity and redundancy for SCFA production in non-industrial populations, variation exists when examining the SCFAs and populations individually. At the genus-level, the pastoral and rural agricultural populations have increased richness of genera encoding genes involved in acetate and butyrate synthesis, while there is similarity across the different lifestyles for genus richness for propionate encoding taxa. Although hunter-gatherers have similar, or lower, genus richness as industrial populations, they have significantly higher diversity at Hill number orders 1 and 2 and Gini-Simpson indices for butyrate and propionate. Additionally, the pastoralists have a generally similar profile to the industrial populations for acetate and propionate Hill number diversity, as well as similarity to the industrial populations in species PD, which may be linked to this pastoralist group having a diet similar to some industrial populations; namely, a diet high in dairy and red meat consumption, coupled with few dietary sources of plant-derived fibers23. This paints a complex picture. Non-industrial populations have a high diversity of genera encoding butyrate synthesis, and butyrate production is spread more evenly across genera in non-industrial populations than in industrial populations. Hunter-gatherers and rural agriculturalists have significantly greater evenness of propionate production, even though they have fewer number of total genera encoding this SCFA. Finally, the richness and evenness of genera encoding acetate is similar between industrial and non-industrial populations. Ecologically, we would expect the industrial populations to be less resilient for production of butyrate and propionate when faced with a shift in taxonomic composition, while non-industrial populations may be only marginally more resilient for acetate production compared to industrial populations. Intriguingly, SCFA relative abundance does not appear to correlate to resilience profile. Acetate and butyrate are significantly more abundant in non-industrial populations but only butyrate shows much stronger resilience profile for non-industrial populations. Additionally, propionate is slightly more abundant in industrial populations, although not significantly, yet our results indicate greater resilience in non-industrial groups for propionate production. This indicates that measuring only total gene, and/or molar, abundance is not enough to make statements about metabolic processes in the human microbiome; rather, ecological approaches are necessary to understand diversity in functional potential of the human microbiome.
The increased species-level alpha diversity in industrial populations initially runs counter to the genus-level results but the genus and species level results ultimately yield similar interpretations after accounting for ecology and ascertainment bias, as discussed below. The substantially higher species richness in industrial populations is striking; however, the differences in PD between industrial and non-industrial populations are not nearly as extreme. This means that the high species richness in the industrial populations is driven by species that are closely phylogenetically related. Indeed, we observed SCFA producing genera found at high abundance in industrial populations (Bacteroides and Clostridium) to have up to nine species encoding SCFAs, while highly abundant non-industrial genera only have one or two species. Therefore, what first appears to indicate high species-level ecological resilience in SCFA production in the industrial populations is actually the result of closely related species performing the same function. It follows that closely related species may be prone to changes in abundance or even elimination after certain types of ecosystem shift events. For example, narrow-spectrum antibiotics33 and exposure to various xenobiotic compounds that lead to variable bacterial metabolic responses34 are events that can affect a limited range of bacteria and lead to shifts in microbial abundance and metabolic activity. While this result has ecological implications, it is also likely the result of historical trends of microbiology research. Bacterial taxa at high abundance in non-industrial gut microbiomes have not been a focus of microbiological isolation and species identification until recently; therefore, we expect more species to be identified from non-industrial gut microbiomes in the future35. Additionally, classification of bacteria into distinct genera and species is undergoing a revolution in the genomic era36 meaning that the high number of species classified to Bacteroides and Clostridium may ultimately be reclassified to different genera. Nevertheless, the fact that we observe a large jump in species richness, but only a minor increase in species PD, in the industrial gut microbiomes suggests that the high industrial species richness is driven by closely related species and therefore, results in the same interpretation as the genus richness results: diversity is high in non-industrial populations.
Ascertainment bias extends to the databases used to identify taxa and genes: fewer genes were identified in non-industrial populations and a smaller proportion of these genes can be linked back to bacteria at every taxonomic level, in non-industrial gut microbiomes. In some cases, such as butyrate synthesis genes, less than 10% of genes are identified to species for non-industrial populations, while over 50% of such identifications were possible for industrial populations. A decreased ability to identify the genus and species encoding SCFA synthesis genes in non-industrial populations means that the ecological metrics underestimate the true ecological diversity of these genes. Moreover, the drop-off in classification from the genus to the species level was significantly greater in non-industrial populations compared to industrial populations. This drop-off means a much lesser ability to identify species compared to genera in non-industrial populations, which helps explain why species diversity was substantially lower in non-industrial populations. Nevertheless, the statistically significant differences observed at the genus-level send a strong signal of the high functional diversity, and potential resilience, of SCFA synthesis genes in non-industrial gut microbiomes.
The metagenome-wide poor performance in terms of gene identification and classifying SCFA genes to genera and species indicates a bias in reference databases that underrepresents diversity in non-industrial gut microbiomes, which is unsurprising. Bias is expected because the vast majority of human gut microbiome studies have used samples from industrial populations. There is an immense challenge in including non-industrial communities in biomedical research, including recruiting research participants, sustaining longitudinal sampling, building culturally appropriate community relationships, and even securing transport of samples35. This has resulted in comparatively few metagenomic studies of human gut microbiomes from non-industrial settings35. Nevertheless, our data demonstrate the extent of this bias and how it can hinder more in-depth study of human gut microbiome health. Given this sizable ascertainment bias favored industrial populations, the non-industrial populations are likely even more diverse, more resilient, than our databases can sufficiently characterize, making our genus-level results even stronger. Without a serious investment to include such populations, the characterization of microbiomes will remain naive to the ecological breadth of the core, healthy, human gut. Imagine studying forest ecology, with only city parks at your disposal. This has been, overwhelmingly, the analogous practice of human microbiome research.
The relative lack of microbiome studies with non-industrial populations means an underrepresentation of not only metagenomic data and genome annotation but also fewer opportunities for cultivation and validation of novel species of bacteria. This ultimately leads to an inequality in the depth to which researchers can describe microbiome samples from non-industrial communities, compared to industrial microbiomes, as diverse groups of novel taxa may be grouped into a single group of “unknown” or “unclassified” bacteria35. Similarly, an incomplete picture of microbial functional potential means that genes may be misidentified or even unannotated completely. Unknown taxa and misidentified genes may be playing key roles in ecological and metabolic processes but researchers are unable to confidently identify them, let alone make statements about their importance in a microbial ecology35. Recent human gut microbiome metagenome studies from diverse populations will undoubtedly improve database representation but the number of studies and metagenomic samples from non-industrial populations still pales in comparison to industrial gut microbiomes26,35,37,38.
Limitations in annotating the full extent of microbial diversity impacts health research. Recently proposed ‘Microbiota Insufficiency Syndrome (MIS)’2 postulates that, while the microbiome has adapted to industrialization, these adaptations are maladaptive to human health. The decreased phylogenetic diversity and loss of specific taxa (e.g. Prevotellaceae, Succinivibrionaceae, and Spirochaetaceae) observed in industrial gut microbiomes may contribute to the increase in non-communicable chronic diseases found at higher prevalence in industrial populations. The root cause of MIS in industrial populations is undoubtedly multifactorial; however, diet is suggested to play a major role2. This syndrome is compelling and we postulate that this insufficiency precisely rests on the stability of functional capacity. Our findings of decreased resilience in industrial populations, as well as species-level diversity driven by a few closely related species, fits in well with MIS. Low resilience in SCFA production may ultimately manifest itself as altered colonocyte function and/or autoimmune disruptions (both symptoms of MIS) due to a decrease in SCFA bioavailability after a group SCFA-producing bacteria were wiped-out during an ecological shift, such as antibiotic or xenobiotic exposure. Similar to MIS, diet is likely to play an important role in SCFA resilience. The non-industrial populations studied in this paper consume much more fiber than industrial populations, on average3,5,14,25,26, and microbial fermentation of dietary fibers is a major source of SCFAs in the human digestive tract39. A diet poor in dietary fiber means less substrate for microbial fermentation and therefore less SCFA production and also higher competition for that fiber, potentially resulting in competitive exclusion and less microbial diversity. Nevertheless, if we are unable to fully characterize and annotate non-industrial gut microbiomes then we will be unable to paint a complete picture of MIS. Currently, we have confidence that there is a wealth of undiscovered resilience in non-industrial gut microbiomes. Once we describe the extent of this diversity/resilience, through increased sampling and focus on partnerships with research institutes in industrializing countries, we will have a more complete picture of MIS and possibly develop therapeutic approaches to combat non-communicable chronic diseases related to the human gut microbiome.
Improved sampling, metabolic profiling, and annotation will not only improve our understanding of SCFA resilience, but it will also permit more detailed picture microbiome wide resilience. Our work shows the value of focusing on specific SCFA genes, due to their importance in human biology and previously reported variation in SCFA molar abundance between industrial and non-industrial populations31,32; however, future work will undoubtedly add to our findings. One avenue for future work is through analyzing SCFA molar concentrations in fecal samples in a longitudinal setting and comparing these results to predicted SCFA resilience from metagenome panels. Unlike genomic data, where we can infer about SCFA production potential via taxonomic diversity, one-time measures of fecal SCFA molar concentrations will not inform about future resilience because SCFA molar concentrations carry no information about which taxa produce each SCFA. Longitudinal SCFA concentration and metagenomic data from non-industrial populations, or animal models, is necessary to inform about SCFA resilience and production in diverse lifestyles. Another avenue for future work is to focus resilience analysis on other microbiome functions of interest, such as resilience of antibiotic resistance genes and amino acid biosynthetic pathways. These valuable studies would be valuable for comparing microbiome resilience dynamics for different functions, with the caveat that there is sufficient genomic annotation data to yield interpretable results.
Lack of sample diversity is not unique to human microbiome research, as human genetics research has been grappling with this very issue for decades. In 2009, 96% of individuals included in human genome-wide association studies (GWAS) claimed European ancestry, as compared to 78% in 201940. Thus, while there have been improvements, GWAS clearly fail to reflect the breadth of human diversity. Incorporating diverse populations in human genome and microbiome research has the potential to greatly benefit the scientific community’s understanding of human biology and develop treatments that are based on human diversity rather than European-ancestry genetics and microbiomes. A key component of increasing representation in genetics and microbiome studies is that these studies are designed as partnerships with minority and/or indigenous communities in a manner that builds both trust between the community and researchers, as well as facilitates the ability for the sample donors to exercise their rights on how data are treated and shared41. More

1.
Molinos, J. G. et al. Climate velocity and the future global redistribution of marine biodiversity. Nat. Clim. Change 6, 83–87 (2016).
Google Scholar 
2.
Pecl, G. T. et al. Biodiversity redistribution under climate change: impacts on ecosystems and human well-being. Science 355, eaai9214 (2017).
PubMed  Article  CAS  Google Scholar 

3.
Brook, B. W. & Alroy, J. Pattern, Process, Inference and Prediction in Extinction Biology (Royal Society, 2017).

4.
Barnosky, A. D. et al. Has the Earth’s sixth mass extinction already arrived? Nature 471, 51–57 (2011).
Article  CAS  Google Scholar 

5.
Thomas, C. D. et al. Extinction risk from climate change. Nature 427, 145–148 (2004).
Google Scholar 

6.
Urban, M. C. Accelerating extinction risk from climate change. Science 348, 571–573 (2015).
CAS  Article  Google Scholar 

7.
Collins, K. S., Edie, S. M., Hunt, G., Roy, K. & Jablonski, D. Extinction risk in extant marine species integrating palaeontological and biodistributional data. Proc. Biol. Sci. 285, 20181698 (2018).
PubMed  PubMed Central  Google Scholar 

8.
Finnegan, S. et al. Paleontological baselines for evaluating extinction risk in the modern oceans. Science 348, 567–570 (2015).
CAS  PubMed  Article  Google Scholar 

9.
van Woesik, R. et al. Hosts of the Plio-Pleistocene past reflect modern-day coral vulnerability. Proc. R. Soc. B 279, 2448–2456 (2012).
PubMed  Article  Google Scholar 

10.
Harnik, P. G. et al. Extinctions in ancient and modern seas. Trends Ecol. Evol. 27, 608–617 (2012).
PubMed  Article  Google Scholar 

11.
Kiessling, W. & Kocsis, Á. T. Adding fossil occupancy trajectories to the assessment of modern extinction risk. Biol. Lett. 12, 20150813 (2016).
PubMed  PubMed Central  Article  Google Scholar 

12.
Calosi, P., Putnam, H. M., Twitchett, R. J. & Vermandele, F. Marine metazoan modern mass extinction: improving predictions by integrating fossil, modern, and physiological data. Annu. Rev. Mar. Sci. 11, 369–390 (2019).
Article  Google Scholar 

13.
Wiens, J. J. & Graham, C. H. Niche conservatism: integrating evolution, ecology, and conservation biology. Annu. Rev. Ecol. Evol. Syst. 36, 519–539 (2005).
Article  Google Scholar 

14.
Beaugrand, G. Theoretical basis for predicting climate-induced abrupt shifts in the oceans. Philos. Trans. R. Soc. B 370, 20130264 (2015).
Article  Google Scholar 

15.
Lord, J. P., Barry, J. P. & Graves, D. Impact of climate change on direct and indirect species interactions. Mar. Ecol. Prog. Ser. 571, 1–11 (2017).
Article  Google Scholar 

16.
Bond, D. P. G. & Grasby, S. E. On the causes of mass extinctions. Palaeogeogr. Palaeoclimatol. Palaeoecol. 478, 3–29 (2017).
Article  Google Scholar 

17.
Penn, J. L., Deutsch, C., Payne, J. L. & Sperling, E. A. Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction. Science 362, eaat1327 (2018).
PubMed  Article  CAS  Google Scholar 

18.
Reddin, C. J., Nätscher, P. S., Kocsis, Á. T., Pörtner, H.-O. & Kiessling, W. Marine clade sensitivities to climate change conform across timescales. Nat. Clim. Change 10, 249–253 (2020).
Article  Google Scholar 

19.
Bolker, B. M. et al. Generalized linear mixed models: a practical guide for ecology and evolution. Trends Ecol. Evol. 24, 127–135 (2009).
PubMed  PubMed Central  Article  Google Scholar 

20.
Burnham, K. P. & Anderson, D. R. Model Selection and Multimodel Inference. A Practical Information-Theoretic Approach 2nd edn (Springer, 2010).

21.
McKinney, M. L. Extinction vulnerability and selectivity: combining ecological and paleontological views. Annu. Rev. Ecol. Syst. 28, 495–516 (1997).
Article  Google Scholar 

22.
Foote, M., Crampton, J. S., Beu, A. G. & Cooper, R. A. On the bidirectional relationship between geographic range and taxonomic duration. Paleobiology 34, 421–433 (2008).
Article  Google Scholar 

23.
Payne, J. L., Truebe, S., Nützel, A. & Chang, E. T. Local and global abundance associated with extinction risk in late Paleozoic and early Mesozoic gastropods. Paleobiology 37, 616–632 (2011).
Article  Google Scholar 

24.
Harnik, P. G. Direct and indirect effects of biological factors on extinction risk in fossil bivalves. Proc. Natl Acad. Sci. USA 108, 13594–13599 (2011).
CAS  PubMed  Article  Google Scholar 

25.
Harnik, P. G., Simpson, C. & Payne, J. L. Long-term differences in extinction risk among the seven forms of rarity. Proc. R. Soc. B 279, 4969–4976 (2012).
PubMed  Article  Google Scholar 

26.
Hopkins, M. J., Simpson, C. & Kiessling, W. Differential niche dynamics among major marine invertebrate clades. Ecol. Lett. 17, 314–323 (2014).
PubMed  Article  Google Scholar 

27.
Svenning, J. ‐C. & Skov, F. Limited filling of the potential range in European tree species. Ecol. Lett. 7, 565–573 (2004).
Article  Google Scholar 

28.
Normand, S. et al. Postglacial migration supplements climate in determining plant species ranges in Europe. Proc. R. Soc. B 278, 3644–3653 (2011).
PubMed  Article  Google Scholar 

29.
Stigall, A. L. When and how do species achieve niche stability over long time scales? Ecography 37, 1123–1132 (2014).
Google Scholar 

30.
Steinbauer, M. J. et al. Biogeographic ranges do not support niche theory in radiating Canary Island plant clades. Glob. Ecol. Biogeogr. 25, 792–804 (2016).
Article  Google Scholar 

31.
Foster, G. L., Hull, P., Lunt, D. J. & Zachos, J. C. Placing our Current ‘Hyperthermal’ in the Context of Rapid Climate Change in our Geological Past (The Royal Society Publishing, 2018).

32.
Leckie, R. M., Bralower, T. J. & Cashman, R. Oceanic anoxic events and plankton evolution: biotic response to tectonic forcing during the mid‐Cretaceous. Paleoceanography 17, 13-1–13-29 (2002).
Article  Google Scholar 

33.
Coxall, H. K. & Pearson, P. N. in Deep Time Perspectives on Climate Change: Marrying the Signal From Computer Models and Biological Proxies Vol. 2 (eds Williams, M. et al.) 351–387 (Geological Society of London, 2007).

34.
Ritterbush, K. A. & Foote, M. Association between geographic range and initial survival of Mesozoic marine animal genera: circumventing the confounding effects of temporal and taxonomic heterogeneity. Paleobiology 43, 209–223 (2017).
Article  Google Scholar 

35.
Stigall, A. L. Analysing links between biogeography, niche stability and speciation: the impact of complex feedbacks on macroevolutionary patterns. Palaeontology 56, 1225–1238 (2013).
Article  Google Scholar 

36.
BouDagher-Fadel, M. K. Evolution and Geological Significance of Larger Benthic Foraminifera (UCL Press, 2018).

37.
Reddin, C. J., Kocsis, Á. T. & Kiessling, W. Marine invertebrate migrations trace climate change over 450 million years. Glob. Ecol. Biogeogr. 27, 704–713 (2018).
Article  Google Scholar 

38.
Valentine, J. W. Temporal bias in extinctions among taxonomic categories. J. Paleontol. 48, 549–552 (1974).

39.
Kocsis, Á. T., Reddin, C. J., Alroy, J. & Kiessling, W. The r package divDyn for quantifying diversity dynamics using fossil sampling data. Methods Ecol. Evol. 10, 735–743 (2019).
Article  Google Scholar 

40.
Gradstein, F. M., Ogg, J. G., Schmitz, M. & Ogg, G. The Geologic Time Scale 2012 (Elsevier, 2012).

41.
Foote, M. Origination and extinction components of taxonomic diversity: general problems. Paleobiology 26, 74–102 (2000).
Article  Google Scholar 

42.
Veizer, J. & Prokoph, A. Temperatures and oxygen isotopic composition of Phanerozoic oceans. Earth Sci. Rev. 146, 92–104 (2015).
CAS  Article  Google Scholar 

43.
Song, H., Wignall, P. B., Song, H., Dai, X. & Chu, D. Seawater temperature and dissolved oxygen over the past 500 million years. J. Earth Sci. 30, 236–243 (2019).
Article  CAS  Google Scholar 

44.
Grossman, E. L. Applying oxygen isotope paleothermometry in deep time. Paleontol. Soc. Pap. 18, 39–68 (2012).
Article  Google Scholar 

45.
Henkes, G. A. et al. Temperature evolution and the oxygen isotope composition of Phanerozoic oceans from carbonate clumped isotope thermometry. Earth Planet. Sci. Lett. 490, 40–50 (2018).
CAS  Article  Google Scholar 

46.
Ryb, U. & Eiler, J. M. Oxygen isotope composition of the Phanerozoic ocean and a possible solution to the dolomite problem. Proc. Natl Acad. Sci. USA 115, 6602–6607 (2018).
CAS  PubMed  Article  Google Scholar 

47.
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2019).

48.
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Soft. https://doi.org/10.18637/jss.v067.i01 (2015).

49.
Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).

50.
Quené, H. & van den Bergh, H. Examples of mixed-effects modeling with crossed random effects and with binomial data. J. Mem. Lang. 59, 413–425 (2008).
Article  Google Scholar 

51.
Malik, W. A., Marco-Llorca, C., Berendzen, K. & Piepho, H.-P. Choice of link and variance function for generalized linear mixed models: a case study with binomial response in proteomics. Commun. Stat. Theory Methods 49, 1–20 (2019).

52.
Zuur, A., Ieno, E. N., Walker, N., Saveliev, A. A. & Smith, G. M. Mixed Effects Models and Extensions in Ecology with R (Springer Science & Business Media, 2009).

53.
Fox, G. A., Negrete-Yankelevich, S. & Sosa, V. J. Ecological Statistics: Contemporary Theory and Application (Oxford Univ. Press, 2015).

54.
Lieberman, M. D. & Cunningham, W. A. Type I and Type II error concerns in fMRI research: re-balancing the scale. Soc. Cogn. Affect. Neurosci. 4, 423–428 (2009).
PubMed  PubMed Central  Article  Google Scholar 

55.
Durbin, J. & Watson, G. S. Testing for serial correlation in least squares regression. Biometrika 58, 1–19 (1971).
Google Scholar 

56.
Nakagawa, S., Johnson, P. C. D. & Schielzeth, H. The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. J. R. Soc. Interface 14, 20170213 (2017).
PubMed  PubMed Central  Article  Google Scholar 

57.
IPCC. Special Report on Global Warming of 1.5 °C (eds Masson-Delmotte, V. et al.) (WMO, 2018). More

This is the first large scale analysis of feline desexing practices in Australia using outcomes documented in the patient medical record. The findings complement those of previous studies that used survey data to analyse the attitudes and opinions of veterinary professionals and owners to desexing31,34,38. The prevalence of desexing among cats in Australia, found to be 83.6%, confirms that desexing rates in Australia are among the highest reported internationally. Survey-based studies have reported that approximately 90% of cats in Australia are desexed, compared with 80% in the USA, and 43% in Italy39,40,41,42,43,44,45. A recent EPR-based study conducted in the UK reported the prevalence of feline desexing as 77%46. While population-wide analyses of desexing status provide a useful snapshot of practices in a region, most do not consider reproductive history.
A clear shift over time towards desexing cats at a younger age was evident here. EAD was 1.76 times more likely to have been carried out among desexed cats born between 2010 to 2017, than in those born between 1995 and 2009. This move towards earlier desexing was apparent in all age groups studied. Despite this trend, EAD had been carried out in only 21.5% of desexed females in the recent period. In fact, only 59.8% of females had been desexed by 6 months of age, which is the traditional recommendation and the most common recommendation reported by vets in Australia31,47. Despite a move towards earlier desexing, opportunities to control reproduction by prepubertal desexing are still being lost.
For an individual female cat, desexing at 6 months or later may be of little consequence, since they may not yet have reached puberty or had access to a mate. A recent survey of cat owners in Australia and New Zealand however found that 66% of cats had outdoor access45. From a population control perspective, eliminating the possibility of pregnancy by adopting EAD as standard has merit. The body of scientific evidence generated specifically to address the short-term and long-term safety of EAD overwhelmingly validates this practice4,17,22,23,24,25,26,27,28,29.
The impact that tighter control of reproduction among owned cats would have on shelter and stray populations is not yet clear. Populations of owned cats (completely reliant on humans) feral cats (living independently of humans) and stray cats (intermediate relationship with humans) do not exist in isolation48. Anthropogenic factors, including the provision of food, abandonment, and failure to curb reproduction, influence cat abundance and movement through these populations. Modelling population dynamics in owned, unowned (stray and feral) and shelter-housed cats holds promise to inform cat management strategies in the future49.
In multivariable models, for cats born 2010–2017, sex, breed, state and socioeconomic indices were all significantly associated with both desexing status and age at surgery. Females were less likely than males to be desexed and, among desexed cats, females were less likely than males to have been desexed at ≤ 4 months, supporting future measures to promote EAD in female cats. The reasons for this difference were not investigated but, conceivably, it may be due to higher fees for desexing females at some practices, or a greater awareness of spraying and roaming behaviours in males than pregancies in young female cats.
Not surprisingly purebred cats were less likely to be desexed than mixed breeds. In contrast, the finding that purebred cats were 2.7 times less likely to undergo EAD was unexpected because breeders commonly request EAD so that progeny for the pet market can be sold without delay50. It is plausible that this result reflects a greater demand in Australia for EAD from the charity and shelter sector, where mixed breed cats predominate, than from breeders. In line with this possibility, recent surveys found 70–80% of cats in Australia and New Zealand are of mixed breed and acquired from shelters, veterinary clinics, friends and as strays40,45,51. The higher odds of EAD in males than females was even greater among purebred cats, a result that may have been influenced by the practice of retaining more entire females than males for breeding.
The breeding season in Australia and New Zealand extends year round with peaks of kittening in spring and summer inferred from shelter admissions9,52,53. Cats born in winter had the lowest odds of being desexed in each age group. One explanation for this finding is that promotion of desexing by veterinary practices and welfare groups is less likely in winter because fewer kittens are born. This seasonal difference is certainly seen in the UK, where the RSPCA conducts desexing campaigns in Autumn to prevent the peak of spring litters54.
State or territory influenced both whether a cat was desexed, and the odds of EAD. Compared with cats in New South Wales (NSW), those in Victoria (VIC) and South Australia (SA) were more likely, and those in Queensland (QLD) less likely to be desexed. Again, compared with NSW, the odds of being desexed at ≤ 4 months were 1.45 greater for cats in VIC and 1.5–2.3 times less for those in QLD, SA and ACT. Desexing is handled inconsistently between Australia’s states and territories. Mandatory desexing legislation exists in ACT (by 3 months of age) and in SA, Tasmania (TAS), WA (by 6 months of age), with some exceptions. No legal requirement to desex cats exists at state level in VIC, QLD, NSW or Northern Territory (NT), although desexing is indirectly incentivized in NSW (by 6 months of age) and VIC (by 3 months of age) where registration is mandatory, and reduced registration fees are applied for desexed cats. No consistent relationship between our findings and state legislation related to desexing cats was identified. In fact, in ACT, where desexing of pet cats at 3 months of age has been a legal requirement since 2007, the second lowest odds of EAD were identified. Most veterinarians practicing in ACT (90%), surveyed 10 years after the legislation was introduced, gave recommendations inconsistent with the legislation and 35% were unaware that desexing by 3 months was mandatory in the ACT34. Whether and how legislation might be an appropriate tool to influence reproduction in owned cats and, indirectly, overpopulation should be further investigated.
Socioeconomic conditions influenced both whether a cat was desexed or not, and the age at desexing. Entire cats were more common in remote, low income and disadvantaged areas. This finding is concerning, given that outdoor access was more likely in non-urban than urban areas in a study of households in Australia and New Zealand45, implying more opportunity to find a mate. In addition, stray cat density correlated positively with socioeconomic deprivation in a New Zealand-based study employing geographically weighted regression analyses55. Together, these findings support the promotion of desexing campaigns in non-urban areas.
Economic indicators such as household income influenced whether a cat was desexed; the odds of being desexed were around 1.2 times greater in the highest compared to the lowest income areas. A similar, but more dramatic effect was reported in a study conducted in the USA where the prevalence of desexing increased from 51.4% to 96.2% as household income increased43. Among desexed cats, EAD was least likely in low income areas, but highest in the most socio-economically disadvantaged areas. Although this might seem paradoxical, IRSD is based on broader indicators of disadvantage than income alone. A UK study, similarly, identified that EAD was most likely in the most deprived regions, and that chances of being desexed by 6 months were more likely in higher income areas56. Possible explanations for these observations include the preferential targeting of areas of greatest disadvantage, rather than those with fewer economic resources, by discount desexing programs promoting EAD, or preferential sourcing of kittens in disadvantaged areas from organizations that routinely practice EAD, such as shelters.
There are limitations to our study that should be considered when interpreting the results. Cats that were either not registered with a veterinary practice, or were registered with a practice that did not contribute to VCA during the study periods were not studied. Therefore actual desexing prevalences are almost certainly lower than the estimates reported here. The study population represents cats that are accessible for desexing and is expected to comprise cats kept as pets, for breeding, owned by shelters, semi-owned cats and others. Provenance and lifestyle were not investigated because we chose not to collect data from the examination text field in VCA because of its low positive predictive value57, and because these data are inconsistently recorded. This precluded the analysis of other variables that may have been related to desexing outcomes such access to outdoors and the number and species of pets. Data collection was not uniform across Australia and variations in sample size, for example between states, may have affected our results. Also it is possible that data for the same cat presenting at more than one practice could be counted more than once, although a previous study using VCA found that  More

1.
Chen M, Zeng G, Zhang J, Xu P, Chen A, Lu L. Global landscape of total organic carbon, nitrogen and phosphorus in lake water. Sci Rep. 2015;5:15043.
CAS  PubMed  PubMed Central  Article  Google Scholar 
2.
Toming K, Kotta J, Uuemaa E, Sobek S, Kutser T, Tranvik LJ. Predicting lake dissolved organic carbon at a global scale. Sci Rep. 2020;10:8471.
CAS  PubMed  PubMed Central  Article  Google Scholar 

3.
Tranvik LJ, Downing JA, Cotner JB, Loiselle SA, Striegl RG, Ballatore TJ, et al. Lakes and reservoirs as regulators of carbon cycling and climate. Limnol Oceanogr. 2009;54:2298–314.
CAS  Article  Google Scholar 

4.
Tranvik LJ, Cole JJ, Prairie YT. The study of carbon in inland waters-from isolated ecosystems to players in the global carbon cycle. Limnol Oceanogr Lett. 2018;3:41–48.
Article  Google Scholar 

5.
Jaffé R, McKnight D, Maie N, Cory R, McDowell WH, Campbell JL. Spatial and temporal variations in DOM composition in ecosystems: the importance of long-term monitoring of optical properties. J Geophys Res Biogeosci. 2008;113:G04032.
Article  CAS  Google Scholar 

6.
Crump BC, Kling GW, Bahr M, Hobbie JE. Bacterioplankton community shifts in an Arctic Lake correlate with seasonal changes in organic matter source. Appl Environ Microbiol. 2003;69:2253–68.
PubMed  PubMed Central  Article  Google Scholar 

7.
Fasching C, Behounek B, Singer GA, Battin TJ. Microbial degradation of terrigenous dissolved organic matter and potential consequences for carbon cycling in brown-water streams. Sci Rep. 2014;4:4981.
CAS  PubMed  PubMed Central  Article  Google Scholar 

8.
Yakimovich KM, Emilson EJS, Carson MA, Tanentzap AJ, Basiliko N, Mykytczuk NCS. Plant litter type dictates microbial communities responsible for greenhouse gas production in amended lake sediments. Front Microbiol. 2018;9:2662.
PubMed  PubMed Central  Article  Google Scholar 

9.
Attermeyer K, Hornick T, Kayler ZE, Bahr A, Zwirnmann E, Grossart HP, et al. Enhanced bacterial decomposition with increasing addition of autochthonous to allochthonous carbon without any effect on bacterial community composition. Biogeosciences. 2014;11:1479–89.
Article  Google Scholar 

10.
Fabian J, Zlatanovic S, Mutz M, Premke K. Fungal-bacterial dynamics and their contribution to terrigenous carbon turnover in relation to organic matter quality. ISME J. 2017;11:415–25.
CAS  PubMed  Article  PubMed Central  Google Scholar 

11.
Quigley LNM, Edwards A, Steen AD, Buchan A. Characterization of the interactive effects of labile and recalcitrant organic matter on microbial growth and metabolism. Front Microbiol. 2019;10:493.
PubMed  PubMed Central  Article  Google Scholar 

12.
Tranvik LJ. Degradation of dissolved organic matter in humic waters by bacteria. In: Hessen DOTLJ, editor. Aquatic Humic Substances. Berlin, Heidelberg: Springer; 1998.

13.
Søndergaard M, Borch NH, Riemann B. Dynamics of biodegradable DOC produced by freshwater plankton communities. Aquat Micro Ecol. 2000;23:73–83.
Article  Google Scholar 

14.
Berg B, McClaugherty C. Initial litter chemical composition. Plant Litter. 2014;3:53–83.
Article  Google Scholar 

15.
Bugg TD, Ahmad M, Hardiman EM, Rahmanpour R. Pathways for degradation of lignin in bacteria and fungi. Nat Prod Rep. 2011;28:1883–96.
CAS  PubMed  Article  PubMed Central  Google Scholar 

16.
Klotzbücher T, Kaiser K, Guggenberger G, Gatzek C, Kalbitz K. A new conceptual model for the fate of lignin in decomposing plant litter. Ecology. 2011;92:1052–62.
PubMed  Article  PubMed Central  Google Scholar 

17.
Grey J, Jones RI, Sleep D. Seasonal changes in the importance of the source of organic matter to the diet of zooplankton in Loch Ness, as indicated by stable isotope analysis. Limnol Oceanogr. 2001;46:505–13.
Article  Google Scholar 

18.
Cole JJ, Carpenter SR, Kitchell JF, Pace ML. Pathways of organic carbon utilization in small lakes: results from a whole-lake 13C addition and coupled model. Limnol Oceanogr. 2002;47:1664–75.
CAS  Article  Google Scholar 

19.
Guenet B, Danger M, Abbadie L, Lacroix G. Priming effect: bridging the gap between terrestrial and aquatic ecology. Ecology. 2010;91:2850–61.
PubMed  Article  PubMed Central  Google Scholar 

20.
Bianchi TS. The role of terrestrially derived organic carbon in the coastal ocean: a changing paradigm and the priming effect. Proc Natl Acad Sci USA. 2011;108:19473–81.
CAS  PubMed  Article  PubMed Central  Google Scholar 

21.
Bengtsson MM, Attermeyer K, Catalán N. Interactive effects on organic matter processing from soils to the ocean: are priming effects relevant in aquatic ecosystems? Hydrobiologia. 2018;822:1–17.
CAS  Article  Google Scholar 

22.
Kuzyakov Y, Friedel JK, Stahr K. Review of mechanisms and quantification of priming effects. Soil Biol Biochem. 2000;32:1485–98.
CAS  Article  Google Scholar 

23.
Bianchi TS, Ward ND. Editorial: the role of priming in terrestrial and aquatic ecosystems. Front Earth Sci. 2019;7:321.
Article  Google Scholar 

24.
Halvorson HM, Francoeur SN, Findlay RH, Kuehn KA. Algal-mediated priming effects on the ecological stoichiometry of leaf litter decomposition: a meta-analysis. Front Earth Sci. 2019;7:76.
Article  Google Scholar 

25.
Kayler ZE, Premke K, Gessler A, Gessner MO, Griebler C, Hilt S, et al. Integrating aquatic and terrestrial perspectives to improve insights into organic matter cycling at the landscape scale. Front Earth Sci. 2019;7:127.
Article  Google Scholar 

26.
Danger M, Cornut J, Chauvet E, Chavez P, Elger A, Lecerf A. Benthic algae stimulate leaf litter decomposition in detritus-based headwater streams: a case of aquatic priming effect? Ecology. 2013;94:1604–13.
PubMed  Article  PubMed Central  Google Scholar 

27.
Guenet B, Danger M, Harrault L, Allard B, Jauset-Alcala M, Bardoux G, et al. Fast mineralization of land-born C in inland waters: first experimental evidences of aquatic priming effect. Hydrobiologia. 2013;721:35–44.
Article  CAS  Google Scholar 

28.
Ward ND, Sawakuchi HO, Richey JE, Keil RG, Bianchi TS. Enhanced aquatic respiration associated with mixing of clearwater tributary and turbid Amazon river waters. Front Earth Sci. 2019;7:101.
Article  Google Scholar 

29.
Bianchi TS, Thornton DCO, Yvon-Lewis SA, King GM, Eglinton TI, Shields MR, et al. Positive priming of terrestrially derived dissolved organic matter in a freshwater microcosm system. Geophys Res Lett. 2015;42:5460–67.
CAS  Article  Google Scholar 

30.
Tang KW, Gladyshev MI, Dubovskaya OP, Kirillin G, Grossart H-P. Zooplankton carcasses and non-predatory mortality in freshwater and inland sea environments. J Plankton Res. 2014;36:597–612.
CAS  Article  Google Scholar 

31.
Cauchie HM, Jaspar-Versali MF, Hoffmann L, Thomé JP. Analysis of the seasonal variation in biochemical composition of Daphnia magna Straus (Crustacea: Branchiopoda: Anomopoda) from an aerated wastewater stabilisation pond. Ann Limnol. 1999;35:223–31.
Article  Google Scholar 

32.
Smirnov NN. Physiology of the Cladocera, 2nd ed., London, Academic Press: Elsevier; 2017.

33.
Dubovskaya OP, Tang KW, Gladyshev MI, Kirillin G, Buseva Z, Kasprzak P, et al. Estimating in situ zooplankton non-predation mortality in an oligo-mesotrophic lake from sediment trap data: caveats and reality check. PLoS ONE. 2015;10:e0131431.
PubMed  PubMed Central  Article  CAS  Google Scholar 

34.
Kirillin G, Grossart H-P, Tang KW. Modeling sinking rate of zooplankton carcasses: effects of stratification and mixing. Limnol Oceanogr. 2012;57:881–94.
Article  Google Scholar 

35.
Tang KW, Hutalle KML, Grossart HP. Microbial abundance, composition and enzymatic activity during decomposition of copepod carcasses. Aquat Micro Ecol. 2006;45:219–27.
Article  Google Scholar 

36.
Tang KW, Bickel SL, Dziallas C, Grossart HP. Microbial activities accompanying decomposition of cladoceran and copepod carcasses under different environmental conditions. Aquat Micro Ecol. 2009;57:89–100.
Article  Google Scholar 

37.
Kolmakova OV, Gladyshev MI, Fonvielle JA, Ganzert L, Hornick T, Grossart HP. Effects of zooplankton carcasses degradation on freshwater bacterial community composition and implications for carbon cycling. Environ Microbiol. 2019;21:34–49.
CAS  PubMed  Article  PubMed Central  Google Scholar 

38.
Corno G, Salka I, Pohlmann K, Hall AR, Grossart HP. Interspecific interactions drive chitin and cellulose degradation by aquatic microorganisms. Aquat Micro Ecol. 2015;76:27–37.
Article  Google Scholar 

39.
Masigol H, Khodaparast SA, Woodhouse JN, Rojas‐Jimenez K, Fonvielle J, Rezakhani F, et al. The contrasting roles of aquatic fungi and oomycetes in the degradation and transformation of polymeric organic matter. Limnol Oceanogr. 2019;64:2662–78.
CAS  Article  Google Scholar 

40.
Gessner MO, Chauvet E. Importance of stream microfungi in controlling breakdown rates of leaf litter. Ecology. 1994;75:1807–17.
Article  Google Scholar 

41.
Grossart HP, Van den Wyngaert S, Kagami M, Wurzbacher C, Cunliffe M, Rojas-Jimenez K. Fungi in aquatic ecosystems. Nat Rev Microbiol. 2019;17:339–54.
CAS  PubMed  Article  PubMed Central  Google Scholar 

42.
Osono T. Functional diversity of ligninolytic fungi associated with leaf litter decomposition. Ecol Res. 2019;35:30–43.
Article  CAS  Google Scholar 

43.
Taube R, Ganzert L, Grossart HP, Gleixner G, Premke K. Organic matter quality structures benthic fatty acid patterns and the abundance of fungi and bacteria in temperate lakes. Sci Total Environ. 2018;610-611:469–81.
CAS  PubMed  Article  PubMed Central  Google Scholar 

44.
Cragg SM, Beckham GT, Bruce NC, Bugg TD, Distel DL, Dupree P, et al. Lignocellulose degradation mechanisms across the Tree of Life. Curr Opin Chem Biol. 2015;29:108–19.
CAS  PubMed  PubMed Central  Article  Google Scholar 

45.
Wilhelm RC, Singh R, Eltis LD, Mohn WW. Bacterial contributions to delignification and lignocellulose degradation in forest soils with metagenomic and quantitative stable isotope probing. ISME J. 2019;13:413–29.
CAS  PubMed  Article  PubMed Central  Google Scholar 

46.
Kuikman PJ, Jansen AG, van Veen JA, Zehnder AJB. Protozoan predation and the turnover of soil organic carbon and nitrogen in the presence of plants. Biol Fertil Soils. 1990;10:22–28.
CAS  Article  Google Scholar 

47.
White DC, Davis WM, Nickels JS, King JD, Bobbie RJ. Determination of the sedimentary microbial biomass by extractible lipid phosphate. Oecologia. 1979;40:51–62.
CAS  PubMed  Article  PubMed Central  Google Scholar 

48.
Romaní AM, Fischer H, Mille-Lindblom C, Tranvik LJ. Interactions of bacteria and fungi on decomposing litter: differential extracellular enzyme activities. Ecology. 2006;87:2559–69.
PubMed  Article  PubMed Central  Google Scholar 

49.
Hutalle-Schmelzer KM, Zwirnmann E, Kruger A, Grossart HP. Changes in pelagic bacteria communities due to leaf litter addition. Micro Ecol. 2010;60:462–75.
Article  Google Scholar 

50.
Smith EJ, Davison W, Hamilton-Taylor J. Methods for preparing synthetic freshwaters. Water Res. 2002;36:1286–96.
CAS  PubMed  Article  PubMed Central  Google Scholar 

51.
Attermeyer K, Premke K, Hornick T, Hilt S, Grossart HP. Ecosystem-level studies of terrestrial carbon reveal contrasting bacterial metabolism in different aquatic habitats. Ecology. 2013;94:2754–66.
PubMed  Article  PubMed Central  Google Scholar 

52.
Halbedel S (2015) Protocol for CO2 sampling in waters by the use of the headspaceequilibration technique, based on the simple gas equation; second update. Protoc Exch. https://assets.researchsquare.com/files/nprot-4275/v1/manuscript.pdf

53.
Cheng W. Measurement of rhizosphere respiration and organic matter decomposition using natural 13C. Plant Soil. 1996;183:263–68.
CAS  Article  Google Scholar 

54.
Taube R, Fabian J, Van den Wyngaert S, Agha R, Baschien C, Gerphagnon M, et al. Potentials and limitations of quantification of fungi in freshwater environments based on PLFA profiles. Fungal Ecol. 2019;41:256–68.
Article  Google Scholar 

55.
Zhang Z, Qu Y, Li S, Feng K, Wang S, Cai W, et al. Soil bacterial quantification approaches coupling with relative abundances reflecting the changes of taxa. Sci Rep. 2017;7:4837.
PubMed  PubMed Central  Article  CAS  Google Scholar 

56.
Mangelsdorf K, Karger C, Zink K-G. Phospholipids as life markers in geological habitats. Hydrocarbons, oils and lipids: diversity, origin, chemistry and fate. 2019. pp. 1–29.

57.
Frostegård Å, Tunlid A, Bååth E. Microbial biomass measured as total lipid phosphate in soils of different organic content. J Microbiol Methods. 1991;14:151–63.
Article  Google Scholar 

58.
Nercessian O, Noyes E, Kalyuzhnaya MG, Lidstrom ME, Chistoserdova L. Bacterial populations active in metabolism of C1 compounds in the sediment of Lake Washington, a freshwater lake. Appl Environ Microbiol. 2005;71:6885–99.
CAS  PubMed  PubMed Central  Article  Google Scholar 

59.
Klindworth A, Pruesse E, Schweer T, Peplies J, Quast C, Horn M, et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 2013;41:e1.
CAS  PubMed  Article  PubMed Central  Google Scholar 

60.
Tedersoo L, Anslan S, Bahram M, Põlme S, Riit T, Liiv I, et al. Shotgun metagenomes and multiple primer pair-barcode combinations of amplicons reveal biases in metabarcoding analyses of fungi. MycoKeys. 2015;10:1–43.
Article  Google Scholar 

61.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJ, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.
CAS  PubMed  PubMed Central  Article  Google Scholar 

62.
Murali A, Bhargava A, Wright ES. IDTAXA: a novel approach for accurate taxonomic classification of microbiome sequences. Microbiome. 2018;6:140.
PubMed  PubMed Central  Article  Google Scholar 

63.
De Caceres M, Legendre P. Associations between species and groups of sites: indices and statistical inference. Ecology. 2009;90:3566–74.
PubMed  Article  PubMed Central  Google Scholar 

64.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinforma. 2008;9:559.
Article  CAS  Google Scholar 

65.
Taipale SJ, Kainz MJ, Brett MT. Diet-switching experiments show rapid accumulation and preferential retention of highly unsaturated fatty acids in Daphnia. Oikos. 2011;120:1674–82.
Article  Google Scholar 

66.
Corno G, Jurgens K. Structural and functional patterns of bacterial communities in response to protist predation along an experimental productivity gradient. Environ Microbiol. 2008;10:2857–71.
PubMed  Article  PubMed Central  Google Scholar 

67.
Fierer N, Lauber CL, Ramirez KS, Zaneveld J, Bradford MA, Knight R. Comparative metagenomic, phylogenetic and physiological analyses of soil microbial communities across nitrogen gradients. ISME J. 2012;6:1007–17.
CAS  PubMed  Article  PubMed Central  Google Scholar 

68.
Song HK, Song W, Kim M, Tripathi BM, Kim H, Jablonski P, et al. Bacterial strategies along nutrient and time gradients, revealed by metagenomic analysis of laboratory microcosms. FEMS Microbiol Ecol. 2017;93:fix114.
Article  CAS  Google Scholar 

69.
Bardgett RD, Kandeler E, Tscherko D, Hobbs PJ, Bezemer TM, Jones TH, et al. Below-ground microbial community development in a high temperature world. Oikos. 1999;85:193–203.
Article  Google Scholar 

70.
Hammel KE, Kapich AN, Jensen KA, Ryan ZC. Reactive oxygen species as agents of wood decay by fungi. Enzym Micro Technol. 2002;30:445–53.
CAS  Article  Google Scholar 

71.
Rojas-Jimenez K, Fonvielle JA, Ma H, Grossart H-P. Transformation of humic substances by the freshwater Ascomycete Cladosporium sp. Limnol Oceanogr. 2017;62:1955–62.
CAS  Article  Google Scholar 

72.
Nierman WC, Feldblyum TV, Laub MT, Paulsen IT, Nelson KE, Eisen JA, et al. Complete genome sequence of Caulobacter crescentus. Proc Natl Acad Sci USA. 2001;98:4136–41.
CAS  PubMed  Article  PubMed Central  Google Scholar 

73.
Zheng W, Lehmann A, Ryo M, Valyi KK, Rillig MC. Growth rate trades off with enzymatic investment in soil filamentous fungi. Sci Rep. 2020;10:11013.
CAS  PubMed  PubMed Central  Article  Google Scholar 

74.
Bärlocher F, Boddy L. Aquatic fungal ecology—How does it differ from terrestrial? Fungal Ecol. 2016;19:5–13.
Article  Google Scholar 

75.
Lange L, Barrett K, Pilgaard B, Gleason F, Tsang A. Enzymes of early-diverging, zoosporic fungi. Appl Microbiol Biotechnol. 2019;103:6885–902.
CAS  PubMed  PubMed Central  Article  Google Scholar 

76.
Janusz G, Pawlik A, Swiderska-Burek U, Polak J, Sulej J, Jarosz-Wilkolazka A, et al. Laccase properties, physiological functions, and evolution. Int J Mol Sci. 2020;21:966.
CAS  PubMed Central  Article  Google Scholar 

77.
Catalán N, Kellerman AM, Peter H, Carmona F, Tranvik LJ. Absence of a priming effect on dissolved organic carbon degradation in lake water. Limnol Oceanogr. 2015;60:159–68.
Article  Google Scholar 

78.
Bengtsson MM, Wagner K, Burns NR, Herberg ER, Wanek W, Kaplan LA, et al. No evidence of aquatic priming effects in hyporheic zone microcosms. Sci Rep. 2014;4:5187.
CAS  PubMed  PubMed Central  Article  Google Scholar 

79.
Tanentzap AJ, Fitch A, Orland C, Emilson EJS, Yakimovich KM, Osterholz H, et al. Chemical and microbial diversity covary in fresh water to influence ecosystem functioning. Proc Natl Acad Sci USA. 2019;116:24689–95.
CAS  PubMed  Article  PubMed Central  Google Scholar 

80.
Orland C, Emilson EJS, Basiliko N, Mykytczuk NCS, Gunn JM, Tanentzap AJ. Microbiome functioning depends on individual and interactive effects of the environment and community structure. ISME J. 2019;13:1–11.
CAS  PubMed  Article  PubMed Central  Google Scholar 

81.
Winder M, Sommer U. Phytoplankton response to a changing climate. Hydrobiologia. 2012;698:5–16.
Article  Google Scholar 

82.
Pothoven SA, Fahnenstiel GL. Spatial and temporal trends in zooplankton assemblages along a nearshore to offshore transect in southeastern Lake Michigan from 2007 to 2012. J Gt Lakes Res. 2015;41:95–103.
Article  Google Scholar 

83.
Selmeczy GB, Abonyi A, Krienitz L, Kasprzak P, Casper P, Telcs A, et al. Old sins have long shadows: climate change weakens efficiency of trophic coupling of phyto- and zooplankton in a deep oligo-mesotrophic lowland lake (Stechlin, Germany)—a causality analysis. Hydrobiologia. 2018;831:101–17.
Article  CAS  Google Scholar  More

This study comprehensively evaluated the risk associations between cause-specific ambulance services, extreme temperatures, and mass concentrations of PM2.5 and its constituents. The significant cold effects on chest pain and headache/dizziness/vertigo/fainting/syncope and heat effects on coma and unconsciousness and lying at public were observed, while the risk of ambulance services of OHCA was elevated in both extreme heat and cold environments. Ambulance services of respiratory distress, lying at public, and OHCA increased as the PM2.5 concentration increased, and the risk was significant at the PM2.5 concentration of 20–60 60 μg/m3 for ambulance services of lying at public and higher than 60 μg/m3 for respiratory distress. After controlling for effects of daily average temperature and PM2.5 concentration, this study still identified the significant effects of sulfate and EC on ambulance services of lying at public and OC on headache/dizziness/vertigo/fainting/syncope as the concentrations of PM2.5 constituents were at 90th percentile.
Limited studies assessed associations between ambulance calls and ambient environment9,13,19,22,23,24,25,26,27. Studies in Emilia-Romagna in Italy23, Brisbane in Australia26, Taiwan19, and Huainan and Luoyang in China22,24, have indicated the numbers of ambulance calls associated with extreme heat; the risks generally increase as the daily temperature exceeds 27 °C19,23,26. However, no consistent finding for cold threshold was identified19,28. Kaohsiung City has a tropical climate (daily temperature ranging from 13.5 °C to 31.5 °C), but it is cooler than cities located near the equator, e.g., Singapore and Manila. Except for ambulance service of OHCA, we found that the significant risks associated with temperature were only identified in environments with extreme temperatures ( 90th percentiles; Fig. 3).
Fine particulate matter (PM2.5) are characterized with a small diameter ( More

Study area
The study area is a one-hectare (ha) subtropical mixed young forest in the age of 20 to 30 years of Hong Kong. It is located at Shek Kong (22.428774, 114.114968), in between the Kadoorie Farm and Botanic Garden (KFBG) and Kadoorie Institute, The University of Hong Kong (HKU) as shown in Fig. 1. This is a forest dynamic plot under ‘The Center for Tropical Forest Science–Forest Global Earth Observatory (CTFS-ForestGEO) (https://forestgeo.si.edu/). It acts as a research base on the forest dynamics, forest biodiversity, carbon sequestration and more, which also provide opportunity for public involvement in scientific research.
Figure 1

The study area: 1 hectare subtropical moist young forest plot, a full demonstration plot of the Forest Global Earth observatory (ForestGEO) project (https://forestgeo.si.edu/), located in Hong Kong in understanding the long-term forest dynamic. The map is generated by Authors using ArcMap version 10.5 (https://desktop.arcgis.com/en/arcmap/).

Full size image

The one-hectare forest plot was demarcated by 25 quadrats in 20 m × 20 m each, which were further sub-divided into sixteen 5 m × 5 m sub-quadrats; delineated with permanent marker poles by the professional surveying team. The forest survey was launched on January 11th, 2012 and completed on September 6th, 2012. A total of 63 species, 10,442 individual trees with 20,888 stems were recorded in the site. The stem locations (UTM WGS84 coordinate system) were recorded in nearest 5 cm, and DBH were measured at 1.3 m at breast height in nearest 1 mm. The dominant species are Litsea rotundifolia, Psychotria asiatica, Ilex asprella and Aporosa dioica, which all are native species and accounted for over 80% of stems with the study area. The forest condition of the study area is shown in Table 2.
Table 2 The forest condition and descriptive statistics.
Full size table

Methods
A three-stage methodological framework was outlined for this study as shown in Figs. 2 and 3. In the first stage, parameters including DBH, wood density and stem location were recorded for field-measured AGB computation. DBH was directly obtained from site, wood density (in g cm−3) was obtained from the global wood density database43,44 or World Agroforestry database45 (http://db.worldagroforestry.org//wd), up to species level. If the species was not recorded in the database, the value was replaced by its genus averaged wood density value. Tree height was derived from the LiDAR data with reference to the recorded x,y location of the stems. The second stage is LiDAR AGB model derivation. Five allometric models (Model 1 to Model 5) were used to estimate the AGB of individual trees based on field-measured parameters. The allometric model with the lowest model error was selected to compute the ‘Field-measured AGB’, which would be used as the dependent variable to develop the LiDAR-derived AGB model. The LiDAR plots metrics were generated in various plot-size (i.e., 10 m radius, 5 m radius and 2.5 m radius) within the study area. Stepwise linear regression was used to select the important and significant predictors in each regression model. Three different model forms were tested (Model I, II, III as discussed in “Stage 2: LiDAR AGB model derivation” section), The LiDAR derived AGB regression models would be evaluated by means of assumption tests and bootstrapping; as well as cross-validation (CV) in the last stage of model evaluation and validation.
Figure 2

The workflow of the study. Stage 1: allometric modeling to compute ground-truth AGB. Five allometric models were compared and the best one is selected to calibrate model in the next stage. Stage 2: LiDAR AGB Model Derivation, by comparing LiDAR plot metrics derived in different plot size (i.e. 10 m, 5 m, 2.5 m radius). Stage 3: Model Evaluation & Validation on the best model selected from Stage 2 (Source: Authors).

Full size image

Figure 3

The graphical abstract of the study (Source: Authors).

Full size image

Stage 1: allometric modeling
20,888 stems from 63 species, were planned to be used to develop the allometric models. Amongst, the 81% of the wood density value were measured up to species level (i.e., 51 species) and 19% were up to genus level (i.e., 12 species). However, as the tree height parameter was retrieved from the LiDAR 1 m Canopy Height Model (CHM), 263 stems found no value and thus the sample size reduced to 20,625 stems. The DBH of stems ranged from 1.0 cm to 57.1 cm, with mean of 2.55 cm and S.D. of 2.38 cm.
The selected five allometric models included two pantropical models by Chave et al.16; two subtropical models by Xu et al.46 and one local model by Nichol and Sarker47 were assessed. The inclusion of the subtropical models from Xu et al.46 is because of the similar forest type (i.e., subtropical moist forest from southern China) with our study area; whereas the pantropical models from Chave et al.16 is due to its large database across the pantropical regions yet not being explored adequately its applicability in subtropical forests.
The five allometric models were compared by One-way ANOVA and the Tukey HSD post-hoc analysis48. Model 1 (Eq. (3)) was the best model proposed by Chave et al.16, comprised of wood density (ρ), diameter-at-breast height (DBH) and height (H) as shown below (AIC = 3130):

$$AGB=0.0673times {left(rho {(DBH)}^{2}Hright)}^{0.976}$$
(3)

However, it is very unlikely to have a precise tree height data in a closed canopy49. Therefore, an alternative allometric equation without height was then proposed by Chave et al.16 and adopted as Model 2 (Eq. (4)) in our study (AIC =  − 4293):to represent the subtropical AGB

$$AGB=expleft[-1.803-0.976E+0.976text{ln}left(rho right)+2.673text{ln}left(DBHright)-0.0299left[{(text{ln}(DBH))}^{2}right]right]$$
(4)

The bioclimatic variable, “E”, which made up of three parameters to account for climatic variations: Temperature seasonality (TS), Long-term Maximum Climatological Water Deficit (CWD) and Precipitation Seasonality (PS). The formula of the bioclimatic variable, “E”, is shown below (Eq. (5)):

$$E={left(0.178 times TS-0.938 times CWD-6.61 times PSright)}^{{10}^{-3}}$$
(5)

The monthly mean temperature (℃), precipitation (mm) and evapotranspiration (mm) data were obtained from the Hong Kong Observatory50 in the period of 1981 to 2010 to compute the three variables. E was then computed as 0.261 and was input into Eq. (4) to derive the ground-truth AGB of individual stems.
The AGB models developed by Xu et al.46 were conducted in the subtropical mixed-species moist forest in southern part of China. Two allometric models from the study were selected to represent the subtropical AGB model. Model 3 (Eq. (6)) assumed tree height is available (by extracting from LiDAR):

$$AGB=text{exp}(-2.334+2.118text{ln}(DBH)+0.5436text{ln}(H)+0.5953text{ln}(rho ))$$
(6)

Meanwhile, an alternative allometric model would be used by assuming tree height was not available, which is the Model 4 (Eq. (7)) of this study:

$$AGB=text{exp}(-1.8226+2.4105text{ln}left(DBHright)+0.5781text{ln}(rho ))$$
(7)

Nichol and Sarker47 developed a local allometric model for Hong Kong by harvesting 75 trees from 15 dominant species of Hong Kong. DBH and H within 50 circular sample plots from a variety of tree stands were measured, which were then used to establish allometric model with field measured AGB. The best allometric model was a model using DBH as the sole parameter and thus selected as the Model 5 (Eq. (8)) of this study:

$$AGB=text{exp}(-1.8226+2.4105text{ln}left(DBHright)+0.5781text{ln}(rho ))$$
(8)

Stage 2: LiDAR AGB model derivation
The airborne LiDAR data used in this study was captured by Optech Gemini ALTM Airborne Laser Terrain Mapper and acquired by The Government of the Hong Kong Special Administrative Region from December 1st, 2010 to January 8th, 2011. A total of 5575 ground truth points were generated, with horizontal accuracy of 0.294 m (95% confidence interval (C.I.)). Vertical accuracy was also assessed against the orthometric heights of Hong Kong Principal Datum (HKPD), the average vertical accuracy was 0.1 m (95% C.I.). Multiple returns were recorded per pulse up to four range measurements (i.e., first, second third and last). The LiDAR acquisition parameters are shown in Table 3.
Table 3 The LiDAR acquisition parameters.
Full size table

Prior to generation of the plot metrics, the LiDAR data was pre-processed by creating the ground TIN by extracting only the ground returns. The extraction and processing were performed in ArcMap 10.5 with LAStools extension (https://rapidlasso.com/lastools/) and the FUSION 3.7 (http://forsys.cfr.washington.edu/fusion/fusionlatest.html). The canopy surface was defined using all non-ground returns with height above 2.0 m. The reason of choosing 2.0 m as the threshold was to avoid canopy returns to be mixed with ground returns51. Moreover, ‘all returns’ were used instead of the ‘first returns’, since the former provided more information on the lower canopies or understory51, while over 50% of returned points in our study area were classified as medium or low vegetation. The ground TIN was subtracted from the canopy surface to compute the normalized height value of each canopy point. The LiDAR metrics were then derived from the normalized height point cloud.
To be compatible with the LiDAR dataset and to facilitate the plot delineation procedure, the entire 1 ha study area was divided into circular plots with three plot sizes: (1) twenty-five 10 m radius plots (16,182 stems), (2) one-hundred 5 m radius plot (15,538 stems) and (3) four-hundred 2.5 m radius plots (15,100 stems) as shown in Fig. 4a–d. Circular plots were considered more favorable than rectangular or square plots, since the periphery-to-area ratio was the smallest and thus minimized the number of edge trees52.
Figure 4

The (a) 1 ha rectangular plot clipped into (b) twenty-five 10 m radius circular plots; (c) one-hundred 5 m radius plots; and (d) four-hundred 2.5 m radius plots respectively; for the LiDAR metrics derivation. The maps are generated by Authors using ArcMap version 10.5 (https://desktop.arcgis.com/en/arcmap/).

Full size image

The 59 plot metrics, under the descriptive, height, intensity and canopy cover categories, were derived by the ‘cloudmetrics’ function in FUSION version 3.7. The LiDAR metrics, and its log-transformed metrics were input into stepwise regression model as independent predictors of AGB. Significant predictors were selected (F  1.0 were removed) into the regression model. Three sets of regression models were generated (Table 4) and the allometric models tended to be linear and normal after logarithmic transformation18,19.
Table 4 The input variables for the three regression models.
Full size table

Observing the normal Q–Q plot (Fig. 5a) of the dependent variable (i.e., AGB) and the Kolmogorov–Smirnov statistic (d = 0.216, p  0.200) indicated a normal distribution (Fig. 5b) and which became the Model II. The further log-transformation into Model III did not indicate significant improvement and thus Model II in various plot sizes were to be further explored.
Figure 5

The normal Q–Q plot, in 10 m radius plot size, of (a) raw AGB and (b) log-AGB. After logarithmic transformation, the AGB (dependent variable) was normalized.

Full size image

Assumption tests including test on normality, homoscedasticity and absence of multi-collinearity were conducted. The Normal P–P plot, scatter plots on residuals, ‘Tolerance’ index and Variance Inflation Factor (VIF) were applied to detect the violation of assumptions. If the ‘Tolerance’ index is smaller than 0.1, or VIF is greater than 10, it indicates a significance chance of collinearity53.
Stage 3: model evaluation and validation
AGB regression models were to be evaluated in this stage. The Model R2, Adjusted R2, Mean-absolute-deviation (MAE) and Root-mean-squared-error (RMSE) were reported to indicate the explanatory power of the model. R2 showed the amount of explained variance by the model, MAE was the average magnitude of error without considering the direction (Eq. (9)). RMSE (Eq. (10)) was the square root of the averaged squared residual and being sensitive to outliers (i.e., larger error) as the errors are squared. The RMSE would be reported in the unit of kg/ha.

$$MAE =frac{sum_{i=1}^{n}left|(widehat{y}-y)right|}{n}$$
(9)

$$RMSE= sqrt[]{frac{{{sum }_{i=1}^{n}(widehat{y}-y)}^{2}}{n}}$$
(10)

(widehat{text{y}}) is the ith estimate for AGB of each plot, y is the ith observation of that plot, divided by the (n), denoting the sample size.
However, as the regression model was built with limited sample plots, it might subject to model overfitting. Bootstrapping was adopted to assess statistical accuracy in terms of the confidence intervals54. Ultimately, it provides a robust estimation of the standard errors, confidence intervals for estimates, including the model regression coefficient, mean and correlation coefficients. The model uncertainty of this study was assessed by 1000 runs of bootstrapping. The 95% “Bias-corrected and accelerated (BCa) confidence interval (C.I.)” on the beta coefficient of model predictors was reported.
Leave-one-out cross validation (LOOCV) was utilized for model validation. The model would be trained by all data points except the one being left out for validation55. The Predicted Residual Error Sum of Squares, PRESS, was the sum of the ‘squared deleted residuals’ (SSDR) of the n−1 observation. Predicted R2 was computed by dividing PRESS by the total sum of square residual (SSTO) expressed in Eq. (11). The RMSE of the CV model was also reported as Eq. (12). The RMSE of the CV model without overfitting shall be approximate to that of the original model.

$${R}_{pred}^{2}=1-frac{PRESS}{SSTO}$$
(11)

$$C{V}_{RMSE}= sqrt[]{frac{PRESS}{d.f.}},$$
(12)

d.f. stands for degree of freedom (i.e., d.f. = 21). These model calibration and validation work were conducted in IBM SPSS Statistics version 24. More

1.
Devendra, C. & Thomas, D. Smallholder farming systems in Asia. Agric. Syst. 71, 17–25 (2002).
Article  Google Scholar 
2.
Herrero, M. et al. Smart investments in sustainable food production: revisiting mixed crop-livestock systems. Science 327, 822–825 (2010).
ADS  CAS  PubMed  Article  Google Scholar 

3.
Wright, I. A. et al. Integrating crops and livestock in subtropical agricultural systems. J. Sci. Food Agric. 92, 1010–1015 (2011).
PubMed  Article  CAS  Google Scholar 

4.
Halstead, P. Pastoralism or household herding? Problems of scale and specialization in early Greek animal husbandry. World Archaeol. 28, 20–42 (1996).
Article  Google Scholar 

5.
Bogaard, A. et al. Crop manuring and intensive land management by Europe’s first farmers. Proc. Natl. Acad. Sci. 110, 12589–12594 (2013).
ADS  CAS  PubMed  Article  Google Scholar 

6.
Altieri, M. A., Nicholls, C. I., Henao, A. & Lana, M. A. Agroecology and the design of climate change-resilient farming systems. Agron. Sustain. Dev. 35, 869–890 (2015).
Article  Google Scholar 

7.
Garrett, R. D. et al. Drivers of decoupling and recoupling of crop and livestock systems at farm and territorial scales. Ecol. Soc. 25, 24 (2020).
Article  Google Scholar 

8.
Verhoeven, J. T. A., Arheimer, B., Yin, C. & Hefting, M. M. Regional and global concerns over wetlands and water quality. Trends Ecol. Evol. 21, 96–103 (2006).
PubMed  Article  Google Scholar 

9.
Liu, J. et al. A high-resolution assessment on global nitrogen flows in cropland. Proc. Natl. Acad. Sci. 107, 8035–8040 (2010).
ADS  CAS  PubMed  Article  Google Scholar 

10.
Macdonald, J. M., & Mcbride, W. D. The transformation of U.S. livestock agriculture: scale, efficiency, and risks (2009).

11.
Gerber, P. J. et al. Tackling climate change through livestock: a global assessment of emissions and mitigation opportunities. Food and Agriculture Organization of the United Nations (2013).

12.
Lin, B. B. Resilience in agriculture through crop diversification: adaptive management for environmental change. Bioscience 61, 183–193 (2011).
Article  Google Scholar 

13.
Gaudin, A. C. M. et al. Increasing crop diversity mitigates weather variations and improves yield stability. PLoS ONE 10, e0113261 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

14.
Peterson, C. A., Eviner, V. T. & Gaudin, A. C. M. Ways forward for resilience research in agroecosystems. Agric. Syst. 162, 19–27 (2018).
Article  Google Scholar 

15.
Bowles, T. M. et al. Long-term evidence shows that crop-rotation diversification increases agricultural resilience to adverse growing conditions in North America. One Earth 2, 284–293 (2020).
Article  Google Scholar 

16.
Lobell, D. B. & Field, C. B. Global scale climate-crop yield relationships and the impacts of recent warming. Environ. Res. Lett. 2, 014002 (2007).
ADS  Article  Google Scholar 

17.
Gornall, J. et al. Implications of climate change for agricultural productivity in the early twenty-first century. Philos. Trans. R. Soc. B 365, 2973–2989 (2010).
Article  Google Scholar 

18.
Osborne, T. M. & Wheeler, T. R. Evidence for a climate signal in trends of global crop yield variability over the past 50 years. Environ. Res. Lett. 8, 024001 (2013).
ADS  Article  Google Scholar 

19.
United Nations. Population division of the department of economic and social affairs of the United Nations: world population prospects. https://population.un.org/wpp/DataQuery/ (2019).

20.
Schmidhuber, J. & Tubiello, F. N. Global food security under climate change. Proc. Natl. Acad. Sci. 104, 19703–19708 (2007).
ADS  CAS  PubMed  Article  Google Scholar 

21.
Bullock, J. M. et al. Resilience and food security: rethinking an ecological concept. J. Ecol. 105, 880–884 (2017).
Article  Google Scholar 

22.
Knapp, S. & van der Heijden, M. G. A. A global meta-analysis of yield stability in organic and conservation agriculture. Nat. Commun. 9, 1–9 (2018).
CAS  Article  Google Scholar 

23.
de Moraes, A. et al. Integrated crop-livestock systems in the Brazilian subtropics. Eur. J. Agron. 57, 4–9 (2014).
Article  Google Scholar 

24.
Niles, M. T., Garrett, R. D. & Walsh, D. Ecological and economic benefits of integrating sheep into viticulture production. Agron. Sustain. Dev. 38, 1–11 (2018).
Article  Google Scholar 

25.
Companhia Nacional de Abastecimento (CONAB). Safra brasileira de grãos: Tabela de levantamento. https://www.conab.gov.br/info-agro/safras/graos (2020).

26.
Empresa Brasileira de Pesquisa Agropecuária (Embrapa). ILPF em números. https://ainfo.cnptia.embrapa.br/digital/bitstream/item/158636/1/2016-cpamt-ilpf-em-numeros.pdf (2016).

27.
Garrett, R. D. et al. Social and ecological analysis of commercial integrated crop livestock systems: current knowledge and remaining uncertainty. Agric. Syst. 155, 136–146 (2017).
Article  Google Scholar 

28.
Bell, L. W. & Moore, A. D. Integrated crop-livestock systems in Australian agriculture: trends, drivers and implications. Agric. Syst. 111, 1–12 (2012).
Article  Google Scholar 

29.
Sulc, R. M. & Franzluebbers, A. J. Exploring integrated crop-livestock systems in different ecoregions of the United States. Eur. J. Agron. 57, 21–30 (2014).
Article  Google Scholar 

30.
Carvalho, P. C. F. et al. Animal production and soil characteristics from integrated crop-livestock systems: toward sustainable intensification. J. Anim. Sci. 96, 3513–3525 (2018).
PubMed  PubMed Central  Article  Google Scholar 

31.
Russelle, M. P., Entz, M. H. & Franzluebbers, A. J. Reconsidering integrated crop-livestock systems in North America. Agron. J. 99, 325–334 (2007).
Article  Google Scholar 

32.
Carvalho, P. C. F. et al. Managing grazing animals to achieve nutrient cycling and soil improvement in no-till integrated systems. Nutr. Cycl. Agroecosyst. 88, 259–273 (2010).
Article  Google Scholar 

33.
Oliveira, C. A. O. et al. Comparison of an integrated crop-livestock system with soybean only: economic and production responses in southern Brazil. Renew. Agric. Food Syst. 29, 230–238 (2013).
Article  Google Scholar 

34.
Ryschawy, J., Choisis, N., Choisis, J. P., Joannon, A. & Gibon, A. Mixed crop-livestock systems: an economic and environmental-friendly way of farming?. Animal 6, 1722–1730 (2012).
CAS  PubMed  Article  Google Scholar 

35.
Peterson, C. A., Bell, L. W., Carvalho, P. C. F. & Gaudin, A. C. M. Resilience of an integrated crop–livestock system to climate change: a simulation analysis of cover crop grazing in southern Brazil. Front. Sustain. Food Syst. 4, 604099 (2020).
Article  Google Scholar 

36.
Chávez, L. F. et al. Diversidade metabólica e atividade microbiana no solo em sistema de integração lavoura-pecuária sob intensidades de pastejo. Pesq. Agropec. Bras. 46, 1254–1261 (2011).
Article  Google Scholar 

37.
Peterson, C. A. et al. Winter grazing does not affect soybean yield despite lower soil water content in a subtropical crop-livestock system. Agron. Sustain. Dev. 39, 26 (2019).
CAS  Article  Google Scholar 

38.
Assmann, J. M. et al. Soil carbon and nitrogen stocks and fractions in a long-term integrated crop-livestock system under no-tillage in southern Brazil. Agric. Ecosyst. Environ. 190, 52–59 (2014).
CAS  Article  Google Scholar 

39.
Peyraud, J. L. & Peeters, A. The role of grassland based production system in the protein security. Grassland Science in Europe – The multiple roles of grassland in the European bioeconomy 21, 29–43 (2016).
Google Scholar 

40.
Harrison, G. W. Stability under environmental stress: resistance, resilience, persistence, and variability. Am. Nat. 113, 659–669 (1979).
MathSciNet  Article  Google Scholar 

41.
Lehman, C. L. & Tilman, D. Biodiversity, stability, and productivity in competitive communities. Am. Nat. 156, 534–552 (2000).
PubMed  Article  Google Scholar 

42.
Isbell, F. et al. Biodiversity increases the resistance of ecosystem productivity to climate extremes. Nature 526, 574–577 (2015).
ADS  CAS  PubMed  Article  Google Scholar 

43.
Lightfoot, C. W. F., Dear, K. B. G. & Mead, R. Intercropping sorghum with cowpea in dryland farming systems in Botswana. II. Comparative stability of alternative cropping systems. Exp. Agric. 23, 435–442 (1987).
Article  Google Scholar 

44.
Li, M., Peterson, C. A., Tautges, N. E., Scow, K. M. & Gaudin, A. C. M. Yields and resilience outcomes of organic, cover crop, and conventional practices in a Mediterranean climate. Sci. Rep. 9, 12283 (2019).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

45.
Nielsen, D. C. & Vigil, M. F. Wheat yield and yield stability of eight dryland crop rotations. Agron. J. 110, 594–601 (2018).
Article  Google Scholar 

46.
Temesgen, T., Keneni, G., Sefera, T. & Jarso, M. Yield stability and relationships among stability parameters in faba bean (Vicia faba L.) genotypes. Crop J. 3, 258–268 (2015).
Article  Google Scholar 

47.
Finlay, K. W. & Wilkinson, G. N. The analysis of adaptation in a plant-breeding programe. Aust. J. Agric. Res. 14, 742–754 (1963).
Article  Google Scholar 

48.
Raun, W. R., Barreto, H. J. & Westerman, R. L. Use of stability analysis for long-term soil fertility experiments. Agron. J. 85, 159–167 (1993).
Article  Google Scholar 

49.
Williams, A. et al. Soil water holding capacity mitigates downside risk and volatility in US rainfed maize: time to invest in soil organic matter?. PLoS ONE 11, e0160974 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

50.
Williams, A. et al. A regionally-adapted implementation of conservation agriculture delivers rapid improvements to soil properties associated with crop yield stability. Sci. Rep. 8, 8467 (2018).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

51.
Tilman, D., Reich, P. B. & Knops, J. M. H. Biodiversity and ecosystem stability in a decade-long grassland experiment. Nature 441, 629–632 (2006).
ADS  CAS  PubMed  Article  Google Scholar 

52.
Craven, D. et al. Multiple facets of biodiversity drive the diversity-stability relationship. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-018-0647-7 (2018).
Article  PubMed  Google Scholar 

53.
Bennett, J. A. et al. Resistance of soil biota and plant growth to disturbance increases with plant diversity. Ecol. Lett. https://doi.org/10.1111/ele.13408 (2019).
Article  PubMed  Google Scholar 

54.
Instituto Nacional de Meteorologia (INMET). Normais climatológicas do Brasil. https://portal.inmet.gov.br/normais (2020).

55.
Soil Survey Staff. Soil Taxonomy: a basic system of soil classification for making and interpreting soil surveys (USDA Natural Resources Conservation Service, 1999).

56.
Comissão de Química e Fertilidade do Solo – RS/SC (CQFS RS/SC). Manual de adubação e calagem para os Estados do Rio Grande do Sul e Santa Catarina (Sociedade Brasileira de Ciência do Solo, 2004).

57.
Barthram, G. T. (1985). Experimental techniques: The HFRO sward stick. In Alcok, M. M. The Hill farming research organization Biennial report 1984/1985, pp. 29–30 (1985).

58.
Mott, G. O., & Lucas, H. L. The design, conduct, and interpretation of grazing trials on cultivated and improved pastures. In Proceedings of the international grassland congress, pp. 1380–1386 (1952).

59.
Klingman, D. L., Miles, S. R. & Mott, G. O. The Cage Method for determining consumption and yield of pasture herbage. Agron. J. 35, 739–746 (1943).
Article  Google Scholar 

60.
Nunes, P. A. A. et al. Grazing intensity determines pasture spatial heterogeneity and productivity in an integrated crop-livestock system. Grassl. Sci. 65, 49–59 (2019).
Article  Google Scholar 

61.
van Zanten, H. H. E., Mollenhorst, H., Klootwijk, C. W., van Middelaar, C. E. & Boer, I. J. M. Global food supply: land use efficiency of livestock systems. Int. J. Life Cycle Assess. 21, 747–758 (2016).
Article  CAS  Google Scholar 

62.
Gil, J. D. B. et al. Tradeoffs in the quest for climate smart agricultural intensification in Mato Grosso, Brazil. Environ. Res. Lett. 13, 064025 (2018).
ADS  Article  Google Scholar 

63.
National Research Council. Growth and Body Reserves. In: Nutrient Requirements of Beef Cattle, pp. 22–39 (NRC, 2016).

64.
USDA. Agricultural research service of the United States Department of Agriculture: FoodData Central. https://fdc.nal.usda.gov/fdc-app.html#/food-details/174270/nutrients (2019).

65.
Banco Central do Brasil. Correção de valores pela caderneta de poupança. https://www.bcb.gov.br (2020).

66.
Agrolink. Cotações dos produtos agropecuários: Bovinos. https://www.agrolink.com.br/cotacoes/historico/rs/boi-gordo-kg-vivo-1kg (2019).

67.
Agrolink. Cotações dos produtos agropecuários: Soja. https://www.agrolink.com.br/cotacoes/historico/rs/soja-em-grao-sc-60kg (2019).

68.
Banco Central do Brasil. Correção de valores pelo Índice Geral de Preços do Mercado (IGP-M/FGV). https://www3.bcb.gov.br/CALCIDADAO/publico/corrigirPorIndice.do?method=corrigirPorIndice (2020).

69.
International Monetary Fund. Exchange rate archives by month. https://www.imf.org/external/np/fin/data/param_rms_mth.aspx (2019).

70.
Companhia Nacional de Abastecimento (CONAB). Planilhas de custos de produção – Séries históricas. https://www.conab.gov.br/info-agro/custos-de-producao/planilhas-de-custo-de-producao/itemlist/category/414-planilhas-de-custos-de-producao-series-historicas (2019).

71.
R Core Team. R: a language and environment for statistical computing (2018).

72.
Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).
Article  Google Scholar 

73.
Lenth, R. emmeans: estimated marginal means, aka least-squares means. R package version 1.3.1. (2018).

74.
Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest package: tests in linear mixed effects models. J. Stat. Softw. 82, 1–26 (2017).
Article  Google Scholar 

75.
Fox, J. & Weisberg, S. An {R} companion to applied regression (2011).

76.
Hothorn, T., Bretz, F. & Westfall, P. Simultaneous inference in general parametric models. Biometrical J. 50, 346–363 (2008).
MathSciNet  MATH  Article  Google Scholar 

77.
de Mendiburu, F. agricolae: statistical procedures for agricultural research. R package version 1.2–8. (2017).

78.
Kunrath, T. R., Carvalho, P. C. F., Cadenazzi, M., Bredemeier, C. & Anghinoni, I. Grazing management in an integrated crop-livestock system: soybean development and grain yield. Rev. Ciência Agronômica 46, 645–653 (2015).
Google Scholar 

79.
Peterson, C. A., Deiss, L. & Gaudin, C. M. Commercial integrated crop-livestock systems achieve comparable crop yields to specialized production systems: a meta-analysis. PLoS ONE 15, e0231840 (2020).
CAS  PubMed  PubMed Central  Article  Google Scholar 

80.
Franzluebbers, A. J. & Stuedemann, J. A. Soil physical responses to cattle grazing cover crops under conventional and no tillage in the Southern Piedmont USA. Soil Tillage Res. 100, 141–153 (2008).
Article  Google Scholar 

81.
Tracy, B. F. & Zhang, Y. Soil compaction, corn yield response, and soil nutrient pool dynamics within an integrated crop-livestock system in Illinois. Crop Sci. 48, 1211–1218 (2008).
CAS  Article  Google Scholar 

82.
Schmitt, J. Nematoides fitoparasitas e de vida livre como bioindicadores de qualidade do solo de um sistema de integração lavoura-pecuária (Universidade Federal de Santa Maria, 2019).

83.
Bronick, C. J. & Lal, R. Soil structure and management: a review. Geoderma 124, 3–22 (2005).
ADS  CAS  Article  Google Scholar 

84.
Ingram, L. J. et al. Grazing impacts on soil carbon and microbial communities in a mixed-grass ecosystem. Soil Sci. Soc. Am. J. 72, 939–948 (2008).
ADS  CAS  Article  Google Scholar 

85.
Rosenfeld, J. S. Functional redundancy in ecology and conservation. Oikos 98, 156–162 (2002).
Article  Google Scholar 

86.
Oliver, T. H. et al. Biodiversity and resilience of ecosystem functions. Trends Ecol. Evol. 30, 673–684 (2015).
PubMed  Article  Google Scholar 

87.
Noy-Meir, I. et al. Stability of grazing systems: an application of predator-prey graphs. J. Ecol. 63, 459–481 (1975).
Article  Google Scholar 

88.
Franzluebbers, A. J. et al. Well-managed grazing systems: a forgotten hero of conservation. J. Soil Water Conserv. 67, 100A-104A (2012).
Article  Google Scholar 

89.
Schuster, M. Z. et al. Grazing intensities affect weed seedling emergence and the seed bank in an integrated crop-livestock system. Agric. Ecosyst. Environ. 232, 232–239 (2016).
Article  Google Scholar 

90.
Kunrath, T. R. et al. Sward height determines pasture production and animal performance in a long-term soybean-beef cattle integrated system. Agric. Syst. 177, 102716 (2020).
Article  Google Scholar 

91.
Mott, G. O. Grazing pressure and the measurement of pasture production. In: Proceedings of the International Grassland Congress, pp. 606–611 (1960).

92.
Maraschin, G. E. et al. Native pasture, forage on offer and animal response. In: Proceedings of the international grassland congress, pp. 26–27 (1997).

93.
de Souza Filho, W. et al. Mitigation of enteric methane emissions through pasture management in integrated crop-livestock systems: trade-offs between animal performance and environmental impacts. J. Clean. Prod. 213, 968–975 (2019).

94.
Soussana, J.-F. & Lemaire, G. Coupling carbon and nitrogen cycles for environmentally sustainable intensification of grasslands and crop-livestock systems. Agric. Ecosyst. Environ. 190, 9–17 (2014).
CAS  Article  Google Scholar 

95.
Food and Agriculture Organization of the United Nations (FAO). Trade and markets: the FAO meat price index. http://www.fao.org/economic/est/est-commodities/meat/en/ (2020).

96.
Center for Advanced Studies on Applied Economics (CEPEA). Agricultural prices: Soybean. https://www.cepea.esalq.usp.br/en/indicator/soybean.aspx (2020).

97.
Center for Advanced Studies on Applied Economics (CEPEA). Agricultural prices: Cattle. https://www.cepea.esalq.usp.br/en/indicator/cattle.aspx (2020). More

1.
Aitken, S. N., Yeaman, S., Holliday, J. A., Wang, T. & Curtis-McLane, S. Adaptation, migration or extirpation: climate change outcomes for tree populations. Evol. Appl. 1, 95–111 (2008).
Article  Google Scholar 
2.
Davis, M. B. & Shaw, R. G. Range shifts and adaptive responses to Quaternary climate change. Science 292, 673–679 (2001).
CAS  Article  Google Scholar 

3.
Ikeda, D. H. et al. Genetically informed ecological niche models improve climate change predictions. Glob. Change Biol. 23, 164–176 (2017).
Article  Google Scholar 

4.
Maguire, K. C., Shinneman, D. J., Potter, K. M. & Hipkins, V. D. Intraspecific niche models for ponderosa pine (Pinus ponderosa) suggest potential variability in population-level response to climate change. Syst. Biol. https://doi.org/10.1093/sysbio/syy017 (2018).

5.
Smith, A. B., Godsoe, W., Rodríguez-Sánchez, F., Wang, H.-H. & Warren, D. Niche estimation above and below the species level. Trends Ecol. Evol. 34, 260–273 (2019).
Article  Google Scholar 

6.
Razgour, O. et al. Considering adaptive genetic variation in climate change vulnerability assessment reduces species range loss projections. Proc. Natl Acad. Sci. USA 116, 10418–10423 (2019).
CAS  Article  Google Scholar 

7.
Wang, T., O’Neill, G. A. & Aitken, S. N. Integrating environmental and genetic effects to predict responses of tree populations to climate. Ecol. Appl. 20, 153–163 (2010).
CAS  Article  Google Scholar 

8.
Radeloff, V. C. et al. The rise of novelty in ecosystems. Ecol. Appl. 25, 2051–2068 (2015).
Article  Google Scholar 

9.
Williams, J. W. & Jackson, S. T. Novel climates, no-analog communities, and ecological surprises. Front. Ecol. Environ. 5, 475–482 (2007).
Article  Google Scholar 

10.
Aitken, S. N. & Whitlock, M. C. Assisted gene flow to facilitate local adaptation to climate change. Annu. Rev. Ecol. Evol. Syst. 44, 367–388 (2013).
Article  Google Scholar 

11.
Vitt, P., Havens, K., Kramer, A. T., Sollenberger, D. & Yates, E. Assisted migration of plants: changes in latitudes, changes in attitudes. Biol. Conserv. 143, 18–27 (2010).
Article  Google Scholar 

12.
Williams, M. I. & Dumroese, R. K. Preparing for climate change: forestry and assisted migration. J. For. 111, 287–297 (2013).
Google Scholar 

13.
Keller, S. R., Levsen, N., Olson, M. S. & Tiffin, P. Local adaptation in the flowering-time gene network of balsam poplar, Populus balsamifera L. Mol. Biol. Evol. 29, 3143–3152 (2012).
CAS  Article  Google Scholar 

14.
Keller, S. R., Chhatre, V. E. & Fitzpatrick, M. C. Influence of range position on locally adaptive gene–environment associations in Populus flowering time genes. J. Hered. 109, 47–58 (2018).
CAS  Article  Google Scholar 

15.
Chuine, I. Why does phenology drive species distribution? Phil. Trans. R. Soc. B 365, 3149–3160 (2010).
Article  Google Scholar 

16.
Morin, X., Viner, D. & Chuine, I. Tree species range shifts at a continental scale: new predictive insights from a process-based model. J. Ecol. 96, 784–794 (2008).
Article  Google Scholar 

17.
Ferrier, S., Manion, G., Elith, J. & Richardson, K. Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment. Divers. Distrib. 13, 252–264 (2007).
Article  Google Scholar 

18.
Ellis, N., Smith, S. J. & Pitcher, C. R. Gradient forests: calculating importance gradients on physical predictors. Ecology 93, 156–168 (2012).
Article  Google Scholar 

19.
Fitzpatrick, M. C. & Keller, S. R. Ecological genomics meets community-level modelling of biodiversity: mapping the genomic landscape of current and future environmental adaptation. Ecol. Lett. 18, 1–16 (2015).
Article  Google Scholar 

20.
Fei, S. et al. Divergence of species responses to climate change. Sci. Adv. 3, e1603055 (2017).
Article  Google Scholar 

21.
VanDerWal, J. et al. Focus on poleward shifts in species’ distribution underestimates the fingerprint of climate change. Nat. Clim. Change 3, 239–243 (2013).
Article  Google Scholar 

22.
Shaw, R. G. From the past to the future: considering the value and limits of evolutionary prediction. Am. Nat. 193, 1–10 (2018).
Article  Google Scholar 

23.
Hampe, A. & Petit, R. J. Conserving biodiversity under climate change: the rear edge matters. Ecol. Lett. 8, 461–467 (2005).
Article  Google Scholar 

24.
Yun, J. et al. Influence of winter precipitation on spring phenology in boreal forests. Glob. Change Biol. 24, 5176–5187 (2018).
Article  Google Scholar 

25.
Fu, Y. H. et al. Unexpected role of winter precipitation in determining heat requirement for spring vegetation green-up at northern middle and high latitudes. Glob. Change Biol. 20, 3743–3755 (2014).
Article  Google Scholar 

26.
Peterson, M. L., Doak, D. F. & Morris, W. F. Incorporating local adaptation into forecasts of species’ distribution and abundance under climate change. Glob. Change Biol. 25, 775–793 (2019).
Article  Google Scholar 

27.
Atkins, K. E. & Travis, J. M. J. Local adaptation and the evolution of species’ ranges under climate change. J. Theor. Biol. 266, 449–457 (2010).
CAS  Article  Google Scholar 

28.
Chen, I.-C., Hill, J. K., Ohlemuller, R., Roy, D. B. & Thomas, C. D. Rapid range shifts of species associated with high levels of climate warming. Science 333, 1024–1026 (2011).
CAS  Article  Google Scholar 

29.
Groom, Q. J. Some poleward movement of British native vascular plants is occurring, but the fingerprint of climate change is not evident. PeerJ 1, e77 (2013).
Article  Google Scholar 

30.
Olson, M. S. et al. The adaptive potential of Populus balsamifera L. to phenology requirements in a warmer global climate. Mol. Ecol. 22, 1214–1230 (2013).
CAS  Article  Google Scholar 

31.
Fitzpatrick, M., Chhatre, V., Soolanayakanahally, R. & Keller, S. Experimental support for genomic prediction of climate maladaptation using the machine learning approach Gradient Forests. Preprint at https://doi.org/10.22541/au.159863198.86187354 (2020).

32.
Blois, J. L., Zarnetske, P. L., Fitzpatrick, M. C. & Finnegan, S. Climate change and the past, present, and future of biotic interactions. Science 341, 499–504 (2013).
CAS  Article  Google Scholar 

33.
Keller, S. R. et al. Climate-driven local adaptation of ecophysiology and phenology in balsam poplar, Populus balsamifera L. (Salicaceae). Am. J. Bot. 98, 99–108 (2011).
Article  Google Scholar 

34.
Alberto, F. J. et al. Potential for evolutionary responses to climate change—evidence from tree populations. Glob. Change Biol. 19, 1645–1661 (2013).
Article  Google Scholar 

35.
Little, E. L. Atlas of United States Trees (US Dept of Agriculture, Forest Service, 1971).

36.
Romero-Lankao, P. et al. in Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part B: Regional Aspects (eds Barros, V. R. et al.) 1439–1498 (Cambridge Univ. Press, 2014).

37.
Fetter, K. C., Gugger, P. F. & Keller, S. R. in Comparative and Evolutionary Genomics of Angiosperm Trees (eds Groover, A. & Cronk, Q.) 303–333 (Springer International, 2017); https://doi.org/10.1007/7397_2016_19

38.
Soolanayakanahally, R. Y., Guy, R. D., Silim, S. N., Drewes, E. C. & Schroeder, W. R. Enhanced assimilation rate and water use efficiency with latitude through increased photosynthetic capacity and internal conductance in balsam poplar (Populus balsamifera L.). Plant Cell Environ. 32, 1821–1832 (2009).
CAS  Article  Google Scholar 

39.
Chhatre, V. E. et al. Climatic niche predicts the landscape structure of locally adaptive standing genetic variation. Preprint at https://doi.org/10.1101/817411 (2019).

40.
Günther, T. & Coop, G. Robust identification of local adaptation from allele frequencies. Genetics 195, 205–220 (2013).
Article  Google Scholar 

41.
Frichot, E., Schoville, S. D., Bouchard, G. & François, O. Testing for associations between loci and environmental gradients using latent factor mixed models. Mol. Biol. Evol. 30, 1687–1699 (2013).
CAS  Article  Google Scholar 

42.
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 

43.
Frichot, E. & François, O. LEA: an R package for landscape and ecological association studies. Methods Ecol. Evol. 6, 925–929 (2015).
Article  Google Scholar 

44.
Goudet, J. & Jombart, T. hierfstat: Estimation and tests of hierarchical F-statistics. R package version 0.04-22 (2015).

45.
Manion, G., Lisk, M., Nieto-Lugilde, D., Mokany, K. & Fitzpatrick, M. gdm: Generalized dissimilarity modeling. R package version 1.3.11 (2018).

46.
Hijmans, R. J. geosphere: Spherical trigonometry. R package version 1.5-10 (2019).

47.
Gougherty, A. V., Chhatre, V. E., Keller, S. R. & Fitzpatrick, M. C. Contemporary range position predicts the range-wide pattern of genetic diversity in balsam poplar (Populus balsamifera L.). J. Biogeogr. 47, 1246–1257 (2020).
Article  Google Scholar 

48.
Vallejos, R., Osorio, F. & Bevilacqua, M. Spatial Relationships Between Two Georeferenced Variables: With Applications in R (Springer, 2018). More