Philippa Kaur

Snir, A. et al. The origin of cultivation and proto-weeds, long before Neolithic farming. PLoS ONE 10(7), e0131422. https://doi.org/10.1371/journal.pone.0131422 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Groman-Yaroslavski, I., Weiss, E. & Nadel, D. Composite sickles and cereal harvesting methods at 23,000-years-old Ohalo II Israel. PLoS ONE 11(11), e0167151. https://doi.org/10.1371/journal.pone.0167151 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Edwards, P. C. A 14000 year-old hunter-gatherer’s toolkit. Antiquity 81(314), 865–876. https://doi.org/10.1017/S0003598X0009596X (2007).Article 

Google Scholar 
Le Dosseur, G. Bone Objects in the Southern Levant from the Thirteenth to the Fourth Millennia. Bulletin du Centre de recherche français à Jérusalem 12, 111–125 (2003).
Google Scholar 
Garrard, A., & Yazbeck, C. The Natufian of Moghr el-Ahwal in the Qadisha valley, northern Lebanon. in Natufian Foragers in the Levant. International Monographs in Prehistory (eds. Bar-Yosef, O. & Valla, F. R.). 17–47. (Michigan, Ann Arbor, 2013).Belfer-Cohen, A. The Natufian in the Levant. Annu. Rev. Anthropol. 20, 167–186. https://doi.org/10.1146/annurev.an.20.100191.001123 (1991).Article 

Google Scholar 
Stordeur, D. Le Natoufien et son évolution à travers les artefacts en os in Natufian Foragers in the Levant. International Monographs in Prehistory (eds. Bar-Yosef, O. & Valla, F. R.). 457–482. (Michigan, Ann Arbor, 2013).Rosen, S. A. Lithics after the Stone Age: a handbook of stone tools from the Levant. (Rowman Altamira, 1997).Anderson, P. C. Prehistory of agriculture: new experimental and ethnographic approaches. (Cotsen Institute of Archaeology Press, 1999).Ibáñez, J. J., González Urquijo, J. E., & Rodríguez, A. The evolution of technology during the PPN in the Middle euphrates. A view from use wear analysis of lithic tools. in Systèmes techniques et communautés du Néolithique Préceramique au Proche Orient. Technical Systems and Near Eastern PPN Communities (eds. Astruc, L., Binder, D. & Briois, F.) 153–165 (Editions APDCA, 2007).Maeda, O., Lucas, L., Silva, F., Tanno, K. I. & Fuller, D. Q. Narrowing the harvest: Increasing sickle investment and the rise of domesticated cereal agriculture in the Fertile Crescent. Quatern. Sci. Rev. 145, 226–237. https://doi.org/10.1016/j.quascirev.2016.05.032 (2016).ADS 
Article 

Google Scholar 
Pichon, F. Exploitation of the cereals during the Pre-pottery Neolithic of Dja’de-el-Mughara: Preliminary results of the functional study of the glossy blades. Quatern. Int. 427, 138–151. https://doi.org/10.1016/j.quaint.2016.01.064 (2017).Article 

Google Scholar 
Borrell, F., & Molist, M. Projectile Points, Sickle Blades and Glossed Points. Tools and Hafting Systems at Tell Halula (Syria) during the 8th millennium cal. BC Paléorient, 33(2), 59–77 (2007). https://doi.org/10.2307/41496812.Douka, K., Efstratiou, N., Hald, M., Henriksen, P. & Karetsou, A. Dating Knossos and the arrival of the earliest Neolithic in the southern Aegean. Antiquity 91(356), 304–321. https://doi.org/10.15184/aqy.2017.29 (2017).Article 

Google Scholar 
Perlès, C. From the Near East to Greece: Let’s reverse the focus. Cultural elements that didn’t transfer. in How did farming reach Europe? (ed. Lichter, C.) 275–290 (Istanbul, Ege Yayınları, 2005).Gijn A.L. van & Wentink K. The role of flint in mediating identities: The microscopic evidence. in Mobilty, meaning & transformations of things, shifting contexts of material culture through time and space. (eds. Hahn, H.P. & Weiss, H.) 120–132 (Oxford, Oxbow Books, 2013).Guilaine, J. The neolithic transition: From the Eastern to the Western Mediterranean. in Times of Neolithic Transition along the Western Mediterranenn. (eds. O., García-Puchol & D. C., Salazar-García) 15–31 (New York, Springer, 2017). https://doi.org/10.1007/978-3-319-52939-4_2.Forenbaher, S. & Miracle, P. T. The spread of farming in the Eastern Adriatic. Antiquity 79(305), 514–528 (2005).Article 

Google Scholar 
Gabriele, M. et al. Long-distance mobility in the North-Western Mediterranean during the Neolithic transition using high resolution pottery sourcing. J. Archaeol. Sci. Rep. 28, 102050. https://doi.org/10.1016/j.jasrep.2019.102050 (2019).Article 

Google Scholar 
Manen, C., Perrin, T., Guilaine, J., Bouby, L., Bréhard, S., Briois, F., Durand, F., Marinval, P. & Vigne, J. D. The Neolithic transition in the western Mediterranean: A complex and non-linear diffusion process—the radiocarbon record revisited. Radiocarbon 61(2), 531–571 (2019). https://doi.org/10.1017/RDC.2018.98Ibáñez, J. J., Clemente Conte, I., Gassin, B., Gibaja, J. F., Gonzáles Urquijo, J. E., Márquez, B., Philibert, S., Rodríguez Rodríguez, A. Harvesting technology during the Neolithic in south-west Europe. in Prehistoric technology 40 years later: functional studies and the Russian legacy (eds. Longo L. & Skakun, N.) 183–95 (Oxford, Archaeopress, 2008).Gibaja, J. F., Ibáñez, J. J., González Urquijo, J. E. Neolithic Sickles in the Iberian Peninsula. in Exploring and Explaining Diversity in Agricultural Technology, EARTH 2 (eds. van Gijn, A., Whittaker, P. & Anderson, P.) 112–118 (Oxford, Oxbow Books, 2014).Mazzucco, N., Capuzzo, G., Petrinelli-Pannocchia, C., Ibáñez, J. J., Gibaja, J. F. Harvesting tools and the spread of the Neolithic into the Central-Western Mediterranean area. Quat. Int. 470(Part B), 511–528 (2018). https://doi.org/10.1016/j.quaint.2017.04.018.Mazzucco, N., Guilbeau, D., Kačar, S., Podrug, E., Forenbaher, S., Radić, D., Moore, A. T. M. The time is ripe for a change. The evolution of harvesting technologies in Central Dalmatia during the Neolithic period (6th millennium cal BC). J. Anthropol. Archaeol. 51, 88–103 (2018). https://doi.org/10.1016/j.jaa.2018.06.003Mazzucco, N. et al. Migration, adaptation, innovation: The spread of Neolithic harvesting technologies in the Mediterranean. PLoS ONE 15(4), e0232455. https://doi.org/10.1371/journal.pone.0232455 (2020).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Fugazzola Delpino, M. A., D’Eugenio, G. & Pessina, A. “La Marmotta” (Anguillara Sabazia, RM): Scavi 1989—un abitato perilacustre di età Neolitica. Bull. Paletnol. Ital. 84, 181–315 (1993).
Google Scholar 
Fugazzola Delpino, M. A., Pessina, A. Le village néolithique submergé de La Marmotta (lac de Bracciano, Rome). in Le Néolithique du Nord-Ouest méditerranéen (ed. Vaquer, J.) 35–38 (Société préhistorique française, Paris, 1999)Fugazzola Delpino, M. A. La Marmotta. in Le ceramiche impresse nel Neolitico antico. Italia e Mediterraneo (eds. Fugazzola, M.A., Pessina, A. & Tiné, V) 373–395 (Istituto Poligrafico e Zecca dello Stato, Roma, 2002).Grantham, G. L. faucille et la faux. Études rurales 151–152, 103–131 (1999).Article 

Google Scholar 
Sigaut, F. Identification des techniques de récolte des graines alimentaires. J. Agric. Trad. Bot. Appl. 25(3), 145–161 (1978).
Google Scholar 
Anderson, P. C., Sigaut, F. Introduction: reasons for variability in harvesting techniques and tools. in Exploring and Explaining Diversity in Agricultural Technology, EARTH 2 (eds. van Gijn, A., Whittaker, P. & Anderson, P.) 85–93 (Oxford, Oxbow Books, 2014).Halstead, P. Two oxen ahead: Pre-mechanized farming in the Mediterranean (John Wiley & Sons, 2014).Book 

Google Scholar 
Fugazzola Delpino, M. A. & Mineo, M. La piroga neolitica di Bracciano (La Marmotta 1). Bull. Paletnol. Ital. 86, 197–266 (1995).
Google Scholar 
Fugazzola Delpino, M. A., Tinazzi, O. Dati di cronologia da un villaggio del Neolitico Antico. Le indagini dendrocronologiche condotte sui legni de La Marmotta (lago di Bracciano-Roma). in Miscellanea in ricordo di Francesco Nicosia, 1–10 (Studia Erudita, Fabrizio Serra Editore, 2010).Salavert, A. et al. Direct dating reveals the early history of opium poppy in western Europe. Sci. Rep. 10, 20263. https://doi.org/10.1038/s41598-020-76924-3 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Ghiselli, L. et al. Nutritional characteristics of ancient Tuscan varieties of Triticum aestivum L. Ital. J. Agron. 11(4), 237–245 (2016).Article 

Google Scholar 
Pichon, F. Une moisson expérimentale de céréales, Séranon (août 2016), ArchéOrient – Le Blog, 14 octobre2016, (2016). https://archeorient.hypotheses.org/6667.Banks, W. E. & Kay, M. High-resolution casts for lithic use-wear analysis. Lithic Technol. 28(1), 27–34. https://doi.org/10.1080/01977261.2003.11721000 (2003).Article 

Google Scholar 
Ibáñez, J. J., Anderson, P. C., Gonzalez-Urquijo, J. & Gibaja, J. Cereal cultivation and domestication as shown by microtexture analysis of sickle gloss through confocal microscopy. J. Archaeol. Sci. 73, 62–81. https://doi.org/10.1016/j.jas.2016.07.011 (2016).Article 

Google Scholar 
Caruso Fermé, L. Modalidades de adquisición y uso del material leñoso entre grupos cazadores-recolectores patagónicos (Argentina). Métodos y técnicas de estudios del material leñoso arqueológico. PhD Dissertation (Universidad Autónoma de Barcelona, Barcelona, 2012).Caruso Fermé, L., Clemente, I., Civalero, M.T. A use-wear analysis of wood technology of patagonian hunter-gatherers. The case of Cerro Casa de Piedra 7, Argentina. J. Archaeol. Sci. 15, 315–321 (2015). https://doi.org/10.1016/j.jas.2015.03.015.Caruso Fermé, L., Aschero, C. Manufacturing and use of the wooden artifacts. A use-wear analysis of wood technology in hunter-gatherer groups (Cerro Casa de Piedra 7 site, Argentina). J. Archaeol. Sci. 31, 102291 (2020). https://doi.org/10.1016/j.quaint.2020.10.067.Schweingruber, F. H. Anatomy of European wood: An atlas for the identification of European trees, shrubs and dwarf shrubs (Paul Haupt, 1990).
Google Scholar 
Rageot, M. et al. Birch bark tar production: Experimental and biomolecular approaches to the study of a common and widely used prehistoric adhesive. J. Archaeol. Method Theory 26, 276–312. https://doi.org/10.1007/s10816-018-9372-4 (2019).Article 

Google Scholar 
Rageot, M. et al. New insights into Early Celtic consumption practices: Organic residue analyses of local and imported pottery from Vix-Mont Lassois. PLoS ONE 14(6), e0218001. https://doi.org/10.1371/journal.pone.0218001 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Arobba, D., Caramiello, R., Martino, G. P. Analisi paleobotaniche di resine dal relitto navale romano del Golfo Dianese. Rivista di Studi Liguri, LXIII-LXIV: 339–355 (1999).Marshall, D. M. Archaeological pollen: extraction from ancient resins. The American Association of Stratigraphic Palynologists. Prog. and Abstr., 38th Ann. Mtg., 34 (2005).Berglund, B. E., Ralska-Jasiewiczowa, M. Pollen analysis and pollen diagrams. in Handbook of Holocene Palaeoecology and Palaeohydrology. (eds. Berglund, B. E.) 455–484 (Chichester, Wiley, 1986).Traverse, A. Paleopalynology. Second Edition, 813 p. (Dordrecht, Springer, 2007).Punt W. (ed.) The Northwest European pollen flora (NEPF), vol. 2 (1980), vol. 3 (1981), vol. 4 (1984) vol. 5 (1988), vol. 6 (1991), vol. 7 (1996), vol. 8 (2003) (Elsevier, Wim Punt, Amsterdam, 1980–2003)Fægri, K. & Iversen, J. Textbook of pollen analysis (John Wiley and Sons, 1989).
Google Scholar 
Moore, P. D., Webb, J. A. & Collinson, M. E. Pollen analysis 2nd edn. (Blackwell, 1991).
Google Scholar 
Beug, H.-J. Leitfaden der Pollenbestimmung für Mitteleuropa und angrenzende Gebiete (Pfeil, 2004).
Google Scholar 
Reille, M. Pollen et spores d’Europe et d’Afrique du Nord. (Marseille, Laboratoire de Botanique Historique et Palynologie, 1992).Katz, O. et al. Rapid phytolith extraction for analysis of phytolith concentrations and assemblages during an excavation: An application at Tell es-Safi/Gath Israel. J. Archaeol. Sci. 37(7), 1557–1563. https://doi.org/10.1016/j.jas.2010.01.016 (2010).Article 

Google Scholar 
Brown, D. A. Prospects and limits of a phytolith key for grasses in the central United States. J. Archaeol. Sci. 11, 345–368. https://doi.org/10.1016/0305-4403(84)90016-5 (1984).Article 

Google Scholar 
Rosen, A. M. Preliminary identification of silica skeletons from Near Eastern archaeological sites: an anatomical approach. in Phytolith Systematics: Emerging Issues, Advances in Archaeological and Museum Science (eds. Rapp, G. Jr. & Mulholland, S. C.) 129–148 (New York, Plenum Press, 1992)Mulholland, S. C., Rapp Jr. G. A morphological classification of grass silica-bodies. in Phytolith Systematics: Emerging Issues, Advances in Archaeological and Museum Science (eds. Rapp, G. Jr. & Mulholland, S. C.) 65–89 (New York, Plenum Press, 1992)Piperno, D. R. Phytoliths: A comprehensive Guide for Archaeologists and Paleoecologists (Altamira Press, 2006).
Google Scholar 
Albert, R. M., & Weiner, S. Study of phytoliths in prehistoric ash layers from Kebara and Tabun caves using a quantitative approach. in Phytoliths: applications in earth sciences and human history, (eds. Meunier, J.D. & Colin, F.) 251–266 (Tokyo, Balkema Publisher, 2001)Albert, R. M. et al. Phytolith-rich layers from the Late Bronze and Iron Ages at Tel Dor (Israel): Mode of formation and archaeological significance. J. Archaeol. Sci. 35(1), 57–75. https://doi.org/10.1016/j.jas.2007.02.015 (2008).Article 

Google Scholar 
Albert, R. M., Ruíz, J. A. & Sans, A. PhytCore ODB: A new tool to improve efficiency in the management and exchange of information on phytoliths. J. Archaeol. Sci. 68, 98–105 (2016).Article 

Google Scholar 
Portillo, M., Kadowaki, S., Nishiaki, Y. & Albert, R. M. Early Neolithic household behavior at Tell Seker al-Aheimar (Upper Khabur, Syria): A comparison to ethnoarchaeological study of phytoliths and dung spherulites. J. Archaeol. Sci. 42, 107–118 (2014).Article 

Google Scholar 
Tsartsidou, G. et al. The phytolith archaeological record: strengths and weaknesses evaluated based on a quantitative modern reference collection from Greece. J. Archaeol. Sci. 34, 1262–1275. https://doi.org/10.1016/j.jas.2006.10.017 (2007).Article 

Google Scholar 
Neumann, K., Strömberg , A. E. C., Ball, T., Albert, R. M., Vrydaghs, L. Scott-Cummings, L. (International Committee for Phytolith Taxonomy ICPT). International Code for Phytolith Nomenclature (ICPN) 2.0. Annals of Botany, 124(2): 189–199 (2019).Anderson, P. C. Insight into plant harvesting and other activities at Hatoula, as revealed by microscopic functional analysis of selected chipped stone tools. Le site de Hatoula en Judée occidental. (eds. Lechevallier, M. & Ronen, A.) 277–293 (Paris, Association Paléorient, 1994)Fugazzola Delpino, M.A. La vita quotidiana del Neolitico. Il sito della Marmotta sul Lago di Bracciano. in Settemila anni fa il primo pane. Ambienti e culture delle società neolitiche (eds. Pessina, A. & Muscio G.) 185–192 (Udine, Museo Friulano di Storia Naturale, 1998–1999)Mineo, M. Monossili d’Europa: costruite anche per le rotte marine? in Ubi minor: le isole minori del Mediterraneo centrale: dal Neolitico ai primi contatti coloniali (eds. Guidi, A., Cazzella, A. & Nomi, F.). Scienze dell’Antichità 22, 453–475 (2016)Helwig, K., Monahan, V. & Poulin, J. The identification of hafting adhesive on a slotted antler point from a southwest Yukon ice patch. Am. Antiq. 73, 279–288. https://doi.org/10.1017/S000273160004227X (2008).Article 

