Philippa Kaur

Birkeland, C. Coral Reefs in the Anthropocene (Springer, 2015).Book 

Google Scholar 
Kleypas, J. A., Mcmanus, J. W. & Meñez, L. A. B. Environmental limits to coral reef development: Where do we draw the line?. Am. Zool. 39(1), 146–159. https://doi.org/10.1093/icb/39.1.146 (1999).Article 

Google Scholar 
Perry, C. T. & Larcombe, P. Marginal and non-reef-building coral environments. Coral Reefs 22, 427–432. https://doi.org/10.1007/s00338-003-0330-5 (2003).Article 

Google Scholar 
Wilkinson, C. R. Global and local threats to coral reef functioning and existence: review and predictions. Mar. Freshw. Res. 50, 867–878. https://doi.org/10.1071/mf99121 (1999).Article 

Google Scholar 
Mies, M. et al. South atlantic coral reefs are major global warming refugia and less susceptible to bleaching. Front. Mar. Sci. 7, 514. https://doi.org/10.3389/fmars.2020.00514 (2020).Article 

Google Scholar 
Leão, Z. M. A. N. et al. Brazilian coral reefsin a period of global change: A synthesis. Braz. J. Oceanogr. 64, 97–116. https://doi.org/10.1590/S1679-875920160916064sp2 (2016).Article 

Google Scholar 
Coker, D. J., Wilson, S. K. & Pratchett, M. S. Importance of live coral habitat for reef fishes. Rev. Fish Biol. Fish. 24, 89–126. https://doi.org/10.1007/s11160-013-9319-5 (2014).Article 

Google Scholar 
Alvarez-Filip, L., Gill, J. A. & Dulvy, N. K. Complex reef architecture supports more small-bodied fishes and longer food chains on Caribbean reefs. Ecosphere 2, 118. https://doi.org/10.1890/ES11-00185.1 (2011).Article 

Google Scholar 
Wilson, S. K., Graham, N. A. J., Pratchett, M. S., Jones, G. P. & Polunin, N. V. C. Multiple disturbances and the global degradation of coral reefs: Are reef fishes at risk or resilient?. Glob. Change Biol. 12, 2220–2234. https://doi.org/10.1111/j.1365-2486.2006.01252.x (2006).ADS 
Article 

Google Scholar 
Sully, S., Burkepile, D. E., Donovan, M. K., Hodgson, G. & van Woesik, R. A global analysis of coral bleaching over the past two decades. Nat. Commun. 10, 1264. https://doi.org/10.1038/s41467-019-09238-2 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Bellwood, D. R., Hughes, T. P., Folke, C. & Nystrom, M. Confronting the coral reef crisis. Nature 429, 827–833. https://doi.org/10.1038/nature02691 (2004).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Hughes, T. P. et al. climate change, human impacts, and the resilience of coral reefs. Science 301, 929–933. https://doi.org/10.1126/science.1085046 (2003).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Holbrook, N. J. et al. Keeping pace with marine heatwaves. Nat. Rev. Earth Environ. 1, 482–493. https://doi.org/10.1038/s43017-020-0068-4 (2020).ADS 
Article 

Google Scholar 
Bleuel, J., Pennino, M. G. & Longo, G. O. Coral distribution and bleaching vulnerability areas in Southwestern Atlantic under ocean warming. Sci. Rep. 11, 12833. https://doi.org/10.1038/s41598-021-92202-2 (2021).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Fontoura, L. et al. The macroecology of reef fish agonistic behaviour. Ecography 43, 1278–1290. https://doi.org/10.1111/ecog.05079 (2020).Article 

Google Scholar 
Inagaki, K. Y., Pennino, M. G., Floeter, S. R., Hay, M. E. & Longo, G. O. Trophic interactions will expand geographically but be less intense as oceans warm. Glob. Change Biol. 26, 6805–6812. https://doi.org/10.1111/gcb.15346 (2020).ADS 
Article 

Google Scholar 
Longo, G. O., Hay, M. E., Ferreira, C. E. L. & Floeter, S. R. Trophic interactions across 61 degrees of latitude in the Western Atlantic. Glob. Ecol. Biogeogr. 28, 107–117. https://doi.org/10.1111/geb.12806 (2019).Article 

Google Scholar 
Pratchett, M. S. et al. Effects of climate-induced coral bleaching on coral-reef fishes: Ecological and economic consequences. Oceanogr. Mar. Biol. Annu. Rev. 46, 251–296. https://doi.org/10.1201/9781420065756.ch6 (2008).Article 

Google Scholar 
Graham, N. A. J. et al. Lag effects in the impacts of mass coral bleaching on coral reef fish, fisheries, and ecosystems. Conserv. Biol. 21, 1291–1300. https://doi.org/10.1111/j.1523-1739.2007.00754.x (2007).Article 
PubMed 

Google Scholar 
Strona, G. et al. Global tropical reef fish richness could decline by around half if corals are lost. Proc. R. Soc. B 288, 20210274. https://doi.org/10.1098/rspb.2021.0274 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
McClenachan, L. Extinction risk in reef fishes 199–207 (Cambridge University Press, 2015).
Google Scholar 
Power, M. E. et al. Challenges in the quest for keystones. Bioscience 46, 609–620. https://doi.org/10.2307/1312990 (1996).Article 

Google Scholar 
Pereira, P. H. C. et al. The influence of multiple factors upon reef fish abundance and species richness in a tropical coral complex. Ichthyol. Res. 61, 375–384. https://doi.org/10.1007/s10228-014-0409-8 (2014).Article 

Google Scholar 
Coni, E. O. C. et al. An evaluation of the use of branching fire-corals (Millepora spp.) as refuge by reef fish in the Abrolhos Bank, eastern Brazil. Environ. Biol. Fish. 96, 45–55. https://doi.org/10.1007/s10641-012-0021-6 (2013).Article 

Google Scholar 
Graham, N. A. J. et al. Extinction vulnerability of coral reef fishes. Ecol. Lett. 14, 341–348. https://doi.org/10.1111/j.1461-0248.2011.01592.x (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Cornwell, W. K., Schwilk, D. W. & Ackerly, D. D. A trait-based test for habitat filtering: convex hull volume. Ecology 87, 1465–1471. https://doi.org/10.1890/0012-9658(2006)87[1465:ATTFHF]2.0.CO;2 (2006).Article 
PubMed 

Google Scholar 
Mouillot, D., Graham, N. A. J., Villéger, S., Mason, N. W. H. & Bellwood, D. R. A functional approach reveals community responses to disturbances. Trends Ecol. Evol. 28(3), 167–177. https://doi.org/10.1016/j.tree.2012.10.004 (2013).Article 
PubMed 

Google Scholar 
Pimiento, C. et al. Functional diversity of marine megafauna in the Anthropocene. Sci. Adv. 6, 7650. https://doi.org/10.1126/sciadv.aay7650 (2020).ADS 
Article 

Google Scholar 
Loiola, M. et al. Structure of marginal coral reef assemblages under different turbidity regime. Mar. Environ. Res. 147, 138–148. https://doi.org/10.1016/j.marenvres.2019.03.013 (2019).CAS 
Article 
PubMed 

Google Scholar 
Aued, A. W. et al. Large-scale patterns of benthic marine communities in the Brazilian Province. PLoS ONE 13, e0198452. https://doi.org/10.1371/journal.pone.0198452 (2018).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Leão, Z. M. A. N., Kikuchi, R. K. P. & Testa, V. Corals and Coral Reefs of Brazil 9–52 (Elsevier Publisher, 2003).
Google Scholar 
Pinheiro, H. T. et al. South-western Atlantic reef fishes: Zoogeographical patterns and ecological drivers reveal a secondary biodiversity centre in the Atlantic Ocean. Divers. Distrib. 24, 951–965. https://doi.org/10.1111/ddi.12729 (2018).Article 

Google Scholar 
Floeter, S. R. et al. Atlantic reef fish biogeography and evolution. J. Biogeogr. 35, 22–47. https://doi.org/10.1111/j.1365-2699.2007.01790.x (2008).Article 

Google Scholar 
Cord, I. et al. Brazilian marine biogeography: A multi-taxa approach for outlining sectorization. Mar. Biol. 169(5), 61. https://doi.org/10.1007/s00227-022-04045-8 (2022).Article 

Google Scholar 
Leal, I. C. S., Araújo, M. E. D., Cunha, S. R. D. & Pereira, P. H. C. The influence of fire-coral colony size and agonistic behaviour of territorial damselfish on associated coral reef fish communities. Mar. Environ. Res. 108, 45–54. https://doi.org/10.1016/j.marenvres.2015.04.009 (2015).CAS 
Article 
PubMed 

Google Scholar 
Kéry, M. & Royle, J. A. Applied hierarchical modeling in ecology: Analysis of distribution abundance and species richness in R and BUGS. In Prelude and Static Models Vol. 1 (eds Kéry, M. & Royle, J. A.) (Academic Press, 2016).MATH 

