Ecology

We observed a clear reduction in the acoustic activity of sperm whales and beaked whales during the period when sonar signals were recorded at Station 5, indicating that whales ceased foraging in this area while military sonars were in use. The acoustic detection rate of sperm whales returned to pre-exercise baseline levels within the days following the CF16 exercise, while the observed reduction in beaked whale acoustic activity was more prolonged. Detection rates of Cuvier’s beaked whale clicks remained low throughout the 8-day period immediately following the exercise, and UMBW clicks were largely absent during this period. This study is observational and limited to showing correlation rather than cause and effect; nonetheless, these results are consistent with previous experimental research on the responses of beaked whales to simulated and real military sonars and suggest that whales were disturbed from normal foraging behaviour and likely displaced from the affected area during the CF16 exercise.The scale and duration of sonar use recorded during this study provides important context for the observed results. Much of the experimental work conducted to date on the responses of beaked whales and other odontocetes to sonar has involved controlled exposure experiments using animal-borne tags to record the fine-scale movements and acoustic behavior of individuals, allowing responses to be examined on the scale of minutes to hours e.g.,7,8,10. Experimental exposures to simulated sonar signals lasting approximately 15–30 min have elicited pronounced avoidance responses in Blainville’s beaked whales7, Cuvier’s beaked whales8, Baird’s beaked whales16, and northern bottlenose whales9,10. Generally, these studies were focused on the onset of the response and did not always assess the duration over which altered behaviour continued. However, the absence of foraging behaviour for several hours following exposure was noted in some cases, and focal animals performed sustained directed movement away from the exposure location during this time, covering distances of up to tens of kilometers10. In broader-scale studies examining responses of Blainville’s beaked whales to real multi-ship naval training operations on the Atlantic Undersea Test and Evaluation Center (AUTEC) in the Bahamas, displacements of up to 68 km were observed, lasting 2–4 days before whales returned to foraging in the area where they were exposed7,28. In the present study, the duration of naval sonar activity recorded during the CF16 exercise was considerably more prolonged, with bouts of sonar continuing for up to 13 consecutive hours and occurring repeatedly over an 8-day period. Although we can only make inference on species-level rather than individual-level responses based on the absence of clicks in our recordings, it is plausible that military sonar activity at this scale led to wide spatial avoidance of the affected area over an extended period.The absence of sperm whale click detections in the Station 5 recordings for 6 consecutive days during the CF16 exercise is notable; few prior studies have demonstrated sustained changes in foraging behaviour or substantial displacement of sperm whales following sonar exposure. Behavioural response studies conducted in northern Norway using controlled experimental exposures showed varying responses by sperm whales, which included changes in orientation and direction of horizontal movement, changes in acoustic behaviour, and altered dive profiles23. Exposure to lower frequency sonar signals in the 1–2 kHz range generally prompted stronger responses, including a reduction in foraging effort or transition from a foraging to non-foraging state, while exposure to higher frequency sonar signals in the range of 6–7 kHz did not appear to trigger changes in foraging behaviour21,29. More recently, Isojunno et al.30 quantified the responses of sperm whales to continuous and pulsed active sonars, and found that sound exposure level was more important than amplitude in predicting a change in foraging effort. We were not able to investigate differential responses to frequency or other sonar characteristics in this study, due to the observational nature of the study and the absence of sperm whale clicks throughout most of the exercise period. Likewise, we cannot exclude the physical presence of ships, aircraft, and submarines in the area or additional types of noise produced during maneuvers as potential factors contributing to the cessation of sperm whale and beaked whale click production and foraging behaviour.The observed changes in acoustic activity were more easily quantified for sperm whales than for beaked whales, due to higher baseline hourly presence of sperm whale clicks in the recordings. Sperm whales produce powerful echolocation clicks throughout their foraging dives, which can be recorded at ranges of 16 km or more31, and a single individual foraging in the vicinity of a hydrophone may be detected continuously throughout multiple dive cycles. Our analysis was based on sperm whale click detections that met a threshold signal-to-noise ratio (SNR), and the results therefore provide a minimum estimate of sperm whale presence in the vicinity of the recorder. Reporting results at the level of hourly presence rather than the number of individual click detections largely mitigated the effects of excluding low-SNR clicks recorded at greater distances from the hydrophone or during higher ambient noise conditions. Likewise, the presence of sperm whales on an hourly time scale is not likely to be substantially underestimated when recordings are collected using a low duty cycle32. By contrast, beaked whales produce echolocation clicks at higher frequencies and lower source levels, with highly directional beam patterns33. These clicks are likely only detected at ranges of up to approximately 4 km when the whale is oriented toward the hydrophone, and at lesser distances when clicks are received off-axis34. As a result, there is greater variability and lower baseline detection rates of beaked whale clicks on fixed passive acoustic recorders, which reduces statistical power to assess temporal changes in acoustic activity. Moreover, the duty-cycled recording schedule used at Station 5 provided only 65 s of high-frequency data 3 times per hour, and the presence of beaked whales is likely to be underestimated by this duty cycle, with potentially greater underestimation of Mesoplodont species compared to Cuvier’s beaked whales35.Continuous recordings were collected at the East Gully and Central Gully recording sites, but included only partial temporal coverage of the exercise period and no pre-exercise baseline data. No comparable recordings were available from these locations in a prior or subsequent year to form a control dataset. As a result, we were not able to use these datasets to assess changes in acoustic activity associated with the CF16 exercise. A slight decrease in hourly presence of northern bottlenose whale clicks in the Central Gully recordings occurred on September 19th–20th, 2016; however, we are aware that an oceanographic research vessel was coincidentally in the area deploying scientific instrumentation in close proximity to the Central Gully recording site on these dates, creating an additional source of potential disturbance. Despite these limitations, we included an analysis of the recordings collected at the East and Central Gully sites for two reasons: first, to provide perspective on the geographic extent over which activities associated with the CF16 exercise occurred; and second, to illustrate the diversity in beaked whale species composition at different locations across the region. Analysis of the recordings for sonar signals revealed that higher levels of sonar activity occurred near the Station 5 recording site than near the East or Central Gully locations. Due to the distance between recording sites and the timing of the sonar signals recorded, it appears that the recorded sonar signals came from multiple source locations over the duration of the exercise. Recordings from Central Gully contained the fewest sonar signals and lowest measured received levels, likely due to the deliberate avoidance of the Gully MPA and surrounding area by exercise participants during CF16. The Gully was established as an MPA in 2004, and is one of three adjacent canyons on the eastern Scotian Shelf currently designated as critical habitat areas for the endangered Scotian Shelf population of northern bottlenose whales36. The Station 5 recording site was located approximately 300 km to the southwest, and experienced higher levels of naval sonar activity during CF16. However, none of the locations were chosen specifically to monitor CF16, and we do not have access to information on the general exercise areas used, specific locations of naval vessels, submarines, or aircraft participating in the CF16 exercise, or the source levels of transmitted sonar signals. Due to the opportunistic nature of the recordings, the received levels of sonar signals measured at Station 5 likely do not represent the highest sound levels introduced into the marine environment during the CF16 exercise.Unlike many areas where behavioural responses to sonar are commonly studied, there are no instrumented naval training ranges off eastern Canada, and cetaceans inhabiting this region are unlikely to be accustomed to regularly hearing naval active sonars. Other than during the CF16 exercise, sonar signals were not noted during a large-scale analysis of cetacean call occurrence and soundscape characterization in 2 years of recordings collected at Station 5 and numerous other passive acoustic monitoring sites off eastern Canada26. Exposure context and familiarity with a signal may be important factors influencing an individual’s response to acoustic disturbance15. Experimental research on Cuvier’s beaked whales near a U.S. naval training range located off southern California demonstrated possible distance-mediated effects of sonar exposure, with more pronounced behavioural responses occurring with closer source proximity, even when received levels from the closer source were likely lower than those from more distant, high-powered sonar transmissions, which did not elicit as strong a response15. The movement and predictability of the sound source as well as the timing and duration of sonar transmissions may also be important factors influencing the behavioural response15. Whales inhabiting waters off southern California are likely habituated to hearing distant sonar due to routine naval training activities occurring on the range. Conversely, Wensveen et al.10 found that northern bottlenose whales in the eastern North Atlantic exhibited similar responses to simulated sonar signals played at various distances up to 28 km, suggesting that they perceived this novel stimuli as a potential threat even from a distance and at relatively low received levels. Bernaldo de Quiros et al.5 hypothesized that beaked whales not regularly exposed to active sonar signals may respond more strongly, both physiologically and behaviourally, which poses a concern for a region where military training activities involving the use of sonar are relatively infrequent, but occur periodically in the form of large-scale exercises involving the extensive use of active sonars and creating significant potential for acoustic disturbance.Behavioural disturbance due to anthropogenic noise may have energetic, health, and fitness consequences for deep-diving odontocete species. Disruption of normal diving patterns creates energetic costs due to the significant investment in each dive and the reduction of time available for prey intake when foraging dives are interrupted. Recent studies on the functional relationship between beaked whales and deep-sea prey resources suggest that certain characteristics of prey, including minimum size and density thresholds, are required for beaked whales to successfully meet their energetic needs12,37. While the distribution and characteristics of deep-sea prey are challenging to study and largely unknown in most regions, considerable environmental heterogeneity may be present, causing the quality of foraging habitat to vary significantly over even small horizontal scales12,37. This patchiness in habitat quality has important implications for behavioural disturbance, as even short-term displacement from high-quality habitat areas can affect the fitness of individuals and potentially lead to population-level consequences13.In addition to the consequences of sublethal disturbance, it is important to note that the likelihood of observing more acute impacts of exposure to naval active sonar, including injuries or fatalities, is extremely low in offshore regions. Individual and mass strandings of beaked whales and other cetaceans associated with military activities have typically been documented on oceanic islands with populated coastlines1,3,6. Factors affecting the probability that cetacean carcasses will wash ashore include buoyancy and decomposition rates in local water conditions, oceanic surface currents, the topography of coastlines, and the location of habitat relative to shore6. Off Nova Scotia, potential beaked whale and sperm whale habitat (consisting of water depths greater than 500 m) is located more than 100 km from the coastline, and injuries or fatalities occurring in deep water habitat in this region are unlikely to result in observed strandings. Stranding incidents involving sperm whales and beaked whales have been reported in Nova Scotia, but the cause of mortality is usually unknown38. Cetacean mortality is highly underestimated even in the aftermath of catastrophic events such as large oil spills39, and a lack of observed injuries or mortalities following offshore military activities should not be construed as evidence that no direct or immediate harm was caused.This study offered a unique opportunity to use existing passive acoustic monitoring (PAM) data to assess disturbance of poorly-known odontocete species during a real-world, large-scale military sonar exercise in a region where military sonar use at this scale is relatively uncommon. Ideally, a PAM study designed to examine disturbance in this context would collect continuous rather than duty-cycled recordings, and include ample baseline data surrounding the period of interest as well as in prior and subsequent years. Additionally, multiple acoustic sensors arranged in a dense array surrounding exercise locations would provide further insight into the spatial context of exposure and patterns of disturbance. Despite the data limitations in the present study, our results demonstrate that changes in odontocete foraging behaviour associated with acute, large-scale disturbance may be evident in PAM data even at low duty cycles. The nature of the observed effect (e.g., temporary disruption of foraging, spatial displacement, or more acute injury or distress) remains unknown, as do the number of individuals affected and the longer-term health and fitness implications. Broader baseline data on species occurrence and an improved understanding of species’ ecology and habitat use in the region are necessary for making informed mitigation decisions, allowing key habitat areas to be avoided, and understanding the impacts of naval active sonar exposure in this region on individuals and populations. More