Google Scholar 
Steigenberger, G. & Herm, C. Natural resins and balsams from an eighteenth-century pharmaceutical collection analysed by gas chromatography/mass spectrometry. Anal. Bioanal. Chem. 401, 1771–1784. https://doi.org/10.1007/s00216-011-5169-y (2011).CAS 
Article 
PubMed 

Google Scholar 
van den Berg, K. J., Boon, J. J., Pastorova, I. & Spetter, L. F. M. Mass spectrometric methodology for the analysis of highly oxidized diterpenoid acids in Old Master paintings. J. Mass Spectrom. 35, 512–533. https://doi.org/10.1002/(SICI)1096-9888(200004)35:4%3c512::AID-JMS963%3e3.0.CO;2-3 (2000).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Behre K. E. Anthropogenic Indicators in Pollen Diagrams, A.A. (Rotterdam, Balkema, 1986).Mercuri, A. M. et al. Anthropogenic Pollen Indicators (API) from archaeological sites as local evidence of human-induced environments in the Italian peninsula. Ann. Bot. 3, 143–153 (2013).
Google Scholar 
Andersen, S.-T., Identification of wild grass and cereal pollen. in Danmarks Geologiske Undersøgelse (ed. Aaby, B.) 69–92 (Geological Survey of Denmark, 1978).Bottema, S. Cereal-type pollen in the Near East as indicators of wild or domestic crops. in Préhistoire de l’agriculture: nouvelles approches expérimentales et ethnographiques (ed. Anderson P. C.) 95–106 (Paris, CRA, 1992). https://doi.org/10.1007/BF00217499.Lagerås, P. Long-term history of land-use and vegetation at Femtingagölen—a small lake in the Småland Uplands, southern Sweden. Veg. Hist. Archaeobot. 5, 215–228 (1996).Article 

Google Scholar 
Joly, C., Barillé, L., Barreau, M., Mancheron, A. & Visset, L. Grain and annulus diameter as criteria for distinguishing pollen grains of cereals from wild grasses. Rev. Palaeobot. Palynol. 146, 221–233. https://doi.org/10.1016/j.revpalbo.2007.04.003 (2007).Article 

Google Scholar 
Punt, W. Umbelliferae. Rev. Palaeobot. Palynol. 42, 155–364 (1984).Article 

Google Scholar 
Ellis, M. B. & Ellis, J. P. Microfungi of Land Plants. An Identification Handbook (London, Croom Helm, 1985) (Figure 1270).Ellis, M. B. & Ellis, J. P. Microfungi of Land Plants. An Identification Handbook (London, Croom Helm, 1985) (Figures 174; 176).Rottoli, M., Pessina, A. Neolithic agriculture in Italy: an update of archaeobotanical data with particular emphasis on northern settlements. in The Origins and Spread of Domestic Plants in Southwest Asia and Europe. (eds. Colledge, S. & Conolly, J.) 141–154 (Routledge, New York, 2016)Gurova, M. Prehistoric sickles in the collection of the National Museum of Archaeology in Sofia. in Southeast Europe and Anatolia in Prehistory: Essays in Honor of Vassil Nikolov on his 65th Anniversary (eds. Bacvarov, K. & Gleser, E.) 159–165 (Bonn, Verlag Dr. Rudolf Habelt GmbH, 2016)Sidéra, I. Nouveaux éléments d’origine proche-orientale dans le Néolithique ancien balkanique. in Analyse de l’industrie osseuse. in Préhistoire d’Anatolie. Genèse de deux mondes (ed. Otte, M.), 215–239 (Liège, ERAUL, 1997)Mellaart, J. Excavations at Hacılar: Fourth preliminary report, 1960. Anat. Stud. Anat. Stud. 11, 39–75 (1961).Article 

Google Scholar 
Nag, P. K., Goswami, A., Ashtekar, S. P. & Pradhan, C. K. Ergonomics in sickle operation. Appl. Ergon. 19(3), 233–239 (1988).CAS 
Article 

Google Scholar 
Astruc, L., Tkaya, M. B. & Torchy, L. D. l’efficacité des faucilles néolithiques au Proche-Orient: approche expérimentale. Bulletin de la Société préhistorique française 109(4), 671–687 (2012).Article 

Google Scholar 
Sigaut, F. Les techniques de récolte des grains : identification, localisation, problèmes d’interprétation. in Rites et rythmes agraires (ed. Cauvin, M.-C.) 31–43 (Lyon, Maison de l’Orient et de la Méditerranée Jean Pouilloux, 1991)Magri, D. Late Quaternary vegetation history at Lagaccione near Lago di Bolsena (central Italy). Rev. Palaeobot. Palynol. 106(3–4), 171–208 (1999).Article 

Google Scholar 
Gale, R., & Cutler, D. F. Plants in archaeology: identification manual of vegetative plant materials used in Europe and the Southern Mediterranean to c. 1500 (Westbury and Royal Botanic Gardens, Kew, 2000).Chabal, L. & Feugère, M. L. Le mobilier organique des puits antiques et autres contextes humides de Lattara. Lattara 18, 137–188 (2005).
Google Scholar 
Chabal, L. (ed.) Quatre puits de l’agglomération routière gallo-romaine d’Ambrussum (Villetelle, Hérault). Supplément. Revue Archéologique de Narbonnaise, 42: 65–71 (2013).Caruso Fermé, L. & Piqué Huerta, R. Landscape and forest exploitation at the ancient Neolithic site of La Draga (Banyoles, Spain). The Holocene, 24(3): 266 (2014).Boschian, G. Il Riparo “Ermanno de Pompeis” presso l’Eremo di San Bartolomeo di Legio. Scavi 1990–1999. in Atti della XXXVI Riunione Scientifica IIPP, Preistoria e Protostoria dell’Abruzzo, Chieti-Celano, 27–30 settembre 2001, 105–116 (IIPP; Firenze, 2003).Radi, G. & Danese, E. L’abitato di Colle Santo Stefano di Ortucchio (L’Aquila). in Atti della XXXVI Riunione Scientifica IIPP, Preistoria e Protostoria dell’Abruzzo, Chieti-Celano, 27–30 settembre 2001, 145–161 (IIPP; Firenze, 2003).De Francesco, A. M., Bocci, M., Crisci, G. M., & Francaviglia, V. Obsidian provenance at several Italian and Corsican archaeological sites using the non-destructive X-ray fluorescence method. in Obsidian and ancient manufactured glass (eds. Liritzis, I., & Stevenson, C. M.) 115–129 (Albuquerque, UNM Press, 2012).Degano, I. et al. Hafting of Middle Paleolithic tools in Latium (central Italy): New data from Fossellone and Sant’Agostino caves. PLoS ONE 14, e0213473 (2019).CAS 
Article 

Google Scholar 
Nardella, F. et al. Chemical investigations of bitumen from Neolithic archaeological excavations in Italy by GC/MS combined with principal component analysis. Anal. Methods 11, 1449–1459. https://doi.org/10.1039/c8ay02429d (2019).CAS 
Article 

Google Scholar 
Rageot, M. et al. Management systems of adhesive materials throughout the Neolithic in the North-West Mediterranean. J. Archaeol. Sci. 126, 105309 (2021).Article 

Google Scholar 
Binder, D., Bourgeois, G., Benoist, F. & Vitry, C. Identification de brai de bouleau (betula) dans le néolithique de Giribaldi (Nice, France) par la spectrométrie de masse. Revue d’Archéométrie 14, 37–42 (1990).Article 

Google Scholar 
Vuorela, I. Relative pollen rain around cultivated fields. Acta Bot. Fenn. 102, 1–27 (1973).
Google Scholar 
Robinson, M. & Hubbard, R. N. L. B. The transport of pollen in the bracts of hulled cereals. J. Archaeol. Sci. 4(2), 197–199. https://doi.org/10.1016/0305-4403(77)90067-X (1977).Article 

Google Scholar 
Hall, V.A., The role of harvesting techniques in the dispersal of pollen grains of Cerealia. Pollen et Spores, XXX, 2, pp. 265–270.Portillo, M., Llergo, Y., Ferrer, A. & Albert, R. M. Tracing microfossil residues of cereal processing in the archaeobotanical record: an experimental approach. Veg. Hist. Archaeobot. 26(1), 59–74. https://doi.org/10.1007/s00334-016-0571-1 (2017).Article 

Google Scholar 
Negri, G. Nuovo erbario figurato (Hoepli ed., Milano, 1981).Paris R. R. & Moyse H. Matière Médicale. Vol 2°, (Masson, Paris. 1976).Bulgarelli, G. & Flamigni, S. Le piante tossiche e velenose (Hoepli ed., Milano, 2010).Les, D. H. Aquatic Dicotyledons of North America: Ecology, Life History, and Systematics (CRC Press, 2017).Book 

Google Scholar 
Curti, L. Herbarium, un’inedita collezione di piante del XVIII secolo conservata presso l’orto Botanico dell’Università di Padova (Offset Invicta S.p.A., Padova, 1992).Rottoli, M. Zafferanone selvatico (Carthamus lanatus) e cardo della Madonna (Silybum marianum), piante raccolte o coltivate nel Neolitico antico a “La Marmotta”? Bollettino di Paletnologia Italiana, 91–92, 47–61 (2000–2001).Rottoli, M. “La Marmotta”, Anguillara Sabazia (RM), scavi 1989. Analisi paletnobotaniche: prime risultanze. Bullettino di Paletnologia Italiana 84, 305–315 (1993).Van Geel, B. Non-pollen palynomorphs. in Tracking Environmental Change Using Lake Sediments: Terrestrial, vol. 3. (ed. Smol, J. P., Birks, H. J. B., Last W. M.) 99–119 (Algal and Siliceous Indicators, New York, 2001)Hawksworth, David L., van Geel, Bas, Wiltshire, Patricia E. J. The enigma of the Diporotheca palynomorph. Rev. Palaeobot. Palynol. 235, 94–98 (2016). https://doi.org/10.1016/j.revpalbo.2016.09.010.Krug, J. C., Benny, G. L., Keller, H. W. Coprophilous fungi. In Biodiversity of Fungi. Inventory and Monitoring Methods (ed. Foster M., Bill, G.) 467–499 (Elsevier Science, Amsterdam, 2004). More

Jiang, D. B., Lang, X. M., Tian, Z. P. & Wang, T. Considerable model-data mismatch in temperature over China during the mid-Holocene: results of PMIP simulations. J. Clim. 25, 4135–4153 (2012).ADS 
Article 

Google Scholar 
Marcott, S. A., Shakun, J. D., Clark, P. U. & Mix, A. C. A reconstruction of regional and global temperature for the past 11300 years. Science 339, 1198–1201 (2013).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Liu, Z. et al. The Holocene temperature conundrum. Proc. Natl Acad. Sci. USA 111, E3501–E3505 (2014).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Marsicek, J., Shuman, B., Bartlein, P., Shafer, S. L. & Brewer, S. Reconciling divergent trends and millennial variations in Holocene temperatures. Nature 554, 92–96 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Affolter, S., Huselmann, A., Fleitmann, D., Edwards, R. L. & Leuenberger, M. Central Europe temperature constrained by speleothem fluid inclusion water isotopes over the past 14,000 years. Sci. Adv. 5, eaav3809 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Bova, S., Rosenthal, Y., Liu, Z. & Yan, M. Seasonal origin of the thermal maxima at the Holocene and the last interglacial. Nature 589, 548–553 (2021).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Meyer, H. et al. Long-term winter warming trend in the Siberian Arctic during the mid- to late Holocene. Nat. Geosci. 8, 122–125 (2015).ADS 
CAS 
Article 

Google Scholar 
Baker, J., Lachniet, M., Chervyatsova, O., Asmerom, Y. & Polyak, V. J. Holocene warming in western continental Eurasia driven by glacial retreat and greenhouse forcing. Nat. Geosci. 10, 430–435 (2017).ADS 
CAS 
Article 

Google Scholar 
Mann, M., Schmidt, G., Miller, S. & LeGrande, A. Potential biases in inferring Holocene temperature trends from long-term borehole information. Geophys. Res. Lett. 36, L05708 (2009).ADS 
Article 

Google Scholar 
Liu, T. S. Loess and the Environment (In Chinese). (China Ocean Press, Beijing, 1985).Rousseau, D. D. & Wu, N. Q. A new molluscan record of the monsoon variability over the past 130 000 yr in the Luochuan loess sequence, China. Geology 25, 275–278 (1997).ADS 
Article 

Google Scholar 
Wu, N. Q., Li, F. J. & Rousseau, D. D. Terrestrial mollusk records from Chinese loess sequences and changes in the East Asian monsoonal environment. J. Asian Earth Sci. 155, 35–48 (2018).ADS 
Article 

Google Scholar 
Qian, L. Q. Climate of Loess Plateau (in Chinese). (China Meteorological Press, Beijing, 1991).Chen, D. & Gao, J. Economic Fauna Sinica of China: Terrestrial Mollusca (in Chinese). (Science Press, Beijing, 1987).Proćków, M., Drvotová, M., Juřičková, L. & Kuźnik-Kowalska, E. Field and laboratory studies on the life-cycle, growth and feeding preference in the hairy snail Trochulus hispidus (L., 1758) (Gastropoda: Pulmonata: Hygromiidae). Biologia 68, 131–141 (2013).Article 

Google Scholar 
Rousseau, D. Climatic transfer function from Quaternary molluscs in European loess deposits. Quat. Res. 36, 195–209 (1991).Article 

Google Scholar 
Rousseau, D., Preece, R. & Limondin-Lozouet, N. British late glacial and Holocene climatic history reconstructed from land snail assemblages. Geology 26, 651–654 (1998).ADS 
Article 

Google Scholar 
Wang, Z. H., Fang, J. Y., Tang, Z. Y. & Lin, X. Patterns, determinants and models of woody plant diversity in China. Proc. R. Soc. B 278, 2122–2132 (2011).PubMed 
Article 

Google Scholar 
Gu, Z. Y., Liu, Z. X., Xu, B. & Wu, N. Q. Stable carbon and oxygen isotopes in land snail carbonate shells from a last glacial loess sequence and their implications of environmental changes (in Chinese). Quat. Sci. 29, 13–22 (2009).CAS 

Google Scholar 
Sun, X. H., Gu, Z. Y. & Xu, B. Oxygen isotopic variations in the shells collected monthly from a live species of land snails at local in Zhenjiang, Jiangsu Province, China (in Chinese). Quat. Sci. 29, 976–980 (2009).CAS 

Google Scholar 
Huang, L., Wu, N., Gu, Z. & Chen, X. Variability of snail growing season at the Chinese Loess Plateau during the last 75 ka. Chin. Sci. Bull. 57, 1036–1045 (2012).CAS 
Article 

Google Scholar 
Dong, Y. J. et al. Paleorecords reveal the increased temporal instability of species diversity under biodiversity loss. Quat. Sci. Rev. 269, 107147 (2021).Article 

Google Scholar 
Horsák, M. Mollusc assemblages in palaeoecological reconstructions: an investigation of their predictive power using transfer function models. Boreas 40, 459–467 (2011).Article 

Google Scholar 
Sümegi, P. & Gulyás, S. Some notes on the interpretation and reliability of malacological proxies in paleotemperature reconstructions from loess- comments to Obreht et al.‘s “A critical reevaluation of paleoclimate proxy records from loess in the Carpathian Basin”. Earth-Sci. Rev. 221, 103675 (2021).Samartin, S. et al. Warm Mediterranean mid-Holocene summers inferred from fossil midge assemblages. Nat. Geosci. 10, 207–212 (2017).ADS 
CAS 
Article 

Google Scholar 
Seppä, H., Birks, H., Odland, A., Poska, A. & Veski, S. A modern pollen-climate calibration set from northern Europe: developing and testing a tool for palaeoclimatological reconstructions. J. Biogeogr. 31, 251–267 (2004).Article 

Google Scholar 
Allen, J. et al. Rapid environmental changes in southern Europe during the last glacial period. Nature 400, 740–743 (1999).ADS 
CAS 
Article 

Google Scholar 
Rioual, P. et al. High-resolution record of climate stability in France during the last interglacial period. Nature 413, 293–296 (2001).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Salonen, J. S. et al. Abrupt high-latitude climate events and decoupled seasonal trends during the Eemian. Nat. Commun. 9, 2851 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Lu, H. et al. Seasonal climatic variation recorded by phytolith assemblages from Baoji loess sequence in central China over the last 150000 a. Sci. China, Ser. D. 26, 629–639 (1996).
Google Scholar 
Telford, R. J. & Birks, H. J. B. A novel method for assessing the statistical significance of quantitative reconstructions inferred from biotic assemblages. Quat. Sci. Rev. 30, 1272–1278 (2011).ADS 
Article 

Google Scholar 
Lu, H. Y., Wu, N. Q., Liu, K. B., Jiang, H. & Liu, T. S. Phytoliths as quantitative indicators for the reconstruction of past environmental conditions in China II: palaeoenvironmental reconstruction in the Loess Plateau. Quat. Sci. Rev. 26, 759–772 (2007).ADS 
Article 

Google Scholar 
Sun, J. M., Diao, G. Y., Wen, Q. Z. & Zhou, H. Y. A preliminary study on quantitative estimate of Palaeoclimate by using geochemical transfer function in the Loess Plateau (In Chinese). Geochimica 28, 265–272 (1999).CAS 

Google Scholar 
Wen, R. et al. Pollen–climate transfer functions intended for temperate eastern Asia. Quat. Int. 311, 3–11 (2013).Article 

Google Scholar 
Xu, Q., Xiao, J., Li, Y., Tian, F. & Nakagawa, T. Pollen-based quantitative reconstruction of Holocene climate changes in the Daihai lake area, inner Mongolia, China. J. Clim. 23, 2856–2868 (2010).ADS 
Article 