Google Scholar 
Hadj-Hammou, J., Mouillot, D. & Graham, N. A. J. Response and effect traits of coral reef fish. Front. Mar. Sci. https://doi.org/10.3389/fmars.2021.640619 (2021).Article 

Google Scholar 
McLean, M. et al. Trait similarity in reef fish faunas across the world’s oceans. PNAS 118(12), e2012318118. https://doi.org/10.1073/pnas.2012318118 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Brandl, S. J. et al. Coral reef ecosystem functioning: eight core processes and the role of biodiversity. Front. Ecol. Environ. 17, 445–454. https://doi.org/10.1002/fee.2088 (2019).Article 

Google Scholar 
Eggertsen, L. et al. Seaweed beds support more juvenile reef fish than seagrass beds in a south-western Atlantic tropical seascape. Estuar. Coast. Shelf S. 196, 97–108. https://doi.org/10.1016/j.ecss.2017.06.041 (2017).ADS 
Article 

Google Scholar 
Mouillot, D. et al. Functional over-redundancy and high functional vulnerability in global fish faunas on tropical reefs. PNAS 111, 13757–13762. https://doi.org/10.1073/pnas.1317625111 (2014).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Briggs, J. C. Marine Zoogeography (McGraw-Hill, 1974).
Google Scholar 
Garcia, G. S., Dias, M. S. & Longo, G. O. Trade-off between number and length of remote videos for rapid assessments of reef fish assemblages. J. Fish Biol. 99(3), 896–904. https://doi.org/10.1111/jfb.14776 (2021).Article 
PubMed 

Google Scholar 
Quimbayo, J. P. et al. Life-history traits, geographical range, and conservation aspects ofreef fishes from the Atlantic and Eastern Pacific. Ecology 102, e03298. https://doi.org/10.1002/ecy.3298 (2021).Article 
PubMed 

Google Scholar 
Katsanevakis, S. et al. Monitoring marine populations and communities: methods dealing with imperfect detectability. Aquat. Biol. 16, 31–52. https://doi.org/10.3354/ab00426 (2012).Article 

Google Scholar 
Villéger, S., Mason, N. W. H. & Mouillot, D. New multidimensional functional diversity indices for a multifaceted framework in functional ecology. Ecology 89, 2290–2301. https://doi.org/10.1890/07-1206.1 (2008).Article 
PubMed 

Google Scholar 
Maire, E., Grenouillet, G., Brosse, S. & Villéger, S. How many dimensions are needed to accurately assess functional diversity? A pragmatic approach for assessing the quality of functional spaces. Glob. Ecol. Biogeogr. 24, 728–740. https://doi.org/10.1111/geb.12299 (2015).Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2021)Kellner, K. jagsUI: A Wrapper Around ‘rjags’ to Streamline ‘JAGS’ Analyses. R package version 1.5.2. https://CRAN.R-project.org/package=jagsUI (2021)Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).Book 

Google Scholar 
Ferreira, C. E. L., Gonçalves, J. E. A. & Coutinho, R. Community structure of fishes and habitat complexity on a tropical rocky shore. Environ. Biol. Fish. 61, 353–369 (2001).Article 

Google Scholar 
Fulton, C. J. et al. Macroalgal meadow habitats support fish and fisheries in diverse tropical seascapes. Fish Fish. 21, 700–717. https://doi.org/10.1111/faf.12455 (2020).Article 

Google Scholar 
Ferreira, L. C. L. et al. Different responses of massive and branching corals to a major heatwave at the largest and richest reef complex in South Atlantic. Mar. Biol. 168, 54. https://doi.org/10.1007/s00227-021-03863-6 (2021).CAS 
Article 

Google Scholar 
Lonzetti, B. C., Vieira, E. A. & Longo, G. O. Ocean warming can help zoanthids outcompete branching hydrocorals. Coral Reefs 41, 175–189. https://doi.org/10.1007/s00338-021-02212-9 (2022).Article 

Google Scholar 
Grillo, A. C., Candido, C. F., Giglio, V. J. & Longo, G. O. Unusual high coral cover in a Southwestern Atlantic subtropical reef. Mar. Biodivers. 51, 77. https://doi.org/10.1007/s12526-021-01221-9 (2021).Article 

Google Scholar 
Matheus, Z. et al. Benthic reef assemblages of the Fernando de Noronha Archipelago, tropical South-west Atlantic: Effects of depth, wave exposure and cross-shelf positioning. PLoS ONE 14(1), e0210664. https://doi.org/10.1371/journal.pone.0210664 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Meirelles, P. M. et al. Baseline assessment of mesophotic reefs of the vitória-trindade seamount chain based on water quality, microbial diversity, benthic cover and fish biomass data. PLoS ONE 10(6), e0130084. https://doi.org/10.1371/journal.pone.0130084 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Ferreira, C. E. L., Floeter, S. R., Gasparini, J. L., Ferreira, B. P. & Joyeux, J. C. Trophic structure patterns of Brazilian reef fishes: A latitudinal comparison. J. Biogeogr. 31, 1093–1106. https://doi.org/10.1111/j.1365-2699.2004.01044.x (2004).Article 

Google Scholar 
Fontoura, L. et al. Climate-driven shift in coral morphological structure predicts decline of juvenile reef fishes. Glob. Change Biol. 26, 557–567. https://doi.org/10.1111/gcb.14911 (2020).ADS 
Article 

Google Scholar 
MacNeil, M. A. et al. Accounting for detectability in reef-fish biodiversity estimates. Mar. Ecol.-Prog. Ser. 367, 249–260. https://doi.org/10.3354/meps07580 (2008).ADS 
Article 

Google Scholar 
Capitani, L., de Araujo, J. N., Vieira, E. A., Angelini, R. & Longo, G. O. Ocean warming will reduce standing biomass in a Tropical Western Atlantic reef ecosystem. Ecosystems 25, 843–857. https://doi.org/10.1007/s10021-021-00691-z (2022).Article 

Google Scholar 
Fogliarini, C. O., Longo, G. O., Francini-Filho, R. B., McClenachan, L. & Bender, M. G. Sailing into the past: Nautical charts reveal changes over 160 years in the largest reef complex in the South Atlantic Ocean. PECON 20(3), 231–239. https://doi.org/10.1007/10.1016/j.pecon.2022.05.003 (2022).Article 

Google Scholar 
Gasparini, J. L., Floeter, S. R., Ferreira, C. E. L. & Sazima, I. Marine ornamental trade in Brazil. Biodivers. Conserv. 14, 2883–2899. https://doi.org/10.1007/s10531-004-0222-1 (2005).Article 

Google Scholar 
Francini-Filho, R. B. et al. Brazil 163–198 (Springer, 2019).
Google Scholar 
Bellwood, D. R., Goatley, C. H. R. & Bellwood, O. The evolution of fishes and corals on reefs: Form, function and interdependence. Biol. Rev. 92, 878–901. https://doi.org/10.1111/brv.12259 (2017).Article 
PubMed 

Google Scholar 
Nunes, L. T. et al. Ecology of Prognathodes obliquus, a butterflyfish endemic to mesophotic ecosystems of St. Peter and St. Paul’s Archipelago. Coral Reefs 38, 955–960. https://doi.org/10.1007/s00338-019-01822-8 (2019).ADS 
Article 

Google Scholar 
Liedke, A. et al. Abundance, diet, foraging and nutritional condition of the banded butterflyfish (Chaetodon striatus) along the western Atlantic. Mar. Biol. 163, 6. https://doi.org/10.1007/s00227-015-2788-4 (2016).CAS 
Article 

Google Scholar  More