This research focuses on the quantitative and qualitative analyses of cyanobacteria and microalgae present in rainfall during the summer phytoplankton bloom season of August–September 2019. In addition, a continuous episode of rainfall over several days was selected to demonstrate the washout process of microorganisms from the air with rain.Quantity of cyanobacteria and microalgae washed out with rain during the growing seasonCurrently, there is a growing number of scientific articles on cyanobacteria and microalgae in the atmosphere8. Unfortunately, there is a reference methodology for efficiently counting the microorganisms present in the air or in rainfall. A popular method for quantifying cyanobacteria and microalgae in the air is to show the number of taxa found in the collected samples after growth6,31,42,43,44,45,46. In this study, a total of 16 taxa of airborne cyanobacteria and microalgae were found in the samples. In the rainwater samples obtained during the summer of 2019, 11 taxa of cyanobacteria and microalgae were distinguished. The green algae in the rainwater samples included Bracteacoccus sp., Oocystis sp., Coenochloris sp., Chlorella sp., and Chlorococcum sp., while the cyanobacteria included Leptolyngbya sp., Pseudanabaena sp., Synechococcus sp., and Synechocystis sp. In addition, Chrysochromulina sp., which belongs to Haptophyta, was observed.Other studies recorded the presence of several to several dozen taxa in the air6,31,42,43,44,45,46. Certainly, a number of factors, starting with atmospheric conditions and ending with physical and chemical parameters of the surrounding waters, influence the diversity of cyanobacteria and microalgae in the atmospheric air. Analyzing global trends, only cyanobacteria have been found in the atmosphere of every region of the world31. However, according to Dillon et al.47, cyanobacteria have been detected in clouds at variable abundances between ~ 1% and 50% of the total microbial community. Xu et al.48 found that cyanobacteria constituted only 1.1% of the total bacterial community in clouds. It needs to be highlighted that there is still a lack of research available to provide this type of information for rainfall samples.For the period from July to September 2019, the results showed that the number of cyanobacteria and microalgae cells present in rainfall varied over time (Fig. 1) and ranged between 100 cells L–1 and 342.2 × 103 cells L–1. From July to the end of August, the cell number was relatively low, ranging from 100 cells L–1 to 28.6 × 103 cells L–1. This variability was related to the change in the biomass of blue green algae in the Gulf of Gdańsk (Table S2; Fig. 1). Therefore, this research also shows the close relationship between the processes taking place in the Baltic Sea and the presence of cyanobacteria and microalgae in the atmosphere. As the biomass of cyanobacteria in the Baltic Sea increased, the number of cyanobacteria and microalgae cells in the rainfall samples also increased (***p  More

1.Román-Palacios, C. & Wiens, J. J. Recent responses to climate change reveal the drivers of species extinction and survival. Proc. Natl Acad. Sci. USA 117, 4211–4217 (2020).
Google Scholar 
2.Fuller, A. et al. Physiological mechanisms in coping with climate change. Physiol. Biochem. Zool. 83, 713–720 (2010).
Google Scholar 
3.Sinervo, B. et al. Erosion of lizard diversity by climate change and altered thermal niches. Science 328, 894–899 (2010).CAS 

Google Scholar 
4.Brawn, J. D., Benson, T. J., Stager, M., Sly, N. D. & Tarwater, C. E. Impacts of changing rainfall regime on the demography of tropical birds. Nat. Clim. Change 7, 133–136 (2016).
Google Scholar 
5.Summers, B. A. Climate change and animal disease. Vet. Pathol. 46, 1185–1186 (2009).CAS 

Google Scholar 
6.Randall, C. J. & van Woesik, R. Contemporary white-band disease in Caribbean corals driven by climate change. Nat. Clim. Change 5, 375–379 (2015).
Google Scholar 
7.Munson, L. et al. Climate extremes promote fatal co-infections during canine distemper epidemics in African lions. PLoS ONE 3, e2545 (2008).
Google Scholar 
8.Rohr, J. R. et al. Frontiers in climate change–disease research. Trends Ecol. Evol. 26, 270–277 (2011).
Google Scholar 
9.Zarnetske, P. L., Skelly, D. K. & Urban, M. C. Biotic multipliers of climate change. Science 336, 1516–1518 (2012).CAS 

Google Scholar 
10.Cohen, J. M., Sauer, E. L., Santiago, O., Spencer, S. & Rohr, J. R. Divergent impacts of warming weather on wildlife disease risk across climates. Science 370, eabb1702 (2020).CAS 

Google Scholar 
11.Cornwallis, C. K. et al. Cooperation facilitates the colonization of harsh environments. Nat. Ecol. Evol. 1, 0057 (2017).
Google Scholar 
12.Koenig, W. D. & Dickinson, J. L. (eds) Cooperative Breeding in Vertebrates: Studies of Ecology, Evolution, and Behavior (Cambridge Univ. Press, 2016).13.Groenewoud, F. & Clutton-Brock, T. Meerkat helpers buffer the detrimental effects of adverse environmental conditions on fecundity, growth and survival. J. Anim. Ecol. 90, 641–652 (2020).
Google Scholar 
14.Langwig, K. E. et al. Sociality, density-dependence and microclimates determine the persistence of populations suffering from a novel fungal disease, white-nose syndrome. Ecol. Lett. 15, 1050–1057 (2012).
Google Scholar 
15.Vicente, J., Delahay, R. J., Walker, N. J. & Cheeseman, C. L. Social organization and movement influence the incidence of bovine tuberculosis in an undisturbed high-density badger Meles meles population. J. Anim. Ecol. 76, 348–360 (2007).CAS 

Google Scholar 
16.Bermejo, M. et al. Ebola outbreak killed 5000 gorillas. Science 314, 1564 (2006).CAS 

Google Scholar 
17.Hanya, G. et al. Mass mortality of Japanese macaques in a western coastal forest of Yakushima. Ecol. Res. 19, 179–188 (2004).
Google Scholar 
18.Angulo, E. et al. Allee effects in social species. J. Anim. Ecol. 87, 47–58 (2018).
Google Scholar 
19.Woodroffe, R., Groom, R. & McNutt, J. W. Hot dogs: high ambient temperatures impact reproductive success in a tropical carnivore. J. Anim. Ecol. 86, 1329–1338 (2017).
Google Scholar 
20.Brandell, E. E., Dobson, A. P., Hudson, P. J., Cross, P. C. & Smith, D. W. A metapopulation model of social group dynamics and disease applied to Yellowstone wolves. Proc. Natl Acad. Sci. USA 118, 33649227 (2021).
Google Scholar 
21.Clutton-Brock, T. H. & Manser, M. in Cooperative Breeding in Vertebrates: Studies of Ecology, Evolution, and Behavior (eds Koenig, W. D. & Dickinson, J. L.) 294–317 (Cambridge Univ. Press, 2016).22.Drewe, J. A. Who infects whom? Social networks and tuberculosis transmission in wild meerkats. Proc. R. Soc. B 277, 633–642 (2010).
Google Scholar 
23.Parsons, S. D. C., Drewe, J. A., van Pittius, N. C. G., Warren, R. M. & van Helden, P. D. Novel cause of tuberculosis in meerkats, South Africa. Emerg. Infect. Dis. 19, 2004–2007 (2013).
Google Scholar 
24.Duncan, C., Manser, M., & Clutton-Brock, T. H. Decline and fall: the causes of group failure in cooperatively breeding meerkats. Ecol. Evol. https://doi.org/10.1002/ece3.7655 (2021).25.Drewe, J. A., Foote, A. K., Sutcliffe, R. L. & Pearce, G. P. Pathology of Mycobacterium bovis infection in wild meerkats (Suricata suricatta). J. Comp. Pathol. 140, 12–24 (2009).CAS 

Google Scholar 
26.van Wilgen, N. J., Goodall, V. & Holness, S. Rising temperatures and changing rainfall patterns in South Africa’s national parks. Aquat. Microb. Ecol. 36, 706–721 (2016).
Google Scholar 
27.Conradie, S. R., Woodborne, S. M., Cunningham, S. J. & McKechnie, A. E. Chronic, sublethal effects of high temperatures will cause severe declines in southern African arid-zone birds during the 21st century. Proc. Natl Acad. Sci. USA 116, 14065–14070 (2019).CAS 

Google Scholar 
28.Fischer, E. M., Beyerle, U. & Knutti, R. Robust spatially aggregated projections of climate extremes. Nat. Clim. Change 3, 1033–1038 (2013).
Google Scholar 
29.Bourne, A. R., Cunningham, S. J., Spottiswoode, C. N. & Ridley, A. R. Hot droughts compromise interannual survival across all group sizes in a cooperatively breeding bird. Ecol. Lett. 23, 1776–1788 (2020).
Google Scholar 
30.Van de Ven, T. M. F. N., Fuller, A. & Clutton‐Brock, T. H. Effects of climate change on pup growth and survival in a cooperative mammal, the meerkat. Funct. Ecol. 34, 194–202 (2020).
Google Scholar 
31.Katale, B. Z. et al. Prevalence and risk factors for infection of bovine tuberculosis in indigenous cattle in the Serengeti ecosystem, Tanzania. BMC Vet. Res. 9, 267 (2013).
Google Scholar 
32.Paniw, M., Maag, N., Cozzi, G., Clutton-Brock, T. & Ozgul, A. Life history responses of meerkats to seasonal changes in extreme environments. Science 363, 631–635 (2019).CAS 

Google Scholar 
33.Dwyer, R. A., Witte, C., Buss, P., Goosen, W. J. & Miller, M. Epidemiology of tuberculosis in multi-host wildlife systems: implications for black (Diceros bicornis) and white (Ceratotherium simum) rhinoceros. Front. Vet. Sci. 7, 580476 (2020).
Google Scholar 
34.Patterson, S., Drewe, J. A., Pfeiffer, D. U. & Clutton-Brock, T. H. Social and environmental factors affect tuberculosis related mortality in wild meerkats. J. Anim. Ecol. 86, 442–450 (2017).
Google Scholar 
35.Dubuc, C. et al. Increased food availability raises eviction rate in a cooperative breeding mammal. Biol. Lett. 13, 20160961 (2017).
Google Scholar 
36.Maag, N., Cozzi, G., Clutton-Brock, T. H. & Ozgul, A. Density‐dependent dispersal strategies in a cooperative breeder. Ecology 99, 1932–1941 (2018).
Google Scholar 
37.Ekernas, L. S. & Cords, M. Social and environmental factors influencing natal dispersal in blue monkeys, Cercopithecus mitis stuhlmanni. Anim. Behav. 73, 1009–1020 (2007).
Google Scholar 
38.Ozgul, A., Bateman, A. W., English, S., Coulson, T. & Clutton-Brock, T. H. Linking body mass and group dynamics in an obligate cooperative breeder. J. Anim. Ecol. 83, 1357–1366 (2014).
Google Scholar 
39.Tomlinson, A. J., Chambers, M. A., Wilson, G. J., McDonald, R. A. & Delahay, R. J. Sex-related heterogeneity in the life-history correlates of Mycobacterium bovis infection in European badgers (Meles meles). Transbound. Emerg. Dis. 60, 37–45 (2013).
Google Scholar 
40.Courchamp, F., Grenfell, B. & Clutton-Brock, T. H. Population dynamics of obligate cooperators. Proc. R. Soc. B 266, 557–563 (1999).
Google Scholar 
41.Lerch, B. A., Nolting, B. C. & Abbott, K. C. Why are demographic Allee effects so rarely seen in social animals? J. Anim. Ecol. 87, 1547–1559 (2018).
Google Scholar 
42.Borg, B. L., Brainerd, S. M., Meier, T. J. & Prugh, L. R. Impacts of breeder loss on social structure, reproduction and population growth in a social canid. J. Anim. Ecol. 84, 177–187 (2015).
Google Scholar 
43.Brown, P. T. & Caldeira, K. Greater future global warming inferred from Earth’s recent energy budget. Nature 552, 45–50 (2017).CAS 

Google Scholar 
44.Zscheischler, J. et al. Future climate risk from compound events. Nat. Clim. Change 8, 469–477 (2018).
Google Scholar 
45.Blois, J. L., Zarnetske, P. L., Fitzpatrick, M. C. & Finnegan, S. Climate change and the past, present, and future of biotic interactions. Science 341, 499–504 (2013).CAS 

Google Scholar 
46.Blackwood, J. C., Streicker, D. G., Altizer, S. & Rohani, P. Resolving the roles of immunity, pathogenesis, and immigration for rabies persistence in vampire bats. Proc. Natl Acad. Sci. USA 110, 20837–20842 (2013).CAS 

Google Scholar 
47.Fenner, A. L., Godfrey, S. S. & Michael Bull, C. Using social networks to deduce whether residents or dispersers spread parasites in a lizard population. J. Anim. Ecol. 80, 835–843 (2011).
Google Scholar 
48.Paniw, M. et al. The myriad of complex demographic responses of terrestrial mammals to climate change and gaps of knowledge: a global analysis. J. Anim. Ecol. 90, 1398–1407 (2021).
Google Scholar 
49.McDonald, J. L. et al. Demographic buffering and compensatory recruitment promotes the persistence of disease in a wildlife population. Ecol. Lett. 19, 443–449 (2016).
Google Scholar 
50.Plowright, R. K., Sokolow, S. H., Gorman, M. E., Daszak, P. & Foley, J. E. Causal inference in disease ecology: investigating ecological drivers of disease emergence. Front. Ecol. Environ. 6, 420–429 (2008).
Google Scholar 
51.Russell, R., DiRenzo, G. V., Szymanski, J., Alger, K. & Grant, E. H. C. Principles and mechanisms of wildlife population persistence in the face of disease. Front. Ecol. Evol. 8, 344 (2020).
Google Scholar 
52.Baudouin, A. et al. Disease avoidance, and breeding group age and size condition the dispersal patterns of western lowland gorilla females. Ecology 100, e02786 (2019).
Google Scholar 
53.Townsend, A. K., Hawley, D. M., Stephenson, J. F. & Williams, K. E. G. Emerging infectious disease and the challenges of social distancing in human and non-human animals. Proc. R. Soc. B 287, 20201039 (2020).CAS 

Google Scholar 
54.Schisler, G. J., Bergersen, E. P. & Walker, P. G. Effects of multiple stressors on morbidity and mortality of fingerling rainbow trout infected with Myxobolus cerebralis. Trans. Am. Fish. Soc. 129, 859–865 (2000).
Google Scholar 
55.Härkönen, T., Harding, K., Rasmussen, T. D., Teilmann, J. & Dietz, R. Age- and sex-specific mortality patterns in an emerging wildlife epidemic: the phocine distemper in European harbour seals. PLoS ONE 2, e887 (2007).
Google Scholar 
56.Clutton-Brock, T. H. et al. Reproduction and survival of suricates (Suricata suricatta) in the southern Kalahari. Afr. J. Ecol. 37, 69–80 (1999).
Google Scholar 
57.Clutton-Brock, T. H., Hodge, S. J. & Flower, T. P. Group size and the suppression of subordinate reproduction in Kalahari meerkats. Anim. Behav. 76, 689–700 (2008).
Google Scholar 
58.Bateman, A. W., Ozgul, A., Coulson, T. & Clutton-Brock, T. H. Density dependence in group dynamics of a highly social mongoose, Suricata suricatta. J. Anim. Ecol. 81, 628–639 (2012).
Google Scholar 
59.Adler, R. F. et al. The Global Precipitation Climatology Project (GPCP) monthly analysis (new version 2.3) and a review of 2017 global precipitation. Atmosphere 9, 138 (2018).
Google Scholar 
60.Moss, R. H. et al. The next generation of scenarios for climate change research and assessment. Nature 463, 747–756 (2010).CAS 

Google Scholar 
61.Parding, K. M. et al. GCMeval – an interactive tool for evaluation and selection of climate model ensembles. Clim. Serv. 18, 100167 (2020).
Google Scholar 
62.Delahay, R. J., Langton, S., Smith, G. C., Clifton-Hadley, R. S. & Cheeseman, C. L. The spatio-temporal distribution of Mycobacterium bovis (bovine tuberculosis) infection in a high-density badger population. J. Anim. Ecol. 69, 428–441 (2000).
Google Scholar 
63.Delahay, R. J. et al. Long-term temporal trends and estimated transmission rates for Mycobacterium bovis infection in an undisturbed high-density badger (Meles meles) population. Epidemiol. Infect. 141, 1445–1456 (2013).CAS 

Google Scholar 
64.Buzdugan, S. N., Chambers, M. A., Delahay, R. J. & Drewe, J. A. Diagnosis of tuberculosis in groups of badgers: an exploration of the impact of trapping efficiency, infection prevalence and the use of multiple tests. Epidemiol. Infect. 144, 1717–1727 (2016).CAS 

Google Scholar 
65.Akaike, H. in Selected Papers of Hirotugu Akaike (eds Parzen, E. et al.) 199–213 (Springer, 1998).66.Wood, S. N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J. R. Stat. Soc. B Stat. 73, 3–36 (2011).
Google Scholar 
67.Grimm, V. et al. The ODD protocol: a review and first update. Ecol. Model. 221, 2760–2768 (2010).
Google Scholar 
68.Wood, S. N. Statistical inference for noisy nonlinear ecological dynamic systems. Nature 466, 1102–1104 (2010).CAS 

Google Scholar 
69.Fronzek, S., Carter, T. R., Räisänen, J., Ruokolainen, L. & Luoto, M. Applying probabilistic projections of climate change with impact models: a case study for sub-Arctic palsa mires in Fennoscandia. Clim. Change 99, 515–534 (2010).
Google Scholar  More

1.Westerling, A. L. R. Increasing western US forest wildfire activity: sensitivity to changes in the timing of spring. Phil. Trans. R. Soc. B 371, 20150178 (2016).PubMed 
PubMed Central 

Google Scholar 
2.Abatzoglou, J. T. & Williams, A. P. Impact of anthropogenic climate change on wildfire across western US forests. Proc. Natl Acad. Sci. USA 113, 11770–11775 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
3.Gonzalez, P. et al. Southwest: Impacts, Risks, and Adaptation in the United States: Fourth National Climate Assessment (U.S. Global Change Research Program, 2018).4.McLauchlan, K. K. et al. Fire as a fundamental ecological process: research advances and frontiers. J. Ecol. 108, 2047–2069 (2020).
Google Scholar 
5.Bowman, D. M. J. S. et al. Fire in the Earth system. Science 324, 481–484 (2009).6.Davis, K. T. et al. Wildfires and climate change push low-elevation forests across a critical climate threshold for tree regeneration. Proc. Natl Acad. Sci. USA 116, 6193–6198 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
7.Stephens, S. L. et al. Drought, tree mortality, and wildfire in forests adapted to frequent fire. Bioscience 68, 77–88 (2018).
Google Scholar 
8.Radeloff, V. C. et al. Rapid growth of the US wildland–urban interface raises wildfire risk. Proc. Natl Acad. Sci. USA 115, 3314–3319 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
9.Syphard, A. D., Keeley, J. E., Pfaff, A. H. & Ferschweiler, K. Human presence diminishes the importance of climate in driving fire activity across the United States. Proc. Natl Acad. Sci. USA 114, 13750–13755 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
10.Mietkiewicz, N. et al. In the line of fire: consequences of human-ignited wildfires to homes in the U.S. (1992–2015). Fire 3, 50 (2020).11.Balch, J. K. et al. Human-started wildfires expand the fire niche across the United States. Proc. Natl Acad. Sci. USA 114, 2946–2951 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
12.McKenzie, D. & Littell, J. S. Climate change and the eco-hydrology of fire: will area burned increase in a warming western USA. Ecol. Appl. 27, 26–36 (2017).PubMed 

Google Scholar 
13.Littell, J. S., Mckenzie, D., Peterson, D. L. & Westerling, A. L. Climate and wildfire area burned in western U.S. ecoprovinces, 1916–2003. Ecol. Appl. 19, 1003–1021 (2009).PubMed 

Google Scholar 
14.Jensen, D. et al. The sensitivity of US wildfire occurrence to pre-season soil moisture conditions across ecosystems. Environ. Res. Lett. 13, 014021 (2018).PubMed 
PubMed Central 

Google Scholar 
15.Vicente-Serrano, S. M., Quiring, S. M., Peña-Gallardo, M., Yuan, S. & Domínguez-Castro, F. A review of environmental droughts: increased risk under global warming? Earth Sci. Rev. 201, 102953 (2020).
Google Scholar 
16.Ficklin, D. L. & Novick, K. A. Historic and projected changes in vapor pressure deficit suggest a continental-scale drying of the United States atmosphere. J. Geophys. Res. 122, 2061–2079 (2017).
Google Scholar 
17.Sarhadi, A., Ausín, M. C., Wiper, M. P., Touma, D. & Diffenbaugh, N. S. Multidimensional risk in a nonstationary climate: joint probability of increasingly severe warm and dry conditions. Sci. Adv. 4, eaau3487 (2018).PubMed 
PubMed Central 

Google Scholar 
18.Abatzoglou, J. T., Williams, A. P., Boschetti, L., Zubkova, M. & Kolden, C. A. Global patterns of interannual climate–fire relationships. Glob. Change Biol. 24, 5164–5175 (2018).
Google Scholar 
19.Williams, A. P. & Abatzoglou, J. T. Recent advances and remaining uncertainties in resolving past and future climate effects on global fire activity. Curr. Clim. Change Rep. 2, 1–14 (2016).
Google Scholar 
20.Bradstock, R. A. A biogeographic model of fire regimes in Australia: current and future implications. Glob. Ecol. Biogeogr. 19, 145–158 (2010).
Google Scholar 
21.Krawchuk, M. A. & Moritz, M. A. Constraints on global fire activity vary across a resource gradient. Ecology 92, 121–132 (2011).PubMed 

Google Scholar 
22.Scarff, F. R. et al. Effects of plant hydraulic traits on the flammability of live fine canopy fuels. Funct. Ecol. 35, 835–846 (2021).23.Ruffault, J., Martin-StPaul, N., Pimont, F. & Dupuy, J. L. How well do meteorological drought indices predict live fuel moisture content (LFMC)? An assessment for wildfire research and operations in Mediterranean ecosystems. Agric. For. Meteorol. 262, 391–401 (2018).
Google Scholar 
24.Pivovaro, A. L. et al. The effect of ecophysiological traits on live fuel moisture content. Fire 2, 28 (2019).25.Nolan, R. H., Hedo, J., Arteaga, C., Sugai, T. & Resco de Dios, V. Physiological drought responses improve predictions of live fuel moisture dynamics in a Mediterranean forest. Agric. For. Meteorol. 263, 417–427 (2018).26.Skelton, R. P., West, A. G. & Dawson, T. E. Predicting plant vulnerability to drought in biodiverse regions using functional traits. Proc. Natl Acad. Sci. USA 112, 5744–5749 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
27.Ma, W. et al. Assessing climate change impacts on live fuel moisture and wildfire risk using a hydrodynamic vegetation model. Biogeosciences 18, 4005–4020 (2021).CAS 

Google Scholar 
28.McColl, K. A. et al. The global distribution and dynamics of surface soil moisture. Nat. Geosci. 10, 100–104 (2017).CAS 

Google Scholar 
29.Chuvieco, E., González, I., Verdú, F., Aguado, I. & Yebra, M. Prediction of fire occurrence from live fuel moisture content measurements in a Mediterranean ecosystem. Int. J. Wildland Fire 18, 430–441 (2009).30.Rao, K., Williams, A. P., Flefil, J. F. & Konings, A. G. SAR-enhanced mapping of live fuel moisture content. Remote Sens. Environ. 245, 111797 (2020).
Google Scholar 
31.Nolan, R. H., Boer, M. M., Resco De Dios, V., Caccamo, G. & Bradstock, R. A. Large-scale, dynamic transformations in fuel moisture drive wildfire activity across southeastern Australia. Geophys. Res. Lett. 43, 4229–4238 (2016).
Google Scholar 
32.Dennison, P. E. & Moritz, M. A. Critical live fuel moisture in chaparral ecosystems: a threshold for fire activity and its relationship to antecedent precipitation. Int. J. Wildland Fire 18, 1021–1027 (2009).
Google Scholar 
33.Tumino, B. J., Duff, T. J., Goodger, J. Q. D. & Cawson, J. G. Plant traits linked to field-scale flammability metrics in prescribed burns in Eucalyptus forest. PLoS ONE 14, e0221403 (2019).34.Rodman, K. C. et al. A trait-based approach to assessing resistance and resilience to wildfire in two iconic North American conifers. J. Ecol. 109, 313–326 (2021).
Google Scholar 
35.Resco de Dios, V. Plant–Fire Interactions (Springer, 2020).36.Hurteau, M. D., Liang, S., Westerling, A. L. R. & Wiedinmyer, C. Vegetation–fire feedback reduces projected area burned under climate change. Sci. Rep. 9, 2838 (2019).37.Littell, J. S., McKenzie, D., Wan, H. Y. & Cushman, S. A. Climate change and future wildfire in the western United States: an ecological approach to nonstationarity. Earths Future 6, 1097–1111 (2018).38.Abatzoglou, J. T. & Kolden, C. A. Relationships between climate and macroscale area burned in the western United States. Int. J. Wildland Fire 22, 1003–1020 (2013).
Google Scholar 
39.Goss, M. et al. Climate change is increasing the likelihood of extreme autumn wildfire conditions across California. Environ. Res. Lett. 15, 094016 (2020).40.Bradshaw, L. S., Deeming, J. E., Burgan, R. E. & Cohen, J. D. The 1978 National Fire-Danger Rating System: Technical Documentation General Technical Report INT-169 (US Department of Agriculture, Forest Service, Intermountain Forest and Range Experiment Station,1984); https://doi.org/10.2737/INT-GTR-16941.Hardy, C. C. & Hardy, C. E. Fire danger rating in the United States of America: an evolution since 1916. Int. J. Wildland Fire 16, 217–231 (2007).42.Rabin, S. S. et al. The Fire Modeling Intercomparison Project (FireMIP), phase 1: experimental and analytical protocols with detailed model descriptions. Geosci. Model Dev. 10, 1175–1197 (2017).CAS 

Google Scholar 
43.Hantson, S. et al. The status and challenge of global fire modelling. Biogeosciences 13, 3359–3375 (2016).
Google Scholar 
44.Anderegg, W. R. L. Spatial and temporal variation in plant hydraulic traits and their relevance for climate change impacts on vegetation. New Phytol. 205, 1008–1014 (2015).PubMed 

Google Scholar 
45.Konings, A. G. & Gentine, P. Global variations in ecosystem-scale isohydricity. Glob. Change Biol. 23, 891–905 (2017).
Google Scholar 
46.Forkel, M. et al. Emergent relationships with respect to burned area in global satellite observations and fire-enabled vegetation models. Biogeosciences 16, 57–76 (2019).
Google Scholar 
47.Brodribb, T. J., Powers, J., Cochard, H. & Choat, B. Hanging by a thread? Forests and drought. Science 368, 261–266 (2020).CAS 
PubMed 

Google Scholar 
48.Trugman, A. T., Anderegg, L. D. L., Shaw, J. D. & Anderegg, W. R. L. Trait velocities reveal that mortality has driven widespread coordinated shifts in forest hydraulic trait composition. Proc. Natl Acad. Sci. USA 117, 8532–8538 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
49.Williams, A. P. et al. Correlations between components of the water balance and burned area reveal new insights for predicting forest fire area in the southwest United States. Int. J. Wildland Fire 24, 14–26 (2015).
Google Scholar 
50.Knapp, P. A. Spatio-temporal patterns of large grassland fires in the Intermountain West U.S.A. Glob. Ecol. Biogeogr. Lett. 7, 259–272 (1998).
Google Scholar 
51.Keeley, J. & Syphard, A. Climate change and future fire regimes: examples from California. Geosciences 6, 37 (2016).
Google Scholar 
52.Badia, A., Serra, P. & Modugno, S. Identifying dynamics of fire ignition probabilities in two representative Mediterranean wildland–urban interface areas. Appl. Geogr. 31, 930–940 (2011).
Google Scholar 
53.Fusco, E. J., Abatzoglou, J. T., Balch, J. K., Finn, J. T. & Bradley, B. A. Quantifying the human influence on fire ignition across the western USA. Ecol. Appl. 26, 2390–2401 (2016).
Google Scholar 
54.Syphard, A. D. et al. Human influence on California fire regimes. Ecol. Appl. 17, 1388–1402 (2007).PubMed 

Google Scholar 
55.Ager, A. A., Finney, M. A., Kerns, B. K. & Maffei, H. Modeling wildfire risk to northern spotted owl (Strix occidentalis caurina) habitat in central Oregon, USA. For. Ecol. Manage. 246, 45–56 (2007).56.Thomas, D., Butry, D., Gilbert, S., Webb, D. & Fung, J. The Costs and Losses of Wildfires: A Literature Survey NIST Special Publication 1215 (NIST, 2017); https://doi.org/10.6028/NIST.SP.121557.Wang, D. et al. Economic footprint of California wildfires in 2018. Nat. Sustain. 4, 252–260 (2021).
Google Scholar 
58.Burke, M. et al. The changing risk and burden of wildfire in the United States. Proc. Natl Acad. Sci. USA 118, e2011048118 (2021).59.García, M., Chuvieco, E., Nieto, H. & Aguado, I. Combining AVHRR and meteorological data for estimating live fuel moisture content. Remote Sens. Environ. 112, 3618–3627 (2008).
Google Scholar 
60.Matthews, S. Dead fuel moisture research: 1991–2012. Int. J. Wildland Fire 23, 78–92 (2014).
Google Scholar 
61.Cohen, J. D. et al. The National Fire-Danger Rating System: Basic Equations Vol. 82 (US Department of Agriculture, Forest Service, Pacific Southwest Forest and Range Experiment Station, 1985).62.Pellizzaro, G., Cesaraccio, C., Duce, P., Ventura, A. & Zara, P. Relationships between seasonal patterns of live fuel moisture and meteorological drought indices for Mediterranean shrubland species. Int. J. Wildland Fire 16, 232–241 (2007).63.Liu, L., Zhang, Y., Wu, S., Li, S. & Qin, D. Water memory effects and their impacts on global vegetation productivity and resilience. Sci. Rep. 8, 2962 (2018).PubMed 
PubMed Central 

Google Scholar 
64.Anderegg, W. R. L. et al. Woody plants optimise stomatal behaviour relative to hydraulic risk. Ecol. Lett. 21, 968–977 (2018).PubMed 

Google Scholar 
65.Meinzer, F. C., Johnson, D. M., Lachenbruch, B., McCulloh, K. A. & Woodruff, D. R. Xylem hydraulic safety margins in woody plants: coordination of stomatal control of xylem tension with hydraulic capacitance. Funct. Ecol. 23, 922–930 (2009).
Google Scholar 
66.National Fuel Moisture Database (United States Forest Service, 2018); https://www.wfas.net/nfmd/public/index.php67.Abatzoglou, J. T. Development of gridded surface meteorological data for ecological applications and modelling. Int. J. Climatol. 33, 121–131 (2011).
Google Scholar 
68.Homer, C. et al. Completion of the 2006 National Land Cover Database for the conterminous United States. Photogramm. Eng. Remote Sens. 77, 858–864 (2011).69.Williams, A. P. et al. Observed impacts of anthropogenic climate change on wildfire in California. Earths Future 7, 892–910 (2019).
Google Scholar 
70.Boschetti, L., Roy, D., Hoffman, A. A. & Humber, M. Collection 5 MODIS Burned Area Product User Guide Version 3.0.1 (NASA EOSDIS Land Processes DAAC, 2013).71.PRISM Climate Data (Prism Climate Group, Oregon State University, accessed 16 December 2020); https://prism.oregonstate.edu72.Simard, M., Pinto, N., Fisher, J. B. & Baccini, A. Mapping forest canopy height globally with spaceborne lidar. J. Geophys. Res. 116, G04021 (2011).73.Fan, Y., Miguez-Macho, G., Jobbágy, E. G., Jackson, R. B. & Otero-Casal, C. Hydrologic regulation of plant rooting depth. Proc. Natl Acad. Sci. USA 114, 10572–10577 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
74.Montzka, C., Herbst, M., Weihermüller, L., Verhoef, A. & Vereecken, H. A global data set of soil hydraulic properties and sub-grid variability of soil water retention and hydraulic conductivity curves. Earth Syst. Sci. Data 9, 529–543 (2017).
Google Scholar 
75.Liu, S. et al. NACP MsTMIP: Unified North American Soil Map (ORNL DAAC, 2014); https://doi.org/10.3334/ornldaac/124276.Pedregosa, F. et al. Scikit-learn: machine learning in Python. J. Mach. Learn. Res. 12, 2825–2830 (2011).
Google Scholar 
77.Martinuzzi, S. et al. The 2010 Wildland–Urban Interface of the Conterminous United States (USDA, 2015).78.Medlyn, B. E. et al. Reconciling the optimal and empirical approaches to modelling stomatal conductance. Glob. Change Biol. 17, 2134–2144 (2011).
Google Scholar  More

1.Le Loeuff, J., Buffetaut, E. & Martin, M. The last stages of dinosaur faunal history in Europe: a succession of Maastrichtian dinosaur assemblages from the Corbières (southern France). Geol. Mag. 131, 625–630 (1994).
Google Scholar 
2.Vila, B., Sellés, A. G. & Brusatte, S. L. Diversity and faunal changes in the latest Cretaceous dinosaur communities of southwestern Europe. Cretac. Res. 57, 552–564 (2016).
Google Scholar 
3.Fondevilla, V. et al. Chronostratigraphic synthesis of the latest Cretaceous dinosaur turnover in south-western Europe. Earth Sci. Rev. 191, 168–189 (2019).
Google Scholar 
4.Sanz, J. L., Powell, J. E., Le Loeuff, J., Martínez, R. & Pereda-Suberbiola, X. Sauropod remains from the Upper Cretaceous of Laño (northcentral Spain). Titanosaur phylogenetic relationships. Est. Mus. Cienc. Nat. Alava 14, 235–255 (1999).
Google Scholar 
5.Garcia, G., Amico, S., Fournier, F., Thouand, E. & Valentin, X. A new titanosaur genus (Dinosauria, Sauropoda) from the Late Cretaceous of southern France and its paleobiogeographic implications. Bull. Soc. Géol. Fr. 181, 269–277 (2010).
Google Scholar 
6.Díez Díaz, V. et al. A new titanosaur (Dinosauria: Sauropoda) from the Upper Cretaceous of Velaux La-Bastide Neuve (southern France). Hist. Biol. https://doi.org/10.1080/08912963.2020.1841184 (2020).7.Le Loeuff, J. Ampelosaurus atacis (nov. gen., nov. sp.), a new titanosaurid (Dinosauria, Sauropoda) from the Late Cretaceous of the Upper Aude Valley (France). C. R. Acad. Sci. II 321, 693–700 (1995).
Google Scholar 
8.Díez Díaz, V. et al. A new titanosaur (Dinosauria, Sauropoda) from the Upper Cretaceous of Lo Hueco (Cuenca, Spain). Cretac. Res. 68, 49–60 (2016).
Google Scholar 
9.Company, J. Bone histology of the titanosaur Lirainosaurus astibiae (Dinosauria: Sauropoda) from the latest Cretaceous of Spain. Naturwissenschaften 98, 67–78 (2011).CAS 
PubMed 

Google Scholar 
10.Klein, N. et al. Modified laminar bone in Ampelosaurus atacis and other titanosaurs (Sauropoda): implications for life history and physiology. PLoS ONE 7, e36907 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
11.Díez Díaz, V. et al. The titanosaurian dinosaur Atsinganosaurus velauciensis (Sauropoda) from the Upper Cretaceous of southern France: new material, phylogenetic affinities, and palaeobiogeographical implications. Cretac. Res. 91, 429–456 (2018).
Google Scholar 
12.Benítez-López, A. et al. The island rule explains consistent patterns of body size evolution in terrestrial vertebrates. Nat. Ecol. Evol. 5, 768–786 (2021).PubMed 

Google Scholar 
13.Benton, M. J. et al. Dinosaurs and the island rule: the dwarfed dinosaurs from Haţeg Island. Palaeogeogr. Palaeoclimatol. Palaeoecol. 293, 438–454 (2010).
Google Scholar 
14.Canudo, J. I. Descripción de un fragmento proximal de fémur de Titanosauridae (Dinosauria, Sauropoda) del Maastrichtiense superior de Serraduy (Huesca). In Proc. XVII Jornadas de la Sociedad Española de Paleontología (eds Meléndez, G. et al.) 255–262 (Sociedad Española de Paleontología y Área y Museo de Paleontología de la Universidad de Zaragoza, 2001).15.Vila, B. et al. The diversity of sauropods and their first taxonomic succession from the latest Cretaceous of south-western Europe: clues to demise and extinction. Palaeogeogr. Palaeoclimatol. Palaeoecol. 350–352, 19–38 (2012).
Google Scholar 
16.Sallam, H. M. et al. New Egyptian sauropod reveals Late Cretaceous dinosaur dispersal between Europe and Africa. Nat. Ecol. Evol. 2, 445–451 (2018).PubMed 

Google Scholar 
17.Buffetaut, E. Archosaurian reptiles with Gondwanan affinities in the Upper Cretaceous of Europe. Terra Nova 1, 69–74 (1989).
Google Scholar 
18.Le Loeuff, J. The Campano-Maastrichtian vertebrate faunas from southern Europe and their relationships with other faunas in the world; palaeobiogeographical implications. Cretac. Res. 12, 93–114 (1991).
Google Scholar 
19.Pereda-Suberbiola, X. Biogeographical affinities of Late Cretaceous continental tetrapods of Europe: a review. Bull. Soc. Geol. Fr. 180, 57–71 (2009).
Google Scholar 
20.Csiki-Sava, Z., Buffetaut, E., Ősi, A., Pereda-Suberbiola, X. & Brusatte, S. L. Island life in the Cretaceous – faunal composition, biogeography, evolution, and extinction of land-living vertebrates on the Late Cretaceous European archipelago. ZooKeys 469, 1–161 (2015).
Google Scholar 
21.Ezcurra, M. D. & Agnolín, F. L. A new global palaeobiogeographical model for the late Mesozoic and early Tertiary. Syst. Biol. 61, 553–566 (2012).PubMed 

Google Scholar 
22.Sellés, A. G. & Vila, B. Re-evaluation of the age of some dinosaur localities from the southern Pyrenees by means of megaloolithid oospecies. J. Iber. Geol. 41, 125–139 (2015).
Google Scholar 
23.Bonaparte, J. F. & Coria, R. A. Un nuevo y gigantesco saurópodo titanosaurio de la Formación Río Limay (Albiano–Cenomaniano) de la Provincia del Neuquén, Argentina. Ameghiniana 30, 217–282 (1993).
Google Scholar 
24.Curry Rogers, K. The postcranial osteology of Rapetosaurus krausei (Sauropoda: Titanosauria) from the Late Cretaceous of Madagascar. J. Vertebr. Paleontol. 29, 1046–1086 (2009).
Google Scholar 
25.Zurriaguz, V. & Powell, J. New contributions to the presacral osteology of Saltasaurus loricatus (Sauropoda, Titanosauria) from the Upper Cretaceous of northern Argentina. Cretac. Res. 54, 283–300 (2015).
Google Scholar 
26.Coria, R. A., Filippi, L. S., Chiappe, L. M., García, R. & Arcucci, A. B. Overosaurus paradasorum gen. et sp. nov., a new sauropod dinosaur (Titanosauria: Lithostrotia) from the Late Cretaceous of Neuquén, Patagonia, Argentina. Zootaxa 3683, 357–376 (2013).PubMed 

Google Scholar 
27.Calvo, J. O., González Riga, B. J. & Porfiri, J. D. A new titanosaur sauropod from the Late Cretaceous of Neuquén, Patagonia, Argentina. Arq. Mus. Nac. 65, 485–504 (2007).
Google Scholar 
28.Jain, S. L. & Bandyopadhyay, S. New titanosaurid (Dinosauria: Sauropoda) from the Late Cretaceous of central India. J. Vertebr. Paleontol. 17, 114–136 (1997).
Google Scholar 
29.Gorscak, E. & O’Connor, P. M. A new African titanosaurian sauropod dinosaur from the middle Cretaceous Galula Formation (Mtuka Member), Rukwa Rift Basin, southwestern Tanzania. PLoS ONE 14, e0211412 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
30.Kellner, A. W. A. & de Azevedo, S. A. K. A new sauropod dinosaur (Titanosauria) from the Late Cretaceous of Brazil. Nat. Sci. Mus. Monogr. 15, 111–142 (1999).
Google Scholar 
31.Novas, F. E. et al. Paleontological discoveries in the Chorrillo Formation (upper Campanian-lower Maastrichtian, Upper Cretaceous), Santa Cruz Province, Patagonia, Argentina. Rev. Mus. Argent. Cienc. Nat. 21, 217–293 (2019).
Google Scholar 
32.Wilson, J. A., D’Emic, M. D., Curry Rogers, K. A., Mohabey, D. M. & Sen, S. Reassessment of the sauropod dinosaur Jainosaurus (=“Antarctosaurus”) septentrionalis from the Upper Cretaceous of India. Contrib. Mus. Paleontol. Univ. Mich. 32, 17–40 (2009).
Google Scholar 
33.Powell, J. E. Revision of South American titanosaurid dinosaurs: palaeobiological, palaeobiogeographical and phylogenetic aspects. Rec. Queen Vic. Mus. 111, 1–173 (2003).
Google Scholar 
34.Smith, J. B. et al. A giant sauropod dinosaur from an Upper Cretaceous mangrove deposit in Egypt. Science 292, 1704–1706 (2001).CAS 
PubMed 

Google Scholar 
35.Otero, A. & Vizcaíno, S. F. Hindlimb musculature and function of Neuquensaurus australis (Sauropoda: Titanosauria). Ameghiniana 45, 333–348 (2008).
Google Scholar 
36.von Huene, F. Los saurisquios y ornitisquios del Cretáceo Argentino. Mus. La Plata 3, 1–196 (1929).
Google Scholar 
37.Mannion, P. D. & Otero, A. A reappraisal of the Late Cretaceous Argentinean sauropod dinosaur Argyrosaurus superbus, with a description of a new titanosaur genus. J. Vertebr. Paleontol. 32, 614–638 (2012).
Google Scholar 
38.Mocho, P., Pérez-García, A., Martín Jiménez, M. & Ortega, F. New remains from the Spanish Cenomanian shed light on the Gondwanan origin of European Early Cretaceous titanosaurs. Cretac. Res. 95, 164–190 (2019).
Google Scholar 
39.Díez Díaz, V., Pereda Suberbiola, X. & Sanz, J. L. Appendicular skeleton and dermal armour of the Late Cretaceous titanosaur Lirainosaurus astibiae (Dinosauria: Sauropoda) from Spain. Palaeontol. Electronica 16, 19A (2013).
Google Scholar 
40.Le Loeuff, J. in Thunder-Lizards: The Sauropodomorph Dinosaurs (eds Tidwell, V. & Carpenter, K.) 115–137 (Indiana Univ. Press, 2005).41.Borsuk-Bialynicka, M. A new camarasaurid sauropod Opisthocoelicaudia skarzynskii gen. n., sp. n. from the Upper Cretaceous of Mongolia. Acta Palaeontol. Pol. 37, 5–63 (1977).
Google Scholar 
42.Filippi, L. S., García, R. A. & Garrido, A. A new sauropod titanosaur from the Plottier Formation (Upper Cretaceous) of Patagonia (Argentina). Geol. Acta 9, 1–12 (2011).
Google Scholar 
43.Salgado, L., Coria, R. A. & Calvo, J. O. Evolution of titanosaurid sauropods. I: phylogenetic analysis based on the postcranial evidence. Ameghiniana 34, 3–32 (1997).
Google Scholar 
44.D’Emic, M. D. The early evolution of titanosauriform sauropod dinosaurs. Zool. J. Linn. Soc. 166, 624–671 (2012).
Google Scholar 
45.Tschopp, E. & Mateus, O. Clavicles, interclavicles, gastralia, and sternal ribs in sauropod dinosaurs: new reports from Diplodocidae and their morphological, functional and evolutionary implications. J. Anat. 222, 321–340 (2013).PubMed 

Google Scholar 
46.Wilson, J. A. Sauropod dinosaur phylogeny: critique and cladistic analysis. Zool. J. Linn. Soc. 136, 217–276 (2002).
Google Scholar 
47.Powell, J. E. in Los Dinosaurios y Su Entorno Biótico (eds Sanz, J. L. & Buscalioni, A. D.) 165–230 (Instituto ‘Juan de Valdés’, 1992).48.Otero, A. The appendicular skeleton of Neuquensaurus, a Late Cretaceous saltasaurine sauropod from Patagonia, Argentina. Acta Palaeontol. Pol. 55, 399–426 (2010).
Google Scholar 
49.Gilmore, C. W. Reptilian Fauna of the North Horn Formation of Central Utah Professional Paper 210-C (USGS Numbered Series, 1946).50.Ullmann, P. V. & Lacovara, K. J. Appendicular osteology of Dreadnoughtus schrani, a giant titanosaurian (Sauropoda, Titanosauria) from the Upper Cretaceous of Patagonia, Argentina. J. Vertebr. Paleontol. 36, e1225303 (2016).
Google Scholar 
51.Poropat, S. F. Carl Wiman’s sauropods: the Uppsala Museum of Evolution’s collection. GFF 135, 104–119 (2013).CAS 

Google Scholar 
52.Cerda, I. A., Salgado, L. & Powell, J. E. Extreme postcranial pneumaticity in sauropod dinosaurs from South America. Paläontol. Z. 86, 441–449 (2012).
Google Scholar 
53.Wilson, J. A. & Carrano, M. T. Titanosaurs and the origin of “wide-gauge” trackways: a biomechanical and systematic perspective on sauropod locomotion. Paleobiol. 25, 252–267 (1999).
Google Scholar 
54.Upchurch, P., Barrett, P. & Dodson, P. in The Dinosauria (eds Weishampel, D. B. et al.) 259–324 (Univ. California Press, 2004).55.Lehman, T. M. & Coulson, A. B. A juvenile specimen of the sauropod dinosaur Alamosaurus sanjuanensis from the Upper Cretaceous of Big Bend National Park, Texas. J. Paleontol. 76, 156–172 (2002).
Google Scholar 
56.Gallina, P. A. & Otero, A. Reassessment of Laplatasaurus araukanicus (Sauropoda: Titanosauria) from the Upper Cretaceous of Patagonia, Argentina. Ameghiniana 52, 487–501 (2015).
Google Scholar 
57.Wilson, J. A. & Upchurch, P. Redescription and reassessment of the phylogenetic affinities of Euhelopus zdanskyi (Dinosauria: Sauropoda) from the Early Cretaceous of China. J. Syst. Palaeontol. 7, 199–239 (2009).
Google Scholar 
58.Sereno, P. C. A rationale for phylogenetic definitions, with application to the higher-level taxonomy of Dinosauria. Neues Jahrb. Geol. Paläontol. Abh. 210, 41–83 (1998).
Google Scholar 
59.Chiappe, L. M. et al. Sauropod dinosaur embryos from the Late Cretaceous of Patagonia. Nature 396, 258–261 (1998).CAS 

Google Scholar 
60.Haq, B. U. Cretaceous eustasy revisited. Glob. Planet. Change 113, 44–58 (2014).
Google Scholar 
61.Gheerbrant, E. & Rage, J.-C. Paleobiogeography of Africa: how distinct from Gondwana and Laurasia? Palaeogeogr. Palaeoclimatol. Palaeoecol. 241, 224–246 (2006).
Google Scholar 
62.Canudo, J. I. et al. What Iberian dinosaurs reveal about the bridge said to exist between Gondwana and Laurasia in the Early Cretaceous. Bull. Soc. Géol. Fr. 180, 5–11 (2009).
Google Scholar 
63.Dal Sasso, C., Pierangelini, G., Famiani, F., Cau, A. & Nicosia, U. First sauropod bones from Italy offer new insights on the radiation of Titanosauria between Africa and Europe. Cretac. Res. 64, 88–109 (2016).
Google Scholar 
64.Stein, K. et al. Small body size and extreme cortical bone remodeling indicate phyletic dwarfism in Magyarosaurus dacus (Sauropoda: Titanosauria). Proc. Natl Acad. Sci. USA 107, 9258–9263 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
65.Botfalvai, G. et al. ‘X’ marks the spot! Sedimentological, geochemical and palaeontological investigations of Upper Cretaceous (Maastrichtian) vertebrate fossil localities from the Vălioara valley (Densuş-Ciula Formation, Hațeg Basin, Romania). Cretac. Res. 123, 104781 (2021).
Google Scholar 
66.Csiki-Sava, Z. et al. The east side story–the Transylvanian latest Cretaceous continental vertebrate record and its implications for understanding Cretaceous–Paleogene boundary events. Cretac. Res. 57, 662–698 (2016).
Google Scholar 
67.Campione, N. E. & Evans, D. C. A universal scaling relationship between body mass and proximal limb bone dimensions in quadrupedal terrestrial tetrapods. BMC Biol. 10, 60 (2012).PubMed 
PubMed Central 

Google Scholar 
68.González Riga, B. J., Lamanna, M. C., Ortiz David, L. D., Calvo, J. O. & Coria, J. P. A gigantic new dinosaur from Argentina and the evolution of the sauropod hind foot. Sci. Rep. 6, 19165 (2016).PubMed 
PubMed Central 

Google Scholar 
69.Seebacher, F. New method to calculate allometric length–mass relationships of dinosaurs. J. Vertebr. Paleontol. 21, 51–60 (2001).
Google Scholar 
70.Mallison, H. & Wings, O. Photogrammetry in paleontology— a practical guide. J. Paleontol. Tech. 12, 1–31 (2014).
Google Scholar 
71.Matthews, N., Noble, T. & Breithaupt, B. H. in Dinosaur Tracks—The Next Steps (eds Falkingham, P. L. et al.) 28–55 (Indiana Univ. Press, 2016).72.Falkingham, P. L. et al. A standard protocol for documenting modern and fossil ichnological data. Palaeontology 61, 469–480 (2018).
Google Scholar 
73.Huelsenbeck, J. P. & Ronquist, F. MRBAYES: Bayesian inference of phylogeny. Bioinformatics 17, 754–755 (2001).CAS 

Google Scholar 
74.Stadler, T., Kühnert, D., Bonhoeffer, S. & Drummond, A. J. Birth–death skyline plot reveals temporal changes of epidemic spread in HIV and hepatitis C virus (HCV). Proc. Natl Acad. Sci. USA 110, 228–233 (2013).CAS 
PubMed 

Google Scholar 
75.Matzke, N. J. Probabilistic historical biogeography: new models for founder-event speciation, imperfect detection, and fossils allow improved accuracy and model-testing. Front. Biogeogr. 5, 242–248 (2013).
Google Scholar 
76.Ogg, J. G. & Hinnov, L. A. in The Geological Time Scale (eds Gradstein, F. M. et al.) 793–853 (Elsevier, 2012).77.Vianey-Liaud, M., Khosla, A. & Garcia, G. Relationships between European and Indian dinosaur eggshells of the oofamily Megaloolithidae. J. Vertebr. Paleontol. 23, 575–585 (2003).
Google Scholar  More

1.Dellwig O, Schnetger B, Brumsack H-J, Grossart H-P, Umlauf L. Dissolved reactive manganese at pelagic redoxclines (part II): hydrodynamic conditions for accumulation. J Mar Syst. 2012;90:31–41.
Google Scholar 
2.Taylor GT, Iabichella M, Ho T, Scranton MI, Thunell RC, Muller-Karger F, et al. Chemoautotrophy in the redox transition zone of the Cariaco Basin: a significant midwater source of organic carbon production. Limnol Oceanogr. 2001;46:148–63.CAS 

Google Scholar 
3.Zopfi J, Ferdelman TG, Jørgensen BB, Teske A, Thamdrup B. Influence of water column dynamics on sulfide oxidation and other major biogeochemical processes in the chemocline of Mariager Fjord (Denmark). Mar Chem. 2001;74:29–51.CAS 

Google Scholar 
4.Trefry JH, Presley BJ, Keeney-Kennicutt WL, Trocine RP. Distribution and chemistry of manganese, iron, and suspended particulates in Orca Basin. Geo-Mar Lett. 1984;4:125–30.
Google Scholar 
5.Dahl TW, Anbar AD, Gordon GW, Rosing MT, Frei R, Canfield DE. The behavior of molybdenum and its isotopes across the chemocline and in the sediments of sulfidic Lake Cadagno, Switzerland. Geochim Cosmochim Acta. 2010;74:144–63.CAS 

Google Scholar 
6.Özsoy E, Ünlüata Ü. Oceanography of the Black Sea: a review of some recent results. Earth-Sci Rev. 1997;42:231–72.
Google Scholar 
7.Wegwerth A, Eckert S, Dellwig O, Schnetger B, Severmann S, Weyer S, et al. Redox evolution during Eemian and Holocene sapropel formation in the Black Sea. Palaeogeogr Palaeoclimatol Palaeoecol. 2018;489:249–60.
Google Scholar 
8.Murray JW, Jannasch HW, Honjo S, Anderson RF, Reeburgh WS, Top Z, et al. Unexpected changes in the oxic/anoxic interface in the Black Sea. Nature. 1989;338:411–3.CAS 

Google Scholar 
9.Schulz-Vogt HN, Pollehne F, Jürgens K, Arz HW, Bahlo R, Dellwig O, et al. Effect of large magnetotactic bacteria with polyphosphate inclusions on the phosphate profile of the suboxic zone in the Black Sea. ISME J. 2019;13:1198–208.CAS 
PubMed 
PubMed Central 

Google Scholar 
10.Dellwig O, Wegwerth A, Schnetger B, Schulz H, Arz HW. Dissimilar behaviors of the geochemical twins W and Mo in hypoxic-euxinic marine basins. Earth-Sci Rev. 2019;193:1–23.CAS 

Google Scholar 
11.Stanev EV, Poulain PM, Grayek S, Johnson KS, Claustre H, Murray JW. Understanding the dynamics of the oxic-anoxic interface in the Black Sea. Geophys Res Lett. 2018;45:864–71.CAS 

Google Scholar 
12.Trouwborst RE. Soluble Mn(III) in suboxic zones. Science. 2006;313:1955–7.CAS 
PubMed 

Google Scholar 
13.Vliet DM, Meijenfeldt FAB, Dutilh BE, Villanueva L, Sinninghe Damsté JS, Stams AJM, et al. The bacterial sulfur cycle in expanding dysoxic and euxinic marine waters. Environ Microbiol. 2021;23:2834–57.PubMed 

Google Scholar 
14.Konovalov SK, Luther GW, Friederich GE, Nuzzio DB, Tebo BM, Murray JW, et al. Lateral injection of oxygen with the Bosporus plume-fingers of oxidizing potential in the Black Sea. Limnol Oceanogr. 2003;48:2369–76.CAS 

Google Scholar 
15.Lewis BL, Landing WM. The biogeochemistry of manganese and iron in the Black Sea. Deep Sea Res A Oceanogr Res Pap. 1991;38:S773–S803.
Google Scholar 
16.Yakushev EV, Pollehne F, Jost G, Kuznetsov I, Schneider B, Umlauf L. Analysis of the water column oxic/anoxic interface in the Black and Baltic seas with a numerical model. Mar Chem. 2007;107:388–410.CAS 

Google Scholar 
17.Gregg MC, Yakushev E. Surface ventilation of the Black Sea’s cold intermediate layer in the middle of the western gyre. Geophys Res Lett. 2005;32:1–4.
Google Scholar 
18.Schnetger B, Dellwig O. Dissolved reactive manganese at pelagic redoxclines (part I): a method for determination based on field experiments. J Mar Syst. 2012;90:23–30.
Google Scholar 
19.Tebo BM, Bargar JR, Clement BG, Dick GJ, Murray KJ, Parker D, et al. Biogenic manganese oxides: Properties and mechanisms of formation. Annu Rev Earth Planet Sci. 2004;32:287–328.CAS 

Google Scholar 
20.Glockzin M, Pollehne F, Dellwig O. Stationary sinking velocity of authigenic manganese oxides at pelagic redoxclines. Mar Chem. 2014;160:67–74.CAS 

Google Scholar 
21.Dellwig O, Leipe T, März C, Glockzin M, Pollehne F, Schnetger B, et al. A new particulate Mn-Fe-P-shuttle at the redoxcline of anoxic basins. Geochim Cosmochim Acta. 2010;74:7100–15.CAS 

Google Scholar 
22.Burdige DJ, Nealson KH. Chemical and microbiological studies of sulfide-mediated manganese reduction. Geomicrobiol J. 1986;4:361–87.CAS 

Google Scholar 
23.Yao W, Millero FJ. The rate of sulfide oxidation by δMnO2 in seawater. Geochim Cosmochim Acta. 1993;57:3359–65.CAS 

Google Scholar 
24.Parks DH, Chuvochina M, Waite DW, Rinke C, Skarshewski A, Chaumeil PA, et al. A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat Biotechnol. 2018;36:996.CAS 
PubMed 

Google Scholar 
25.Henkel JV, Dellwig O, Pollehne F, Herlemann DPR, Leipe T, Schulz-Vogt HN. A bacterial isolate from the Black Sea oxidizes sulfide with manganese(IV) oxide. Proc Natl Acad Sci USA. 2019;116:12153–5.CAS 
PubMed 
PubMed Central 

Google Scholar 
26.Henkel JV, Vogts A, Werner J, Neu TR, Spröer C, Bunk B, et al. Candidatus Sulfurimonas marisnigri sp. nov. and Candidatus Sulfurimonas baltica sp. nov., thiotrophic manganese oxide reducing chemolithoautotrophs of the class Campylobacteria isolated from the pelagic redoxclines of the Black Sea and the Baltic Sea. Syst Appl Microbiol. 2021;44:1–11.27.Grote J, Jost G, Labrenz M, Herndl GJ, Jürgens K. Epsilonproteobacteria represent the major portion of chemoautotrophic bacteria in sulfidic waters of pelagic redoxclines of the Baltic and Black Seas. Appl Environ Microbiol. 2008;74:7546–51.CAS 
PubMed 
PubMed Central 

Google Scholar 
28.Pernthaler A, Pernthaler J, Amann R. Fluorescence in situ hybridization and catalyzed reporter deposition for the identification of marine bacteria. Appl Environ Microbiol. 2002;68:3094–101.CAS 
PubMed 
PubMed Central 

Google Scholar 
29.Sekar R, Pernthaler A, Pernthaler J, Warnecke F, Posch T, Amann R. An improved protocol for quantification of freshwater Actinobacteria by fluorescence in situ hybridization. Appl Environ Microbiol. 2003;69:2928–35.CAS 
PubMed 
PubMed Central 

Google Scholar 
30.Grote J, Labrenz M, Pfeiffer B, Jost G, Jürgens K. Quantitative distributions of Epsilonproteobacteria and a Sulfurimonas subgroup in pelagic redoxclines of the central Baltic Sea. Appl Environ Microbiol. 2007;73:7155–61.CAS 
PubMed 
PubMed Central 

Google Scholar 
31.Daims H, Bruhl 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. 1999;22:434–44.CAS 
PubMed 

Google Scholar 
32.Wallner G, Amann R, Beisker W. Optimizing fluorescent in situ hybridization with rRNA-targeted oligonucleotide probes for flow cytometric identification of microorganisms. Cytometry. 1993;11:136–43.
Google Scholar 
33.Glöckner FO, Yilmaz P, Quast C, Gerken J, Beccati A, Ciuprina A, et al. 25 years of serving the community with ribosomal RNA gene reference databases and tools. J Biotechnol. 2017;261:169–76.PubMed 

Google Scholar 
34.Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:590–6.
Google Scholar 
35.Konstantinidis KT, Tiedje JM. Towards a genome-based taxonomy for prokaryotes. J Bacteriol. 2005;187:6258–64.CAS 
PubMed 
PubMed Central 

Google Scholar 
36.Buchfink B, Reuter K, Drost H-G. Sensitive protein alignments at tree-of-life scale using DIAMOND. Nat Methods. 2021;18:366–8.CAS 
PubMed 
PubMed Central 

Google Scholar 
37.Von Meijenfeldt FAB, Arkhipova K, Cambuy DD, Coutinho FH, Dutilh BE. Robust taxonomic classification of uncharted microbial sequences and bins with CAT and BAT. Genome Biol. 2019;20:1–14.
Google Scholar 
38.Schulz HD. Conceptual models and computer models. In: Schulz HD, Zabel M, editors. Marine geochemistry. Springer: Berlin, Heidelberg; 2006. p. 513–47.39.Diepenbroek M, Glöckner FO, Grobe P, Güntsch A, Huber R, König-Ries B, et al. Towards an integrated biodiversity and ecological research data management and archiving platform: the German federation for the curation of biological data (GFBio). In: Plödereder E, Grunske L, Schneider E, Ull D, editors. Informatik 2014. Bonn: Gesellschaft für Informatik e.V.; 2014.p. 1711–21.40.Yilmaz P, Kottmann R, Field D, Knight R, Cole JR, Amaral-Zettler L, et al. Minimum information about a marker gene sequence (MIMARKS) and minimum information about any (x) sequence (MIxS) specifications. Nat Biotechnol. 2011;29:415–20.CAS 
PubMed 
PubMed Central 

Google Scholar 
41.Revsbech NP, Thamdrup B, Dalsgaard T, Canfield DE. Construction of STOX oxygen sensors and their application for determination of O2 concentrations in oxygen minimum zones. Methods Enzymol. 2011;486:325–41.CAS 
PubMed 

Google Scholar 
42.Dahl C. A biochemical view on the biological sulfur cycle. In: Environmental technologies to treat sulphur pollution: principles and engineering. IWA Publishing: London; 2020;2:55–96.43.Murray JW, Yakushev EV. Past and present water column anoxia. Past and present water column anoxia. Dordrecht: Springer Netherlands; 2006.44.Schulz HD. Quantification of early diagenesis: dissolved constituents in pore water and signals in the solid phase. In: Schulz HD, Zabel M, editors. Marine geochemistry. Berlin/Heidelberg: Springer-Verlag; 2006. p. 73–124.45.Tebo BM. Manganese(II) oxidation in the suboxic zone of the Black Sea. Deep Res A. 1991;38:883–905.
Google Scholar 
46.Konovalov S, Samodurov A, Oguz T, Ivanov L. Parameterization of iron and manganese cycling in the Black Sea suboxic and anoxic environment. Deep Res Part I Oceanogr Res Pap. 2004;51:2027–45.CAS 

Google Scholar 
47.Lahme S, Callbeck CM, Eland LE, Wipat A, Enning D, Head IM, et al. Comparison of sulfide-oxidizing Sulfurimonas strains reveals a new mode of thiosulfate formation in subsurface environments. Environ Microbiol. 2020;22:1784–1800.CAS 
PubMed 

Google Scholar 
48.Grote J, Schott T, Bruckner CG, Glockner FO, Jost G, Teeling H, et al. Genome and physiology of a model Epsilonproteobacterium responsible for sulfide detoxification in marine oxygen depletion zones. Proc Natl Acad Sci USA. 2012;109:506–10.CAS 
PubMed 

Google Scholar 
49.Sievert SM, Scott KM, Klotz MG, Chain PSG, Hauser LJ, Hemp J, et al. Genome of the Epsilonproteobacterial chemolithoautotroph Sulfurimonas denitrificans. Appl Environ Microbiol. 2008;74:1145–56.CAS 
PubMed 

Google Scholar 
50.Friedrich CG, Bardischewsky F, Rother D, Quentmeier A, Fischer J. Prokaryotic sulfur oxidation. Curr Opin Microbiol. 2005;8:253–9.CAS 
PubMed 

Google Scholar 
51.Götz F, Pjevac P, Markert S, McNichol J, Becher D, Schweder T, et al. Transcriptomic and proteomic insight into the mechanism of cyclooctasulfur- versus thiosulfate-oxidation by the chemolithoautotroph Sulfurimonas denitrificans. Environ Microbiol. 2019;21:244–58.PubMed 

Google Scholar 
52.Pjevac P, Meier DV, Markert S, Hentschker C, Schweder T, Becher D, et al. Metaproteogenomic profiling of microbial communities colonizing actively venting hydrothermal chimneys. Front Microbiol. 2018;9:1–12.
Google Scholar 
53.Meier DV, Pjevac P, Bach W, Hourdez S, Girguis PR, Vidoudez C, et al. Niche partitioning of diverse sulfur-oxidizing bacteria at hydrothermal vents. ISME J. 2017;11:1545–58.CAS 
PubMed 
PubMed Central 

Google Scholar 
54.Wang S, Jiang L, Hu Q, Liu X, Yang S, Shao Z. Elemental sulfur reduction by a deep‐sea hydrothermal vent Campylobacterium Sulfurimonas sp. NW10. Environ Microbiol. 2021;23:965–79.CAS 
PubMed 

Google Scholar 
55.Yao W, Millero FH. Oxidation of hydrogen sulfide by Mn(IV) and Fe(III) (hydr)oxides in seawater. Mar Chem. 1996;52:1–16.CAS 

Google Scholar 
56.Herszage J, dos Santos Afonso M. Mechanism of hydrogen sulfide oxidation by manganese(IV) oxide in aqueous solutions. Langmuir. 2003;19:9684–92.CAS 

Google Scholar 
57.Glazer BT, Luther GW, Konovalov SK, Friederich GE, Nuzzio DB, Trouwborst RE, et al. Documenting the suboxic zone of the Black Sea via high-resolution real-time redox profiling. Deep Res II Top Stud Oceanogr. 2006;53:1740–55.
Google Scholar 
58.Jørgensen BB, Fossing H, Wirsen CO, Jannasch HW. Sulfide oxidation in the anoxic Black Sea chemocline. Deep Sea Res A Oceanogr Res Pap. 1991;38:1083–103.
Google Scholar 
59.Yiǧiterhan O, Murray JW. Trace metal composition of particulate matter of the Danube River and Turkish rivers draining into the Black Sea. Mar Chem. 2008;111:63–76.
Google Scholar 
60.Brewer PG, Spencer DW. Distribution of some trace elements in Black Sea and their flux between dissolved and particulate phases: water. In: The Black Sea–Geology, Chemistry, and Biology. AAPG Special Volumes. AAPG; 1974;137–43.61.Fuchsman CA, Kirkpatrick JB, Brazelton WJ, Murray JW, Staley JT. Metabolic strategies of free-living and aggregate-associated bacterial communities inferred from biologic and chemical profiles in the Black Sea suboxic zone. FEMS Microbiol Ecol. 2011;78:586–603.CAS 
PubMed 

Google Scholar 
62.Kelly DP. Biochemistry of the chemolithotrophic oxidation of inorganic sulphur. Philos Trans R Soc Lond B Biol Sci. 1982;298:499–528.CAS 
PubMed 

Google Scholar 
63.Kirkpatrick JB, Fuchsman CA, Yakushev EV, Egorov AV, Staley JT, Murray JW. Dark N2 fixation: nifH expression in the redoxcline of the Black Sea. Aquat Micro Ecol. 2018;82:43–58.
Google Scholar 
64.Glaubitz S, Kießlich K, Meeske C, Labrenz M, Jürgens K. SUP05 Dominates the gammaproteobacterial sulfur oxidizer assemblages in pelagic redoxclines of the central baltic and black seas. Appl Environ Microbiol. 2013;79:2767–76.CAS 
PubMed 
PubMed Central 

Google Scholar 
65.Shah V, Chang BX, Morris RM. Cultivation of a chemoautotroph from the SUP05 clade of marine bacteria that produces nitrite and consumes ammonium. ISME J. 2017;11:263–71.CAS 
PubMed 

Google Scholar 
66.Rogge A, Vogts A, Voss M, Jürgens K, Jost G, Labrenz M. Success of chemolithoautotrophic SUP05 and Sulfurimonas GD17 cells in pelagic Baltic Sea redox zones is facilitated by their lifestyles as K- and r -strategists. Environ Microbiol. 2017;19:2495–506.CAS 
PubMed 

Google Scholar 
67.Overmann J, Cypionka H, Pfennig N. An extremely low-light-adapted phototrophic sulfur bacterium from the Black Sea. Limnol Oceanogr. 1992;37:150–5.CAS 

Google Scholar 
68.Jensen MM, Kuypers MMM, Lavik G, Thamdrup B. Rates and regulation of anaerobic ammonium oxidation and denitrification in the Black Sea. Limnol Oceanogr. 2008;53:23–36.CAS 

Google Scholar 
69.Hannig M, Lavik G, Kuypers MMM, Woebken D, Martens-Habbena W, Jürgens K. Shift from denitrification to anammox after inflow events in the central Baltic Sea. Limnol Oceanogr. 2007;52:1336–45.CAS 

Google Scholar 
70.Engström P, Dalsgaard T, Hulth S, Aller RC. Anaerobic ammonium oxidation by nitrite (anammox): Implications for N2 production in coastal marine sediments. Geochim Cosmochim Acta. 2005;69:2057–65.
Google Scholar 
71.Dapena-Mora A, Fernández I, Campos JL, Mosquera-Corral A, Méndez R, Jetten MSM. Evaluation of activity and inhibition effects on Anammox process by batch tests based on the nitrogen gas production. Enzym Micro Technol. 2007;40:859–65.CAS 

Google Scholar 
72.Havig JR, McCormick ML, Hamilton TL, Kump LR. The behavior of biologically important trace elements across the oxic/euxinic transition of meromictic Fayetteville Green Lake, New York, USA. Geochim Cosmochim Acta. 2015;165:389–406.CAS 

Google Scholar 
73.Jürgens K, Taylor GT. Microbial ecology and biogeochemistry of oxygen-deficient water columns. Microbial Ecology of the Ocean, 3rd ed. Hoboken: Wiley; 2018. p. 231–88.74.Jost G, Martens-Habbena W, Pollehne F, Schnetger B, Labrenz M. Anaerobic sulfur oxidation in the absence of nitrate dominates microbial chemoautotrophy beneath the pelagic chemocline of the eastern Gotland Basin, Baltic Sea. FEMS Microbiol Ecol. 2010;71:226–36.CAS 
PubMed 

Google Scholar 
75.Aller RC, Rude PD. Complete oxidation of solid phase sulfides by manganese and bacteria in anoxic marine sediments. Geochim Cosmochim Acta. 1988;52:751–65.CAS 

Google Scholar 
76.King GM. Effects of added manganic and ferric oxides on sulfate reduction and sulfide oxidation in intertidal sediments. FEMS Microbiol Ecol. 1990;73:131–8.CAS 

Google Scholar  More

A theoretical framework of interacting bacterial strainsOur mathematical model was motivated by experimental assays used to establish colony biofilms where the founding inoculum is placed on the surface of solidified nutrient agar. Within the inoculum footprint, individual (or small clusters of) bacteria settle at random locations and grow over time into a mature structured macroscale community (Fig. 1A). In the mathematical model, all the founding cells are assumed to have identical properties. However, to track the dynamics of biofilm growth we divided the founding cells into two groups, denoted by ({B}_{1}) (shown in magenta) and ({B}_{2}) (shown in green) (Fig. 1B). Note that we refer to ({B}_{1}) and ({B}_{2}) as strains for brevity, even though they represent two isogenic cell lineages that express different fluorescent proteins in a single-strain biofilm (Fig. 1A). In our theoretical framework, biofilm dynamics were reduced to the fundamental processes of local growth and spatial spread (more details below), which provided a species-independent representation of dual-strain biofilm growth. Suitably nondimensionalised (see Section S3), the model is given by$$frac{partial {B}_{1}}{partial t}=nabla cdot left({Id}left(1-left({B}_{1}+{B}_{2}right)right){nabla B}_{1}right)+{B}_{1}left(1-left({B}_{1}+{B}_{2}right)right),$$$$frac{partial {B}_{2}}{partial t}=nabla cdot left({Id}left(1-left({B}_{1}+{B}_{2}right)right)nabla {B}_{2}right)+{B}_{2}left(1-({B}_{1}+{B}_{2})right),$$
(1)
where, the variables ({0le B}_{1}left({{{{{boldsymbol{x}}}}}},tright),{B}_{2}left({{{{{boldsymbol{x}}}}}},tright)le 1) denote the scaled densities of each strain, respectively at time (t, > ,0) (one nondimensional time unit corresponding to approx. 2.9 h) and at spatial position ({{{{{boldsymbol{x}}}}}}in Omega) (one nondimensional space unit corresponding to approx. 0.15 mm). The spatial domain (Omega ={{{{{{boldsymbol{x}}}}}}in {{mathbb{R}}}^{2}:{||}{{{{{boldsymbol{x}}}}}}{||}le R}) is a two-dimensional disk, representing the biofilm growth medium (Fig. 1C). This simplification provided a significant reduction in computational cost and was motivated by an analysis of a previously published data set, in which we found a two-order of magnitude difference between biofilm diameter and biofilm thickness in B. subtilis NCIB 3610 [27]. The model is therefore unable to explicitly resolve density distributions along the vertical axis, for example, layering of subpopulation caused by gradients in environmental conditions [28,29,30] or topographical features such as ‘wrinkles’ [31]. However, it is fully capable of capturing overlap between subpopulations that are below the environmental carrying capacity and thus can track spatio-temporal coexistence. Moreover, as we show below, we find strong agreement between data obtained from two-dimensional in silico biofilms and data gathered from laboratory grown biofilms, which further supports the model simplification.Fig. 1: Experimental and modelling set-up.A An example of the experimental assay. Founder cells carry either a constitutively produced copy of GFP (green) or mTagBFP (magenta). The bacteria were mixed in a 1:1 ratio and images taken after 24 h and 72 h of incubation. The number of founder cells was approx. 10 CFUs. The scalebars are 5 mm long. B An example realisation of the mathematical model. In the right-hand plots green and magenta are used to differentiate two subsets of the initial patches ((t=0), top) and their subsequent development ((t=25), bottom). Black areas indicate the computational domain, (varOmega). The plot of initial condition is a blow-up of the centre of the whole domain. The scalebars represent 7 nondimensional space units. C Schematic of model initial condition. Initial populations (filled coloured circles) are placed in ({varOmega }_{0}), a small subdomain of the whole computational domain (varOmega) (both centred at the origin (O)).Full size imageThe initial conditions of the theoretical framework were motivated by the random positions at which bacteria settle on the agar within the inoculum footprint (Fig. 1A). In our theoretical framework, we represented the experimental inoculum footprint by a small disk ({Omega }_{0}=left{{{{{{boldsymbol{x}}}}}}in Omega :{||}{{{{{boldsymbol{x}}}}}}{||} ; < ; {R}_{0}right}) in the centre of the computational domain (Fig. 1C). We modelled the random deposition of bacteria by randomly placing ‘microcolonies’ within ({Omega }_{0}) at nodes of a triangulated spatial mesh of linear geometric order, used in the application of a finite element method to numerically solve the model equations (Fig. 1B, C). Each initial microcolony was assumed to only contain one strain and to be at carrying capacity (i.e., ({B}_{1}=1) or ({B}_{2}=1) within each microcolony). Unless otherwise stated, we used an even number ((N)) of initial microcolonies and assigned exactly (N/2) to each strain at random. At spatial locations other than the assigned microcolonies, both densities were set to zero.The size of a spatial mesh element used in the model (approx. (0.008{m}{m}^{2}) in experimental parameters) was much larger than that of a single bacterial cell. This means that the initial conditions represented the experimental assays shortly after inoculation (typically after 24 h of incubation), at which time each bacterium (or small cluster of bacteria) had formed a distinct, spatially separated microcolony. Hence, the number of in silico microcolonies, (N,) represented the number of bacteria used in the initial inoculum. Resolving the initial data at this spatial scale allowed analysis for founder densities (0le Nle 824). Using a selected set of values from that range was sufficient to capture clear trends (see below). The range covers biologically relevant founder densities, which generate mature colony biofilms with broadly similar morphologies (Supplementary Fig. S1). Additionally, to verify whether the observed trends could be extrapolated to (N ; > ; 824), we represented high founder densities by piecewise spatially homogeneous initial conditions ({B}_{1}={B}_{2}=0.5) in ({Omega }_{0}) and ({B}_{1}={B}_{2}=0) otherwise.The strains were assumed to grow logistically, with growth being limited by the total population, which could not exceed unity (after nondimensionalisation). Moreover, spatial propagation was described by diffusion as is common [32]. However, in our model, we employed a diffusion coefficient that decreased with increasing population size. This density dependence prevented merging of initially separated founding patches in the model and was invoked to capture experimental observations that indicated such colonies abut rather than merge on meeting [33, 34]. The indicator function ({Id}=1) if ({B}_{1}+{B}_{2}le 1) and ({Id}=0) otherwise guaranteed nonnegativity of the diffusion coefficients; this constrained the model to the physically relevant case and moreover ensured numerical stability during simulation.Finally, we defined the competitive outcome score (for ({B}_{1})) of the interaction to be the relative mass of strain ({B}_{1}) i.e., ({B}_{1}^{Omega }/({B}_{1}^{Omega }+{B}_{2}^{Omega })) at the chosen end point ((t=T)) of our model simulation, where$${B}_{i}^{Omega }:={int }_{Omega }{B}_{i}({{{{{boldsymbol{x}}}}}},T){{{{{rm{d}}}}}}{{{{{boldsymbol{x}}}}}},,i=1,2.$$The competitive outcome score lies in the interval (left[{{{{mathrm{0,1}}}}}right]) with the value 0.5 signifying a 1:1 ratio between the strains. Note that we could swap the indices without loss of generality to equivalently define the competitive outcome to be the relative mass of strain(,{B}_{2}) at the chosen end point.Low founder densities yield large variability in competitive outcomesIn the absence of spatial dynamics, the mathematical model predicted that the ratio between both strains would always remain constant (left(frac{d}{{dt}}big(frac{{B}_{1}}{{B}_{2}}big)=0right)) and therefore that the competitive outcome would be determined by the initial ratio. To test whether such a relationship continued to hold in the full, spatially extended system, we examined data from simulations over a test range of initial founding cell densities. The initial strain ratio was selected to be 1:1 for each test.Model simulations using homogeneous initial conditions (representing high founder densities) consistently resulted in a competitive outcome score of 0.5 (i.e., strains in 1:1 ratio) with the strains remaining homogeneously distributed in space across the colony (Fig. 2A, Supplementary Movie S1). By contrast, independent model realisations using a specified number of microcolonies placed at randomly chosen locations representing low (({N}=6)) and intermediate (({N}=824)) founder densities, revealed significant variation in competitive outcome (Fig. 2B, C, Supplementary Movies S2 and S3). To explore this observed variability in more detail, we employed a Monte Carlo approach. For each fixed founder density (N) within the selected set, 1000 independent model realisations were conducted. Data from these simulations revealed that the competitive outcome score for each founder density was normally distributed with mean 0.5. The standard deviation was relatively large for low founder densities ((N={{{{mathrm{4,6,8,10}}}}})) and decreased with further increases in (N) (Fig. 2D). (Note the small standard deviation for (N=2); see supplementary information for a discussion of this special case). Finally, our model predicted significant changes in the spatial organisation of the two strains within the biofilm in response to changing founder density, consistent with previous studies [14]. For high founder densities, isogenic in silico strains were predicted to coexist homogenously (Fig. 2A). However, as the founder density was decreased (decreasing (N)), homogeneous coexistence was gradually replaced by the formation of spatial sectors dominated by one strain or the other. Full segregation occurred for low founder densities (Fig. 2B, C).Fig. 2: Spatial structure and variability in competitive outcome depend on founder density.A–C Example model realisations for different founder densities. All plots show the system’s initial conditions ((t=0)) and the outcomes after 25 time units. Plots visualising the systems’ states at (t=0) show a blow-up of the subdomain ({varOmega }_{0}); plots visualising outcomes at (t=25) show the full computational domain (varOmega) (black background). The scalebars are seven unit lengths long. A The outcome of simulations initialised with piecewise spatially homogeneous populations representing high founder density. The ‘Merged’ image channel shows both strains (grey colour corresponds to overlap); the ({B}_{1})(green) and ({B}_{2}) (magenta) channels only show single strain filters of the plot. B The range of outcomes observed for low founder density (number of initial cell patches ({N}=6)). C The range of outcomes for intermediate founder densities ((N=824)). In (B, C) only the ‘Merged’ channel is shown. D Variability in competitive outcome increases with decreasing founder density. Each boxplot contains data from 1000 model realisations. Blue and red boxplots correspond to the founder densities in B and C, respectively.Full size imageAccess to free space determines competitive outcomeNext, we attempted to uncover the mechanism(s) by which low founder densities drive variability in competitive outcome. Motivated by [14], we first tested whether the initial separation between initial microcolonies of different types was the simple determinant. We did not find this to be the case for isogenic strain pairings in the mathematical model (Supplementary Fig. S2).As an alternative, we hypothesised that a microcolony surrounded by others may have little impact on competitive outcome as its contribution to biofilm growth would be ultimately limited. On the other hand, microcolonies located close to the boundary of the biofilm inoculum would be free to expand radially and thus could make a more significant contribution to the competitive outcome (for an example timelapse video see Movie S3). Hence, we explored whether competitive outcome was correlated to a strain’s potential for radial expansion beyond the inoculum. To do so, we assumed the potential for radial expansion to be solely determined by the geographical locations of a strain’s initial microcolonies. We then defined an appropriate score for this potential as follows. First, a circle was drawn that enclosed the initial microcolonies. Second, each point on the circle was associated with the nearest microcolony and assigned to that strain. Third, the total arc length on the circle associated with each strain was computed. Finally, the access to free space score (AFS score) for strain ({B}_{1}), denoted AFS1, was then computed as the ratio of the total arc length associated with ({B}_{1}) to the circumference of the circle. Therefore, (0le {{{{{rm{AF}}}}}}{{{{{{rm{S}}}}}}}_{1}le 1) quantified strain ({B}_{1})’s hypothesised potential to contribute to radial biofilm expansion. It is straightforward to confirm that the AFS score for strain ({B}_{2}), ({{{{{rm{AF}}}}}}{{{{{{rm{S}}}}}}}_{2}=1-{{{{{rm{AF}}}}}}{{{{{{rm{S}}}}}}}_{1}). See Section S4.2 and Supplementary Figs. S3 and S4 for a mathematically rigorous definition of the AFS score.We explored the utility of the AFS score using (N=6) and (N=824) as representatives of low and intermediate founder cell densities, respectively. We increased the number of model realisations to 5000 for each of the selected values of N to ensure improved accuracy of our data analysis. The AFS score was then calculated for each of the 10,000 initial conditions (see examples Fig. 3A, B). On completion of each simulation, the corresponding competitive outcome score was computed. Analysis of these model data confirmed that the AFS score accurately predicts competitive outcome: for each fixed founder density, the AFS score unfolds the variation shown in Fig. 2D, yielding a positive, linear relationship between AFS1 and competitive outcome for ({B}_{1}) (Fig. 3C, D). For each of the selected values of (N), initial configurations of microcolonies with a low AFS1 score predictably generated a low competitive outcome for ({B}_{1}). Correspondingly, initial configurations with a high AFS1 score predictably generated a high competitive outcome for ({B}_{1}). The slope of this linear relationship provided a deterministic quantification of the variability of competitive outcomes for a given founder density (cf. Fig. 3C, D, Supplemental text).Fig. 3: Access to free space determines competitive outcome.A, B Example model realisations for different founder densities. All plots show system initial conditions ((t=0)) with the reference circle used to compute the AFS score (the circle is rescaled for visualisation purposes) and outcomes after 25 time units. The founder densities are (N=824) and (N=6) in A and B, respectively. Plots visualising system states at (t=0) show a blow-up of the subdomain ({varOmega }_{0}); plots visualising outcomes at (t=25) show the full computational domain (varOmega) (black background). The scalebars are seven unit lengths long. C, D The relation between the AFS score ({AF}{S}_{1}), and competitive outcome is shown for intermediate founder density ((N=824)) and low founder density ((N=6)) in C and D, respectively. Data were obtained from 5000 model realisations and cover the continuum of ({AF}{S}_{1}). The observed probability density function for AFS is shown (circular markers); along with the density function of a fitted normal distribution ((mu approx 0.5,sigma approx 0.10) in C, (mu approx 0.5,sigma approx 0.16) in D) (solid line). E The relation between the standard deviations of the AFS score ({AF}{S}_{1}) and the competitive outcome. Each data point (circle) represents a different founder density and contains information from 1000 model realisations.Full size imageWe subsequently established that the predictive power of the AFS score was maintained across the range of founder densities considered in the model. Additionally, the variation in the AFS score was shown to decrease with increasing founder density (cf. Fig. 3C, D). Further, we revealed strong correlation between variation in AFS score and variation in competitive outcome (Fig. 3E). Therefore, for increasing founder density, the observed decrease in variation in competitive outcome can be directly attributed to the decrease in variation in the AFS score.Dual strain single-isolate biofilm assays confirm modelling hypothesesNext, we aimed to test the hypotheses put forward by the mathematical model. We selected an isogenic pair of Bacillus subtilis strains derived from isolate NCIB 3610 that constitutively produced the green fluorescent protein GFP (NRS6942, shown in green, Table S1) and the blue fluorescent protein mTagBFP (NRS6932, shown in magenta, Tables S1 and S2), respectively. In line with the modelling assumption, the isolates were mixed in a 1:1 ratio at a defined initial cell density (we used an OD600 of 1) and this cell culture was serially diluted prior to inoculating the colony biofilms (Section S7). Thus, biofilms were inoculated using ~106 CFUs and dilutions in 10-fold increments to order 1 CFU. For each founder density, 12 technical replicates were performed to provide a meaningful sample size, and the experiment was repeated on three independent occasions. We used a non-destructive colony biofilm image analysis approach, to measure the relative mass (and hence the competitive outcome) of the two isogenic strains at 24 h, 48 h, 72 h after inoculation (see Section S10). We confirmed that the output from the image analysis correlated well with data generated by disruption of the colony biofilm and analysis of the relative strain proportions determined using single cells analysis by flow cytometry (Fig. 4A) (see also [35]). The mTagBFP labelled strain consistently performed marginally worse than the GFP labelled competitor at high founder densities in co-culture, which suggests some impact on competitive fitness (Fig. 4B, C). To allow comparison with results from the mathematical model, we denoted the mTagBFP (NRS6932, shown in magenta) and GFP (NRS6942, shown in green) strains as ({B}_{1}) and ({B}_{2}), respectively, with associate AFS scores AFS1 and AFS2Fig. 4: Experimental data confirm modelling hypotheses.A Comparison of image analysis with flow cytometry. A scatter plot comparing measurements of relative density of the mTagBFP-labelled strain obtained from image analysis and flow cytometry is shown. Each data point corresponds to one biofilm, which was imaged before being analysed by flow cytometry. The data contains measurements taken from all strain pairs, all founder densities, and all time points. The solid blue line shows the identity (x=y), with the coefficient of determination being ({R}^{2}=0.91). B Example images of single-strain biofilms consisting of GFP (green(,{B}_{1})) and mTagBFP (magenta, ({B}_{2})) labelled copies of 3610. Taken after 72 h of incubation and shown for two different founder densities (scalebar 5 mm). C Strain density data. Competitive outcome measurements taken after 24 h, 48 h and 72 h of biofilm incubation. Plotted are technical repeats from one biological repeat of the experiment. The full data set is presented in Fig. S5A. D Example visualisations of AFS score calculations. Three example biofilms images at 24 h (left), 48 h (middle) and 72 h (right). The strains are as described in B. Images at 24 h show the reference circle used for the AFS1 score. E The relationship between AFS1 and competitive outcome for ({B}_{1}). AFS was calculated from images taken at 24 h, and competitive outcome for ({B}_{1}) after 48 h (left, (n=30)) and 72 h (right, (n=25)). The linear correlation coefficient (rho) is indicated.Full size imageOur experimental analysis proved consistent with the model predictions. High founder densities resulted in a broadly homogenous distribution of both strains over the footprint of the biofilm, while low founder densities led to a high degree of spatial segregation of the strains within the mature biofilm (Fig. 4B, see also [14]). Additionally, analysis of experimental data confirmed that variability in competitive outcome increased with decreasing founder density (Fig. 4B, C, Supplementary Fig. S5A). For founder densities equivalent to (sim)103 to (sim)106 CFUs, the competitive outcome was consistent across each set of technical replicates. By contrast, for founder densities equivalent to (sim)1 to (sim)102 CFUs, the competitive outcome was variable across each set of technical replicates. We noted that variability in competitive outcome, at all initial founder densities, was marginally amplified over time.We assumed the process of repeated dilution and selection of the inoculum volume may not guarantee an exact cell count and/or initial strain ratio of 1:1 at lower founder densities. Indeed, for low founder densities after 24 hrs incubation, we observed inconsistencies in the number and ratio of CFUs deposited (Supplementary Fig. S5B). We therefore considered whether these inconsistencies in the biofilm inocula contributed to the observed variability in competitive outcome. To explore this in more detail, we first implemented a combinatorial ‘cell picking’ model that mathematically simulated the process of selecting the small inoculum volume from a larger cell culture (see Section S4.3). This process identified a threshold of ({sim} {10}^{2}) CFUs below which variability in cell number and/or strain ratio could measurably deviate from their intended values in our experimental assay. Above this threshold, the combinatorial argument predicted limited deviation from the intended values (Supplementary Fig. S6A). Coupling these theoretical predictions with our experimental observations (Supplementary Fig. S5B), we concluded that any observed variability in competitive outcome cannot be a consequence of a measurable deviation in the inoculum composition for colony biofilms founded with (sim {10}^{2}) CFUs or higher.We next wanted to determine whether the predictive power of the AFS score could be used to connect experimental initial configurations of the bacteria with the observed competitive outcome. To do this accurately, we required that the founding bacteria remained spatially separated as small colonies until an image was taken at 24 h (the earliest imaging time-point, see Fig. 4D). Therefore, we only used founder densities lower than 102 CFUs. However, the above noted inconsistencies in initial strain ratios and cell counts at these densities raised the question of whether AFS could still accurately predict competitive outcome. To test this, we repeated our Monte Carlo simulations of (1) in which the number of initial microcolonies for each strain was drawn using the combinatorial cell picking model, rather than being a fixed number and in a 1:1 ratio. Analysing the resulting simulation data for model (1) confirmed that the predictive power of the AFS score was robust to any ‘naturally-occurring’ variation in the initial strain ratio (Supplementary Fig. S6B). Correspondingly, our analysis of the experimental data revealed a strong correlation between a strain’s AFS score and the competitive outcome measured at 48 h and 72 h after incubation (Fig. 4E).A modelling framework for non-isogenic strainsWe have established that for isogenic strains, the initial configuration of founding bacteria determines the competitive outcome in a ‘race for space’ and that the AFS score can accurately predict which strain will dominate. A natural question that follows is what would happen if this race for space was influenced by antagonistic interactions such as killing or growth inhibition. Therefore, we considered the effect of introducing a local (e.g., contact-dependent or short-range non-contact dependent) antagonistic mechanism that causes a reduction in strain net growth. In an extension of our theoretical framework (1), constants describing the ratios between the strains’ maximum growth rates in the absence of competition ((r)), diffusion coefficients ((d)) and competition coefficients ((c)) were introduced to allow for the possibility of differences in strain properties. This resulted in the following system obtained after a suitable nondimensionalisation (see Section S3):$$frac{partial {B}_{1}}{partial t}=nabla cdot left({Id}left(1-frac{{B}_{1}+{B}_{2}}{k}right){nabla B}_{1}right)+{B}_{1}left(1-frac{{B}_{1}+{B}_{2}}{k}right)-{B}_{1}{B}_{2},$$$$frac{partial {B}_{2}}{partial t}=nabla cdot left({Id}cdot dleft(1-frac{{B}_{1}+{B}_{2}}{k}right)nabla {B}_{2}right)+{{rB}}_{2}left(1-frac{{B}_{1}+{B}_{2}}{k}right)-c{B}_{1}{B}_{2}.$$
(2)
Here, the indicator function ({Id}=1) if ({B}_{1}+{B}_{2}le k) and ({Id}=0) otherwise, where k is the nondimensional carrying capacity. To start, strains were assumed to possess identical growth dynamics in the absence of competitors (i.e., r (=1,{d}=1)), but to significantly differ in their ability to negatively impact the competitor strain. For the simulations we set (c=0.2) representing a five-fold difference in competition strength, with ({B}_{2}) being the more effective competitor. A linear stability analysis of model [4] confirmed that in this case and for a homogeneous initial distribution of the strains in a 1:1 ratio, ({B}_{2}) wins the interaction. For this reason, we therefore refer to ({B}_{2}) as the (intrinsically) stronger strain and to ({B}_{1}) as the (intrinsically) weaker strain in the following.The assumption of identical growth dynamics allowed us to focus on the impact of antagonistic interactions on competitive outcome. We anticipated that this assumption was unlikely to hold for non-isogenic strains in experimental settings and therefore we examined (as will be discussed later) the impact of changes to the parameters (r,{d}) and (c). Subsequently, we showed the effect of such parameter variation to be limited.Spatial segregation induced by low founder densities enables coexistenceIn the context of local antagonistic interactions, low founder densities were expected to offer protection for the weaker strain by driving spatial segregation and the formation of enclaves. Test simulations supported this hypothesis. Model realisations with high (spatially uniform initial conditions) and intermediate ((N=824)) founder densities consistently led to competitive exclusion of the weaker strain (Fig. 5A, B, Supplementary Movies S4 and S5), while model realisations with low founder densities ((N=6)) resulted in coexistence with the strains being spatially segregated (Fig. 5C). Once established during early stages of the model simulation, spatial segregation was conserved. However, the stronger strain continually invaded its competitor’s clusters along strain-to-strain interfaces and eventually took over the biofilm centre. Simultaneously, the weaker strain enlarged its sectors due to unimpeded growth on the biofilm edge. Coexistence, as measured by competitive outcome was achieved by a balance of these processes (Supplementary Movie S6).Fig. 5: Modelling data for a non-isogenic strain pair with local antagonistic interactions.A–C Example model realisations for high (A), intermediate (B) and low (C) founder density are shown. A the Merged image channel shows both strains (grey colour corresponds to overlap), the ({B}_{1}) and ({B}_{2}) channels only show single strain filters of the plot. In B, C only the Merged channel is shown. Plots visualising system states at (t=0) show a blow-up of the subdomain ({varOmega }_{0}) and the circles used to calculate the AFS scores around the initial conditions are not to scale. Plots visualising outcomes at (t=25) show the full computational domain (varOmega) (black background). The scalebars are seven unit lengths long. D The relation between founder density and competitive outcome. Each boxplot contains data from 1000 model realisations. E The relation between the AFS score ({AF}{S}_{1}), and competitive outcome for one fixed founder density ((N=6)). Data were obtained from 5000 model realisations and covers the continuum of ({AF}{S}_{1}). The observed probability density function for AFS is shown (circular markers); the density function of a fitted normal distribution ((mu approx 0.5,sigma approx 0.16)) as a solid line.Full size imageLow founder densities generated significant variation in competitive outcome (Fig. 5C). In particular, outcomes were observed for which the weaker strain ({B}_{1}) coexisted with, and could even outperform, the stronger strain ({B}_{2}). To better understand the impact of founder density, we performed Monte Carlo simulations with 1000 independent model realisations for each founder density (N) in our test range. Data from these simulations revealed both the mean and variation of competitive outcome for the weaker strain increased with decreasing founder density (Fig. 5D).Access to free space determines competitive outcome for low founder densitiesThe mathematical model consistently predicted competitive exclusion of the weaker strain at intermediate and high founder densities (Fig. 5A, B). Hence, in these cases, the AFS score no longer provided a meaningful predictor of competitive outcome. Rather, the model predicted the outcome to be dominated by the local antagonisms. However, as detailed above, low founder densities ((N) = 6) resulted in a highly variable competitive outcome and therefore we explored whether the AFS score remained an accurate predictor in this case. The simulation data confirmed that for this fixed number (N), the AFS score remained capable of accurately unfolding the observed variation in competitive outcome (Fig. 5E). Thus, initial strain configurations with a low AFS1 predictably generated a low competitive outcome for ({B}_{1}). The reciprocal was also maintained where initial strain configurations with high AFS1 predictably generated high competitive outcome for ({B}_{1}). As for isogenic strains, this relationship was found to be linear with the slope providing a measure of the deterministic range of competitive outcomes for a given founder density. The relationship between AFS and competitive outcome was again shown to be robust to natural variation in the initial strain ratio inherent in low founding cell densities (Supplementary Fig. S6C).Our mathematical model predicted that coexistence remained possible over a range of maximum growth rates, (r) (within a two-fold difference between dimensional strain growth rates in the absence of competition), diffusion coefficients, (d) (within a three-fold difference between dimensional diffusion coefficients), and most surprisingly, any values of the competition coefficient, (c) (Section S6 and Supplementary Fig. S7A–C). In particular, we showed that a strain required extreme competition efficiency ((c) very large) in order to compensate for being slower in growth ((d,r ; > ; 1)) (Supplementary Fig. S7D). Finally, the predictive power of the AFS score was preserved over the parameter range tested (Supplementary Fig. S7E, F).Dual-isolate biofilm assays – selection of a competition partnerTo experimentally test our model predictions, we needed to identify a suitable partner for NCIB 3610. We chose a Bacillus subtilis strain called NRS6153 (hereafter 6153). This selection was made because (i) 6153 is a genetically competent wild type strain with no known auxotrophies [36]); (ii) in liquid culture conditions the generation times of the two strains are within ~1.5-fold of each other (Fig. 6A); (iii) under biofilm conditions, single strain biofilms of both isolates have footprint sizes that are within (sim)2-fold of each other (Fig. 6B); (iv) across a broad range of founder densities, the competitive outcome of an isogenic pairing of 6153 isolates in a colony biofilm is broadly similar to that of an isogenic pairing of 3610 strains, albeit with more variability in the competitive outcome at the 72-h time point for high founder densities (cf. Fig. 4C (Supplementary Fig. S5A) and Fig. 6C (Supplementary Fig. S8A)); (v) when a colony biofilm is founded at high density with marked strains of 3610 and 6153 starting at an initial 1:1 ratio, 6153 is consistently outcompeted by 3610 (and hence defines 3610 as the stronger strain in the context of this study) (Fig. 6D); and (vi) using an antibiosis halo formation assay, interrogation of the interaction between 3610 and 6153 showed no evidence of contact-independent growth inhibition (Fig. 6E). In combination, these data allow us to infer that the mode of competition during co-culture in the colony biofilm is locally antagonistic.Fig. 6: Selection of a competitive strain.A Growth curves of 3610 (black) and 6153 (grey) in MSgg cultures at 30 °C. The three lines shown for each isolate represent separate biological repeats. B Biofilm footprint area of single-strain 3610 and 6153 biofilms. Data from 18 and 16 biofilms are shown for the 24 h and 48 h timepoint, respectively. C Competitive outcome data from colony biofilm assays of isogenic 6153 biofilms are shown after 24 h, 48 h and 72 h of incubation. Plotted are the technical repeats from one biological repeat. The full data set is presented in Supplementary Fig. S8A. D Flow cytometry data of mixed biofilms grown for 24, 48, and 72 h at 30 °C on MSgg media. Isolate names followed by ‘g’ represent strains constitutively producing  GFP, (green on the graph). Isolate names followed by ‘b’ indicate strains constitutively producing mTagBFP, (magenta on the graph). Three biological and three technical replicates were performed for each strain mix and timepoint and all data points are shown. The error bars represent the mean standard deviation. E Halo formation assays on MSgg agar plates at 24 h of growth. Strains producing mTagBFP (magenta) and GFP (green) are shown.Full size imageDual-isolate biofilm assays confirm modelling hypothesesWe performed dual strain biofilm assays competing 3610 and 6153 over a wide range of founder densities. These competitive assays confirmed the modelling prediction that in biofilms inoculated at low founder densities, coexistence within a non-isogenic strain pair is enabled by spatial segregation (Fig. 7A). Under such conditions, the intrinsically weaker strain (6153) formed spatial sectors and thus was able to coexist with the stronger strain (3610) through spatial segregation (Fig. 7A, B). In contrast, and again as predicted by the mathematical model (and reported during the selection of strain 6153 as a competition partner), for biofilms inoculated at high founder density, 3610 competitively excluded 6153 (Fig. 7A, B, Supplementary Fig. S8B). Finally, a computation of AFS scores based on images taken after 24 h of incubation showed strong correlation between a strain’s AFS score and its competitive outcome after both 48 h and 72 h of incubation for both 6153 alone and when in co-culture with 3610 (Supplementary Figs. S9 and 7C).Fig. 7: Experimental data for a non-isogenic strain pair with local antagonistic interactions.A Example dual-strain biofilms (3610 labelled with GFP (green), 6153 labelled with mTagBFP (magenta)). Images taken after 72 h of incubation for two different founder densities. Scalebars as in Fig. 2. B Competitive outcome data for 3610 in the 3610/6153 pair after 24 h, 48 h and 72 h of biofilm incubation. Plotted are technical repeats from one biological repeat of the experiment. The full data set is presented in Supplementary Fig. S8B. C The relationship between AFS and competitive outcome for 6153. AFS1 was calculated based on images taken after 24 h of biofilm incubation, and competitive outcome after 48 h (top, ({n}=22)) and 72 h (bottom, (n=17)).Full size image More