Google Scholar 
Li, J. et al. Quantitative Holocene climatic reconstructions for the lower Yangtze region of China. Clim. Dyn. 50, 1101–1113 (2018).Article 

Google Scholar 
Nakagawa, T., Tarasov, P. E., Nishida, K., Gotanda, K. & Yasuda, Y. Quantitative pollen-based climate reconstruction in central Japan: application to surface and Late Quaternary spectra. Quat. Sci. Rev. 21, 2099–2113 (2002).ADS 
Article 

Google Scholar 
Chen, M.-T. et al. Dynamic millennial-scale climate changes in the northwestern Pacific over the past 40,000 years. Geophys. Res. Lett. 37, L23603 (2010).ADS 

Google Scholar 
Sun, Y., Oppo, D. W., Xiang, R., Liu, W. & Gao, S. Last deglaciation in the Okinawa Trough: Subtropical northwest Pacific link to Northern Hemisphere and tropical climate. Paleoceanography 20, PA4005 (2005).ADS 
Article 

Google Scholar 
de Garidel-Thoron, T. et al. A multiproxy assessment of the western equatorial Pacific hydrography during the last 30 kyr. Paleoceanography 22, PA3204 (2007).ADS 

Google Scholar 
Chen, F., Duan, Y. & Hou, J. An 88 ka temperature record from a subtropical lake on the southeastern margin of the Tibetan Plateau (third pole): new insights and future perspectives. Sci. Bull. 66, 1056–1057 (2021).Article 

Google Scholar 
James, R. P. & Arguez, A. On the estimation of daily climatological temperature variance. J. Atmos. Ocean. Tech. 32, 2297–2304 (2015).Article 

Google Scholar 
Mearns, L. O., Rosenzweig, C. & Goldberg, R. Mean and variance change in climate scenarios: methods, agricultural applications, and measures of uncertainty. Climatic Change 35, 367–396 (1997).Article 

Google Scholar 
Laskar, J., Fienga, A., Gastineau, M. & Manche, H. La2010: a new orbital solution for the long-term motion of the Earth. Astron. Astrophys. 532, A89 (2011).ADS 
MATH 
Article 

Google Scholar 
Bader, J. et al. Global temperature modes shed light on the Holocene temperature conundrum. Nat. Commun. 11, 4726 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Liu, Y. et al. A possible role of dust in resolving the Holocene temperature conundrum. Sci. Rep. 8, 4434 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Thompson, A. J., Zhu, J., Poulsen, C. J., Tierney, J. E. & Skinner, C. B. Northern Hemisphere vegetation change drives a Holocene thermal maximum. Sci. Adv. 8, eabj6535 (2022).PubMed 
PubMed Central 
Article 

Google Scholar 
Lambert, F. et al. The role of mineral-dust aerosols in polar temperature amplification. Nat. Clim. Chang. 3, 487–491 (2013).ADS 
Article 

Google Scholar 
Xu, Y., Wang, H., Liao, H. & Jiang, D. Simulation of the direct radiative effect of mineral dust aerosol on the climate at the last glacial maximum. J. Clim. 24, 843–858 (2011).ADS 
Article 

Google Scholar 
Wohlfahrt, J., Harrison, S. P. & Braconnot, P. Synergistic feedbacks between ocean and vegetation on mid- and high-latitude climates during the mid-Holocene. Clim. Dynam. 22, 223–238 (2004).ADS 
Article 

Google Scholar 
Jahn, A., Claussen, M., Ganopolski, A. & Brovkin, V. Quantifying the effect of vegetation dynamics on the climate of the Last Glacial Maximum. Clim. Past 1, 1–7 (2005).Article 

Google Scholar 
Braconnot, P., Joussaume, S., Marti, O. & de Noblet, N. Synergistic feedbacks from ocean and vegetation on the African monsoon response to mid-Holocene insolation. Geophys. Res. Lett. 26, 2481–2484 (1999).ADS 
Article 

Google Scholar 
Braconnot, P. et al. Evaluation of climate models using palaeoclimatic data. Nat. Clim. Change 2, 417–424 (2012).ADS 
Article 

Google Scholar 
Zhang, X. & Chen, F. Non-trivial role of internal climate feedback on interglacial temperature evolution. Nature 600, E1–E3 (2021).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Li, F. J. et al. Quantitative distribution and calculation of ecological amplitude of land snail Metodontia in the Chinese Loess Plateau and adjacent regions (In Chinese with English abstract). Quat. Sci. 36, 564–574 (2016).
Google Scholar 
Dong, Y. J. et al. Influence of monsoonal water-energy dynamics on terrestrial mollusk species-diversity gradients in northern China. Sci. Total Environ. 676, 206–214 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Dong, Y. J. et al. Anthropogenic modification of soil communities in northern China for at least two millennia: Evidence from a quantitative mollusk approach. Quat. Sci. Rev. 248, 106579 (2020).Article 

Google Scholar 
Cameron, R. A. D. & Pokryszko, B. M. Estimating the species richness and composition of land mollusc communities: Problems, consequences and practical advice. J. Conchol. 38, 529–547 (2005).
Google Scholar 
Dong, Y., Wu, N., Li, F., Huang, L. & Wen, W. Time-transgressive nature of the magnetic susceptibility record across the Chinese Loess Plateau at the Pleistocene/Holocene transition. PLoS One 10, e0133541 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Ter Braak, C. J. F. & Smilauer, P. CANOCO Reference Manual and CanoDraw for Windows User’s Guide: Software for Canonical Community Ordination (Version 4.5). (Microcomputer Power, New York, 2002).Ter Braak, C. J. F. & Smilauer, P. CANOCO Reference Manual and User’s Guide: Software for Ordination (Version 5.0). (Microcomputer Power, New York, 2012).Ter Braak, C. J. F. & Juggins, S. Weighted averaging partial least squares regression (WA-PLS): an improved method for reconstructing environmental variables from species assemblages. Hydrobiologia 269, 485–502 (1993).Article 

Google Scholar 
Birks, H. J. B., Lotter, A. F., Juggins, S. & Smol, J. P. Tracking Environmental Change Using Lake Sediments Volume 5: Data Handling and Numerical Techniques. p. 123–141. (Springer, London, 2012).Ter Braak, C. J. F. Canonical community ordination. Part I: Basic theory and linear methods. Ecoscience 1, 127–140 (1994).Article 

Google Scholar 
Juggins, S. C2 data analysis (version 1.7.4). (Newcastle University, Newcastle, 2011).Juggins, S. Rioja: Analysis of Quaternary Science Data. R package version 0.9-21 http://cran.r-project.org/package=rioja (2017).Simpson, G. L. & Oksanen, J. Analogue: Analogue matching and Modern Analogue. Technique Transfer Function Models. R package version 0.17-4 https://cran.r-project.org/package=analogue (2020).Telford, R. J. palaeoSig: Significance Tests of Quantitative Palaeoenvironmental Reconstructions. R package version 2.0-3 http://cran.r-project.org/package=palaeoSig (2019).Hansen, J., Sato, M. & Ruedy, R. Perception of climate change. Proc. Natl Acad. Sci. USA 109, E2415–E2423 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Olonscheck, D., Schurer, A. P., Lücke, L. & Hegerl, G. C. Large-scale emergence of regional changes in year-to-year temperature variability by the end of the 21st century. Nat. Commun. 12, 7237 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Liu, Y., Ren, G., Kang, H. & Sun, X. A significant bias of Tmax and Tmin average temperature and its trend. J. Appl. Meteorol. Clim. 58, 2235–2246 (2019).ADS 
Article 

Google Scholar 
Parey, S., Dacunha-Castelle, D. & Hoang, T. T. H. Mean and variance evolutions of the hot and cold temperatures in Europe. Clim. Dyn. 34, 345–359 (2010).Article 

Google Scholar 
Dong, Y. SeaTemCon_R code for “The Holocene temperature conundrum answered by mollusk records from East Asia”. Zenodo https://doi.org/10.5281/zenodo.6426798 (2022).Article 

Google Scholar 
Dong, Y. Data repository for “The Holocene temperature conundrum answered by mollusk records from East Asia”. Zenodo https://doi.org/10.5281/zenodo.6426911 (2022).Article 

Google Scholar  More