Animals, housing, and treatmentsAll procedures involving animals were approved by the University of Florida Institutional Animal Care and Use Committee (Protocol #201709925). All methods were performance in accordance with the relevant guidelines and regulations, and permission and informed consent was obtained from the University of Florida (owners) for the use of the steers in this experiment.The experiment was carried out during July and August of 2017 at the Feed Efficiency Facility of University of Florida, North Florida Research and Education Center, located in Marianna, Florida (30°52′N, 85°11″W, 35 m asl). Both ‘Argentine’ bahiagrass and ‘Florigraze’ rhizoma peanut hays were obtained from commercial producers. The hay bales were stored in enclosed barns throughout the duration of the experiment.Twenty-five Brahman × Angus crossbred steers (Bos sp.) were utilized (average BW = 341 ± 17 kg, approx. 16 months of age). The steers were grazing bermudagrass (Cynodon dactylon) pastures, a C4 grass, prior to the start of the study. The day prior to the start of the experiment (e.g. day-1), steers brought to working facilities, where they remained 16 h off feed and water, in order to obtain shrunk bodyweights. On day 0 of the experiment, steers were weighed, blocked by bodyweight, and allocated to five treatments (5 steers per treatment) and housed in grouped pens. Hay intake was recorded utilizing GrowSafe© systems (GrowSafe Systems Ltd., Calgary, AB, Canada), which utilize radio frequency identification to record feed intake by weight change measured to the nearest gram. Water was available ad libitum. Forage treatments were offered ad libitum by providing sufficient hay to maintain full feed troughs throughout each day of the experiment. Treatments were five proportions of ‘Florigraze’ rhizoma peanut hay in ‘Argentine’ bahiagrass hay: (1) 100% bahiagrass hay (0% RP); (2) 25% rhizoma peanut hay + 75% bahiagrass hay (25% RP); (3) 50% rhizoma peanut hay + 50% bahiagrass hay (50% RP); (4) 75% rhizoma peanut hay + 25% bahiagrass hay (75% RP); (5) 100% rhizoma peanut hay (100% RP). Diet chemical composition is presented in Table 1. All treatment proportions were weighed and mixed on as-fed basis. Mixing of diets was done manually; no hay mixers or choppers were used, to minimize leaf shatter.Sample collectionSteers were housed for 32 days and sampling occurred on 0, 8, 16, 24, and 32 days after initiation of treatment diets; exception was for feces, which were collected on d-1 given steers were fasted on d-0 of the experiment. The hay mixtures offered to the steers were collected (10 samples of each diet) and analyzed for nutritive value (Table 1), at the start of the experiment. All sampling occurred between 0700 and 1000 h on each of the sampling days.Fecal samples were collected directly from the rectum and placed in quart-sized plastic bags to avoid contamination. The feces were frozen at −20 °C. All fecal samples were thawed, dried at 55 °C for 72 h, and ground to pass a 2-mm stainless steel screen using a Wiley Mill (Model 4, Thomas-Wiley Laboratory Mill, Thomas Scientific, Swedesboro, NJ, USA). Samples were then ball milled using a Mixer Mill MM400 (Retsch GmbH, Haan, Germany) at 25 Hz for 9 min.Blood was obtained through jugular venipuncture using 10-mL K2 EDTA vials (Becton Dickinson and Company, Franklin Lakes, NJ, USA), and stored in ice until centrifugation. All tubes were centrifuged at 714 G for 20 min using an Allegra X-22R centrifuge (Beckman Coulter, Brea, CA, USA). A 10-mL sample of plasma was collected and placed in a separate glass vial, the remaining plasma, white blood cell, and platelet fractions were discarded. The remaining RBC pellet was re-suspended with 9 vol. 0.9% NaCl solution and mixed at room temperature for 15 min at 2 Hz orbital shaker. The tubes were then centrifuged at 714 G for 20 min. The saline solution from the centrifuged tubes was discarded after centrifugation. The rinse procedure was repeated two more times for a total of three rinses. After the third rinse procedure, a 500-µL sample was removed, frozen at −20 °C, and subsequently freeze-dried for isotopic analyses.Hair clippings were obtained from an area of 20 × 20 cm on the left hindquarter, utilizing veterinary hair clippers (Sunbeam-Oster Inc., Boca Raton, FL, USA). Hair clippings were collected, placed in nylon bags (Ankom Technology, Macedon, NY, USA), and frozen for subsequent analysis. Clippings were always collected in the same location from each animal in order to ensure new hair growth would be analyzed. All hair samples were cleaned using soapy water and defatted following protocols for other keratin-based tissues 31,34. Each sample was sonicated twice for 30 min in a methanol and chloroform solution (2:1, v/v), rinsed with distilled water, and oven dried overnight at 60 °C. Each hair sample was ball milled using a Mixer Mill MM400 (Retsch GmbH, Haan, Germany) at 25 Hz for 9 min.CalculationsAfter processing, all samples were analyzed for total C and N using a CHNS analyzer through the Dumas dry combustion method (Vario MicroCube, Elementar Americas Inc., Ronkonkoma, NY, USA) coupled to an isotope ratio mass spectrometer (IsoPrime 100, Elementar, Elementar Americas Inc., Ronkonkoma, NY, USA). Certified standards of L-glutamic acid (USGS40, USGS41; United States Geological Survey) were used for isotope ratio mass spectrometer calibration. Isotope ratios were as follows: δ13C of −26.39, + 37.63‰, and δ15N of −4.52, 47.57‰ for USGS40 and USGS41, respectively. Calibration of the IRMS was conducted according to Cook, et al. 35, with an accuracy of ≤ 0.06 ‰ for 15N and ≤ 0.08 ‰ for 13C.The isotope ratio for 13C/12C was calculated as:$$delta^{{{13}}} {text{C}} = , left( {^{{{13}}} {text{C}}/^{{{12}}} {text{C}}_{{{text{sample}}}} {-}^{{{13}}} {text{C}}/^{{{12}}} {text{C}}_{{{text{reference}}}} } right)/ , left( {^{{{13}}} {text{C}}/^{{{12}}} {text{C}}_{{{text{reference}}}} times { 1}000} right)$$
(1)

where δ13C is the C isotope ratio of the sample relative to Pee Dee Belemnite (PDB) standard (‰), 13C/12Csample is the C isotope ratio of the sample, and 13C/12Creference is the C isotope ratio of PDB standard 5. The isotope ratio for 15N/14N was calculated as:$$delta^{{{15}}} {text{N}} = , left( {^{{{15}}} {text{N}}/^{{{14}}} {text{N}}_{{{text{sample}}}} -^{{{15}}} {text{N}}/^{{{14}}} {text{N}}_{{{text{reference}}}} } right)/left( {^{{{15}}} {text{N}}/^{{{14}}} {text{N}}_{{{text{reference}}}} times { 1}000} right)$$
(2)
where δ15N is the N isotope ratio of the sample relative to atmospheric nitrogen (‰), 15N/14Nsample is the N isotope ratio of the sample, and 15N/14Nreference is the N isotope ratio of atmospheric N (standard) 5. The fraction factor (Δ) in this study was considered to be the difference between the diet isotopic composition (δ13C or δ15N) and that of the given sample 5.The dietary proportion of rhizoma peanut hay was back-calculated using δ13C and δ15N of each plant on a DM basis 3. This method is advantageous in that it does not require further tissue processing and facilitates implementation at the field scale. The proportion of rhizoma peanut was estimated using Eq. (3), described by Jones et al. 3:$$%RP=100-left{100 times frac{A-C}{B-C}right}$$
(3)
where %RP is the proportion of RP in the diet, A is the δ13C or δ15N of the sample obtained, B is the δ13C or δ15N of bahiagrass, and C is the δ13C or δ15N of RP.Statistical analysisAll response variables were analyzed using linear mixed model procedures as implemented in SAS PROC GLIMMIX (SAS/STAT 15.1, SAS Institute). Individual animals were considered the experimental unit. Treatment, collection day, and their interaction were considered fixed effects, and block was considered a random effect in the model. The data were analyzed as repeated measures, considering collection day as the repeated measure. The best covariance matrix was selected according to the lowest AICC fit statistic. Least squares treatment means were compared through pairwise t test using the PDIFF option of the LSMEAN statement in the aforementioned procedure. Based on the recommendations by Milliken and Johnson 36 and Saville 37, no adjustment for multiple comparisons was made. Orthogonal polynomial contrasts (linear and quadratic effects) were used to test effects of RP inclusion on isotopic responses. The slice option was used when the treatment × collection day interaction was significant (P ≤ 0.05), using collection day as the factor, to test treatment effects across collection days. Significance was declared at P ≤ 0.05. The contrast statement was used to test for linear or quadratic effects. Regression analyses for the dietary predictions were conducted using PROC REG from SAS.Predictions of dietary proportions based on Eq. (3) were generated for 16 subgroups reflecting combinations of isotope (13C or 15N), day (8 or 32), and sample type (feces, plasm, RBC, or hair). Analyses comparing predicted versus actual diet proportions included several components. First, we computed the concordance correlation coefficient (CCC) following the recommendations from Crawford, et al. 38. The CCC is a measure of agreement that encompasses both precision and accuracy, along with corresponding 95% bias accelerated and corrected (BCa) bootstrap confidence intervals. For comparative purposes we calculated the Pearson correlation coefficient which only reflects precision. Both correlation coefficients range from −1.0 to 1.0 and we interpreted values ≥ 0.80 as indicating strong positive agreement/correlation. Next, we regressed the actual percentages on the predicted percentages using linear regression. Perfect prediction corresponds to the estimated regression line having an intercept of zero and a slope of 1.0. We then partitioned the mean square error (MSE) of the predicted (from Eq. (3), not the above linear regression) versus actual percentages as described in Rice and Cochran 39. This partitioning entails calculating the proportion of MSE attributable to three sources of error: the difference in mean predicted and actual values (mean component, denoted “MC”), the error resulting from the slope of the above linear regression deviating from 1.0 (slope component, denoted “SC”), and random error (random component, denoted “RC”). The results from the above analyses were examined to identify subgroups whose predictions were sufficiently good to warrant hypothesis testing. In this context “good” means that the predicted percentages were strongly correlated with the actual percentages and the magnitudes of the predicted percentages were similar to the actual percentages. The objective of the hypothesis testing was to formally evaluate whether dietary proportions could be directly predicted from Eq. (3) (in contrast to generating predictions using the equation from regressing actual dietary percentages on the predicted percentages from Eq. (3)). Paired two one-sided test (TOST) equivalence tests were conducted for the selected subgroups with α = 0.0540. These tests are formulated such that the null hypothesis is “non-equivalence” and the alternative hypothesis is “equivalence”. An input parameter to the test is the equivalence region, a range for which we consider the mean actual minus predicted difference to be unimportant (“equivalent”) from a practical standpoint. For each equivalence test we also computed the 90% confidence interval for the mean actual minus predicted difference which we refer to as the “minimum equivalence region”. Based on the structure of TOST equivalence tests, to reject the null hypothesis at the 0.05 level, the equivalence region specified for the test must completely contain the minimum equivalence region. For example, if the pre-specified equivalence region is (−15%, 15%) and the computed minimum equivalence region is (−16%, −6%) the null hypothesis would not be rejected. Finally, we re-ran all of the analyses described above for the selected subgroups where, prior to analysis, predicted percentages outside of the valid range were assigned the appropriate boundary value (i.e., predicted percentages  100% were assigned a value of 100%). More

The pathosystemPlantago lanceolata L. is a perennial monoecious ribwort plantain that reproduces both clonally via the production side rosettes, and sexually via wind pollination. Seeds drop close to the mother plant and usually form a long-term seed bank47. Podospharea plantaginis (Castagne; U. Braun and S. Takamatsu) (Erysiphales, Ascomycota) is an obligate biotrophic powdery mildew that infects only P. lanceolata and requires living host tissue through its life cycle48. It completes its life cycle as localized lesions on host leaves, only the haustorial feeding roots penetrating the leaf tissue to feed nutrients from its host. Infection causes significant stress for host plant and may increase the host mortality31. The interaction between P. lanceolata and P. plantaginis is strain-specific, whereby the same host genotype may be susceptible to some pathogen genotypes while being resistant to others49. The putative resistance mechanism includes two steps. First, resistance occurs when the host plant first recognizes the attacking pathogen and blocks its growth. When the first step fails and infection takes place, the host may mitigate infection development. Both resistance traits vary among host genotypes49.Approximately 4000 P. lanceolata populations form a network covering an area of 50 × 70 km in the Åland Islands, SW of Finland. Disease incidence (0/1) in these populations has been recorded systematically every year in early September since 2001 by approximately 40 field assistants, who record the occurrence of the fungus P. plantaginis in the local P. lanceolata populations30. At this time, disease symptoms are conspicuous as infected plants are covered by white mycelia and conidia. The coverage (m2) of P. lanceolata in the meadows was recorded between 2001 and 2008 and is used as an estimate of host population size. In the field survey two technicians estimate Plantago population size by visually estimating how much ground/other vegetation P. lanceolata foliage covers (m2) in each meadow. The proportion of P. lanceolata plants in each population suffering from drought is also estimated annually in the survey. Data on average rainfall (mm) in July and August was estimated separately for each population using detailed radar-measured rainfall (obtained by Finnish Meteorological Institute) and it was available for years 2001–2008.Host population connectivity (SH)27 for each local population i was computed with the formula that takes into account the area of host coverage (m²) of all host populations surveyed, denoted with (Aj), and their spatial location compared to other host populations. We assume that the distribution of dispersal distances from a location are described by negative exponential distribution. Under this assumption, the following formula (1) quantifies for a focal population i, the effect of all other host populations, taking into account their population sizes and how strongly they are connected through immigration to it:$${S}_{i}^{H}=mathop{sum}limits_{jne i}{{{{{rm{e}}}}}}^{{-alpha d}_{{ij}}}sqrt{{A}_{j}}.$$
(1)
here, dij is the Euclidian distance between populations i and j and 1/α equals the mean dispersal distance, which was set to be two kilometres based on results from a previous study16.The annual survey data has demonstrated that P. plantaginis infects annually 2–16% of all host populations and persists as a highly dynamic metapopulation through extinctions and re-colonizations of local populations16. The number of host populations has remained relatively stable over the study period49. The first visible symptoms of P. plantaginis infection appear in late June as white-greyish lesions consisting of mycelium supporting the dispersal spores (conidia) that are carried by wind to the same or new host plants. Six to eight clonally produced generations follow one another in rapid succession, often leading to local epidemic with substantial proportion of the infected hosts by late summer within the host local population. Podosphaera plantaginis produces resting structures, chasmothecia, that appear towards the end of growing season in August–September31. Between 20% and 90% of the local pathogen populations go extinct during the winter, and thus the recolonization events play an important role in the persistence of the pathogen regionally16.Inoculation assay: Effect of connectivity and disease history on phenotypic disease resistanceHost and pathogen material for the experimentTo examine whether the diversity and level of resistance vary among host populations depending on their degree of connectivity (SH) and disease history, we selected 20 P. lanceolata populations for an inoculation assay. These populations occur in different locations in the host network, and were selected based on their connectivity values (S H of selected populations was 37–110 in isolated and 237–336 in highly connected category, Fig. 1). We did not include host populations in the intermediate connectivity category that was used in the population dynamic analyses in the inoculation assay due to logistic constraints. Podosphaera plantaginis is an obligate biotrophic pathogen that requires living host tissue throughout its life cycle, and obtaining sufficient inoculum for experiments is extremely time and space consuming. In both isolated and highly connected categories, half of the populations (IDs 193, 260, 311, 313, 337, 507, 1821, 1999, 2818 and 5206) were healthy during the study years 2001–2014, while half of the populations (IDs 271, 294, 309, 321, 490, 609, 1553, 1556, 1676 and 1847) were infected by P. plantaginis for several years during the same period. We collected P. lanceolata seeds from randomly selected ten individual plants around the patch area from each host population in August 2014.To acquire inoculum for the assay, we collected the pathogen strains as infected leaves, one leaf from ten plant individuals from four additional host populations (IDs 3301, 4684, 1784, and 3108) in August 2014. None of the pathogen populations were same as the sampled host populations and hence, the strains used in the assay all represent allopatric combinations. Both host and pathogen populations selected for the study were separated by at least two kilometres. The collected leaves supporting infection were placed in Petri dishes on moist filter paper and stored at room temperature until later use.Seeds from ten mother plants from each population were sown in 2:1 mixture of potting soil and sand, and grown in greenhouse conditions at 20 ± 2 °C (day) and 16 ± 2 °C (night) with 16:8 L:D photoperiod. Due to the low germination rate of collected seeds, population 260 (isolated and healthy population) was excluded from the study. Seedlings of ten different mother plants were randomly selected among the germinated plants for each population (n = 190), and grown in individual pots until the plants were eight weeks old.The pathogen strains were purified through three cycles of single colony inoculations and maintained on live, susceptible leaves on Petri dishes in a growth chamber 20 ± 2 °C with 16:8 L:D photoperiod. Every two weeks, the strains were transferred to fresh P. lanceolata leaves. Purified powdery mildew strains (M1–M4), one representing each allopatric population (3301, 4684, 1784 and 3108), were used for the inoculation assay. To produce enough sporulating fungal material, repeated cycles of inoculations were performed before the assay.Inoculation assay quantifying host resistance phenotypesIn order to study how the phenotypic resistance of hosts varies depending on population connectivity and infection history, we scored the resistance of 190 host genotypes, ten individuals from each study populations (n = 19), in an inoculation assay. Here, one detached leaf from each plant was exposed to a single pathogen strain (M1–M4) by brushing spores gently with a fine paintbrush onto the leaf. Leaves were placed on moist filter paper in Petri dishes and kept in a growth chamber at 20 ± 2 with a 16/8D photoperiod. All the inoculations were repeated on two individual Petri plates, leading to 760 host genotype—pathogen genotype combinations and a total of 1520 inoculations (19 populations * 10 plant genotypes * 4 pathogen strains * 2 replicates). We then observed and scored the pathogen infection on day 12 post inoculation, under dissecting microscope. The resulting plant phenotypic response was scored as 0 = susceptible (infection) when mycelium and conidia were observed on the leaf surface, and as 1 = resistance (no infection), when no developing lesions could be detected under a dissecting microscope. A genotype was defined resistant only if both inoculated replicates showed similar response (1), and susceptible if one or both replicates became infected (0).Statistical analysesBayesian spatio-temporal INLA model of the changes in host population sizeTo study how the pathogen infection influences on host population growth, we analyzed the relative change in host population size (m2) (defined as population size (t) − population size (t−1))/population size (t−1)) between consecutive years utilizing data from 2001 to 2008 in response to pathogen presence-absence status at t−1 (Supplementary Table 2). To assess whether this depends on host population connectivity, we estimated the separate effects of pathogen presence/absence in the previous year for connectivity categories—high-, low, and intermediate—that were based on the 0.2 and 0.8 quantiles of the host-connectivity values (Fig. 1A and Supplementary Figs. 1, 2). This allowed us to directly assess and compare the effect of the pathogen on host population growth in the extreme categories between isolated and highly connected host populations which were represented in the sampling for the inoculation study (Fig. 2).As covariates, we included the proportion (0–100%) of dry host plants measured each year within each local population as well as data on the amount of rainfall at the summer months (June, July, and August) obtained from the satellite images, as these were suggested be relevant for this pathosystem in an earlier analysis16. Observations where the change in host population size, or the host population coverage had absolute values larger than their 0.99 quantiles in the whole data, were regarded as outliers and omitted from the analysis. Before the analyses, all the continuous covariates were scaled and centred, and the categorical variables were transformed into binary variables.The relative changes in local host population size between consecutive years was analyzed by a Bayesian spatio-temporal statistical model that simultaneously considers the effects of a set of biologically meaningful predictors. The linear predictor thus consists of two parts (2,3):$$beta {X}_{t}+{z}_{t}{A}_{t}$$
(2)
where (beta) represents the correlation coefficients corresponding to the effects of environmental covariates, ({z}_{t}) corresponds to the spatiotemporal random effect, and ({X}_{t}) and ({A}_{t}) project these to the observation locations. For ({z}_{t}) we assume that the observations from a location in consecutive time points (t−1) and t are described by 1st order autoregressive process:$${z}_{t}=varphi {z}_{t-1}+{w}_{t}$$
(3)
where ({w}_{t}) corresponds to spatially structured zero-mean random noise, for which a Matern covariance function is assumed. Statistical inference then targets jointly the covariate effects (beta), the temporal autocorrelation (varphi), and the hyperparameters describing the spatial autocorrelation in wt. From these the overall variance, as well as spatial range—a distance after which spatial autocorrelation ceases to be significant—can be inferred (Supplementary Fig. 3). For more detailed description of the structure of the statistical model and how to do efficient inference with it using R-INLA, we refer to refs. 16,50.Identification of resistance phenotypesThe phenotype composition of each study population was defined by individual plant responses to the four pathogen strains, where each response could be “susceptible = 0” or “resistant = 1”. For example, a phenotype “1111” refers to a plant resistant to all four pathogen strains. The diversity of distinct resistance phenotypes within populations was estimated using the Shannon diversity index as implemented in the vegan software package51. The Shannon diversity index for all four study groups was then analyzed using a linear model with class predictors population type (well-connected or isolated), infection history (healthy or infected), and their interaction.Analysis of population connectivity and infection history effects on host resistanceTo test whether host population resistance varied depending on connectivity (SH) and infection history, we analyzed the inoculation responses (0 = susceptible, 1=resistant) of each host-pathogen combination by using a logit mixed-effect model in the lme4 package52. The model included the binomial dependent variable (resistance-susceptible; 1/0), and class predictors population type (well-connected or isolated), infection history (healthy or infected), mildew strain (M1, M2, M3, and M4) and their interactions. Plant individual and population were defined as random effects, with plant genotype (sample) hierarchically nested under population. Model fit was assessed using chi-square tests on the log-likelihood values to compare different models and significant interactions, and the best model was selected based on AIC-values. P-values for regression coefficients were obtained by using the car package53. We ran all the analyses in R software54.The metapopulation modelWe model the ecological and co-evolutionary dynamics of host and pathogen metapopulations to understand key features of the experimental system that impact on the qualitative patterns observed. The structure and parameters in our model are therefore not estimated using experimental data, but rather are chosen to cover a range of possibilities (e.g., low vs high transmission rates, variation in trade-off shapes for fitness costs). We construct the metapopulations in two stages to account for relatively well and poorly connected demes. All demes are identical in quality (i.e., no differences in intrinsic birth or death rates between demes) and only differ in their connectivity. Our metapopulation consists of an outer network of 20 demes, equally spaced around the unit square (0.2 units apart), and a 7×7 inner lattice of demes at a minimum distance of 0.2 units from the outer network (Fig. 3A), giving a total of 69 demes. Demes that are separated by a Euclidean distance of at most 0.2 are then connected to each other. This means that populations near the centre of the metapopulation are highly connected, while those on the boundary of the metapopulation are poorly connected. This also has the effect of making connections between well and poorly connected demes assortative (i.e., well/poorly connected demes tend to be connected to well/poorly connected demes). We relax the assumption of assortativity in a second type of network by randomly reassigning connections between demes, while maintaining the same degree distribution. (i.e., the probability of two demes being connected is proportionate to their degree). While well connected demes still have more connections to other well connected demes than to poorly connected demes, they are not more likely to be connected to a well connected deme than by chance based on the degree distribution. In both types of network structure, we classify a deme as well-connected if it is in the top 20% of the degree distribution and poorly connected if it is in the bottom 20%.We model the genetics using a multilocus gene-for-gene framework with haploid host and pathogen genotypes characterized by (L) biallelic loci, where 0 and 1 represent the presence and absence, respectively, of resistance and infectivity alleles. Host genotype (i) and pathogen genotype (j) are represented by binary strings: ({x}_{i}^{1}{x}_{i}^{2}ldots {x}_{i}^{L}) and ({y}_{j}^{1}{y}_{j}^{2}ldots {y}_{j}^{L}). Resistance acts multiplicatively such that the probability of host (i) being infected when challenged by pathogen (j) is ({Q}_{{ij}}={sigma }^{{d}_{{ij}}}), where (sigma) is the reduction in infectivity per effective resistance allele and ({d}_{{ij}}={sum }_{k=1}^{L}{x}_{i}^{k}big(1-{y}_{j}^{k}big)) is the number of effective resistance alleles (i.e., the number of loci where hosts have a resistance allele but pathogens do not have a corresponding infectivity allele). Hosts and pathogens with more resistance or infectivity alleles are assumed to pay higher fitness costs, ({c}_{H}left(iright)) eq. (4) and ({c}_{P}left(jright)) eq. (5) with:$${c}_{H}left(iright)={c}_{H}^{1}left(frac{1-{{{{{rm{e}}}}}}^{frac{{c}_{H}^{2}}{L}{sum }_{k=1}^{L}{x}_{i}^{k}}}{1-{{{{{rm{e}}}}}}^{{c}_{H}^{2}}}right)$$
(4)
and$${c}_{P}left(jright)={c}_{P}^{1}left(frac{1-{{{{{rm{e}}}}}}^{frac{{c}_{P}^{2}}{L}{sum }_{k=1}^{L}{y}_{j}^{k}}}{1-{{{{{rm{e}}}}}}^{{c}_{P}^{2}}}right)$$
(5)
where (0 , < , {c}_{H}^{1},; {c}_{P}^{1},le, 1) control the overall strength of the costs (i.e., the maximum proportional reduction in reproduction (hosts) or transmission rate (pathogens)) and ({c}_{H}^{2},; {c}_{P}^{2}in {{mathbb{R}}}_{ne 0}) control the shape of the trade-off. When ({c}_{H}^{2},; {c}_{P}^{2}, < , 0) the costs decelerate (increasing returns) and when ({c}_{H}^{2},; {c}_{P}^{2}, > , 0) the costs accelerate the costs accelerate (decreasing returns) (Supplementary Fig. 4). This formulation, therefore, allows for a wide-range of trade-off shapes that may occur in nature.The dynamics of the (finite) host and pathogen populations are modelled stochastically using the tau-leap method with a fixed step size of (tau=1). For population (p), the mean host birth rate at time (t) for host (i) (6) is$${B}_{i}^{p}left(tright)=left(aleft(1-{c}_{H}left(iright)right)-q{N}_{p}left(tright)right){S}_{i}^{p}left(tright)$$
(6)
where (a) is the maximum per-capita birth rate, (q) is the strength of density-dependent competition on births, ({N}_{p}left(tright)={S}_{i}^{p}left(tright)+{I}_{icirc }^{p}left(tright)) is the local host population size, ({S}_{i}^{p}left(tright)) and ({I}_{icirc }^{p}left(tright)={sum }_{j=1}^{n}{I}_{{ij}}^{p}left(tright)) are the local sizes of susceptible and infected individuals of genotype (i), and ({I}_{{ij}}^{p}left(tright)) is the local size of hosts of genotype (i) infected by pathogen (j). Host mutations occur at an average rate of ({mu }_{H}) per loci (limited to at most one mutation per time step), so that the mean number of mutations from host type (i) to ({i}^{{prime} }) is ({mu }_{H}{m}_{i{i}^{{prime} }}{B}_{i}^{p}left(tright)), where ({m}_{i{i}^{{prime} }}=1) if genotypes (i) and ({i}^{{prime} }) differ at exactly one locus, and is 0 otherwise.The mean local mortalities for susceptible and infected individuals are (b{S}_{i}^{p}left(tright)) and (left(b+alpha right){I}_{{ij}}^{p}left(tright)), respectively, where (b) is the natural mortality rate and (alpha) is the disease-associated mortality rate. The average number of infected hosts that recover is (gamma {I}_{{ij}}^{p}left(tright)), where (gamma) is the recovery rate.The mean number of new local infections of susceptible host type (i) by pathogen (j) eq. (7) is:$${INF}_{{ij}}^{p}left(tright)=beta left(1-{c}_{P}left(jright)right){Q}_{{ij}}{S}_{i}^{p}left(tright){Y}_{j}^{p}left(tright)$$
(7)
where (beta) is the baseline transmission rate and ({Y}_{j}^{p}left(tright)) is the local number of pathogen propagules following mutation and dispersal. Pathogen mutations occur in a similar manner to host mutations, with mutations from type (j) to ({j}^{{prime} }) occurring at rate ({mu }_{P}{m}_{j{j}^{{prime} }}{I}_{circ j}^{p}left(tright)) where ({mu }_{P}) is the mutation rate per loci (limited to at most one mutation per timestep) and ({I}_{circ j}^{p}left(tright)={sum }_{i=1}^{n}{I}_{{ij}}^{p}left(tright)) is the local number of pathogen (j.) Following mutation, the local number of pathogens of type (j) eq. (8) is:$${W}_{j}^{p}left(tright)={I}_{circ j}^{p}left(tright)left(1-{mu }_{P}Lright)+{mu }_{P}{m}_{j{j}^{{prime} }}{I}_{circ j}^{p}left(tright)$$
(8)
Pathogen dispersal occurs following mutation at a rate of (rho) between connected demes, given by the adjacency matrix ({G}_{{pr}}), with ({G}_{varSigma p}) the total number of connections for deme (p). The mean local number of pathogen propagules following mutation and dispersal eq. (9) is therefore:$${Y}_{j}^{p}left(tright)={W}_{j}^{p}left(tright)left(1-rho {G}_{varSigma p}right)+rho mathop {sum }limits_{r=1}^{{M}_{varSigma }}{G}_{{pr}}{W}_{j}^{r}left(tright)$$
(9)
We focus our parameter sweep on: (i) the structure of the network (assortative or random connections); (ii) the strength (left({c}_{H}^{1},; {c}_{P}^{1}right)) and shape (left({c}_{H}^{2},; {c}_{P}^{2}right)) of the trade-offs; (iii) the transmission rate (left(beta right)); and (iv) the dispersal rate (left(rho right)), fixing the remaining parameters as described in Supplementary Table 1 (preliminary investigations suggested they had less of an impact on the qualitative outcome) and conducting 100 simulations per parameter set. For each simulation we initially seed all populations with the most susceptible host type and place the least infective pathogen type in one of the well-connected populations to minimize the risk of early extinction. We then solve the dynamics for 10,000 time steps (preliminary investigations indicated this was a sufficient period for the metapopulations to reach a quasi-equilibrium in terms of overall resistance). We calculate the average level of resistance (proportion of loci with a resistance allele) between time steps 4001 and 5000 (transient dynamics) and over the final 1000 time steps (long-term dynamics) for well and poorly connected demes, categorized according to whether the disease is present in (infected) or absent from (uninfected) the local population at a given time point and discarding simulations where the pathogen is driven globally extinct.We compare the mean level of resistance in infected/uninfected poorly/well-connected populations across all simulations to the empirical results. We say that a simulation is a qualitative ‘match’ for the empirical findings if: (i) in poorly connected demes, the infected populations are on average at least 5% more resistant than uninfected populations; and (ii) in well-connected demes, the uninfected populations are on average at least 5% more resistant than infected populations. In other words, if ({R}_{{CS}}) is the mean resistance for a population with connectivity (C) ((C=W) and (C=P) for well and poorly connected demes, respectively) and infection status (S) ((S=U) and (S=I) for uninfected and infected populations, respectively), then a parameter set is a qualitative ‘match’ for the empirical findings if ({R}_{{WU}} > 1.05{R}_{{WI}}) and (1.05{R}_{{PI}}, > , 1.05{R}_{{PU}}). If these criteria are not met, then the parameter set is a qualitative ‘mismatch’ for the empirical findings. The model is not intended to be a replica of an empirical metapopulation, but rather is used to reveal the key factors which lead to qualitatively similar distributions of resistance and disease incidences observed in the study of the Åland islands. Hence, the purpose of the model is to determine which biological factors are likely to be crucial to the patterns observed herein.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