Experimental site, location and climateFive years’ field experimentation on ICM was started in 2014–15 at the ICAR-Indian Agricultural Research Institute (28°35′ N latitude, 77°12′ E longitude, 229 m MSL), New Delhi, India. The study site comes under the ‘Trans IGPs’, being semi-arid with an average annual rainfall of 650 mm, of which ~ 80% occurs in July–September (south-west monsoon). The mean max. / min. air temperature ranges between 20-40ºC and 4-28ºC, respectively. The five years (2014–2019) weather data were recorded from the observatory adjoining to the experimental field, and presented in Supplementary Table 1. Before start of the experiment, a rainy season Sesbania was grown in 2014 to ensure the uniform fertility across the blocks. Initial soil samples (0.0–0.15 m depth) were collected in October 2014 after incorporating the Sesbania residues in soil. The soil samples were processed for the chemical analysis. The study site had a pH of 7.9 (1:2.5 soil and water ratio)68, 3.8 g kg−1 soil organic-C69, 94.1 kg ha−1 KMnO4 oxidizable N70, 97 µg g−1 soil microbial biomass carbon71, 51.3 μg PNP g−1 soil h−1 alkaline phosphatase72, 53.0 μg TPF g−1 soil d−1 dehydrogenase73, and 13.5 μg NH4-N g−1 soil h−1urease74.Description of different ICM modulesThe eight ICM modules were tested, comprising of four conventional tillage (CT)-based (ICM1-4) and four conservation agriculture (CA)-based (ICM5-8) modules, replicated thrice in a complete randomized block design with the plot size of 60 m2 (15 m × 4.5 m) (Table 4). The crop residues were completely removed in the CT-based modules (ICM1-4), while in the ICM5-8 modules, in-situ wheat (~ 3 Mg ha−1 on dry weight basis)) and maize (~ 5 Mg ha−1, on dry weight basis) residues were retained on the soil surface during all the seasons of crops cultivation (Footnote Table 4, Fig. 6a,b).Table 4 Description of integrated crop management (ICM) modules adopted in maize and wheat crops during the five yearsˈ fixed plot experimentation.Full size tableIn the ICM1-4 modules, the field preparation was carried out by sequential tillage operations, such as, deep ploughing using the disc harrow, cultivator/rotavator twice (0.15–0.20 m), followed by levelling in each season. In the ICM3-4, the raised beds of 0.70 m bed width (bed top 0.40 m and furrow 0.30 m) were formed during each cropping cycle using the tractor mounted bed planter, and simultaneously wheat sowing was done (Fig. 6c). In the case of maize, ridges (0.67 m length) were prepared using the ridge maker. In the CA-based ICM5-8 modules, the tillage operations, such as, seed and fertilizer placement were restricted to the crop row-zone in maize and wheat both. In the ICM7&8, the permanent raised beds (0.67 m mid-furrow to mid-furrow, 0.37 m wide flat tops, and 0.15 m furrow depth), were prepared (Fig. 6d). However, these beds were reshaped using the disc coulter at the end of each cropping cycle without disturbing the surface residues. The sowing was accomplished using the raised bed multi-crop planter.Cultural operations and the fertilizer applicationDuring every season, the maize (cv. PMH 1) was sown in the first week of July using 20 kg seed ha−1. The wheat (cv. HD 2967) crop was sown in the first fortnight of November using the seed-cum fertilizer drill (ICM1-2), bed planter (ICM3-4) and zero-till seed drill (ICM5-8) at 100 kg seed ha−1. The chemical fertilizers (N, P and K) were applied as per the modules described in the footnote of Table 4. At sowing, the full doses of phosphorous (P) and potassium (K) were applied using the di-ammonium phosphate (DAP) and muriate of potash (MOP), and the nitrogen (N) supplied through DAP. The remaining N was top-dressed through urea in two equal splits after the first irrigation and tasseling / silking stages in maize, and crown root initiation and tillering stages of wheat. In the modules receiving ¾ fertilizers (ICM2,4,6,8), the seeds were treated with the NPK liquid bio-fertilizer (LBFs) (diluted 250 ml formulation 2.5 L of water ha−1), and an arbuscular mycorrhiza (AMF) was broadcasted at 12 kg ha−1 as has been described by75. This LBFs had the microbial consortia of N-fixer (Azotobacter chroococcum), P (Pseudomonas) and K (Bacillus decolorationis) solubilizers, procured from the commercial biofertilizer production unit of the Microbiology Division, ICAR-Indian Agricultural Research Institute, New Delhi (Patentee: ICAR, Govt. of India). Weeds were managed by integrating the pre- and post-emergence herbicides, and their combinations along with the hand weeding-mulching, as mentioned in the concerned modules (Footnote Table 4). However, in the CA-based modules (ICM5-8), the non-selective herbicide glyphosate (1 kg ha−1) was used 10 days before the sowing. The need-based integrated insect-pests and disease management practices were followed uniformly across the modules.Soil sampling and analysisBefore start of the experiment, the soil sampling was done from 0.0–0.15 m depth. Afterwards, five random samples from each module from 0.0–0.30 m soil depth were collected at the flowering stage of 5th season wheat. These samples were taken from the three soil depths (0.0 to 0.05, 0.05–0.15 and 0.150–0.30 m) using the core sampler. The ground, air-dried soil samples, passed through a 0.2 mm sieve were used for the determination of the Walkley and Black organic carbon (SOC), as described by76. For the soil biological properties, the soil samples were processed, and stored at 5ºC for 18–24 h, then analyzed the soil microbial biomass carbon (SMBC), dehydrogenase (SDH), alkaline phosphate (SAP) and the urease (URE) activities.The soil microbial biomass carbon (SMBC)The SMBC was measured using the fumigation extraction method as proposed by71. The pre-weighed samples from the respective soil depths were fumigated with the ethanol-free chloroform for the 24 h. Separately, a non–fumigated set was also maintained. Further, 0.5 M K2SO4 (soil: extractant 1:4) was added, and kept on a reciprocal shaker for 30 min. and then filtered through a Whatman No. 42 filter paper. OC of the filtrate was measured through the dichromate digestion, followed by the back titration with 0.05 N ferrous ammonium sulphate. The SMBC was then calculated using the equation:$${text{S}}_{{{text{MBC}}}} = {text{EC }} times { 2}.{64}$$where, EC = (Corg in fumigated soil – Corg in non-fumigated soil), and expressed in µg C g−1 soil.The dehydrogenase activity (SDH)The SDH activity (μg TPF g−1 soil d−1) was assessed using the method of73. The soil sample (~ 6 g) was saturated with 1.0 ml freshly prepared 3% triphenyltetrazolium chloride (TTC), and then incubated for 24 h under the dark. Later on, the methanol was added to stop the enzyme activity, and the absorbance of the filtered aliquot was read at 485 nm.The alkaline phosphatase activity (SAP)The APA activity was estimated in 1.0 g soil saturated with 4 ml of the modified universal buffer (MUB) along with 1 ml of p-nitrophenol phosphate followed by incubation at 37 °C for 1 h. After incubation, 1 ml of 0.5 M CaCl2 and 4 mL of NaOH were added and the contents filtered through Whatman No. 1 filter paper. The amount of p-nitrophenol in the sample was determined at 400 nm72 and the enzyme activity was expressed as µg p-NP g−1 soil h−1.The urease activityUrease activity was measured using 10 g soil suspended in 2.5 ml of urea solution (0.5%). After incubating for a day at 37 °C, 50 ml of 1 M KCl solution was added. This was kept on a shaker for 30 min and the aliquot was filtered through Whatman No. 1 filter paper. To the filtrate (10 ml), 5 ml of sodium salicylate and 2 ml of 0.1% sodium dichloro-isocyanide solution were added and the green color developed was measured at 690 nm74. These values are reported as µg NH4-N g−1 soil h−1.Water application and productivityIn experimental modules, water was given through the controlled border irrigation method. The current meter was fixed in the main lined rectangular channel, and the water velocity was measured. To get the flow discharge, then multiplied with area of cross section of the channel. The following formulae were used to calculate the applied irrigation water quantity and depth3:$${text{Irrigation water applied }}left( {text{L}} right) , = {text{ F }} times {text{ t (i)}}$$$${text{Depth }}left( {{text{mm}}} right) , = {text{ L}} div {text{A}}/{ 1}000$$where, F is flow rate (m3 s−1), t is time (s) taken in each irrigation in each module and A is area (m2).The effective precipitation (EP, difference between total rainfall and the actual evapotranspiration) was calculated, and then EP was added to the irrigation water applied to calculate the total water applied in each module. Across the maize and wheat modules (ICM1-8), irrigations were given at the critical growth stages, such as, knee high and silking / tasseling (maize) and crown root formation, maximum tillering, flowering, heading / milking (wheat) stages, and after long dry spell (≥ 10-days).On the basis of the soil water depletion pattern (at the depth of 0.60 m), in each season, 3–6 irrigations were given to maize, while wheat received 5–8 irrigations per season or crop including the pre-sowing irrigation. The rainfall data were obtained from the meteorological observatory located in the adjoining field. The water productivity (kg grains ha−1 mm−1 of water) was measured as per the equation given below:$${text{Water productivity }} = {text{ economic yield }}left( {{text{kg ha}}^{{ – {1}}} } right)/{text{ total water applied }}left( {{text{mm}}} right)$$Additionally, the systems water productivity (SWP) was also estimated by adding the water productivity (WP) of both maize and wheat crops grown under the MWR.Yield measurementsIn each season, the maize and wheat crops were harvested during the months of October and April, respectively, leaving 0.75 m border rows from all the corners of each module. The crops were harvested from the net sampling area (6 m × 3 m, 18 m2) located at the center of each plot. Maize crop was harvested manually and the wheat by using the plot combine harvester. All the harvested produce was sun dried before threshing and the grain and straw / stover yields were weighed separately. The stover/straw yields were measured by subtracting the grain weight from the total biomass. To compare the total (system) productivity of the different ICM modules, the system yield was computed, taking maize as the base crop, i.e., the maize equivalent yield (MGEY) using the equation20:$${text{M}}_{{{text{GEY}}}} left( {{text{Mg ha}}^{{ – {1}}} } right) , = {text{ Ym }} + , left{ {left( {{text{Yw }} times {text{ Pw}}} right) , div {text{ Pm}}} right}$$where, Ym = maize grain yield (Mg ha−1), Yw = wheat grain yield (Mg ha−1), Pm = price of maize grain (US$ Mg−1) and Pw = price of wheat grain (US$ Mg−1).Farm economicsUnder different ICM modules, the variable production costs and economic returns were worked out based on the prevailing market prices for the respective years. The production costs included the cost of various inputs, such as, rental value of land, seeds, pesticides, LBFs / consortia, AMF, labor, and machinery; tillage / sowing operations, irrigation, mineral fertilizers, plant protection, harvesting, and threshing etc. The costs for the crops’ residues were also considered. The system total returns were computed by adding the economic worth of the individual crop, however, the net returns were the differences between the total returns to the variable production costs of the respective module. The Govt. of India’s minimum support prices (MSP) were considered for the conversion of grain yield to the economic returns (profits) during the respective years. Further, the system net returns (SNR) were worked out by summing the net income from both maize and the wheat in Indian rupees (INR), and then converted to the US$, based on the exchange rates for different years.Sustainable yield index (SYI)77,78described the SYI as a quantitative measure of the sustainability of agricultural rotation/practice. The sustainability could be interpreted using the standard deviation (σ) values, where the lower values of the σ indicate the greater sustainability and vice-versa. Total crop productivity of maize and wheat under the different ICM modules was computed based on the five years’ mean yield data. SYI was calculated using equation78.$${text{S}}_{{{text{YI}}}} = , left( {{-}{overline{text{Y}}}_{{{text{a }}{-}}} sigma_{{text{n}}} {-}_{{1}} } right) , /{text{ Y}}^{{{-}{1}}}_{{text{m}}}$$where, –ȳa is the average yield of the crops across the years under the specific management practice, σn–1 is the standard deviation and Y–1 m is the maximum yield obtained under the set of an ICM module.Statistical analysisThe GLM procedure of the SAS 9.4 (SAS Institute, 2003, Cary, NC) was used for the statistical analysis of all the data obtained from different ICM modules to analyze the variance (ANOVA) under the randomized block design79. Tukey’s honest significant difference test was employed to compare the mean effect of the treatments at p = 0.05.Authors have confirmed that all the plant studies were carried out in accordance with relevant national, international or institutional guidelines. More