Model equationsOur host-parasite mathematical model involves the following host population components: ‘susceptible’ hosts denoted by (S), and hosts infected by k distinct types of parasites ((k=1,2,…,n)), the corresponding population numbers of infected hosts are denoted by (I_{i_1,i_2,…,i_k}), where each index (i_j) can take a value from 1, …, n (to avoid repeated counting of the same infection configuration, we require throughout the paper that (i_1 More

Smith, D. C., Simon, M., Alldredge, A. L. & Azam, F. Intense hydrolytic enzyme activity on marine aggregates and implications for rapid particle dissolution. Nature 359, 139–142. https://doi.org/10.1038/359139a0 (1992).ADS 
CAS 
Article 

Google Scholar 
Alldredge, A. L. & Gotschalk, C. C. Direct observations of the mass flocculation of diatom blooms: Characteristics, settling velocities and formation of diatom aggregates. Deep Sea Res. A 36, 159–171. https://doi.org/10.1016/0198-0149(89)90131-3 (1989).ADS 
CAS 
Article 

Google Scholar 
Jackson, G. A. A model of the formation of marine algal flocs by physical coagulation processes. Deep Sea Res. A 37, 1197–1211. https://doi.org/10.1016/0198-0149(90)90038-w (1990).ADS 
CAS 
Article 

Google Scholar 
Kiørboe, T., Lundsgaard, C., Olesen, M. & Hansen, J. L. S. Aggregation and sedimentation processes during a spring phytoplankton bloom: A field experiment to test coagulation theory. J. Mar. Res. 52, 297–323. https://doi.org/10.1357/0022240943077145 (1994).Article 

Google Scholar 
Jackson, G. Coagulation Theory and Models of Oceanic Plankton Aggregation (CRC Press, 2005).
Google Scholar 
Grossart, H. P., Kiorboe, T., Tang, K. & Ploug, H. Bacterial colonization of particles: Growth and interactions. Appl. Environ. Microb. 69, 3500–3509. https://doi.org/10.1128/aem.69.6.3500-3509.2003 (2003).ADS 
CAS 
Article 

Google Scholar 
Kiorboe, T., Tang, K., Grossart, H. P. & Ploug, H. Dynamics of microbial communities on marine snow aggregates: Colonization, growth, detachment, and grazing mortality of attached bacteria. Appl. Environ. Microbiol. 69, 3036–3047. https://doi.org/10.1128/AEM.69.6.3036 (2003).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Martin, J. H., Knauer, G. A., Karl, D. M. & Broenkow, W. W. VERTEX: Carbon cycling in the northeast pacific. Deep Sea Res. A 34, 267–285. https://doi.org/10.1016/0198-0149(87)90086-0 (1987).ADS 
CAS 
Article 

Google Scholar 
Buesseler, K. O. et al. VERTIGO (vertical transport in the global ocean): A study of particle sources and flux attenuation in the North Pacific. Deep Sea Res. II 55, 1522–1539. https://doi.org/10.1016/j.dsr2.2008.04.024 (2008).ADS 
Article 

Google Scholar 
Grossart, H. P., Tang, K. W., Kiorboe, T. & Ploug, H. Comparison of cell-specific activity between free-living and attached bacteria using isolates and natural assemblages. FEMS Microbiol. Lett. 266, 194–200. https://doi.org/10.1111/j.1574-6968.2006.00520.x (2007).CAS 
Article 
PubMed 

Google Scholar 
Martinez, J., Smith, D. C., Steward, G. F. & Azam, F. Variability in ectohydrolytic enzyme activities of pelagic marine bacteria and its significance for substrate processing in the sea. Aquat. Microb. Ecol. 10, 223–230. https://doi.org/10.3354/ame010223 (1996).Article 

Google Scholar 
Kellogg, C. T. E. et al. Evidence for microbial attenuation of particle flux in the Amundsen Gulf and Beaufort Sea: Elevated hydrolytic enzyme activity on sinking aggregates. Polar Biol. 34, 2007–2023. https://doi.org/10.1007/s00300-011-1015-0 (2011).Article 

Google Scholar 
Jiao, N. et al. Microbial production of recalcitrant dissolved organic matter: Long-term carbon storage in the global ocean. Nat. Rev. Microbiol. 8, 593–599. https://doi.org/10.1038/nrmicro2386 (2010).CAS 
Article 
PubMed 

Google Scholar 
Jiao, N. & Zheng, Q. The microbial carbon pump: From genes to ecosystems. Appl. Environ. Microbiol. 77, 7439–7444. https://doi.org/10.1128/AEM.05640-11 (2011).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Buchan, A., LeCleir, G. R., Gulvik, C. A. & Gonzalez, J. M. Master recyclers: Features and functions of bacteria associated with phytoplankton blooms. Nat. Rev. Microbiol. 12, 686–698. https://doi.org/10.1038/nrmicro3326 (2014).CAS 
Article 
PubMed 

Google Scholar 
Smriga, S., Fernandez, V. I., Mitchell, J. G. & Stocker, R. Chemotaxis toward phytoplankton drives organic matter partitioning among marine bacteria. Proc. Natl. Acad. Sci. USA 113, 1576–1581. https://doi.org/10.1073/pnas.1512307113 (2016).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Secchi, E. et al. The effect of flow on swimming bacteria controls the initial colonization of curved surfaces. Nat. Commun. 11, 2851. https://doi.org/10.1038/s41467-020-16620-y (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Acinas, S. G., Antón, J. & Rodríguez-Valera, F. Diversity of free-living and attached bacteria in offshore Western Mediterranean Waters as depicted by analysis of genes encoding 16S rRNA. Appl. Environ. Microb. 65, 514–522 (1999).ADS 
CAS 
Article 

Google Scholar 
Grossart, H. P., Levold, F., Allgaier, M., Simon, M. & Brinkhoff, T. Marine diatom species harbour distinct bacterial communities. Environ. Microbiol. 7, 860–873. https://doi.org/10.1111/j.1462-2920.2005.00759.x (2005).CAS 
Article 
PubMed 

Google Scholar 
Mestre, M. et al. Sinking particles promote vertical connectivity in the ocean microbiome. Proc. Natl. Acad. Sci. USA 115, E6799–E6807. https://doi.org/10.1073/pnas.1802470115 (2018).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Rieck, A., Herlemann, D. P., Jurgens, K. & Grossart, H. P. Particle-associated differ from free-living bacteria in surface waters of the Baltic Sea. Front. Microbiol. 6, 1297. https://doi.org/10.3389/fmicb.2015.01297 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Ziervogel, K., Steen, A. D. & Arnosti, C. Changes in the spectrum and rates of extracellular enzyme activities in seawater following aggregate formation. Biogeosciences 7, 1007–1015. https://doi.org/10.5194/bg-7-1007-2010 (2010).ADS 
CAS 
Article 

Google Scholar 
Stocker, R., Seymour, J. R., Samadani, A., Hunt, D. E. & Polz, M. F. Rapid chemotactic response enables marine bacteria to exploit ephemeral microscale nutrient patches. Proc. Natl. Acad. Sci. USA 105, 4209–4214. https://doi.org/10.1073/pnas.0709765105 (2008).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Lopez-Perez, M. et al. Genomes of surface isolates of Alteromonas macleodii: The life of a widespread marine opportunistic copiotroph. Sci. Rep. 2, 696. https://doi.org/10.1038/srep00696 (2012).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Thiele, S., Fuchs, B. M., Amann, R. & Iversen, M. H. Colonization in the photic zone and subsequent changes during sinking determine bacterial community composition in marine snow. Appl. Environ. Microbiol. 81, 1463–1471. https://doi.org/10.1128/AEM.02570-14 (2015).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Bachmann, J. et al. Environmental drivers of free-living vs particle-attached bacterial community composition in the mauritania upwelling system. Front. Microbiol. 9, 2836. https://doi.org/10.3389/fmicb.2018.02836 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Kirchman, D. The ecology of Cytophaga-Flavobacteria in aquatic environments. FEMS Microbiol. Ecol. 39, 91–100. https://doi.org/10.1016/s0168-6496(01)00206-9 (2002).CAS 
Article 
PubMed 

Google Scholar 
Bizic-Ionescu, M. et al. Comparison of bacterial communities on limnic versus coastal marine particles reveals profound differences in colonization. Environ. Microbiol. 17, 3500–3514. https://doi.org/10.1111/1462-2920.12466 (2015).CAS 
Article 
PubMed 

Google Scholar 
Zhao, Z., Baltar, F. & Herndl, G. J. Linking extracellular enzymes to phylogeny indicates a predominantly particle-associated lifestyle of deep-sea prokaryotes. Sci. Adv. 6, 4354. https://doi.org/10.1126/sciadv.aaz4354 (2020).ADS 
CAS 
Article 

Google Scholar 
Baumas, C. M. J. et al. Mesopelagic microbial carbon production correlates with diversity across different marine particle fractions. ISME J. 15, 1695–1708. https://doi.org/10.1038/s41396-020-00880-z (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Ploug, H., Grossart, H. P., Azam, F. & Jørgensen, B. B. Photosynthesis, respiration, and carbon turnover in sinking marine snow from surface waters of Southern California Bight: Implications for the carbon cycle in the ocean. Mar. Ecol. Prog. Ser. 179, 1–11. https://doi.org/10.3354/meps179001 (1999).ADS 
CAS 
Article 

Google Scholar 
Ploug, H. & Grossart, H.-P. Bacterial growth and grazing on diatom aggregates: Respiratory carbon turnover as a function of aggregate size and sinking velocity. Limnol. Oceanogr. 45, 1467–1475. https://doi.org/10.4319/lo.2000.45.7.1467 (2000).ADS 
CAS 
Article 

Google Scholar 
Ebrahimi, A., Schwartzman, J. & Cordero, O. X. Cooperation and spatial self-organization determine rate and efficiency of particulate organic matter degradation in marine bacteria. Proc. Natl. Acad. Sci. USA 116, 23309–23316. https://doi.org/10.1073/pnas.1908512116 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Grossart, H.-P. & Ploug, H. Microbial degradation of organic carbon and nitrogen on diatom aggregates. Limnol. Oceanogr. 46, 267–277. https://doi.org/10.4319/lo.2001.46.2.0267 (2001).ADS 
CAS 
Article 

Google Scholar 
Datta, M. S., Sliwerska, E., Gore, J., Polz, M. F. & Cordero, O. X. Microbial interactions lead to rapid micro-scale successions on model marine particles. Nat. Commun. 7, 11965. https://doi.org/10.1038/ncomms11965 (2016).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Kiorboe, T., Grossart, H. P., Ploug, H. & Tang, K. Mechanisms and rates of bacterial colonization of sinking aggregates. Appl. Environ. Microbiol. 68, 3996–4006. https://doi.org/10.1128/AEM.68.8.3996-4006.2002 (2002).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Vaqué, D., Duarte, C. M. & Marrasé, C. Influence of algal population dynamics on phytoplankton colonization by bacteria: Evidence from two diatom species. Mar. Ecol. Prog. Ser. 65, 201–203. https://doi.org/10.3354/meps065201 (1990).ADS 
Article 

Google Scholar 
Grossart, H.-P. & Ploug, H. Bacterial production and growth efficiencies: Direct measurements on riverine aggregates. Limnol. Oceanogr. 45, 436–445. https://doi.org/10.4319/lo.2000.45.2.0436 (2000).ADS 
CAS 
Article 

Google Scholar 
Duhamel, S. et al. Growth and specific P-uptake rates of bacterial and phytoplanktonic communities in the Southeast Pacific (BIOSOPE cruise). Biogeosciences 4, 941–956. https://doi.org/10.5194/bg-4-941-2007 (2007).ADS 
Article 

Google Scholar 
Kirchman, D. L. Growth rates of microbes in the oceans. Annu. Rev. Mar. Sci. 8, 285–309. https://doi.org/10.1146/annurev-marine-122414-033938 (2016).ADS 
Article 

Google Scholar 
Brumley, D. R. et al. Cutting through the noise: Bacterial chemotaxis in marine microenvironments. Front. Mar. Sci. 7, 527. https://doi.org/10.3389/fmars.2020.00527 (2020).Article 

Google Scholar 
Thomas, T. et al. Analysis of the Pseudoalteromonas tunicata genome reveals properties of a surface-associated life style in the marine environment. PLoS ONE 3, e3252. https://doi.org/10.1371/journal.pone.0003252 (2008).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Varbanets, L. D. et al. The black sea bacteria-producers of hydrolytic enzymes. Mikrobiol. Z. 73, 9–15 (2011).CAS 
PubMed 

Google Scholar 
Sapp, M. et al. Species-specific bacterial communities in the phycosphere of microalgae?. Microb. Ecol. 53, 683–699. https://doi.org/10.1007/s00248-006-9162-5 (2007).Article 
PubMed 

Google Scholar 
Sarmento, H. & Gasol, J. M. Use of phytoplankton-derived dissolved organic carbon by different types of bacterioplankton. Environ. Microbiol. 14, 2348–2360. https://doi.org/10.1111/j.1462-2920.2012.02787.x (2012).CAS 
Article 
PubMed 

Google Scholar 
Gram, L., Grossart, H. P., Schlingloff, A. & Kiorboe, T. Possible quorum sensing in marine snow bacteria: Production of acylated homoserine lactones by Roseobacter strains isolated from marine snow. Appl. Environ. Microbiol. 68, 4111–4116. https://doi.org/10.1128/AEM.68.8.4111 (2002).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Arandia-Gorostidi, N. et al. Warming the phycosphere: Differential effect of temperature on the use of diatom-derived carbon by two copiotrophic bacterial taxa. Environ. Microbiol. 22, 1381–1396. https://doi.org/10.1111/1462-2920.14954 (2020).CAS 
Article 
PubMed 

Google Scholar 
Sarmento, H., Morana, C. & Gasol, J. M. Bacterioplankton niche partitioning in the use of phytoplankton-derived dissolved organic carbon: Quantity is more important than quality. ISME J 10, 2582–2592. https://doi.org/10.1038/ismej.2016.66 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Grossart, H. P. & Simon, M. Bacterial colonization and microbial decomposition of limnetic organic aggregates (lake snow). Aquat. Microb. Ecol. 15, 127–140. https://doi.org/10.3354/ame015127 (1998).Article 

Google Scholar 
Kiørboe, T. & Jackson, G. A. Marine snow, organic solute plumes, and optimal chemosensory behavior of bacteria. Limnol. Oceanogr. 46, 1309–1318. https://doi.org/10.4319/lo.2001.46.6.1309 (2001).ADS 
Article 

Google Scholar 
Chakraborty, S. et al. Quantifying nitrogen fixation by heterotrophic bacteria in sinking marine particles. Nat. Commun. 12, 4085. https://doi.org/10.1038/s41467-021-23875-6 (2021).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Hygum, B. H., Petersen, J. W. & Søndergaard, M. Dissolved organic carbon released by zooplankton grazing activity-a high-quality substrate pool for bacteria. J. Plankton Res. 19, 97–111. https://doi.org/10.1093/plankt/19.1.97 (1997).CAS 
Article 

Google Scholar 
Suttle, C. A. Marine viruses–major players in the global ecosystem. Nat. Rev. Microbiol. 5, 801–812. https://doi.org/10.1038/nrmicro1750 (2007).CAS 
Article 
PubMed 

Google Scholar 
Bizic-Ionescu, M., Ionescu, D. & Grossart, H. P. Organic particles: Heterogeneous hubs for microbial interactions in aquatic ecosystems. Front. Microbiol. 9, 2569. https://doi.org/10.3389/fmicb.2018.02569 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Arandia-Gorostidi, N., Weber, P. K., Alonso-Saez, L., Moran, X. A. & Mayali, X. Elevated temperature increases carbon and nitrogen fluxes between phytoplankton and heterotrophic bacteria through physical attachment. ISME J. 11, 641–650. https://doi.org/10.1038/ismej.2016.156 (2017).CAS 
Article 
PubMed 

Google Scholar 
Worrich, A. et al. Mycelium-mediated transfer of water and nutrients stimulates bacterial activity in dry and oligotrophic environments. Nat. Commun. 8(1), 15472. https://doi.org/10.1038/ncomms15472 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Iversen, M. H. & Ploug, H. Ballast minerals and the sinking carbon flux in the ocean: Carbon-specific respiration rates and sinking velocity of marine snow aggregates. Biogeosciences 7, 2613–2624. https://doi.org/10.5194/bg-7-2613-2010 (2010).ADS 
CAS 
Article 

Google Scholar 
Baltar, F., Arístegui, J., Gasol, J. M., Sintes, E. & Herndl, G. J. Evidence of prokaryotic metabolism on suspended particulate organic matter in the dark waters of the subtropical North Atlantic. Limnol. Oceanogr. 54, 182–193. https://doi.org/10.4319/lo.2009.54.1.0182 (2009).ADS 
CAS 
Article 

Google Scholar 
Schneider, B., Schlitzer, R., Fischer, G. & Nöthig, E.-M. Depth-dependent elemental compositions of particulate organic matter (POM) in the ocean. Glob. Biogeochem. Cycles https://doi.org/10.1029/2002gb001871 (2003).Article 

Google Scholar 
Jannasch, H. W. & Wirsen, C. O. Microbial activities in undecompressed and decompressed deep-seawater samples. Appl. Environ. Microbiol. 43, 1116–1124. https://doi.org/10.1128/AEM.43.5.1116-1124.1982 (1982).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Tamburini, C., Garcin, J., Ragot, M. & Bianchi, A. Biopolymer hydrolysis and bacterial production under ambient hydrostatic pressure through a 2000m water column in the NW Mediterranean. Deep Sea Res. II(49), 2109–2123. https://doi.org/10.1016/s0967-0645(02)00030-9 (2002).ADS 
Article 

Google Scholar 
Iversen, M. H. & Ploug, H. Temperature effects on carbon-specific respiration rate and sinking velocity of diatom aggregates: Potential implications for deep ocean export processes. Biogeosciences 10, 4073–4085. https://doi.org/10.5194/bg-10-4073-2013 (2013).ADS 
Article 

Google Scholar 
Guillard, R. R. & Ryther, J. H. Studies of marine planktonic diatoms I Cyclotella nana Hustedt, and Detonula confervacea (cleve) Gran. Can. J. Microbiol. 8, 229–239. https://doi.org/10.1139/m62-029 (1962).CAS 
Article 
PubMed 

Google Scholar 
Pernthaler, A., Pernthaler, J. & Amann, R. Fluorescence in situ hybridization and catalyzed reporter deposition for the identification of marine bacteria. Appl. Environ. Microbiol. 68, 3094–3101 (2002).ADS 
CAS 
Article 

Google Scholar 
Amann, R. I., Krumholz, L. & Stahl, D. A. Fluorescent-oligonucleotide probing of whole cells for determinative, phylogenetic, and environmental studies in microbiology. J. Bacteriol. 172, 762–770 (1990).CAS 
Article 

Google Scholar 
Daims, H., Brühl, A., Amann, R., Schleifer, K. & Wagner, M. The domain-specific probe EUB338 is insufficient for the detection of all bacteria: Development and evaluation of a more comprehensive probe set. Syst. Appl. Microbiol. 22, 11 (1999).Article 

Google Scholar 
Eilers, H., Pernthaler, J., Glockner, F. O. & Amann, R. Culturability and in situ abundance of pelagic bacteria from the North Sea. Appl. Environ. Microbiol. 66, 3044–3051 (2000).ADS 
CAS 
Article 

Google Scholar 
Manz, W., Amann, R., Vancanneyt, M., Schleifer, K.-H. & Ludwig, W. Application of a suite of 16S rRNA-specific oligonucleotide probes designed to investigate bacteria of the phylum cytophaga-flavobacter-bacteroides in the natural environment. Microbiology 142, 1097–1106. https://doi.org/10.1099/13500872-142-5-1097 (1996).CAS 
Article 
PubMed 

Google Scholar 
Amann, R. I., Ludwig, W. & Schleifer, K. H. Phylogenetic identification and in situ detection of individual microbial cells without cultivation. Microbiol. Rev. 59, 143–169. https://doi.org/10.1128/mr.59.1.143-169.1995 (1995).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Amann, R. I. et al. Combination of 16S rRNA-targeted oligonucleotide probes with flow cytometry for analyzing mixed microbial populations. Appl. Environ. Microbiol. 56, 1919–1925. https://doi.org/10.1128/AEM.56.6.1919-1925.1990 (1990).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Musat, N. et al. A single-cell view on the ecophysiology of anaerobic phototrophic bacteria. Proc. Natl. Acad. Sci. USA 105, 17861–17866. https://doi.org/10.1073/pnas.0809329105 (2008).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Polerecky, L. et al. Look@NanoSIMS: A tool for the analysis of nanoSIMS data in environmental microbiology. Environ. Microbiol. 14, 1009–1023. https://doi.org/10.1111/j.1462-2920.2011.02681.x (2012).CAS 
Article 
PubMed 

Google Scholar 
Musat, N. et al. The effect of FISH and CARD-FISH on the isotopic composition of (13)C- and (15)N-labeled Pseudomonas putida cells measured by nanoSIMS. Syst. Appl. Microbiol. 37, 267–276. https://doi.org/10.1016/j.syapm.2014.02.002 (2014).CAS 
Article 
PubMed 

Google Scholar 
Meyer, N. R., Fortney, J. L. & Dekas, A. E. NanoSIMS sample preparation decreases isotope enrichment: Magnitude, variability and implications for single-cell rates of microbial activity. Environ. Microbiol. https://doi.org/10.1111/1462-2920.15264 (2020).Article 
PubMed 

Google Scholar  More

Study subjectsSubjects in our study are individually recognized wild capuchins found in and around the Lomas Barbudal Biological Reserve in Guanacaste, Costa Rica. This population has been under observation since 1990 (Perry 2012; Perry et al. 2012), including near continuous observation from January 2002 through March 2020.Data collectionWe use proximity data on subjects collected during group scan sampling between January 2001 and March 2020 (Altmann 1974). Included in scans are the identity of the subject, and the identity of other individuals within approximately 4 meters of them. Scans have been collected on all individuals in study groups since 2002, and on all adults and subadults since 2001. Scans are taken opportunistically, without regard to time of day. At least 10 min separate consecutive scans of the same individual to reduce the non-independence of scans taken close in time.Data in this manuscript were collected by 124 observers, with an average of 7.1 data collectors per month. Observers typically work in teams of two to three and rotate across different groups to reduce potential observer bias. Observers also rotate across observer teams to avoid observer drift in coding, since observer teams could potentially start to code behaviors differently from each other in the absence of overlap in observer composition.Initial pedigree constructionOf the 376 individuals in our behavioral dataset, 280 (74.5%) were first seen within three months of their births, and we could confidently assign maternity to them based on demographic (pregnancies) and behavioral data (primary nursing) even prior to genotyping. Of the remaining individuals, 41 (10.9%) were males of unknown origin that immigrated into our study population, while the rest were natal to our study groups but were first seen as older infants ( >3 months), juveniles, or (sub)adults (14.6%) and required genotyping to assign/confirm maternity. Paternity was assigned based on genetic information when possible (but see Non-genotyped individuals).In total, 287 subjects (76.3%) had two assigned parents, 37 had one assigned parent (9.8%), and 52 (13.8%) had no assigned parent based on demographic, behavioral, and/or genetic parentage information. Most individuals with no assigned parents were immigrant males (78.9%).GenotypingInformation on genetic parentage assignment (at up to 18 microsatellite loci) in our study population is available from previously published work (1996–2005 (Muniz et al. 2006), 2005–2012 (Godoy et al. 2016b)). Partial genotypes (up to 14 loci) have been generated for individuals in this study which more recently entered the study population through birth or immigration (n = 91, 2012–2020) (See SI File 1). Briefly, DNA was extracted primarily from non-invasively collected fecal samples, and occasionally from tissue samples obtained from deceased individuals, then amplified at up to 18 autosomal tetranucleotide microsatellite loci (Muniz and Vigilant 2008) using either a 1-step or 2-step PCR protocol (Arandjelovic et al. 2009). There were no significant deviations from Hardy-Weinberg equilibrium, and no evidence of linkage disequilibrium between loci was found (Muniz 2008).DNA samples were run at a minimum in triplicate, but additional PCRs were performed on low quality samples (e.g., with low quantities of DNA). Genotypes at each of the loci were assigned to be heterozygous when each allele was seen at least twice in independent PCRs, and assigned as homozygous when the allele was seen in at least three independent PCRs in absence of a second allele.Amplicons were analyzed using an ABI PRISM3100 automated sequencer and GeneMapper Software (Applied Biosystems, Foster City, CA, USA). Likelihood-based parentage assignments were performed using CERVUS 2.0 or 3.0 (Marshall et al. 1998; Kalinowski et al. 2007). The average exclusionary power of the 18 microsatellites was 0.9888 for the first parent and 0.9998 for second parent (Muniz et al. 2006).Individuals with unknown parents (e.g., immigrant males, founders) were genotyped twice (i.e., using two independent DNA samples) following the procedures described above to guard against sample mix up. Known mother-offspring pairs were confirmed by ascertaining the absence of Mendelian mismatches across all loci for the pair, though one mismatch was allowed to account for null alleles, mutations, and genotyping errors. We detected one null allele in the population in 19 individuals and traced it back to a male who was either the father or grandfather of those individuals (Muniz et al. 2006; Godoy et al. 2016b).Candidate males for paternity assignment were chosen based on group membership around the time of an infant’s conception (typically 1–10 males). In cases when conceptions occurred prior to the habituation of a study group, we used the identities of all adult males present when the group was first observed. Candidate mothers were similarly chosen for individuals that were first seen as older infants, juveniles, or (sub)adults. For individuals born post-group habituation, CERVUS has always assigned paternity from the pool of potential candidate fathers. Parent-offspring pairs and trios were allowed one mismatch (excluding those at the locus with the known null allele).Pedigree updatingNon-genotyped individualsDuring stable tenures, alpha males in our population sire approximately 73% of infants born in their groups, including 90% of offspring born to unrelated females (Godoy et al. 2016a). There is strong evidence of inbreeding avoidance between alpha males and their female descendants, with relatedness to females as the primary factor constraining alpha male monopolization of paternity within groups (Muniz et al. 2006, 2010; Godoy et al. 2016a, 2016b; Wikberg et al. 2017, 2018). We used this information to update our pedigree, filling missing father information with the identity of the alpha male around the time of a non-genotyped individual’s conception, but only if their mother was not the daughter or granddaughter of the alpha male (i.e., with inbreeding avoidance). This approach allowed us to assign presumed paternity to 21 non-genotyped individuals (5.6% of subjects) who were natal to our study groups.Individuals with missing or incomplete parentageOut of the original four study groups (from which fissions led to eight additional study groups), we lacked parentage information (i.e., neither parent was sampled) for 12 individuals first seen at the time of habituation. We had incomplete parentage on an additional 11 adults (i.e., only one parent was sampled). We used the software program COLONY version 2.0.6.7 to look for evidence of whether these individuals were related to each other at the level of full sibling (Jones and Wang 2010). We also looked for potential full sibling pairs among the non-natal immigrant males in the population, since co-migrant males are typically kin (Perry 2012; Wikberg et al. 2014, 2018). We assigned five full sibling pairs among co-migrant males, and four full sibling pairs among natal founders. For any remaining co-migrant males and natal founder pairs that were not assigned as full siblings, we assumed these to be either paternal (migrants) or maternal (natal) half siblings, as is typical in this study population (Perry 2012). These assignments are likely to have some error. However, based on what we know about kinship in capuchins, it would introduce more error to assume that these pairs are unrelated.We pruned our modified population pedigree using the R package pedantics version 1.01 (Morrissey and Wilson 2010), to include only individuals that were linked to the subjects in our behavioral dataset. The reduction in missing data can improve convergence and mixing of models (Hadfield 2010). The pruned pedigree contained 419 individuals, with 353 maternities, 354 paternities, 209 full sibships, 413 maternal half sibships, and 1496 paternal half sibships. Maximum pedigree depth was six generations (mean = 3.03).Sociality measures (response variables)We generated two related proximity-based measures of sociality—(1) whether an individual was seen in proximity of another monkey (within ~4 meters) during a scan (i.e., they were not alone), and (2) the number of partners an individual has nearby (within ~4 meters) during a scan. The former is measure of the propensity of an individual to be social versus alone, while the latter is more indicative of the gregariousness of an individual. These two phenotypes are not independent, as they are generated from the same data (Fig. 1a).Fig. 1: Distribution of sociality, sampling, group size, and alpha tenure length.The scatterplot in a shows the proportion of scans per individual per month where the subjects were recorded in proximity of others on the x-axis, and the average number of social partners per scan per month for subjects on the y-axis. The sizes of the circles in a are proportional to sample size (range: 5–317 scans per data point). The figure in b shows the number of calendar years of data sampling per subject (range: 1–20), c variation in group size, and d the number of calendar years represented by different alpha tenures in the dataset. Note that d does not represent the full diversity of alpha tenure lengths in the population, only within the dataset: some tenure lengths are left-truncated as data from 1990–2000 are not included in this dataset. Figure produced in R using ggplot2 version 3.3.5 (Wickham 2016) and cowplot version 1.1.1 (Wilke 2020). The capuchin image was generated in R using sketcher version 0.1.3 (Tsuda 2020) based on an image taken by Nicholas Schleissmann.Full size imageWe compiled the scans of individuals by month (mean: 31.9, range: 5–317 scans per month) so that we had counts of (1) the total number of scans where an individual was social and (2) the total number of partners an individual had. With these counts we could look at the (1) proportion of time spent social (versus alone) and (2) the average number of partners an individual had, while still preserving information about sampling density (number of scans).To be included in any month, subjects needed to have at least five scan samples in that period. As we are interested in the repeatability of our measures of social behavior, subjects had to have at least six months of data to be included.We excluded dependent infants (less than one year of age) as potential social partners of their mothers. We also excluded these dependent infants as subjects, since an infant is expected to be in close proximity of its mother, particularly during the first half of their first year of life (Godoy 2010; Perry 2012). Including data from infants would likely introduce upward bias to heritability estimates, because mothers and their dependent offspring (whom share high relatedness) would often be in close proximity of each other, and their measures of proximity to others would thus also be highly correlated.On average, subjects spent just over half of their sampled time within approximately four meters of another monkey (mean: 0.539, standard deviation: 0.193) and had approximately one social partner per scan (mean: 1.057, standard deviation: 0.619) (Fig. 1a). Our dataset consisted of 22,138 monthly sociality scores on 376 subjects generated from 641 140 scans (mean: 56.5 months per subject, range: 6–184 months per subject). Almost all subjects (99.7%, i.e., all but one) were represented by data across more than one calendar year (25, 50, 75% quantiles: 4, 7, 10 different years of data collection) (Fig. 1b).Fixed effectsWe included age (as a cubic function) and sex in our models, as well as their interaction to account for differences in how male and female capuchins sexually mature and age. Age in our dataset was right-skewed with higher representation at younger ages (mean: 9.3 years, standard deviation: 6.9) (Fig. 2). To put the ages in developmental context, mean age at first live birth is around 6.3 years for females in this population, though females can begin reproducing in their 5th year (Perry et al. 2012). Males as young as six years old have been known to sire offspring (Godoy et al. 2016b), but males tend to not reach full adult size until their 10th year (Jack et al. 2014).Fig. 2: Sociality as a function of age and sex.Circles represent individual monthly data. The sizes of the circles are proportional to sample size (range: 5–317 shows per data point). Circles in a represent the proportion of time individuals were seen in proximity of others (not alone) per month, while in b represent the average number of partners for individuals per month. Solid lines represent estimated sociality scores based on age and sex, with all other fixed effects set to the mean. The two x-axes represent age as z-scores and in years. Figure produced in R using ggplot2 version 3.3.5 (Wickham 2016).Full size imageSeasonal environmental changes, such as in food abundance, or temperature and rain, can lead to changes in how individuals cluster near others, for example, because of how food resources become distributed in the environment. For example, in black-crested gibbons (Nomascus concolor), group averages of dyadic proximity have been documented to decrease from the dry season to wet season, with increased average group proximity during cold months and lowered proximity during warm months (Guan et al. 2013). We account for seasonal variation by modelling monthly changes as a sine wave, through inclusion of the sine and cosine functions of a transformed month variable (See SI File 1 for further details).Central America is a region of ENSO-related precipitation, where the El Niño-Southern Oscillation (ENSO) has an impact on large scale patterns of temperature and precipitation (Ropelewski and Halpert 1987). Bimonthly rainfall anomalies are linked with both the warm El Niño and cool La Niña phases in a neighboring tropical dry forest in Costa Rica, where long-term monitoring of wild white-faced capuchins has shown declines in reproductive output associated with El Niño-like conditions (Campos et al. 2015). To account for the large-scale influence of ENSO on group dynamics, we included a factor variable for three different ENSO phases (Average/Neutral, Cool/La Niña, and Warm/El Niño). We used the bi-monthly Multivariate El Niño/Southern Oscillation (ENSO) index (MEI.v2) obtained from the Physical Sciences Laboratory of the National Oceanic and Atmospheric Administration (https://psl.noaa.gov/enso/mei/, retrieved: 2021-11-06) to determine the different phases. MEI.v2 is a composite index of five different variables (sea level pressure, sea surface temperature, surface zonal winds, surface meridional winds, and Outgoing Longwave Radiation) used to create a time series of ENSO conditions from 1979 to present (Zhang et al. 2019). Warm phases correspond to MEI.v2 values of 0.5 or higher, while cool phases correspond to values of −0.5 or lower.Demographic differences between groups and within groups across time can also lead to variation in behavior. For example, group size has been found to correlate with the amount of time that individuals spend grooming in various primate species (Dunbar 1991; Lehmann et al. 2007). Group size is also associated with higher sociality measures such as both the number of strong and weak ties that individuals form in diverse clades of primates (Schülke et al. 2022). We attempt to account for variation that arises from such demographic differences by including group size (mean: 24.7, standard deviation: 7.9) (Fig. 1c) as a fixed effect.In our models, group size and cubic age were centered and scaled to a mean of zero and a standard deviation of one.Random effectsAll models include the identity of the subject (VID, n = 376) as a random factor, as well as subject identity nested within year (VID:Year, n = 3150), the identity of each subject’s mother (VM, n = 142), maternal identity nested within group of residence within year of data collection (VM:GroupAlpha:Year, n = 2085), and a special variable known as the animal term to account for additive genetic variance (VA). These components contribute to long- and/or short-term repeatability of individuals. All models also include year of data collection (VYear, n = 20), month nested within year (VMonth:Year, n = 224), and the identity of each subject’s group of residence both across years (VGroupAlpha, n = 56) and within years (VGroupAlpha:Year, n = 200).VID in the models (since the models also additionally estimate VM and VA) can be thought of as estimating the “permanent environment variance” (i.e., VPE) of an individual, which is the “individual-specific variation in environmental conditions that permanently affect the phenotype (e.g. early-life conditions)” (Dingemanse et al. 2010). VID:Year captures the variance explained by the repeated sampling of the same individuals within a particular year. We use it to estimate the proportion of the phenotypic variance due to similarity in the trait within individuals from data taken closer in time (within the same year). During such a relatively short period, individuals are more likely to be stable in important social traits such as kin availability, dominance rank for adults, and maternal dominance rank for infants and young juveniles.VM estimates the variance explained by maternal effects (m2), specifically similarity between maternal siblings. Maternal identities were not available for all subjects, namely 11 immigrant males of unknown origin who were not assigned by COLONY as having a full sibling. We created unique dummy codes for their maternal identities, so that no two of these individuals shared the same mother. We additionally nested maternal identities (VM:GroupAlpha:Year) to account for similarity between maternal siblings residing in the same group in the same year. Such a nested structure might capture potential upward biases on heritability due to maternal kin biases in spatial association among siblings residing in the same group.We estimate h2 in our models by fitting a random effects term (VA), referred to as the animal term, which in the R package MCMCglmm links to the identities of individuals in our population pedigree (Hadfield 2010; see below for details on the implementation of the models in MCMCglmm). Inclusion of the animal term provides our models with an additive genetic variance component based on the estimated coefficients of relatedness between individuals in our pedigree. In short, if animals that share more alleles are also more like each other in their behavior, then variation in the behavior may well be due to genetic variation in the population (under the assumption that phenotypic similarity is not due to a shared environment, or is adequately controlled for by fixed and random effects in the model).VYear and VMonth:Year were included in order to account for temporal variation in sociality scores not captured by the fixed effects of seasonality or ENSO phase. These could arise from, for example, observer drift in coding (i.e., measurement error) or prevailing environmental conditions (e.g., drought) that could lead to changes to how individuals cluster near others. There were 218 unique observer combinations across the 224 months represented in the dataset, so VMonth:Year should also capture variance due to any differences between observer teams, though we cannot separate out the unique influence of observers.VGroupAlpha represents variance arising from the different alpha tenures within groups in our study population. VGroupAlpha captures both variance due to group of residence effects and the additional influence of alpha tenures within those groups. In capuchins, alpha males are ‘keystone’ individuals, whose influence is disproportionate relative to that of others in the population, and thus play important roles in establishing group dynamics (Jack and Fedigan 2018). Including group of residence, as defined by alpha tenure, is also important because it helps to account for the higher relatedness within groups within alpha tenures which results from high male reproductive skew toward alpha males. At Lomas Barbudal, males can remain in their alpha position for upwards of 18 years. Alpha tenures in this dataset spanned one to 14 years (Fig. 1d), so we additionally nested the identity of alpha males per group within years (VGroupAlpha:Year) so as to separate the within-year and across-year influences of group of residence.Statistical methodsWe ran analyses in R 4.1.2 (R Core Team 2021), using a Bayesian method with the R package MCMCglmm version 2.32 (Hadfield 2010). Data and code used to run all models is provided in the Supplementary Information.For our binary response variable (social versus alone), which was pooled into monthly units, we fit models with a binomial distribution and logit link function (family = “multinomial2”), with the number of scans each individual was documented social (‘successes’) versus the number of times alone (‘failures’).For our other response variable (number of partners), which was also pooled into monthly units, we fit models with a Poisson distribution (family = “poisson”), with the total number (sum) of partners per month. We included the natural log of the number of scans per month as a fixed effect to account for sampling effort. We set a strong prior for the log of sampling effort so that the rate at which events occurred was 1 (i.e., we could look at average number of partners per scan).We used a parameter-expanded prior (V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000) and two inverse Wishart priors (V = 1, nu = 0.002; V = 1, nu = 0.02) for the G structures in our models (i.e., random effects variance components). The prior on the residual variance component was set to one for both the binomial and Poisson models. Estimates for variance components were robust against the choice of prior (SI Fig. 3). We therefore only report findings from models run with parameter-expanded priors in the main text.Pilot runs (thin = 10, burnin = 3000, nitt = 13,000) indicated that autocorrelation values would remain high for some variance components in models run with parameter-expanded priors, even with large thinning intervals. We therefore increased the number of iterations to guarantee effective sample sizes of at least 1000, but ideally closer to 4000. All models were run with a long burn-in period of at least 10,000 iterations.We ran multiple chains (n = 4) of each model and assessed convergence of the chains visually (SI Files 2a-b), as well as through the Gelman-Rubin criterion implemented via the ‘gelman.diag’ function from the coda package in R (version 0.19-4) (Plummer et al. 2006). Scale reduction factors were below 1.02, signifying good convergence. We used Heidelberger and Welch’s convergence diagnostic test for stationarity to check convergence of each chain using the ‘heidel.diag’ function from the coda package. Results are presented from the first chain of each model.Reduced modelsInclusion of fixed effects can potentially have an impact on the estimates of variance components in models because total phenotypic variance (VP) is estimated (and partitioned among the different random effects) after conditioning on the fixed effects. Heritability estimates, for example, can be higher because the variance explained by the fixed effects structure (VFE) is not included in VP, thus making the relative contribution of VA to VP larger compared to the same model without fixed effects (Wilson 2008). Conversely, not adequately controlling for relevant fixed effects that contribute to phenotypic variance among and within individuals may potentially lead to an underestimation of VA and associated heritability (h2).We ran multiple reduced versions of our models to look at the impact of fixed effects on our variance components. We began with an intercept-only version (i.e., no fixed effects), then built-up complexity by adding in versions with the properties of the individuals first (age, sex), then properties of the group (group size), and subsequently environmental properties (seasonality, ENSO phases). Outputs for these reduced models are provided in the Supplementary Information (SI Table 2, SI Table 3).We provide the deviance information criterion (DIC) values for models (automatically generated by the MCMCglmm package). DIC is a generalization for multi-level models of the Akaike Information Criterion (AIC); and as in AIC, lower DIC values indicate better fit.Transformations from unobserved latent scale to observed data scaleOutputs from our MCMCglmm models were on the unobserved latent scale. We used the R package QGglmm (version 0.7.4) to additionally compute parameters of interest on the observed data scale (de Villemereuil et al. 2016; de Villemereuil 2018). We used the functions ‘QCicc’ to compute Intra-Class Correlation (ICC) coefficients and ‘QGparams’ to compute additive genetic variance and thus narrow-sense heritability (h2) on the observed data scale. We implemented the ‘QGparams’ and ‘QGicc’ functions with parameters model = ‘binomN.logit’ and n.obs = 32 (the average number of scans per subject per month in our dataset) for the binomial model and model = ‘Poisson.log’ for the Poisson model. The choice of value for n.obs is somewhat arbitrary, and we show the consequences for changes in values of this parameter (i.e., higher estimates with increasing values of n.obs) in SI Fig. 4.Closed form solutions in QGglmm are not available for integrating over posterior distributions generated from binomial models with logit link functions (de Villemereuil 2016). Consequently, using the ‘QGicc’ function is particularly slow. We therefore estimate ICCs from our binomial models using a random subset of the posterior (n = 1000 iterations).The code used for transforming the MCMCglmm outputs from the latent scale to the original data scale are available online (see DATA AVAILABILITY).Repeatability and the proportion of variance explained by variance componentsTotal phenotypic variance (VP) was the sum of estimates from all variance components and residual variance in a model (VP = VID + VID:Year + VM + VM:GroupAlpha:Year + VA + VGroupAlpha + VGroupAlpha:Year + VMonth:Year + VYear + Vresidual). The proportion of variance explained by each variance component was calculated by including its estimate in the numerator while including total phenotypic variance in the denominator. So, for example the proportion of variance explained by year of data collection was calculated as (left( {frac{{V_{Year}}}{{V_P}}} right)).Long-term repeatability was calculated with the sum of VID, VM, and VA in the numerator. Short-term repeatability was calculated similarly but with inclusion of within-series variances (VID + VM + VA + VID:Year + VM:GroupAlpha:Year) in the numerator to capture additional consistency in among-individual differences resulting from greater environmental similarity within a time series (i.e., year).We report posterior modes and 95% Highest Posterior Density intervals (i.e., 95HPDI in square brackets). Unless mentioned otherwise, we present results on the unobserved latent scale, and without the variance from the fixed effects (VFE) incorporated into VP. For completeness, estimates with VFE included in VP and transformations to the observed data scale are also provided in SI Table 3. More

The top-brand nonglutinous rice variety ‘Koshihikari’, which has a high palatability, is extremely susceptible to blast. Therefore, farmers apply fungicides over four times during the rice production season. As Koshihikari is sold by the Niigata brand, it has been traditionally viewed as having a high eating quality in Japan, and because of this, both farmers and consumers have requested that the multiline variety KO-BL be tested to determine if it is equivalent to Koshihikari before its introduction. Trials comparing Koshihikari and KO-BL were carried out in 2003 and 2004 in 594 and 622 fields covering 236 and 315 ha, respectively. These trials evaluated plant homogeneity, eating quality, and blast suppression using fewer fungicidal sprays. Following favorable results, in 2005, all Koshihikari were converted to KO-BL multiline variety covering an area of 94,000 ha. In addition, seed use and cultivation were restricted to the Niigata area to distinguish KO-BL from Koshihikari grown in other prefectures.Seed production and mixture processes are managed with precision by each prefectural official member (Fig. 1a). Original isogenic lines (ILs) were separately produced from the original stock in the original strain fields by the Niigata prefectural government. Using a precise mixture machine, the mixture of four ILs was then blended by weight in 2000 kg volumes, all multiplied by ten (giving a total volume of 20 t). Original production fields and commercial fields all used blended seeds that had been authorized by seed production farmers and commercial farmers in the 2003 and 2004 trials. Thus, it takes two years for seed production at the original strain field followed by the original production field for the preparation of commercial fields; thus, the seed mixture composition needs to be determined at least two years before introduction. Susceptible and resistant (effective) ILs were mixed at a ratio of 3:7 from 2005 to 2019 (Fig. 1b, Supplementary Table S1). Susceptible ILs, possessing Pia and Pii genes, were always mixed at a ratio of 1:2, but the composition of resistant ILs, containing Pita-2, Piz, Pib, Piz-t, and Pit genes, was changed every two to three years to avoid the breakdown of resistance6. These changes were determined by annually monitoring blast race distributions.Figure 1Representative seed production flow from original stock to commercial field and history of Koshihikari BL composition from 2005 to 2019 in Niigata Prefecture. (a) S1–S2, susceptible KO-BL; R1-R2, resistant KO-BL. Seeds obtained from original stock field at the Niigata Agricultural Research Institute. Seeds obtained from the original strain field and the original production field at both designated farmers’ fields. Commercial field (general farmers field) used for KO-BL production. Each field requires a year for seed production. (b) Pia and Pii, susceptible; Pita-2, Piz, Pib, Piz-t, and Pit, resistant. The proportion of susceptible KO-BLs and resistant KO-BLs was consistently 3:7 across years.Full size imageIn Niigata Prefecture, the predominant 5 blast races distributed from 1994 to 2004 were 001.0 (virulent to Koshihikari [Pik-s]), 003.0 (virulent to Pik-s and Pia), 005.0 (virulent to Pik-s and Pii), 007.0 (virulent to Pik-s, Pia, and Pii), and 037.1 (virulent to Pik-s, Pia, Pii, and Pik) (Fig. 2a, Supplementary Table S2). Because all the 5 races were virulent to Koshihikari, which had been widely cultivated in Niigata area during the years, there were no drastic race changes. In addition, genetic variations in blast resistance indicated that Koshihikari also harbored the Pish gene, and that the Pia, Pii, and Pik genes were also dominant in the Hokuriku region, including Niigata Prefecture21. Virulent blast races against the resistance genes Pish, Pia, Pii, Pi3, Pi5(t), Pik, Pik-s, and Pi19(t) were dominantly distributed in Niigata Prefecture22. These reports confirmed that Koshihikari had been susceptible to dominant blast races before KO-BL introduction.Figure 2Blast race change during the 1994–2019 period in Niigata Prefecture and the worst-case simulation of blast race dynamics in KO-BL during the 2005–2019 (years 1–15) period. Races and virulences are shown in Table 1. (a) A red line indicates the year (2005) when KO-BL was introduced. Races 007.0 and 037.1 became dominant after the introduction. (b) Actual races and their rates in 2004 and annual KO-BL compositions from 2005 to 2019 were set in the simulation. Parameters set in the simulation were as follows: maximum lesion number in a year, 10,000,000; weather condition, 10 (favorable); virulent mutation rate, 10–5; overwintering probability, 0.01; number of simulated years, 15; and number of simulation trials, 1000. The 1000 trial results for the lesion number increase in each race were averaged in each year and transformed into rates to show race dynamics. All simulation results are shown in Supplementary Table 6 in Supplementary information 2. The races 007.0 and 037.1 were also dominant until year 15 (correspond to 2019). Both actual and simulated race dynamics showed no outbreaks of the resistant composition of KO-BL.Full size imageIn the 2005 release year of KO-BL, the predominant blast races, 001.0 (virulent to Pik-s) and 003.0 (virulent to Pik-s and Pia), drastically decreased in distribution from 41.8% to 22.3% and 27.6% to 17.3%, respectively (Fig. 2a, Supplementary Table S2). Interestingly, races 001.0 and 003.0 rapidly decreased by 5.4% and 1.3% in 2006, respectively, even though especially Pia, which can be infected by the race 003.0, was used in the KO-BL composition. Because all ILs in the composition of KO-BL were resistant to race 001.0, and race 003.0 was only virulent to Pia, which made up 10% of the annual KO-BL composition (Table 1). In contrast, races 007.0 (virulent to Pik-s, Pia, and Pii) and 037.1 (virulent to Pik-s, Pia, Pii, and Pik) dominated from 2005 to 2019. The higher rate of race 007.0 detection was affected by 30% of the ILs composing the annual KO-BL were susceptible. The second highest rate of race 037.1 detection was affected by a number of factors: the high susceptibility of a minor cultivar that had Pii and Pik, the mosaic configuration of fields typical in Niigata, and the air-borne spread of race 037.1. To maintain consensus on KO-BL cultivation based on total blast suppression in Niigata, rarely detected races virulent to resistant ILs in commercial fields are strictly supervised by the prefectural government to avoid unnecessary confusion in Niigata residents.Table 1 Susceptible or resistant reaction of Koshihikari and KO-BL against blast races.Full size tableIn 2008, to mathematically support KO-BL composition changes, we developed a simulation software to estimate long-term blast race dynamics in multilines using a plant‒pathogen coevolution system23. The model calculated the persistence of resistant ILs to determine the optimal timing of changes to multiline variety compositions. To simulate race dynamics in KO-BL, we set five currently investigated races, 001.0 (virulent to Pik-s), 003.0 (virulent to Pik-s and Pia), 005.0 (virulent to Pik-s and Pii), 007.0 (virulent to Pik-s, Pia, and Pii), and 037.1 (virulent to Pik-s, Pia, Pii, and Pik), and their rates in 2004, as well as five emerging races, 043.0 (virulent to Pik-s, Pia, and Piz), 303.0 (virulent to Pik-s, Pia, and Pita-2), 003.2 (virulent to Pik-s, Pia, and Pib), 403.0 (virulent to Pik-s, Pia, and Piz-t), and 003.4 (virulent to Pik-s, Pia, and Pit) (see Fig. 2b, Supplementary Table S3) against five newly introduced respective resistant KO-BLs (see Fig. 1b, Supplementary Table S1) and the annual KO-BL compositions from 2005 to 2019. The worst case (severe epidemic) simulation result (Fig. 2b, Supplementary Tables S3 and S6) showed that race 007.0 (virulent to susceptible Pik-s, Pia and Pii) became the predominant race (77.4%), and race 037.1 (virulent to Pik-s, Pia, Pii, and Pik) remained at a low frequency (21.6%) until the fifteenth year (corresponding to 2019). In addition, super-race virulent to all KO-BLs did not emerge in this simulation. These suppression of outbreaks of newly emerged virulent races, including super-race on resistant KO-BL was apparently affected by 2–3 years of change in resistant KO-BL composition, and total suppression of blast occurrence decreasing the blast population. These results indicated that almost all the epidemics analyzed reflected actual race dynamics without affecting other minor races from other susceptible cultivars grown in Niigata, especially up to 2011. Thus, our decision support system provides an evaluation of KO-BL persistence and indicates the KO-BL composition changes needed for blast race population control in large areas. In addition, our simulation model may be useful for evaluating future KO-BL composition changes.Blast occurrence drastically decreased after 2005 (Fig. 3a, Supplementary Table S4). The average occurrence of leaf and panicle blast was 46.1% and 52.9% during the 1995–2004 period and 9.5% and 9.6% during the 2005–2019 period, respectively. This resulted in a blast suppression effect by 70% of the resistant composition in KO-BL. Current seed production fields are rarely contaminated with virulent races against resistant KO-BLs. This suggests that seed sanitation contributes to the suppression of virulent pathogen epidemics in multilines. In addition, induced resistance24,25 may have no effect on the practical use of multilines. Rice plants were found to induce a resistance response when inoculated with avirulent races of blast (those that stimulate protective responses to virulent race attacks). As the detection of several races in one area is rare and blast occurrence tends to be low, conditions that induce resistance in field situations do not occur. Fungicide applications to control blast in KO-BL and other minor cultivars decreased by approximately one-third during the 2005–2019 period compared with 2004 (Fig. 3b, Supplementary Table S5). Thus, the commercial scale use of crop diversity is clearly effective for the environmentally friendly control of airborne diseases.Figure 3Leaf and panicle blast occurrence from 1994 to 2019 and blast control area from 2004 to 2019 in Niigata Prefecture. (a) A red line indicates the year (2005) when KO-BL was introduced. (b) Gross fungicide spray area decreased by approximately one-third during the 2005–2019 period compared with 2004.Full size imageThe optimum long-term solution for pathogen population control using genetic diversity includes multilines. Blast occurrence in KO-BL introduced in Niigata, and the theoretical value of blast suppression in KO-BL tested at small scales, were reduced by approximately 10% compared to that of monoculture plots26,27,28. Thirty percent of susceptible ILs in KO-BL have the potential to improve compatible races with susceptible ILs and become predominant in large areas. This would contribute to the suppression of rapid increases in new virulent races emerging in the blast population. To maintain consensus on KO-BL cultivation based on total blast suppression in Niigata, rarely detected races virulent to resistant KO-BLs in commercial fields are strictly monitored by the prefectural government. Educating Niigata farmers ensures the long-term use of KO-BL. In fact, lower blast occurrence has been attributed to careful KO-BL cultivation and seed management.The implementation of genetically diversified homogeneous seed mixtures, rotations with resistant KO-BL, restricted KO-BL cultivation, and pathogen monitoring allowed rice quality to be maintained, diseases to be suppressed, and environmentally sound agriculture to be economically viable in Niigata. Collaboration among prefectural officers, farmers, and consumers in Niigata has resulted in safer rice production with good agricultural practices (GAPs) that meet sustainable development goals (SDGs). In addition, DNA tests differentiate KO-BL from the original Koshihikari for buyers, thereby prohibiting illegal distribution. Multiline varieties have been used in small areas in two different prefectures. For example, in Miyagi pref., Sasanishiki BL consisted of Pik, Pik-m, and Piz at ratios of 4:3:3 and 3:3:4 in 1995 and 1996, respectively. This composition was changed to Pik, Pik-m, Piz, and Piz-t at a ratio of 1:1:4:4 from 1997 to 2007 to prevent an increase in race 037.1 (virulent to the BL: Pik and Pik-m). In addition, an equal mixture of seven BLs (Pib, Pik, Pik-m, Piz, Piz-t, Pita, and Pita-2) was cultivated in 300 ha areas (maximum 4000 ha) from 2008 to 2014 without any outbreaks observed. In Toyama pref., the Koshihikari Toyama BL, which consists of resistant ILs, Pita-2, Pib, and Pik-p at a ratio of 4:4:2, was cultivated in an area of 300 ha and required a 50% reduction in chemical inputs from 2003 up to the present. Our model also calculated a greater than 50-year persistence in terms of the small area effect in both prefectural cases. This result depends on an insufficient pathogen population increase in virulent mutations against resistant ILs (data not shown). In this way, the practical use of a multiline provides control without the need for as much fungicide with or without a periodic change in IL composition. Our results demonstrate that the management of crop and pathogen coevolution can control diseases at large scales and, thereby, contribute to global food security. More

Data collectionStudy sitesWe selected four study sites in the Upper Gunnison Basin in Colorado, USA. As we intended for these sites to capture a large amount of variation in environmental conditions, the primary criteria we used for site selection was elevation. Mirroring patterns of mountain weather observed worldwide, elevation is correlated with a multitude of environmental factors in the Rocky Mountains. Predominant trends include a negative correlation between altitude and temperature, a positive correlation between altitude and precipitation34, and transitions in plant community composition (from ‘mountain shrub’ to ‘montane’ to ‘subalpine’ to ‘alpine’) with increasing altitude35.Our sites span a large elevation gradient of approximately 1000 m that captures significant variation in environmental conditions. These sites should not be interpreted as a direct ‘elevation for climate substitution’, especially considering the significant changes in weather at each site that occur over the course of our observation periods, which ran from June to July. The lowest site, ‘cement creek’ (CC; approximately 38.82156° N, 106.86893° W), falls within a sage brush meadow on the boundary between the mountain shrub and montane vegetation zones at 2440 m elevation. The second lowest site, ‘bus turnaround’ (BT; approximately 38.97130° N, 106.99595° W), is situated in an open subalpine meadow at the base of Gothic Mountain at ~ 2940 m. The second highest site, ‘gothic mountain’ (GM; approximately 38.97969° N, 107.01937° W), is also in the subalpine zone on a steep hillside within a clearing in an evergreen forest on the lower slopes of Gothic Mountain at 3220 m. The highest site, ‘high meadow’ (HM; approximately 38.96779° N, 107.02184° W), is on an exposed meadow at the upper fringes of the subalpine zone on the shoulder of Gothic Mountain at 3,410 m.We began our observations at each site by mapping the flax population (see “Population mapping” sub-section below). This occurred on 6/15/2020 for CC, 6/17/2020 for BT, 6/16/2020 for GM, and 6/18/2020 for HM. We concluded all observations at CC on 7/27/2020, at BT on 7/29/2020, at GM on 7/28/2020, and at HM on 7/10/2020 (aside from spore trap collection which concluded on 7/21/2020).Weather monitoringWe deployed environmental sensors to record longitudinal weather data at each of our sites. We recorded temperature, humidity, rainfall, and wind speed and direction (at ~ 1 m above ground) every five minutes (see Supplementary Text 1 for logger specifications). We used temperature and humidity data to calculate absolute humidity (Supplementary Text 2). Due to logistical constraints, rainfall and wind speed/direction loggers were deployed later than temperature and humidity loggers.We extracted weather metrics for a given observation period from this data as follows: We calculated the mean, maximum, and minimum temperature as the mean, maximum, or minimum of temperature readings. Likewise, we calculated mean absolute humidity as the mean of absolute humidity records (calculated from temperature and relative humidity every five minutes). Daily mean rainfall was calculated as of mean rainfall records multiplied by the number of 5 min increments in a day (i.e. 288).Future climate dataTo enable forecasting of the effects of climate change on within and between host transmission processes, we compiled projected future weather data. We extracted downscaled climate and hydrological data generated using the Localized Constructed Analogs (LOCA) method and the CESM1(CAM5) model from the “Downscaled CMIP3 and CMIP5 Climate and Hydrology Projections” archive36,37,38,39,40. Specifically, we extracted values for daily mean temperature, relative humidity, and rainfall from the CESM1(CAM5) hydrology model, and values for daily maximum and minimum temperature from the CESM1(CAM5) climate model for both the RCP4.5 and RCP8.5 emissions scenarios41. Temporally, this data spanned June through July in 2020, 2030, 2040, 2045, 2050, 2060, and 2070. Spatially, this data represents projections for the area spanning (38.9375, 39.0) latitude and (− 107.0, − 106.9375) longitude. This area is the finest grained spatial extent including the BT site for which projections were available.Population mappingIn order to uncover drivers of population level epidemiological dynamics, it is necessary to characterize the ‘landscape of susceptibility’ in the host population. For spatially structured infectious disease systems such as flax rust, this involves documenting not only the infection status of each individual and any likely covariates of contact rate or susceptibility (e.g. size), but also the spatial arrangement of individuals. To accomplish this, we established 10 by 20 m transects at each of our field sites. We recorded the coordinates of each transect corner and the compass bearings associated with the transect axes (so that wind direction could be translated to the transect coordinate system) using a Garmin GPS unit (part #: 010-01735-10). At the beginning of the field season, we mapped the location, height, and infection status (healthy or diseased) of each plant within each transect. Plants with heights less than five centimeters were denoted as seedlings. While other metrics of size such as stem count could be used to form a more complete picture of plant size, we chose to use height alone as a proxy for size as it was the only measure that we could feasibly record for the many hundreds of plants within each transect. At the CC site, we were unable to map all healthy plants in certain areas of the transect due to time constraints, but we did record data for all diseased plants in these regions.Epidemiological surveysAfter documenting the initial epidemiological conditions, we tracked the spatiotemporal spread of flax rust within each transect. Approximately once per week, we visually inspected each flax individual in the transect for signs of disease, recording the location of any newly infected plants so they could be matched to a previously uninfected plant. We also recorded the height of all newly infected plants. We conducted these epidemiological surveys at least seven times at each site. Plants identified as infected in the initial population mapping and in these subsequent epidemiological surveys were marked with flags to ensure that they would not be repeatedly recorded as newly infected. Although we did not map the entirety of the CC transect, we covered the entire 10 by 20 m area in all epidemiological surveys.Diseased focal plantsWithin each transect, we designated a subset of diseased plants as ‘focal diseased plants’ and used them to make detailed measurements of within-host disease spread processes. These processes included changes in infection intensity, which we measured at the plant scale, as well as pustule growth, which we measured at the sub-plant scale using changes in pustule area and as a proxy.All initially infected plants were designated as diseased focal plants, and we continued to give this designation to newly diseased plants (except in some cases, diseased plants with height of approximately 5 cm or less, as these plants were too delicate to be manipulated) until we obtained at least 25 infected focal diseased plants per site. This occurred on 7/20/2020 at CC, on 7/1/2020 at BT, and on 7/9/2020 at GM and HM. We inserted metal tags into the ground at the base of diseased focal plants so they could be quickly and reliably identified. To facilitate longitudinal measurements of disease spread at sub-plant scales, we marked up to three infected stems on each focal diseased plant with pieces of colored flagging tape tied around the base of the stem. These stems were chosen haphazardly. On each of these stems (when possible), we marked the tip of one infected leaf with black ink. We preferentially marked leaves with one or a few pustules. Approximately once every three days, we recorded detailed measurements of plant height and within-host disease spread processes.To measure pustule growth, we photographed each marked leaf with a millimeter ruler in frame using a Cannon EOS Rebel T7i DSLR camera fitted with a Canon EF-S 35 mm f/2.8 Macro IS STM lens and a Cannon MR-14EX II Macro Ring Lite. While achieving consistent photographic angles and lighting conditions is impossible in a field setting, we attempted to keep the leaf and the scale in the same plane while maintaining a perpendicular shooting angle. We also adjusted camera settings so that pustule boundaries could be clearly distinguished in images. We used ImageJ software42 to measure the maximum and minimum diameters of each photographed pustule, and calculated pustule area assuming pustules to be elliptical. To enable the size of individual pustules to be tracked across observations, each measured pustule in each image was labeled so that it could be re-identified in subsequent images based on its location on the leaf and position relative to other pustules. We omitted data for any pustule that could not be confidently identified or whose borders became indistinct from other pustules. We marked a new leaf when we were unable to find the previously marked leaf, when the leaf was accidentally removed during data collection, or when the condition of the leaf deteriorated significantly. We suspect that either the presence of ink or our handling of the plant contributed significantly to the rate at which leaf condition deteriorated, but nevertheless we were able to observe the growth of many pustules over many weeks.To measure infection intensity, we first counted the number of infected stems. Next, for each marked stem, we recorded the length of infected tissue as the length of the stem segment spanning the lowest (relative to the ground) pustule to the highest pustule. Finally, we counted the number of pustules present on a haphazardly selected leaf in the middle of the region of infected tissue on each marked stem. We did not use the same leaf at each observation as the middle of the region of infected tissue changed due to disease spread. Using these measurements, we calculated an infection intensity metric as the product of the number of infected stems, the average length of infected tissue, and the average number of pustules present on a leaf in the middle of the region of infected tissue.Healthy focal plantsAs we wished to investigate patterns of growth in both healthy and infected plants, we designated certain uninfected plants as healthy focal plants, and measured their height whenever we measured the height of focal diseased plants. To control for potential effects of stem marking on the growth of diseased focal plants, we made these same markings on healthy focal plants. Healthy focal plants were selected haphazardly, and we added new healthy focal plants over time to keep the number of diseased and healthy focal plants approximately equal. At some sites, the number of healthy focal plants fell towards the end of the observation period as many healthy focal plants became infected. When this occurred, we stopped recording data for the plant, and in some cases designated it as a diseased focal plant.Spore depositionTo measure the distribution of spore deposition from an infected plant, and how it relates to the infection intensity of the source plant and wind patterns, we deployed spore traps in arrays around a subset of focal diseased plants. We chose focal disease plants that were as removed as possible from other diseased plants to minimize the number of spores originating from other plants that would be caught in the spore traps. These spore traps consisted of a ~ 2 cm2 section of Scotch Permanent Clear Mounting Tape (part #: MT76272-5) affixed to plastic backing with double sided tape. These spore traps were secured into the ground with ~ 5 cm nails. We deployed these traps at distances of either {~ 5 cm} (traps at ~ 5 cm were placed as close to the plant as possible), {~ 5 and 25 cm}, {~ 5, 25, and 50 cm}, or {~ 5, 25, 50, and 100 cm} in each of the four directions corresponding to the axes of the transect in which they were located. The traps were left in the field for approximately one week, and then collected. The sections of mounting tape were transferred to microscope slides and sealed in place using clear packing tape. We then used a light microscope to count the number of spores present in a known area. We counted spores in an area of 0.8 cm2 for most spore traps that contained few or no spores, but we counted a minimum area of 0.1 cm2 for all spore traps, including those with greater than 10,000 spores per cm2. Melampsora lini spores were identified visually using a slide prepared with spores as a reference for spore size, shape, and color. The spores matched descriptions in the literature43. We did not observe any rust fungi other than M. lini in or around our transects, so we can be reasonably certain that our counts represent only M. lini spores from infected focal diseased plants.Statistical analysesPlant growth and within-host disease spreadTo infer how different weather factors affect plant growth, pustule growth, and changes in infection intensity, we fit generalized additive models18. We used longitudinal observations of plant height in healthy and diseased focal plants to fit the model of plant growth, and longitudinal observations of within-host disease spread at various scales in focal diseased plants to fit all other models. Because the observation periods associated with our measurements of these processes varied, we formatted the response variables in change-per-day units: For our analyses of plant growth, pustule growth, and infection intensity progression we used change in plant height per day, change in pustule area per day, and change in infection intensity per day respectively as the response variables.To infer the effects of weather variables, we included smooth terms for these factors as predictors in each model. For each observation of change in plant height, pustule area, or infection intensity, we determined the start and end timepoints of the observation period to extract the corresponding weather data. For observations of plant height, we used 12:00:00 on the day of the observation as the start/end timepoints of the observation period. For observations of pustule area, we used photograph timestamps as the start and end timepoints. For observations of infection intensity, we used the average of the timestamps of photographs taken of pustules on that plant’s leaves as the start and/or end timepoints if such photographs were taken, because all data for an individual diseased focal plant was generally collected at once. If such photographs were not taken, we used the mean timestamp of all photographs taken at that plant’s site on the observation date as the start and/or end timepoint. To focus our analysis on fine-grained relationships between weather and within-host disease spread, we discarded data corresponding to observation windows of eight or more days. According to this criterion, 13 of 609 data points were discarded in the analysis of plant growth, 216 of 3535 data points were discarded in the analysis of pustule growth, and 25 of 338 data points were discarded in the analysis of infection intensity progression. Using the start and end timepoints of the observation period as bounds, we extracted the following variables from the weather data: temperature (mean, maximum, and minimum), absolute humidity, and mean daily rainfall.To account for the nested nature of our data, we included smooth terms for study site and plant identity as random effects (accomplished by setting bs = ”re” in the s() function in mgcv) in each model. In the case of the analyses of pustule growth, plant identity refers to the plant on which the pustules were observed. We also included a term in each model that accounted for the value of the previous observation to capture any allometric effects. In the model of change in plant height per day, this term was a full tensor smooth product of last observed height and infection intensity. In the model of change in pustule area per day, this term was a smooth of last observed pustule area. In the model of change in infection intensity per day, this term was a full tensor smooth product of the base 10 logarithm of the last observed infection intensity and last observed maximum plant height.We fit these models using the gam() function in the mgcv R package and the restricted maximum likelihood parameter estimation method18,44. We used gaussian response distributions with identity link functions. The dimensionality of all basis functions was set to the default value determined by mgcv. The one exception to this involved the marginal basis for log10 infection intensity (a part of the tensor term) in the model of change in infection intensity per day. The dimensionality of this marginal basis was increased to 15 after model diagnostics indicated that the default basis dimension was insufficient. To implement variable selection so that insignificant predictor terms would be effectively removed from the models, we used the double penalty approach of Mara and Wood45 (via the ‘select = TRUE’ option in the mgcv gam() function).The coefficients of the fit models can be used to infer the effects of weather factors on plant growth and within host disease spread. To translate these inferences into predictions about how climate change might affect plant growth and within-host disease spread, we simulated trajectories of plant growth, pustule area, and infection intensity under various climate conditions. For each class of simulation, we started by defining a hypothetical starting state. For simulations of plant growth, this was a 10 cm tall uninfected plant (the plant was assumed to remain uninfected for the entire simulation). For simulations of pustule growth, this was a pustule with area 0.1 cm2. For simulations of infection intensity, this was a 25 cm tall plant with infection intensity 1.0. Next, we simulated 100 trajectories from the model posteriors using observed weather data for each site and a step size of seven days. In these simulations, we restricted plant height, pustule area, and infection intensity to remain greater or equal to 5 cm, 0 mm2, and 0.1, respectively. We included the random effect of site (but not plant identity) when making these predictions. After confirming that these simulation results qualitatively recapitulated observed patterns, we simulated sets of 100 trajectories for future climate data sets. We again excluded the random effect of plant identity when making these projections. We used the same set of 100 random number seeds for each set of 100 trajectory simulations.While performing these simulations, we extracted weather variables as follows: When extracting weather variables from observed data sets, we defined the observation period as 00:00:00 on the first day in a step to 00:00:00 on the last day in a step (e.g. midnight on July 1st to midnight on July 8th). We extracted weather metrics as described in the “Weather data” subsection of the “Data collection” section above. Applying these same methods to the projected weather data is not possible, as these data sets only give only daily mean values. Instead, we calculated mean temperature as the mean of daily mean temperature values, maximum temperature as the maximum of daily maximum temperature values, and minimum temperature as the minimum of daily minimum temperature values. Similarly, we used daily mean temperature and relative humidity to calculate daily mean absolute humidity, and calculated mean absolute humidity as the mean of these values. We applied this same procedure to calculate mean rainfall.Spore dispersalTo infer how spore deposition is related to source plant infection intensity, wind speed, and wind direction, we fit a tilted gaussian plume model22 (TGPM) to our spore deposition data. The equation specifying the TGPM is as follows:$$begin{aligned} & dleft(I,H,s,X,Yright)=frac{I k {W}_{s}}{2 pi s {sigma }^{2} }{e}^{frac{{-Y}^{2}}{2{sigma }^{2}}-frac{{(H-frac{{W}_{s }X}{s})}^{2}}{2{sigma }^{2}}} \ & for quad X >0 0 \ & for quad Xle 0.end{aligned}$$This equation describes the concentration of spores deposited at a given point ((d)) as a function of wind speed ((s)), the infection intensity of the source plant ((I)), a constant relating (I) to the source concentration of spores ((k)), and the coordinates of the point (X,Y), relative to the source. The coordinate system has the source at the origin, the X-axis parallel to the wind direction (with wind flowing in a positive direction), the Y-axis perpendicular to the X-axis on the plane defined by the ground (assumed to be flat), and the vertical Z-axis orthogonal to the X and Y axes. The shape of the three dimensional spore plum emanating from the source is defined by (s), along with constants specifying the falling velocity of spores (({W}_{s})), the height of the source ((H)), and the standard deviation of spore dispersion (({sigma }^{2})) in the horizontal and vertical directions. Following Levine and Okubu22, we calculate ({sigma }^{2}) as a function of a diffusion coefficient, (A), assuming that the variance in spore distribution increases linearly with time. We also assume that spores diffuse along the Y and Z axes at equal rates.$${sigma }^{2}=frac{2AX}{s}.$$We fit the parameters (k, {W}_{s},) and (sigma ) by minimizing the sum of squared differences between model predictions and the spore deposition data we collected from the spore traps arranged around diseased focal plants. We assumed that all spores deposited on a spore trap originated from the associated diseased focal plant, and that spore deposition occurred between noon on the day of deployment and ceased at noon two days post-deployment as the spore traps lost their adhesive properties and became saturated with dust. When predicting spore deposition using the TGPM, we set (I) equal to the infection intensity of the source diseased focal plant on the day of spore trap deployment and (H) equal to half of the maximum height of the source diseased focal plant at the day of spore trap deployment. As wind speed and direction were recorded every 5 min, we calculated total spore deposition as the sum of predicted deposition across five minute time windows. In each individual estimate, we set (s) equal to wind speed and calculated (X) and (Y) from wind direction and the location of the spore trap relative to the focal diseased plant (see Supplementary Text 3).TransmissionUsing the fitted TGPM we next sought to connect patterns of spore dispersal to transmission. The first step in this process was to align data on the spatiotemporal spread of disease with the locations of all healthy and infected plants (including non-focal plants) within each transect so that the relative locations of these could be determined (see Supplementary Text 4). Next, we built a dataset (‘corrected plant height data’) describing the height of all plants on the dates of epidemiological surveys. We preferentially used actual measurements from the initial population survey, epidemiological surveys, or focal plant measurements where they existed. For the majority of plants, height was not recorded on the date of a given epidemiological survey. In these cases, we fore- or hind-casted height from a previous or future record of plant height (whichever was closer) using the fitted generalized additive model of plant growth and observed weather data spanning the period between the date of interest and the closest record. We included the random effect of study site when making predictions, but not the random effect of plant identity as we wished to forecast heights of both focal and non-focal plants and only focal plant data was used to fit the model of plant growth.After completing these steps, we next built a ‘transmission data set’ connecting weather conditions, predicted spore deposition, plant height, and infection occurrence across the time windows between each epidemiological survey. We included separate entries for each healthy plant (focal and non-focal) during each time window. We limited the scope of this data to the periods for each transect in which all infected plants (that were not seedlings) were tracked as diseased focal plants. Using noon on the date of one epidemiological survey as the beginning of the time window, and noon on the date of the next epidemiological survey as the end of the time window, we extracted the same weather variables as described in “Plant growth and within-host disease spread”. Plant height values were extracted from the ‘corrected plant height’ data set. We determined that a plant became infected if a plant with identical coordinates (after data alignment) was recorded as newly infected in the epidemiological survey that took place at the end of the observation period. Total spore deposition for each healthy plant was calculated as the sum of predicted spore deposition originating from each of the plants in the same transect that was infected at the beginning of the time window. To predict the spore deposition experienced by a healthy plant that originated from an individual diseased plant, we used the TGMP(.) We extracted the height of the diseased plant at the beginning of the time window and set (H) to half of this value. Likewise, we extracted the infection intensity of the diseased plant at the beginning of the time window, and set I equal to this value. For diseased seedlings that were not tracked as diseased focal plants, we set I = 0.1, as most were observed to be lightly infected. There were 14 instances of missing infection intensity data for a diseased focal plant on a certain date. In these cases, we fore- or hind-casted infection intensity for that plant using the generalized additive model of infection intensity progression fit in “Plant growth and within-host disease spread” (while including the random effect for study site, but not for plant identity), the closest observation of infection intensity for that plant, and the weather data spanning the time period from the missing observation to the closest measurement. We set the values of (k, {W}_{s},) and (sigma ) to those fit in the “spore deposition” section above. As wind speed and direction were recorded every five minutes, we summed spore deposition over 5 min windows that matched the resolution of the wind speed and direction data. For each of these windows, we predicted spore deposition after setting the values of X and Y based on the relative location of the healthy and diseased plants and the direction of the wind (see Supplementary Information).Using this data set, we modeled infection outcomes in a generalized additive model framework. We assumed a binomial response distribution, coding the outcome of an infection as 1 and the outcome of no infection as 0. To account for variation in the duration of observation period, we used a complementary log–log link function, and included the log time as an offset in the model46. We included smooth terms for all weather variables as predictors, along with smooth terms study site and plant identity as random effects (accomplished by setting bs = ”re” in the s() function in mgcv). To infer the connection between spore deposition and infection, we also included a full tensor smooth product between the base 10 logarithm of mean total spore deposition per day and maximum plant height. When fitting the tensor, we fixed the smoothing parameter for the marginal smooth for log10 total spore deposition to (1times {10}^{10}). This effectively sets the shape of this marginal smooth to be linear, preventing any overfitting involving non-monotonic effects of spore deposition. Height was included as a factor in the tensor for two reasons. First, we calculated spore deposition for a point. In reality, the spore deposition experienced by a plant depends on the size of the ‘target’ that it presents. As such, larger plants are presumably challenged by a greater number of spores than a smaller plant for the same level of predicted spore deposition. Secondly, plant size could be correlated with quantitative resistance, because leaf age can influence pustule development20. To achieve computational feasibility, we fit this model using the bam() function in the mgcv R package and the fast restricted maximum likelihood parameter estimation method (“fREML”). We again implemented variable selection so that insignificant predictor terms would be effectively removed from the models via the double penalty approach.To infer the effects of climate change and to illustrate the relationships between spore deposition, plant height, and infection risk, we predicted odds of infection using the estimated parameters of the fitted GAM (Fig. 5; Supplementary Fig. S14) and various weather data sets. In all of these projections, we included the random effect of site but excluded the random effect of plant identity.Epidemiological modelModel descriptionTo predict how the population level epidemiological dynamics of flax-rust might be affected in various climate change scenarios, we constructed a spatio-temporal epidemic model. This model tracks the infection status, infection intensity, and height of plants over time as infection spreads within the modeled population. There are three components required to simulate the model: (1) initial conditions describing the locations and starting states of all plants, (2) weather data, and (3) models of spore dispersal, transmission, infection intensity progression, and plant growth.The initial conditions for individual simulations of the epidemic model were in all cases constructed to represent either the BT or GM study populations on the date of an early epidemiological survey. For BT, this date was chosen as the first which had complete weather data (spanning 00:000:00 to 24:00:00), which was 6/24/2020. For GM, this date was set to 7/2/2020. Weather data was available for earlier dates, but a sharp increase in prevalence occurred during the prior week, perhaps due to prior missed observations of new infections. The model consistently underestimated observed prevalence when initialized prior to this increase in prevalence. We did not use the HM study population to initialize the model as no epidemiological surveys occurred on a date with complete weather data. Likewise, we do not use the CC study population to initialize the model because not all regions of the transect were mapped (see “Population mapping” sub-section in “Data collection” above).The construction of initial conditions to represent a given site proceeded as follows. The locations of all plants were taken from the ‘corrected plant height’ data set. Plants observed to be infected on the date of the first epidemiological survey were initialized as infected, and all other plants were initialized as uninfected. Following the first epidemiological survey, all infected plants (except for several that were approximately 5 cm or less in height) were designated diseased focal plants and their infection intensity was recorded. We used these measurements as the initial infection intensities of infected plants. We fore- or hind-casted infection intensities for any missing records using the method described in the transmission sub-section of the statistical analyses section above. Plant heights were taken from the ‘corrected plant height’ data set.The weather data component of the model is used to predict spore dispersal, transmission, infection intensity progression, and plant growth in conjunction with the fitted statistical models of these processes. As such, it needs to describe mean weather conditions over the period spanning a simulation step. Additionally, because the TGPM is calibrated to make predictions over 5 min intervals, the weather data must also describe average wind speed and direction over five minute windows.We used the fitted statistical model relating weather to the processes spore dispersal (TGPM), transmission, infection intensity progression, and plant growth to simulate these processes in the epidemiological model. These models contain random effects for both site and plant identity. The random effect of plant identity was excluded when making predictions about infection intensity progression and plant growth as random effects were only fit for focal plants. We included random effects of plant identity when making predictions using the transmission model because random effects were fit for all plants in the modeled populations.Using these components, the model simulates changes in infection status, infection intensity and plant height in time steps of 7 days. The simulation for each time step proceeds as follows:

(1)

First, the model predicts the number of spores deposited at the location of each healthy plant over the course of the entire time step. As all infected plants are assumed to be sources of spore deposition, the model predicts total spore deposition as the sum of spore deposition originating from each infected plant. The spore deposition originating from a single infected plant is calculated as the sum of predictions over 5 min windows spanning the time step. Individual predictions were generated using the TGPM with I is set to the infection intensity of the source plant at the beginning of the time step, H set to 0.5 times the height of the source plant at the beginning of the time step, s set to average wind speed, and X and Y set according to locations of the source and target plants and wind direction (see Supplementary Information). The parameters (k, {W}_{s},) and (sigma ) were set to the values fit to spore trap data.

(2)

Next, for each healthy plant, the model infers the probability of that plant becoming infected from predicted spore deposition at that plant’s location and mean weather metrics using the transmission model. Following this, a random number in [0,1] is drawn, and if it is equal or lower to the probability of infection that plant becomes infected. All newly infected plants are set to have an infection intensity of 1.0.

(3)

After simulating the transmission process, the model next simulates changes in infection intensity. For each plant, we simulated a prediction of the change in infection intensity per day by drawing randomly from the posterior of the infection intensity progression model. We used the infection intensity and height of that plant at the beginning of the time step along with mean weather metrics spanning the time step as the model inputs. We multiplied the predicted change in infection intensity per day by the length of the time step and added it to the plant’s infection intensity at the beginning of the time step to calculate the new infection intensity of the plant.

(4)

Finally, we used a similar process to simulate change in plant height. For each plant, we simulated a prediction of the change in height per day by drawing randomly from the posterior of the plant growth model. Again, we used the height of that plant at the beginning of the time step along with mean weather metrics spanning the time step as model inputs. To calculate the new height of the plant, we multiplied the predicted change in height per day by the length of the time step and added this number to the plant’s height at the beginning of the time step.

There are two sources of stochasticity within the epidemiological model. The first of these sources pertains to the probability of infection. Although the predicted probability of infection for a certain combination of predictor values is always the same, we compare this probability to a random number to determine if infection occurs. The second source of stochasticity pertains to predicted changes in infection intensity and plant height. We predict rates of infection intensity change and plant growth by drawing randomly from the posterior distributions of the corresponding fitted models. These sources of stochasticity allow us to understand the variation in epidemiological model predictions via repeated simulation.SimulationTo validate the epidemiological model, we simulated 10 epidemiological trajectories for the BT and GM sites and compared the results to our observations of flax-rust spread to ensure a qualitative match. For each of these sites, we initialized the model using data corresponding to the site as described above and used observed weather (from the same site) data to run the simulation. When extracting weather variables from observed weather data sets, we defined the observation period as 00:00:00 on the first day in a step to 00:00:00 on the last day in a step (e.g. noon on July 1st to noon on July 8th). We included the random effects of site and plant identity when making predictions from the transmission model, and the random effects of site when making predictions from the infection intensity and plant growth models.To infer the effects of climate change of flax-rust epidemiology, we again initialized the model using data corresponding to the for BT and GM sites, and then simulated 10 epidemiological trajectories for each set of weather data (RCP4.5 and 8.5 emissions scenarios and the years 2020, 2030, 2040, 2045, 2050, 2060 and 2070). We extracted weather metrics as described in the “Weather data” subsection of the “Data collection” section above. We used observed records of wind speed and direction to simulate spore deposition using the TGPM as the projected data sets lack these measures. As before we included the random effects of site and plant identity when making predictions from the transmission model, and the random effects of site when making predictions from the infection intensity and plant growth models. For each set of 10 model simulations corresponding to a unique site and weather data set combination, we used the same 10 random number seeds. More

Blanton, L. S. The rickettsioses: A practical update. Infect. Dis. Clin. North Am. 33, 213–229 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Parola, P. et al. Update on tick-borne rickettsioses around the world: A geographic approach. Clin. Microbiol. Rev. 26, 657–702 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
Robinson, M. T., Satjanadumrong, J., Hughes, T., Stenos, J. & Blacksell, S. D. Diagnosis of spotted fever group Rickettsia infections: The Asian perspective. Epidemiol. Infect. 147, e286 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Graves, S. & Stenos, J. Rickettsioses in Australia. Ann. N. Y. Acad. Sci. 1166, 151–155 (2009).ADS 
PubMed 
Article 

Google Scholar 
Niang, M. et al. Prevalence of antibodies to Rickettsia conorii, Ricketsia africae, Rickettsia typhi and Coxiella burnetii in Mauritania. Eur. J. Epidemiol. 14, 817–818 (1998).CAS 
PubMed 
Article 

Google Scholar 
Parola, P. Tick-borne rickettsial diseases: Emerging risks in Europe. Comp. Immunol. Microbiol. Infect. Dis. 27, 297–304 (2004).PubMed 
Article 

Google Scholar 
Nanayakkara, D. M., Rajapakse, R. P. V. J., Wickramasinghe, S. & Kularatne, S. A. M. Serological evidence for exposure of dogs to Rickettsia conorii, Rickettsia typhi, and Orientia tsutsugamushi in Sri Lanka. Vector Borne Zoon. Dis. Larchmt. N 13, 545–549 (2013).Article 

Google Scholar 
Brown, L. D. & Macaluso, K. R. Rickettsia felis, an emerging flea-borne rickettsiosis. Curr. Trop. Med. Rep. 3, 27–39 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Newton, P. N. et al. A prospective, open-label, randomized trial of doxycycline versus azithromycin for the treatment of uncomplicated murine typhus. Clin. Infect. Dis. 68, 738–747 (2019).CAS 
PubMed 
Article 

Google Scholar 
Vallee, J. et al. Contrasting spatial distribution and risk factors for past infection with scrub typhus and murine typhus in Vientiane City, Lao PDR. 4 (2010).Akram, S. M., Jamil, R. T. & Gossman, W. G. Rickettsia Akari (2021).Dong, X., El Karkouri, K., Robert, C., Raoult, D. & Fournier, P.-E. Genome sequence of Rickettsia australis, the agent of Queensland tick typhus. J. Bacteriol. 194, 5129 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Fournier, P.-E. & Raoult, D. Current knowledge on phylogeny and taxonomy of Rickettsia spp. Ann. N. Y. Acad. Sci. 1166, 1–11 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Legendre, K. P. & Macaluso, K. R. Rickettsia felis: A review of transmission mechanisms of an emerging pathogen. Trop. Med. Infect. Dis. 2, E64 (2017).PubMed 
Article 

Google Scholar 
Murray, G. G. R., Weinert, L. A., Rhule, E. L. & Welch, J. J. The phylogeny of rickettsia using different evolutionary signatures: How tree-like is bacterial evolution?. Syst. Biol. 65, 265–279 (2016).PubMed 
Article 

Google Scholar 
Shpynov, S. N., Fournier, P., Pozdnichenko, N. N., Gumenuk, A. S. & Skiba, A. A. New approaches in the systematics of rickettsiae. New Microbes New Infect. 23, 93–102 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Shpynov, S. et al. Detection of a rickettsia closely related to Rickettsia aeschlimannii, ‘Rickettsia heilongjiangensis’, Rickettsia sp. strain RpA4, and Ehrlichia muris in ticks collected in Russia and Kazakhstan. J. Clin. Microbiol. 42, 2221–2223 (2004).PubMed 
PubMed Central 
Article 

Google Scholar 
Aung, A. K., Spelman, D. W., Murray, R. J. & Graves, S. Review article: Rickettsial infections in Southeast Asia: Implications for local populace and febrile returned travelers. Am. J. Trop. Med. Hyg. 91, 451–460 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Rodkvamtook, W. et al. Scrub typhus outbreak in Chonburi Province, Central Thailand, 2013. Emerg. Infect. Dis. 24, 361–365 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Robinson, M. T., Vongphayloth, K., Hertz, J. C., Brey, P. & Newton, P. N. Tick-transmitted human infections in Asia. Microbiol. Aust. 39, 203–206 (2018).Article 

Google Scholar 
Bartoshevic, E. To the issue of rickettsioses. Health Care Kazakhstan 3, 20–24 (1952) (in Russian).
Google Scholar 
Kereyev, N. Human natural focal diseases in Kazakhstan. Alma-ata (1961) (in Russian).Arkhangelskiy, D. Experimental study of tick-borne rickettsial pathogen in Almaty region. In Collection of Scientific Papers of the Institute of Microbiology and Virologoy Vol 4. Physiology and ecology of micro-organisms. Almta-ata 176–85 (1961) (in Russian).Kyraubayev, K. et al. Study of Dermacentor marginatus ticks for Rickettsiae in Central Kazakhstan. Proc. ASM (2014).Shpynov, S. et al. Detection and identification of spotted fever group Rickettsiae in dermacentor ticks from Russia and Central Kazakhstan. Eur. J. Clin. Microbiol. Infect. Dis. 20, 903–905 (2001).CAS 
PubMed 
Article 

Google Scholar 
Shpynov, S., Rudakov, N. & Yastrebov, V. Identification of new genotypes of rickettsia tick-borne spotted fever group in the south of the Ural, Siberia, Far East and Kazakhstan. Epidemiol. Infect. Dis. 1, 23–27 (2005).
Google Scholar 
Hay, J. et al. Biosurveillance in Central Asia: Successes and challenges of tick-borne disease research in Kazakhstan and Kyrgyzstan. Front. Public Health 4, 4 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Yegemberdiyeva, R. & Shapieva, Z. Clinical and epidemiological characteristic of tick-borne rickettsiosis in Kazakhstan. Abstract Book of the International Conference on Zoonoses. Ulaanbaatar 48–51 (2008).Rudakov, N. V., Shpynov, S. N., Samoilenko, I. E. & Tankibaev, M. A. Ecology and epidemiology of spotted fever group Rickettsiae and new data from their study in Russia and Kazakhstan. Ann. N. Y. Acad. Sci. 990, 12–24 (2003).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Sansyzbayev, Y. et al. Survey for Rickettsiae within fleas of Great Gerbils, Almaty Oblast, Kazakhstan. Vector Borne Zoon. Dis. Larchmt. N 17, 172–178 (2017).Article 

Google Scholar 
Kazakhstan Scientific Practical Center of Sanitary Epidemiological Expertise and Monitoring. Almaty. Epidemiological situation of infectious diseases in the Republic of Kazakhstan from 2016. Annual Report (2016) (in Russian).CDC. https://www.cdc.gov/vhf/omsk/index.html (2022).Turebekov, N. et al. Prevalence of Rickettsia species in ticks including identification of unknown species in two regions in Kazakhstan. Parasit. Vectors 12, 1–16 (2019).Article 

Google Scholar 
Jones, K. E. et al. Global trends in emerging infectious diseases. Nature 451, 990–993 (2008).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Tomassone, L. et al. Neglected vector-borne zoonoses in Europe: Into the wild. Vet. Parasitol. 251, 17–26 (2018).PubMed 
Article 

Google Scholar 
Schex, S., Dobler, G. & Riehm, J. Rickettsia spp. in wild small mammals in Lower Bavaria, South-Eastern Germany. Vector Borne Zoon. Dis. 11, 493–502 (2011).Article 

Google Scholar 
Tukhanova, N. et al. Molecular characterisation and phylogeny of Tula virus in Kazakhstan. Viruses 14, 1258 (2022).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Wölfel, R., Essbauer, S. & Dobler, G. Diagnostics of tick-borne rickettsioses in Germany: A modern concept for a neglected disease. Int. J. Med. Microbiol. 298, 368–374 (2008).Article 
CAS 

Google Scholar 
Fournier, P. E., Roux, V. & Raoult, D. Phylogenetic analysis of spotted fever group Rickettsiae by study of the outer surface protein rOmpA. Int. J. Syst. Bacteriol. 48(Pt 3), 839–849 (1998).CAS 
PubMed 
Article 

Google Scholar 
Jado, I. et al. Molecular method for identification of Rickettsia species in clinical and environmental samples. J. Clin. Microbiol. 44, 4572–4576 (2006).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hall, T. A. BioEdit a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp. Ser. 41, 95–98 (1999).CAS 

Google Scholar 
Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35, 1547–1549 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Turebekov, N. et al. Occurrence of anti-Rickettsia spp. antibodies in hospitalized patients with undifferentiated febrile illness in the southern region of Kazakhstan. Am. J. Trop. Med. Hyg. 104, 2000–2008 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
SPC SEEM. Kazakhstan Scientific Practical Center of Sanitary Epidemiological Expertise and Monitoring, Almaty, Kazakhstan (2021).Yamamoto, Y. PCR in diagnosis of infection: detection of bacteria in cerebrospinal fluids. Clin. Vaccine Immunol. 9, 508–514 (2002).CAS 
Article 

Google Scholar 
Turebekov, N. et al. Prevalence of Rickettsia species in ticks including identification of unknown species in two regions in Kazakhstan. Parasit. Vectors 12, 197 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Gajda, E. et al. Spotted fever Rickettsiae in wild-living rodents from south-western Poland. Parasit. Vectors 10, 413 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Essbauer, S., Hofmann, M., Kleinemeier, C., Wölfel, S. & Matthee, S. Rickettsia diversity in southern Africa: A small mammal perspective. Ticks Tick-Borne Dis. 9, 288–301 (2018).PubMed 
Article 

Google Scholar 
Weinert, L. A., Werren, J. H., Aebi, A., Stone, G. N. & Jiggins, F. M. Evolution and diversity of Rickettsia bacteria. BMC Biol. 7, 6 (2009).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
El Karkouri, K., Ghigo, E., Raoult, D. & Fournier, P.-E. Genomic evolution and adaptation of arthropod-associated Rickettsia. Sci. Rep. 12, 3807 (2022).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Zemtsova, G. E., Montgomery, M. & Levin, M. L. Relative sensitivity of conventional and real-time PCR assays for detection of SFG Rickettsia in blood and tissue samples from laboratory animals. PLoS One 10, e0116658 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Burri, C., Schumann, O., Schumann, C. & Gern, L. Are Apodemus spp. mice and Myodes glareolus reservoirs for Borrelia miyamotoi, Candidatus Neoehrlichia mikurensis, Rickettsia helvetica, R. monacensis and Anaplasma phagocytophilum?. Ticks Tick-Borne Dis. 5, 245–251 (2014).CAS 
PubMed 
Article 

Google Scholar 
Tadin, A. et al. Molecular survey of zoonotic agents in rodents and other small mammals in Croatia. Am. J. Trop. Med. Hyg. 94, 466–473 (2015).PubMed 
Article 
CAS 

Google Scholar 
Karbowiak, G., Biernat, B., Stańczak, J., Szewczyk, T. & Werszko, J. The role of particular tick developmental stages in the circulation of tick-borne pathogens affecting humans in Central Europe. 3. Rickettsiae. Ann. Parasitol. 62, 89–100 (2016).PubMed 

Google Scholar 
Zemtsova, G., Killmaster, L. F., Mumcuoglu, K. Y. & Levin, M. L. Co-feeding as a route for transmission of Rickettsia conorii israelensis between Rhipicephalus sanguineus ticks. Exp. Appl. Acarol. 52, 383–392 (2010).CAS 
PubMed 
Article 

Google Scholar 
Rehácek, J., Urvölgyi, J., Kocianová, E. & Jedlicka, L. Susceptibility of some species of rodents to Rickettsiae. Folia Parasitol. (Praha) 39, 265–284 (1992).
Google Scholar 
Rehácek, J., Zupancicová, M., Kovácová, E., Urvölgyi, J. & Brezina, R. Study of rickettsioses in Slovakia. III. Experimental infection of Apodemus flavicollis Melch. by Rickettsiae of the spotted fever (SF) group isolated in Slovakia. J. Hyg. Epidemiol. Microbiol. Immunol. 21, 306–313 (1976).PubMed 

Google Scholar 
Biernat, B., Stańczak, J., Michalik, J., Sikora, B. & Wierzbicka, A. Prevalence of infection with Rickettsia helvetica in Ixodes ricinus ticks feeding on non-rickettsiemic rodent hosts in sylvatic habitats of west-central Poland. Ticks Tick-Borne Dis. 7, 135–141 (2016).PubMed 
Article 

Google Scholar 
Stańczak, J. et al. Prevalence of infection with Rickettsia helvetica in feeding ticks and their hosts in western Poland. Clin. Microbiol. Infect. Off. Publ. Eur. Soc. Clin. Microbiol. Infect. Dis. 15(Suppl 2), 328–329 (2009).
Google Scholar 
Barandika, J. F. et al. Tick-borne zoonotic bacteria in wild and domestic small mammals in northern Spain. Appl. Environ. Microbiol. 73, 6166–6171 (2007).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Spitalská, E., Boldis, V., Kostanová, Z., Kocianová, E. & Stefanidesová, K. Incidence of various tick-borne microorganisms in rodents and ticks of central Slovakia. Acta Virol. 52, 175–179 (2008).PubMed 

Google Scholar 
Guo, L.-P. et al. Rickettsia raoultii in Haemaphysalis erinacei from marbled polecats, China-Kazakhstan border. Parasit. Vectors 8, 461 (2015).PubMed 
PubMed Central 
Article 

Google Scholar  More