To compare soil viral community composition within and across terrestrial habitats on a regional scale, viromes were generated from 34 near-surface (top 15 cm) soil samples, with a total of 30 viromes included in downstream ecological analyses (see Supplementary Methods). The analyzed viromes were collected from four distinct habitats (wetlands, grasslands, chaparral shrublands, and woodlands, each with 7, 14, 4, and 5 viromes, respectively) across five field sites (Fig. S1 for sampling scheme, Table S1 for soil properties). Following quality filtering, the 30 viromes generated an average of 72,950,833 reads and 416 contigs ≥10 Kbp per virome (Table S2). Wetland viromes yielded more contigs ≥10 Kbp than viromes from other habitats, both in total and on average per virome (Table S2). We used VIBRANT to identify 3490 viral contigs in our assemblies, which were clustered into 3,432 viral operational taxonomic units (vOTUs), defined as ≥10 Kbp viral contigs sharing ≥ 95% average nucleotide identity over 85% contig length [17]. Constrained analysis of principal coordinates (CAP analysis) revealed strong clustering by habitat rather than by site, implying that, where environmental parameters are substantially different, environmental conditions are stronger drivers of viral community composition than geographic distance (Fig. S2).Multiple lines of evidence suggest that wetter soils harbored greater viral diversity than drier soils. We recovered the most vOTUs from wetlands, both in total (56% of all vOTUs were from wetlands) and per virome (on average, 307 vOTUs were recovered per wetland virome, compared to 116 from all habitats) (Fig. 1A). Unsurprisingly, wetlands had significantly greater moisture content than other habitats (Fig. 1B; ANOVA followed by Tukey multiple comparisons of means, p 100 Km distances here. Taken together, we propose that soil viral communities often display high heterogeneity within and among habitats, presumably due to a combination of host adaptations and microdiversity, dispersal limitation, and fluctuating environmental conditions over space and time. More

Vega, F. E., Rosenquist, E. & Collins, W. Global project needed to tackle coffee crisis. Nature 425, 343 (2003).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Craparo, A. C. W., Van Asten, P. J. A., Läderach, P., Jassogne, L. T. P. & Grab, S. W. Coffea arabica yields decline in Tanzania due to climate change: global implications. Agric. For. Meteorol. 207, 1–10 (2015).ADS 
Article 

Google Scholar 
Davis, A. P. et al. High extinction risk for wild coffee species and implications for coffee sector sustainability. Sci. Adv. 5, eaav3473 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Davis, A. P., Gole, T. W., Baena, S. & Moat, J. The impact of climate change on indigenous arabica coffee (Coffea arabica): predicting future trends and identifying priorities. PLoS ONE 7, e47981 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Davis, A. P., Mieulet, D., Moat, J., Sarmu, D. & Haggar, J. Arabica-like flavour in a heat-tolerant wild coffee species. Nat. Plants 7, 413–418 (2021).CAS 
PubMed 
Article 

Google Scholar 
Moat, J., Gole, T. W. & Davis, A. P. Least concern to endangered: applying climate change projections profoundly influences the extinction risk assessment for wild Arabica coffee. Global Change Biol. 25, 390–403 (2019).ADS 
Article 

Google Scholar 
Moat, J. et al. Resilience potential of the Ethiopian coffee sector under climate change. Nat. Plants 3, 17081 (2017).PubMed 
Article 

Google Scholar 
Kath, J. et al. Not so robust: Robusta coffee production is highly sensitive to temperature. Global Change Biol. 26, 3677–3688 (2020).ADS 
Article 

Google Scholar 
Liu, L. et al. Soil moisture dominates dryness stress on ecosystem production globally. Nat. Commun. 11, 1–9 (2020).ADS 
CAS 

Google Scholar 
Grossiord, C. et al. Plant responses to rising vapor pressure deficit. New Phytol. 226, 1550–1566 (2020).PubMed 
Article 

Google Scholar 
IPCC Climate Change 2022: Impacts, Adaptation, and Vulnerability (eds. Pörtner, H.-O. et al.) (Cambridge Univ. Press, 2022).Burke, M. et al. Higher temperatures increase suicide rates in the United States and Mexico. Nat. Clim. Change 8, 723–729 (2018).ADS 
Article 

Google Scholar 
Burke, M., Hsiang, S. M. & Miguel, E. Global non-linear effect of temperature on economic production. Nature 527, 235–239 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Duffy, K. A. et al. How close are we to the temperature tipping point of the terrestrial biosphere? Sci. Adv. 7, eaay1052 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Scheffer, M. et al. Early-warning signals for critical transitions. Nature 461, 53–59 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Schneider, S. H. Abrupt non-linear climate change, irreversibility and surprise. Global Environ. Change 14, 245–258 (2004).Article 

Google Scholar 
Lenton, T. M. Early warning of climate tipping points. Nat. Clim. Change 1, 201–209 (2011).ADS 
Article 

Google Scholar 
Lenton, T. M. et al. Climate tipping points—too risky to bet against. Nature. 575, 592–595 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Lobell, D. B., Bänziger, M., Magorokosho, C. & Vivek, B. Nonlinear heat effects on African maize as evidenced by historical yield trials. Nat. Clim. Change 1, 42–45 (2011).ADS 
Article 

Google Scholar 
Lobell, D. B., Deines, J. M. & Tommaso, S. D. Changes in the drought sensitivity of US maize yields. Nat. Food 1, 729–735 (2020).Article 

Google Scholar 
Lobell, D. B. et al. Greater sensitivity to drought accompanies maize yield increase in the US Midwest. Science 344, 516–519 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Rigden, A., Mueller, N., Holbrook, N., Pillai, N. & Huybers, P. Combined influence of soil moisture and atmospheric evaporative demand is important for accurately predicting US maize yields. Nat. Food 1, 127–133 (2020).Article 

Google Scholar 
Schlenker, W. & Roberts, M. J. Nonlinear temperature effects indicate severe damages to US crop yields under climate change. Proc. Natl Acad. Sci. USA 106, 15594–15598 (2009).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
McDowell, N. G. et al. Mechanisms of woody-plant mortality under rising drought, CO2 and vapour pressure deficit. Nat. Rev. Earth Environ. 3, 294–308 (2022).ADS 
CAS 
Article 

Google Scholar 
Sinclair, T. R. et al. Limited-transpiration response to high vapor pressure deficit in crop species. Plant Sci. 260, 109–118 (2017).CAS 
PubMed 
Article 

Google Scholar 
López, J., Way, D. A. & Sadok, W. Systemic effects of rising atmospheric vapor pressure deficit on plant physiology and productivity. Global Change Biol. 27, 1704–1720 (2021).ADS 
Article 

Google Scholar 
McDowell, N. G. & Allen, C. D. Darcy’s law predicts widespread forest mortality under climate warming. Nat. Clim. Change 5, 669–672 (2015).ADS 
Article 

Google Scholar 
Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Sci. Data 5, 170191 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
You, L., Wood, S., Wood-Sichra, U. & Wu, W. Generating global crop distribution maps: from census to grid. Agric. Syst. 127, 53–60 (2014).Article 

Google Scholar 
Fong, Y., Huang, Y., Gilbert, P. B. & Permar, S. R. chngpt: threshold regression model estimation and inference. BMC Bioinformatics 18, 1–7 (2017).Article 

Google Scholar 
Qin, Y. et al. Agricultural risks from changing snowmelt. Nat. Clim. Change 10, 459–465 (2020).ADS 
Article 

Google Scholar 
Forster, P. M., Maycock, A. C., McKenna, C. M. & Smith, C. J. Latest climate models confirm need for urgent mitigation. Nat. Clim. Change 10, 7–10 (2020).ADS 
Article 

Google Scholar 
Forster, P. M. et al. Projections of when temperature change will exceed 2 °C above pre-industrial levels. Nat. Clim. Change 10, 407–412 (2011).
Google Scholar 
Joshi, M., Hawkins, E., Sutton, R., Lowe, J. & Frame, D. Projections of when temperature change will exceed 2 °C above pre-industrial levels. Nat. Clim. Change 1, 407–412 (2011).ADS 
Article 

Google Scholar 
IPCC, 2021: Summary for Policymakers. In Climate Change 2021: The Physical Science Basis (eds Masson-Delmotte, V. et al.) (Cambridge Univ. Press, in press).Lobell, D. B. et al. The critical role of extreme heat for maize production in the United States. Nat. Clim. Change 3, 497–501 (2013).ADS 

Google Scholar 
Sinclair, T. R., Hammer, G. L. & Van Oosterom, E. J. Potential yield and water-use efficiency benefits in sorghum from limited maximum transpiration rate. Funct. Plant Biol. 32, 945–952 (2005).PubMed 
Article 

Google Scholar 
Martins, M. Q. et al. Protective response mechanisms to heat stress in interaction with high [CO2] conditions in Coffea spp. Front. Plant Sci. 7, 947 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Rodrigues, W. P. et al. Long‐term elevated air [CO2] strengthens photosynthetic functioning and mitigates the impact of supra‐optimal temperatures in tropical Coffea arabica and C. canephora species. Global Change Biol. 22, 415–431 (2016).ADS 
Article 

Google Scholar 
Ghini, R. et al. Coffee growth, pest and yield responses to free-air CO2 enrichment. Clim. Change 132, 307–320 (2015).ADS 
Article 

Google Scholar 
Rakocevic, M. et al. The vegetative growth assists to reproductive responses of Arabic coffee trees in a long-term FACE experiment. Plant Growth Regul. 91, 305–316 (2020).CAS 
Article 

Google Scholar 
Hammer, G. L. et al. Designing crops for adaptation to the drought and high‐temperature risks anticipated in future climates. Crop Sci. 60, 605–621 (2020).Article 

Google Scholar 
Gennari, P., Rosero-Moncayo, J. & Tubiello, F. N. The FAO contribution to monitoring SDGs for food and agriculture. Nat. Plants 5, 1196–1197 (2019).PubMed 
Article 

Google Scholar 
Lesk, C., Rowhani, P. & Ramankutty, N. Influence of extreme weather disasters on global crop production. Nature 529, 84–87 (2016).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Ortiz-Bobea, A., Ault, T. R., Carrillo, C. M., Chambers, R. G. & Lobell, D. B. Anthropogenic climate change has slowed global agricultural productivity growth. Nat. Clim. Change 11, 306–312 (2021).ADS 
Article 

Google Scholar 
Davis, A. P. et al. Hot coffee: the identity, climate profiles, agronomy, and beverage characteristics of Coffea racemosa and C. zanguebariae. Front. Sustain. Food Syst. 5, 740137 (2021).Article 

Google Scholar 
Sarmiento-Soler, A. et al. Disentangling effects of altitude and shade cover on coffee fruit dynamics and vegetative growth in smallholder coffee systems. Agric. Ecosyst. Environ. 326, 107786 (2022).CAS 
Article 

Google Scholar 
Wood, S. N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J. R. Stat. Soc. B 73, 3–36 (2011).MathSciNet 
MATH 
Article 

Google Scholar 
Barton, K. MuMIn: multi-model inference. R-Forge http://r-forge.r-project.org/projects/mumin/ (2009).R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing https://www.r-project.org/ (2021).Harrison, X. A. et al. A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ 6, e4794 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Najafi, E., Devineni, N., Khanbilvardi, R. M. & Kogan, F. Understanding the changes in global crop yields through changes in climate and technology. Earths Future 6, 410–427 (2018).ADS 
CAS 
Article 

Google Scholar 
Ovalle-Rivera, O. et al. Assessing the accuracy and robustness of a process-based model for coffee agroforestry systems in Central America. Agrofor. Syst. 94, 2033–2051 (2020).Article 

Google Scholar 
Varma, S. & Simon, R. Bias in error estimation when using cross-validation for model selection. BMC Bioinformatics 7, 1–8 (2006).Article 

Google Scholar 
Yuan, W. et al. Increased atmospheric vapor pressure deficit reduces global vegetation growth. Sci. Adv. 5, eaax1396 (2019).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Son, H. & Fong, Y. Fast grid search and bootstrap-based inference for continuous two-phase polynomial regression models. Environmetrics 32, e2664 (2021).MathSciNet 
Article 

Google Scholar 
Wintgens, J. N. et al. Coffee: Growing, Processing, Sustainable Production. A Guidebook for Growers, Processors, Traders, and Researchers (Wiley, 2004). More

Study site and field data collectionPermanent 1 m2 photoquadrats were established on Heron Reef in 1962/63, using 9 mm diameter mild steel (rebar) pegs, which were replaced over time. From the 1990’s, replacement pegs were stainless steel for greater longevity. Four sites were established, the protected (south) crest, inner flat, exposed (north) crest and exposed pools. Co-ordinates for each site are presented in Table 1, the layout shown in Fig. 2, and sites have been well described previously5,6. At each census, a 1 m2 frame divided into a 5 × 5 grid using string was placed over the pegs, and the quadrat photographed from directly above at low tide. From 1963 until 2003, a 35 mm camera and colour slide film were used. The camera was attached to a tripod affixed to the 1 m2 frame, and captured around 2/3 of the quadrat. The frame (and camera) were then rotated 180 degrees to capture the remainder of the quadrat. After 2003, a hand-held digital camera was used, with the entire quadrat being captured in a single image. Concurrent with each census, mud maps of each quadrat were hand drawn in the field, and all colonies identified in situ by someone with expertise in coral taxonomy.Table 1 Coordinates of the study sites on Heron Island Reef (WGS84).Full size tableFig. 2Quadrat layouts for each of the four sites respectively, noting that the north crest and north ridge have been treated as a single north crest site in previous publications. Underlining indicates original 1962/63 quadrats. Other quadrats were added in or after 2008, as indicated in the text. Contiguous quadrats are pictured bordering each other. Spacing between separate quadrats or groups of quadrats is not shown to scale. Note that up until 2005, NRNW was known as NR. The acronyms in each quadrat represent its name.Full size imageAt the protected (south) crest, a set of six contiguous quadrats were established in 1963 in a 2 × 3 arrangement parallel to the waterline, and about 420 m southeast of the island. This site is exposed at low tide, and was photographed once all water had drained off it. Images of quadrats A, C & E (the shoreward row) from 1963 to 2012 have been fully processed, and the data have been through QA/QC. Data after 2012 exist as images only. These quadrats form the basis of previous analyses1,4,5,6 for this site. Photographs are available for quadrats B, D & F, but apart from 2003–2010, have not been processed. In 2010, an additional two quadrats were established either side of the original six, leading to a 2 × 5 arrangement. Again, only imagery is available for these additional quadrats.At the inner flat, two pairs of contiguous quadrats were established in 1962, 44 m apart, about 70 m south of the island. This site is covered by ~10 cm of water at low tide, so could only be photographed on a still day. Imagery for this site is only available to 2012, after which the marker stakes appear to have been removed in a cleanup of the area. Images for one quadrat in each pair have been processed, but have not been subject to full QA/QC.At the exposed (north) crest main site, a set of four contiguous quadrats was established about 1100 m northeast of the island in 1963. An additional single quadrat (north ridge) was established 326 m to the east. Images from 1963 to 2012 have been fully processed, and the data have been through QA/QC. Data after 2012 exist as images only. In 2005, the single north ridge quadrat was expanded to 4 m2, and in 2008, both subsites were expanded to six quadrats in a 2 × 3 arrangement. These additional quadrats have been digitised up to 2012, but have not been through full QA/QC.The exposed pools are two individual quadrats about 5 m apart about 30 m north of the eastern (north ridge) exposed crest site. These are on the edge of a natural pool, and range from ~5–50 cm deep at low tide, and so could only be photographed on a calm day. Imagery for this site is only available until 2005, after which the marker stakes could not be relocated. Images from 1963 to 1998 have been processed, but have not been through full QA/QC.Retrieval of coral composition data from the photoquadratsProcessing of the images involved scanning the colour slides to produce digital images, and then orthorectifying each image to a 1 m2 basemap in ArcGIS (ESRI Ltd). The corners of the frame, and the holes for the string grid, were used as control points for the orthorectification. For images that originated as colour slides, each half of the quadrat was individually orthorectified to the same basemap, producing a single image of the entire quadrat (see Fig. 3). While contiguous quadrats were orthorectified individually, they were done so against a basemap containing all quadrats in the group, meaning that the resulting images can be easily merged to create a single image of the group. The outlines of all visible coral colonies ( >~1 cm2), and other benthic organisms such as algae and clams, were then digitised in ArcGIS to create a single shapefile for each quadrat for each year. Each colony was represented as an individual feature within the shapefile, and was assigned a unique colony number and species based on the mud maps drawn in the field. Colony numbers were consistent across years, allowing individual colonies to be tracked over time. If a colony underwent fission, the original colony number was retained for each, with the addition of a unique identifier after a decimal point. For example, if colony 35 split in two, the resultant colonies were identified as 35.1 and 35.2. If 35.2 later split again, the resultant colonies were identified as 35.2.1 and 35.2.2. If the colony overlapped the edge of the quadrat, only the area within the quadrat was digitised, and a flag was applied to indicate that only part of the colony was included (edgestatus = 1 in the data). Upon completion of digitisation, ArcGIS was used to calculate the area and perimeter of all colonies. While multiple census were conducted in 1963, 1971 and 1983, only a single census in each year has been processed. There are currently no plans to undertake further digitisation or QA/QC of this data set.Fig. 3Example orthorectified and stitched (prior to 2001) images from the NCNE quadrat, showing the effects of a cyclone that removed all colonies in 1972, and slow recovery over subsequent decades.Full size image More

Characteristics of LVD dataset and assembled vOTUsThe land use virome dataset LVD was derived from 2.6 billion paired clean reads of sequences across 50 viromes of 25 samples with five types of land uses (Supplementary Data 2). A total of 6,442,065 contigs ( >1500 bp) were yielded, of which 764,466 (11.8%) contigs were identified as putative viral genomes through VIBRANT. Subsequently, putative false positive viral genomes were removed (see Methods section), and 27,951 and 48,936 bona fide viral genomes were retained from the 25 intracellular VLPs (iVLPs) and 25 extracellular VLPs (eVLPs) viromes, respectively. These genomes were clustered into 25,941 and 45,152 vOTUs for iVLPs and eVLPs viromes, respectively, in which the iVLPs and eVLPs viromes shared 11,467 (19.2%) vOTUs. Subsequently, they were merged and dereplicated, resulting in 59,626 vOTUs (Supplementary Data 3) for the following analysis. A total of 8112 (13.6%) vOTUs genomes were classified as complete, in which the median length of all and circular vOTUs were 25,183 bp and 45,511 bp, respectively (Supplementary Fig. 4).To explore the taxonomic affiliation of vOTUs in family and genus-level, a gene-sharing network consist of 59,626 vOTUs genomes from this study and 3502 reference phage genomes (from NCBI Viral RefSeq version 201) revealed 6009 VCs comprising of 37,224 vOTUs, of which 34,417 vOTUs were from LVD, besides 2794 singletons (2653 from LVD dataset), 16,056 outliers (15,833 from LVD) and 8492 overlaps (8061 from LVD) were detected (Supplementary Data 4). Of these, only 157 VCs contained genomes from both the RefSeq and LVD dataset (1864 viral genomes) (Supplementary Data 4). Most of VCs (1837, 30.4%) included only two members.At the family level, most of vOTUs were classified into Siphoviridae (712 by vConTACT2 and 29,671 (50.9%) by Demovir, tailed dsDNA), Podoviridae (610 by vConTACT2 and 9923 (17.6 %) by Demovir, tailed dsDNA), Myoviridae (485 by vConTACT2 and 5445 (9.9%) by Demovir, tailed dsDNA), Tectiviridae (50 by vConTACT2 and 10 (0.10%) by Demovir, non-tailed dsDNA) (Fig. 1). Besides, the Eukaryotic viruses Herpesviridae (159 by Demovir, 0.26%, dsDNA), Phycodnaviridae (120 (0.20%) by Demovir, dsDNA); the Virophage Family Lavidaviridae (15 (0.03%) by Demovir) were detected as well, but a majority of vOTUs were unclassified in genus-level.Fig. 1: The taxonomic assignment of LVD.Pie charts showing the affiliation of 56,870 vOTUs at family level assigned by script Demovir (a). and the affiliation of 1864 vOTUs at family level assigned by package vConTACT2 (b). Source data are provided in the Source Data file.Full size imageViral community structures differ across land use typesBray–Curtis dissimilarity of viral communities (median 0.9951) showed strong heterogeneity of viral communities among different sites (Fig. 2a). While, the Bray–Curtis dissimilarity (median: 0.5109) between paired viral communities of iVLPs and eVLPs from each site have a significant lower heterogeneity than inter-sites (Wilcox.test, p  0.05; Fig. 2b). Therefore, the paired iVLPs and eVLPs viromes from each site were merged for subsequently viral community analysis.Fig. 2: The macrodiversity of soil viral communities.a Boxplot showing Bray–Curtis dissimilarity of viral communities of intra-sites (between the corresponding community of iVLPs and eVLPs, n = 25) and inter-sites (between different sample sites, n = 300). The minima, maxima, center, bounds of box and whiskers in boxplots from bottom to top represented percentile 0, 10, 25, 50, 75, 90, and 100, respectively, the difference between different zones was tested using the two-sided Wilcox.test, ****p  More