More stories

  • in

    Using past interglacial temperature maxima to explore transgressions in modern Maldivian coral and Amphistegina bleaching thresholds

    Study site and target foraminiferal speciesThe Maldivian archipelago is a partially drowned carbonate platform within the central, equatorial Indian Ocean. It consists of two rows of north–south orientated atolls, which encompass an Inner Sea. The lowermost neritic carbonate unit sits upon volcanic bedrock and has been dated back to the Eocene19 with continuous drift deposition, within the Inner Sea, starting ~ 12.9 Ma at the establishment of the modern South Asian Monsoon (SAM)35,36. This seasonally reversing, major climatic system has an impact on both the regional precipitation patterns as well as physiochemical oceanographic properties (Fig. 1). The summer southwest SAM brings warm, wet conditions to the Indian subcontinent, as well as higher saline surface waters from the Arabian Sea into the Maldives region. In comparison, the winter northeast SAM results in cool, dry continental conditions and transports lower salinity water from the Bay of Bengal into the central, equatorial Indian Ocean. As a result, the Maldives seasonal salinity depth profiles can vary significantly, yet due to its tropical location the seasonal sea water temperatures are relatively stable.Three symbiont-bearing foraminiferal species are used in this study: Amphistegina lessonii, Globigerinoides ruber (white) and Trilobatus sacculifer (with sac-like final chamber):Amphistegina lessonii is a larger benthic, symbiont-bearing (diatoms) foraminiferal species. It has a shallow depth range (0–50 m)37,38,39 and is globally abundant in tropical coral reef, benthic foraminiferal shoal and general carbonate shelf settings40. Similarly to corals, amphisteginids have been shown to bleach under high temperatures/high irradiance levels with the new development of the Amphistegina Bleaching Index (ABI) as an indicator of photo-inhibitory stress in coral reef settings41,42. From ~ 30 °C this species starts showing signs of thermal stress, with bleaching and mortality reported for temperatures  > 31 °C11,12.Globigerinoides ruber (w) hosts dinoflagellate endosymbionts and is the most common planktonic foraminiferal species in tropical-subtropical waters13 state that while G. ruber (w) is generally considered one of the shallowest-dwelling species, its depth distribution does vary in relation to regional ecological conditions. It has a particular relation to the nutricline depth in less turbid, oligotrophic conditions43 which has been confirmed for the Maldives28. It is omnivorous, however in comparison to other omnivorous, symbiont-bearing species, it has demonstrated an elevated adaptation for consuming phytoplankton protein over zooplankton protein13. From culture experiments, it has a broad temperature (14–31 °C) and salinity (22–49 PSU) tolerance, and has been reported as the most tolerant species to low sea surface salinity (SSS)13. This species occurs year-round and has a fortnightly reproduction13.Trilobatus sacculifer is a planktonic foraminiferal species abundant in tropical-subtropical surface waters and as such is extensively used in paleo-reconstructions. It hosts dinoflagellate endosymbionts yet is omnivorous, feeding predominantly on calanoid copepods13. It is a euryhaline species, with a broad salinity (24–47 PSU) and temperature (14–32 °C) tolerance. Similarly to G. ruber (w), this species occurs year-round and has a monthly reproduction on a synodic lunar cycle13. While a shallow dwelling species, it is generally reported to live slightly deeper in the water column, in comparison to G. ruber (w)28,30,44.SamplingAll planktonic foraminiferal specimens (G. ruber (w) and T. sacculifer (w/s)) for the geochemical analysis (δ18Oc and Mg/Ca) originate from the International Ocean Discovery Program (IODP) Expedition 359, Site U1467 (4° 51.0274′ N, 73° 17.0223′ E) drilled in 2015 within the Inner Sea of the Maldivian archipelago at a water depth of 487 m19. The age model for these samples was adopted from a previous study45 which is based on the correlation of their long-term (0–1800 kyr) Site 359-U1467 C. mabahethi and G. ruber (w) δ18Oc records to the stacked reference curve of46. Recent surface sediment samples (mudline A and B: representing the sample from the sediment/water interface), as well as three samples from the peak of MIS9e (U1467C, 2H6, 0–1 cm; U1467C, 2H6, 15–16 cm; U1467C, 2H6, 18–19 cm) and MIS11c (U1467B, 3H2, 147–148 cm; U1467B, 3H3, 9–10 cm; U1467B, 3H3, 12–13 cm) were analysed19,28 (sample locations are shown on Fig. 3). The mudline is identified as Recent, likely representing the last few hundred years, based on the presence of Rose Bengal (1 g/L) stained ostracods and benthic foraminifera. The study by45 has verified that diagenetic influences, within this shallow, carbonate environment, are not a concern for foraminiferal geochemical compositions over the investigated time-interval (MIS1-11).Rose Bengal stained A. lessonii specimens were obtained from modern surface rubble samples collected by hand, at 10 m water depth, during the 2015 International Union for Conservation of Nature (IUCN) REGENERATE cruise47 (Supplementary Table 6). Samples were collected from the reefs of two islands, Maayafushi and Rasdhoo, both located within the central part of the Maldivian archipelago. As the foraminifera shells were stained pink, it implies they were living at the time of collection. These specimens were used for stable isotopic analysis and their reconstructed temperatures represent modern (a cumulative signal encompassing their lifespan of four to twelve months48) conditions (Supplementary Tables 5–6). A full explanation of the Rose Bengal protein stain for foraminifera is detailed in49.δ18Oc stable isotopic analysisAll samples were initially washed using a 32 μm sieve to remove the finer clay and silt fractions. Subsequently, they were air dried and sieved into discrete sizes for foraminiferal picking. To ensure enough calcite for the measurements, all specimens for Individual Foraminifera Analysis (IFA) for both G. ruber (w) and T. sacculifer (w/s) (n = 632) were picked from the 355–400 μm size fraction. In addition, traditional whole-shell (pooled) measurements for G. ruber (w) (n = 24) were conducted on specimens from the 212–400 μm fraction (2–5 pooled specimens). Trilobatus sacculifer (w/s) traditional whole-shell analysis (n = 21) was measured on specimens (2 pooled specimens) from the 300–355 μm fraction. The majority of these pooled measurements are obtained from28,45,50,51 (Supplementary Table 1). Amphistegina lessonii measurements were run on single specimens  > 250 μm in size. Prior to stable isotopic analysis, all shells were briefly cleaned (1–2 s) by ultrasonication in Milli-Q water to remove any adhering particles. All stable isotopic measurements were conducted at the School of GeoSciences at the University of Edinburgh on a Thermo Electron Delta + Advantage mass spectrometer integrated with a Kiel carbonate III automated extraction line. Samples were reacted with 100% phosphoric acid (H3PO4) at 90 °C for 15 min, with the evolved CO2 gas collected in a liquid nitrogen coldfinger and analysed compared to a reference gas. All samples are corrected using an internal laboratory standard and expressed as parts per mil (‰) relative to Vienna Pee Dee Belemnite (VPDB). Replicate measurements of the standards give the instrument an analytical precision (1σ) of ~ 0.05 ‰ for δ18O and δ13C.Mg/Ca analysisThe Mg/Ca data is obtained from28,45,50,51 (Supplementary Table 1). Each G. ruber (w) Mg/Ca analysis (n = 17; 212–250 μm in size) was conducted on 30 pooled specimens by inductively coupled plasma optical emission spectrometry (ICP-OES) on a Thermoscientific iCap 6300 (dual viewing) at the Institute of Geosciences of the Goethe-University of Frankfurt. All samples were initially cleaned (1–2 s) by ultrasonication in Milli-Q water and then the standard oxidative cleaning protocol of52 followed to prevent clay mineral contamination. The final centrifuged sample solution was diluted with an yttrium solution (1 mg/l) prior to measurement to allow for the correction of matrix effects. In addition, before each analysis five calibration solutions were measured to allow for intensity ratio calibrations. All element/Ca measurements were standardized using an internal consistency standard (ECRM 752–1, 3.761 mmol/mol Mg/Ca). Furthermore, the elements Al, Fe, and Mn were screened and blanks periodically run to monitor for further signs of contamination during the analyses.Establishment of present and past seawater temperaturesPrior to temperature calculations, we test the IFA distributions for normality using the Shapiro‐Wilk test and the Fisher–Pearson coefficient of skewness with bootstrap confidence intervals, to define the skewness of the datasets53 (Supplementary Table 3). The Recent G. ruber (w) and T. sacculifer (w/s) and MIS11c T. sacculifer are normally distributed. In the case of both MIS9e datasets and the MIS11c G. ruber population, the null hypothesis that the data are normally distributed (p ≤ 0.05) is rejected (Supplementary Table 3). Considering bioturbation within the sediment record is a possibility, we use two methods to identify and remove outliers in the IFA datasets. Firstly, the inter-quartile range (IQR) is used for each δ18Oc dataset, which defines a measurement as an outlier if it falls outside the range [Q1 − 1.5 (Q3 − Q1), Q3 + 1.5 (Q3 − Q1)], with IQR = Q3 − Q1 and Q3 and Q1 representing the third and first quartile of the dataset20. But if there is considerable reworking, the IQR method would not necessarily identify reworked glacial measurements (highest δ18Oc values) within the interglacial samples. As such, the Recent IFA datasets, which are both normally distributed, are used to further set a rudimentary cut-off point for the highest δ18Oc (= lowest temperatures) value to expect during past interglacial minima periods for both G. ruber (w) and T. sacculifer (w/s) (this is discussed further in the Supplementary Materials, Supplementary Figs. 1–3).There are innumerable analytical techniques (e.g., traditional mass spectrometry, secondary-ion mass spectrometry, laser ablation inductively coupled plasma mass spectrometry), proxies (Mg/Ca, δ18O, clumped isotopes, TEX86, Uk’37) as well as target medians (e.g., calcitic shells of foraminifera, aragonitic coral skeletons, ice, lipids, alkenones) which are used in marine paleo-temperature reconstructions. Furthermore, different methods exist in the literature to calculate temperature estimates using both planktonic foraminiferal δ18Oc and Mg/Ca measurements with innumerable species-specific δ18O-temperature and Mg/Ca-temperature equations reported20,23,30,54,55,56. Moreover, due to the exponential nature of the Mg/Ca-temperature equations, if inappropriately applied, offsets in the upper temperature range are exacerbated. Additional considerations are species-specific offsets and differential geochemical compositions within the shell (e.g., high versus low Mg banding, gametogenic calcite). Trilobatus sacculifer gametogenic calcite has been reported to be significantly enriched in Mg in comparison to the rest of the shell57. As T. sacculifer specimens selected for use in this study underwent reproduction, indicated by the presence of a sac-like final chamber58, we can expect their Mg/Ca ratios to be biased. As such, to avoid overestimates we chose to use only G ruber (w, pooled) Mg/Ca and δ18Oc data to calculate representative δ18Osw values for each time interval, for use with both the G. ruber (w) and T. sacculifer (w/s) δ18Oc IFA datasets. Considering both planktonic species are considered as shallow-dwellers with similar living depths and an affinity for the DCM, the utilisation of common δ18Osw values is applicable13,28,30.The G. ruber Mg/Ca-temperature Eq. (1) from55 (temperature calibration range: ~ 22–27 °C), similarly applied in the Maldivian study of28, was used in this study:$$Mg/Ca=0.34left(pm 0.08right)mathrm{exp}(0.102left(pm 0.010right)*T)$$
    (1)
    The applied δ18O-temperature species-specific equations (Eqs. 2 and 3) were previously utilised in the local study by28. Both the G. ruber (Eq. 2) and T. sacculifer (Eq. 3) equations are from the Indian Ocean study of59 (temperature calibration range: ~ 20–31 °C):$$T=12.75-5({delta }^{18}{O}_{c}-{delta }^{18}{O}_{sw})$$
    (2)
    $$T=11.95-5.26({delta }^{18}{O}_{c}-{delta }^{18}{O}_{sw})$$
    (3)
    Using the above equations, the range in temperature estimates are obtained as follows (Fig. 4):

    1.

    The mean G. ruber (w) Mg/Ca measurements are used together with Eq. (1) to calculate a temperature estimate for each time point (Supplementary Table 1). Since the Mg/Ca calcification temperatures are based on 30 pooled specimens, they are considered to reflect mean calcification temperatures.

    2.

    The Mg/Ca derived temperature estimates are then used together with the mean traditional (pooled) G. ruber (w) δ18Oc data and Eq. (2) to calculate representative δ18Osw values for each time point (Supplementary Table 2). As these are calculated from pooled samples, they are considered to mirror mean δ18Osw values for both the Recent and fossil populations.

    3.

    The G. ruber (w) and T. sacculifer (w/s) IFA datasets are then used, together with the relevant species-specific δ18O-temperature equations and δ18Osw values, to calculate the spread in temperature estimates (Fig. 4, Supplementary Tables 3–4).

    Trilobatus sacculifer (w/s) data from the glacial maxima of MIS12 are included in the study to illustrate the applicability of the IFA method, however, as they do not contribute to the discussion on bleaching thresholds, they are discussed further in the Supplementary Materials (Supplementary Figs. 1, 3).Finally, the temperature estimates for the shallow-dwelling symbiont-bearing benthic A. lessonii are obtained using the genus-specific δ18O-temperature equation of60 (Eq. 4) (Supplementary Tables 5–6).$$T=16.3-4.24({delta }^{18}{O}_{c}-{delta }^{18}{O}_{sw})$$
    (4)
    Considering the benthic specimens were deemed living at the time of collection (Rose Bengal stained), a mean regional surface (0 m) δ18Osw value (0.49 ‰) is used together with the δ18Oc data in the calculations (Supplementary Tables 5–6). More

  • in

    Ecological recipes for selecting community function

    1.Doolittle, W. F. & Zhaxybayeva, O. BioScience 60, 102–112 (2010).Article 

    Google Scholar 
    2.Liautaud, K., van Nes, E. H., Barbier, M., Scheffer, M. & Loreau, M. Ecol. Lett. 22, 1243–1252 (2019).PubMed 
    PubMed Central 

    Google Scholar 
    3.Morris, A., Meyer, K. & Bohannan, B. Phil. Trans. R. Soc. B Biol. Sci. 375, 20190244 (2020).CAS 
    Article 

    Google Scholar 
    4.De Monte, S. & Rainey, P. B. J. Biosci. 39, 237–248 (2014).PubMed 

    Google Scholar 
    5.Chang, C.-Y. et al. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-021-01457-5 (2021).6.Goldford, J. E. et al. Science 361, 469–474 (2018).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    7.Machado, D. et al. Nat. Ecol. Evol. 5, 195–203 (2021).PubMed 
    Article 

    Google Scholar 
    8.Ratzke, C., Barrere, J. & Gore, J. Nat. Ecol. Evol. 4, 376–383 (2020).PubMed 
    Article 

    Google Scholar 
    9.Wright, E. S. & Vetsigian, K. H. Nat. Commun. 7, 11274 (2016).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    10.Szathmáry, E. & Demeter, L. J. Theor. Biol. 128, 463–486 (1987).PubMed 
    Article 

    Google Scholar 
    11.Xie, L., Yuan, A. E. & Shou, W. PLoS Biol. 17, e3000295 (2019).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    12.van Vliet, S. & Doebeli, M. Proc. Natl Acad. Sci. USA 116, 20591–20597 (2019).PubMed 
    Article 

    Google Scholar 
    13.Doulcier, G., Lambert, A., De Monte, S. & Rainey, P. B. eLife 9, e53433 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    14.Black, A. J., Bourrat, P. & Rainey, P. B. Nat. Ecol. Evol. 4, 426–436 (2020).PubMed 
    Article 

    Google Scholar 
    15.Williams, H. T. P. & Lenton, T. M. Proc. Natl Acad. Sci. USA 104, 8918–8923 (2007).CAS 
    PubMed 
    Article 

    Google Scholar 
    16.Barbier, M., Arnoldi, J.-F., Bunin, G. & Loreau, M. Proc. Natl Acad. Sci. USA 115, 2156–2161 (2018).CAS 
    PubMed 
    Article 

    Google Scholar 
    17.Pearce, M. T., Agarwala, A. & Fisher, D. S. Proc. Natl Acad. Sci. USA 117, 14572–14583 (2020).CAS 
    PubMed 
    Article 

    Google Scholar 
    18.Roy, F., Barbier, M., Biroli, G. & Bunin, G. PLoS Comput. Biol. 16, e1007827 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar  More

  • in

    Faster monitoring of the invasive alien species (IAS) Dreissena polymorpha in river basins through isothermal amplification

    1.EU. No 1143/2014 of the European Parliament and of the Council of 22 October 2014 on the prevention and management of the introduction and spread of invasive alien species. Off. J. Eur. Union 2014, 35–55 (2014).
    Google Scholar 
    2.Ludyanskiy, M. L., McDonald, D. & MacNeill, D. Impact of the Zebra Mussel, a Bivalve Invader. Bioscience 43, 533–544 (1993).Article 

    Google Scholar 
    3.Lalaguna, C. D. & Marco, A. A. The zebra mussel invasion in Spain and navigation rules. Aquat. Invasions 3, 315–324 (2008).Article 

    Google Scholar 
    4.Rajagopal, S. et al. Origin of spanish invasion by the zebra mussel, dreissena polymorpha (pallas, 1771) revealed by amplified fragment length polymorphism (AFLP) fingerprinting. Biol. Invasions 11, 2147–2159 (2009).Article 

    Google Scholar 
    5.CABI Invasive species Compendium: Dreissena polymorpha (zebra mussel) 2017. https://www.cabi.org/isc/datasheet/85295.6.Marescaux, J. & Van Doninck, K. Using DNA barcoding to differentiate invasive Dreissena species (Mollusca, Bivalvia). Zookeys 365, 235–244 (2013).Article 

    Google Scholar 
    7.Benson, A. J., Raikow, D., Larson, J., Fusaro, A., & Bogdanoff, A. K. Dreissena polymorpha (Pallas, 1771): U.S. Geological Survey, Nonindigenous Aquatic Species Database. Gainesville, FL (2017).8.Minchin, D., Lucy, F. & Sullivan, M. Zebra Mussel: Impacts and Spread. Invasive Aquat. Species Eur. Distrib. Impacts Manag. 135–146. © 2002 Kluwer Acad. Publ. Dordreicht, Netherlands. 135–148 (2002) doi:https://doi.org/10.1007/978-94-015-9956-6_15.9.Montero Melendez, J. Control of invasive alien species in Guadalquivir river basin. EURO-RIOB 2017. https://www.riob.org/es/node/404510.Molloy, D. P., Karatayev, A., Burlakova, L. E., Kurandina, D. P. & Laruelle, F. Natural enemies of zebra mussels: predators, parasites, and ecological competitors. Rev. Fish. Sci. 5, 27–97 (1997).Article 

    Google Scholar 
    11.Nalepa, T. F. & Schloesser, D. W. Zebra Mussels Biology, Impacts, and Control (Lewis Publishers, Boca Raton, 1993).
    Google Scholar 
    12.Birnbaum, C. NOBANIS – Invasive Alien Species Fact Sheet – Dreissena polymorpha. Accessed 2 Dec 2019. https://www.nobanis.org (Online Database of the European Network on Invasive Alien Species – NOBANIS, 2011).13.Lowe, S., Browne, M., Boudjelas, S., De Poorter, M. 100 of the World’s Worst Invasive Alien Species A selection from the Global Invasive Species
    Database. (The Invasive Species Specialist Group (ISSG), 2000).14.Glomski, L. M. Zebra Mussel Chemical Control Guide. US Army Corps of Engineers: Waterways Experiment Station. https://erdclibrary.erdc.dren.mil/jspui/bitstream/11681/6966/1/ERDC-EL-TR-15-9.pdf (2015). 15.Boelman, S. F., Neilson, F. M., Dardeau, E. A. & Cross, T. Zebra mussel (Dreissena polymorpha) control handbook for facility operators, first edition . US Army Corps of Engineers: Waterways Experiment Station. https://hdl.handle.net/11681/2966 (1997).16.Durán, C., Lanao, M., Anadón, A. & Touyá, V. Management in practice management strategies for the zebra mussel invasion in the Ebro River basin. Aquat. Invasions 5, 309–316 (2010).Article 

    Google Scholar 
    17.Bij de Vaate, A. Rajagopal, S. & van der Velde, G. The zebra mussel in Europe: summary and synthesis. in The Zebra Mussel in Europe (ed van der Velde, G.
    et al.) 415–421 (Backhuys Publishers, 2010).18.Herder, J. et al. Environmental DNA—a review of the possible applications for the detection of (invasive) species. Report 2013-104. Accessed 3 May 2019. https://www.researchgate.net/publication/283267157_Environmental_DNA_-_a_review_of_the_possible_applications_for_the_detection_of_invasive_species#fullTextFileContent (Stichting RAVON, Nijmegen, 2014).19.Xiong, W., Li, H. & Zhan, A. Early detection of invasive species in marine ecosystems using high-throughput sequencing: technical challenges and possible solutions. Mar. Biol. 163, 1–12 (2016).CAS 
    Article 

    Google Scholar 
    20.Harvey, C. T., Qureshi, S. A. & MacIsaac, H. J. Detection of a colonizing, aquatic, non-indigenous species. Divers. Distrib. 15, 429–437 (2009).Article 

    Google Scholar 
    21.Jerde, C. L., Mahon, A. R., Chadderton, W. L. & Lodge, D. M. ‘Sight-unseen’ detection of rare aquatic species using environmental DNA. Conserv. Lett. 4, 150–157 (2011).Article 

    Google Scholar 
    22.Dejean, T. et al. Improved detection of an alien invasive species through environmental DNA barcoding: the example of the American bullfrog Lithobates catesbeianus. J. Appl. Ecol. 49, 953–959 (2012).Article 

    Google Scholar 
    23.Jerde, C. L. et al. Detection of Asian carp DNA as part of a Great Lakes basin-wide surveillance program. Can. J. Fish. Aquat. Sci. 70, 522–526 (2013).CAS 
    Article 

    Google Scholar 
    24.Laramie, M. B., Pilliod, D. S. & Goldberg, C. S. Characterizing the distribution of an endangered salmonid using environmental DNA analysis. Biol. Conserv. 183, 29–37 (2015).Article 

    Google Scholar 
    25.Takahara, T., Minamoto, T. & Doi, H. Using environmental DNA to estimate the distribution of an invasive fish species in ponds. PLoS ONE 8, e56584 (2013).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    26.Gingera, T. D., Bajno, R., Docker, M. F. & Reist, J. D. Environmental DNA as a detection tool for zebra mussels Dreissena polymorpha (Pallas, 1771) at the forefront of an invasion event in Lake Winnipeg, Manitoba, Canada. Manag. Biol. Invasions 8, 287–300 (2017).Article 

    Google Scholar 
    27.Darling, J. A. & Mahon, A. R. From molecules to management: adopting DNA-based methods for monitoring biological invasions in aquatic environments. Environ. Res. 111, 978–988 (2011).CAS 
    PubMed 
    Article 

    Google Scholar 
    28.Kaprou, G. D. et al. Miniaturized devices for isothermal DNA amplification addressing DNA diagnostics. Microsyst. Technol. 22, 1529–1534 (2016).CAS 
    Article 

    Google Scholar 
    29.Mori, Y. & Notomi, T. Loop-mediated isothermal amplification (LAMP): a rapid, accurate, and cost-effective diagnostic method for infectious diseases. J. Infect. Chemother. 15, 62–69 (2009).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    30.Fang, X., Liu, Y., Kong, J. & Jiang, X. Loop-mediated isothermal amplification integrated on microfluidic chips for point-of-care quantitative detection of pathogens. Anal. Chem. 82, 3002–3006 (2010).CAS 
    PubMed 
    Article 

    Google Scholar 
    31.Rafati, A. & Gill, P. Microfluidic method for rapid turbidimetric detection of the DNA of Mycobacterium tuberculosis using loop-mediated isothermal amplification in capillary tubes. Microchim. Acta 182, 523–530 (2014).Article 
    CAS 

    Google Scholar 
    32.Notomi, T. et al. Loop-mediated isothermal amplification of DNA. Nucleic Acids Res. 28, e63 (2000).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    33.Mori, Y., Kitao, M., Tomita, N. & Notomi, T. Real-time turbidimetry of LAMP reaction for quantifying template DNA. J. Biochem. Biophys. Methods 59, 145–157 (2004).CAS 
    PubMed 
    Article 

    Google Scholar 
    34.Tomita, N., Mori, Y., Kanda, H. & Notomi, T. Loop-mediated isothermal amplification (LAMP) of gene sequences and simple visual detection of products. Nat. Protoc. 3, 877–882 (2008).CAS 
    PubMed 
    Article 

    Google Scholar 
    35.Garrido-Maestu, A., Fuciños, P., Azinheiro, S., Carvalho, J. & Prado, M. Systematic loop-mediated isothermal amplification assays for rapid detection and characterization of Salmonella spp., Enteritidis and Typhimurium in food samples. Food Control 80, 297–306 (2017).CAS 
    Article 

    Google Scholar 
    36.Verkaar, E., Nijman, I., Boutaga, K. & Lenstra, J. Differentiation of cattle species in beef by PCR-RFLP of mitochondrial and satellite DNA. Meat Sci. 60, 365–369 (2002).CAS 
    PubMed 
    Article 

    Google Scholar 
    37.Wilson-Wilde, L., Norman, J. & Robertson, J. et al. Current issues in species identification for forensic science and the validity of using the cytochrome oxidase I
    (COI) gene. Forensic Sci. Med. Pathol. 6, 233–241. https://doi.org/10.1007/s12024-010-9172-y (2010). CAS 
    PubMed 
    Article 

    Google Scholar 
    38.Hebert, P. D. N., Cywinska, A., Ball, S. L. & Jeremy, R. Biological identifications through DNA barcodes. Proc. Biol. Sci. 270(1512), 313–321 (2003).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    39.Staroscik, A. Copy number calculator for realtime PCR. http://www.scienceprimer.com/copy-number-calculator-for-realtime-pcr (2012).40.Ficetola, G. F., Miaud, C., Pompanon, F. & Taberlet, P. Species detection using environmental DNA from water samples. Biol. Lett. 4, 423–425 (2008).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    41.Zanoli, L. M. & Spoto, G. Isothermal amplification methods for the detection of nucleic acids in microfluidic devices. Biosensors 3, 18–43 (2013).CAS 
    PubMed 
    Article 

    Google Scholar 
    42.Egan, S. P. et al. Rapid molecular detection of invasive species in ballast and harbor water by integrating environmental DNA and light transmission spectroscopy. Environ. Sci. Technol. 49, 4113–4121 (2015).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    43.Xia, Z. et al. Early detection of a highly invasive bivalve based on environmental DNA (eDNA). Biol. Invasions 20, 437–447 (2018).Article 

    Google Scholar 
    44.Williams, M. R. et al. Isothermal amplification of environmental DNA (eDNA) for direct field-based monitoring and laboratory confirmation of Dreissena sp. PLoS ONE 12, e0186462 (2017).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    45.Frackman, B. S., Kobs, G., Simpson, D., Storts, D. & Corporation, P. Betaine and DMSO: enhancing agents for PCR. Promega Notes 65, 9–12 (1998).
    Google Scholar 
    46.Wang, D.-G., Brewster, J., Paul, M. & Tomasula, P. Two methods for increased specificity and sensitivity in loop-mediated isothermal amplification. Molecules 20, 6048–6059 (2015).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    47.Jantz, B. & Neumann, D. Growth and reproductive cycle of the zebra mussel in the River Rhine as studied in a river bypass. Oecologia 114, 213–225 (1998).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    48.Grigorovich, I. A., Kelly, J. R., Darling, J. A. & West, C. W. The quagga mussel invades the Lake Superior Basin. J. Great Lakes Res. 34, 342–350 (2008).Article 

    Google Scholar 
    49.Mahon, A. R. et al. Molecular detection of invasive species in heterogeneous mixtures using a microfluidic carbon nanotube platform. PLoS ONE 6, 1–5 (2011).Article 
    CAS 

    Google Scholar 
    50.PrimerExplorer. LAMP Primer Designing Software (Fujitsu Ltd, Tokyo, 2005).
    Google Scholar 
    51.Untergasser, A. et al. Primer3-new capabilities and interfaces. Nucleic Acids Res. 40, e115–e115 (2012).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    52.Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar  More

  • in

    Engineering complex communities by directed evolution

    1.Mueller, U. G. & Sachs, J. L. Engineering microbiomes to improve plant and animal health. Trends Microbiol. 23, 606–617 (2015).CAS 
    Article 

    Google Scholar 
    2.Gilbert, E. S., Walker, A. W. & Keasling, J. D. A constructed microbial consortium for biodegradation of the organophosphorus insecticide parathion. Appl. Microbiol. Biotechnol. 61, 77–81 (2003).CAS 
    Article 

    Google Scholar 
    3.Yoshida, S., Ogawa, N., Fujii, T. & Tsushima, S. Enhanced biofilm formation and 3-chlorobenzoate degrading activity by the bacterial consortium of Burkholderia sp. NK8 and Pseudomonas aeruginosa PAO1. J. Appl. Microbiol. 106, 790–800 (2009).CAS 
    Article 

    Google Scholar 
    4.Piccardi, P., Vessman, B. & Mitri, S. Toxicity drives facilitation between 4 bacterial species. Proc. Natl Acad. Sci. USA 116, 15979–15984 (2019).CAS 
    Article 

    Google Scholar 
    5.Herrera Paredes, S. et al. Design of synthetic bacterial communities for predictable plant phenotypes. PLoS Biol. 16, e2003962 (2018).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    6.Minty, J. J. et al. Design and characterization of synthetic fungal-bacterial consortia for direct production of isobutanol from cellulosic biomass. Proc. Natl Acad. Sci. USA 110, 14592–14597 (2013).CAS 
    Article 

    Google Scholar 
    7.Jiang, Y., Dong, W., Xin, F. & Jiang, M. Designing synthetic microbial consortia for biofuel production. Trends Biotechnol. 38, 828–831 (2020).CAS 
    Article 

    Google Scholar 
    8.Eng, A. & Borenstein, E. Microbial community design: methods, applications, and opportunities. Curr. Opin. Biotechnol. 58, 117–128 (2019).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    9.Fredrickson, J. K. Ecological communities by design. Science 348, 1425–1427 (2015).CAS 
    Article 

    Google Scholar 
    10.Sanchez-Gorostiaga, A., Bajić, D., Osborne, M. L., Poyatos, J. F. & Sanchez, A. High-order interactions distort the functional landscape of microbial consortia. PLoS Biol. 17, e3000550 (2019).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    11.Senay, Y., John, G., Knutie, S. A. & Brandon Ogbunugafor, C. Deconstructing higher-order interactions in the microbiota: a theoretical examination. Preprint at bioRxiv https://doi.org/10.1101/647156 (2019).12.Gould, A. L. et al. Microbiome interactions shape host fitness. Proc. Natl Acad. Sci. USA 115, E11951–E11960 (2018).CAS 
    Article 

    Google Scholar 
    13.Mickalide, H. & Kuehn, S. Higher-order interaction between species inhibits bacterial invasion of a phototroph-predator microbial community. Cell Syst. 9, 521–533.e10 (2019).CAS 
    Article 

    Google Scholar 
    14.Sanchez, A. Defining higher-order interactions in synthetic ecology: lessons from physics and quantitative genetics. Cell Syst. 9, 519–520 (2019).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    15.Guo, X. & Boedicker, J. Q. The contribution of high-order metabolic interactions to the global activity of a four-species microbial community. PLoS Comput. Biol. 12, e1005079 (2016).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    16.Sundarraman, D. et al. Higher-order interactions dampen pairwise competition in the zebrafish gut microbiome. mBio 11, e01667-20 (2020).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    17.Goldman, R. P. & Brown, S. P. Making sense of microbial consortia using ecology and evolution. Trends Biotechnol. 27, 3–4 (2009).CAS 
    Article 

    Google Scholar 
    18.Brenner, K., You, L. & Arnold, F. H. Response to Goldman and Brown: Making sense of microbial consortia using ecology and evolution. Trends Biotechnol. 27, 4 (2009).CAS 
    Article 

    Google Scholar 
    19.Escalante, A. E., Rebolleda-Gómez, M., Benítez, M. & Travisano, M. Ecological perspectives on synthetic biology: insights from microbial population biology. Front. Microbiol. 6, 143 (2015).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    20.Gilmore, S. P. et al. Top-down enrichment guides in formation of synthetic microbial consortia for biomass degradation. ACS Synth. Biol. 8, 2174–2185 (2019).CAS 
    Article 

    Google Scholar 
    21.Cortes-Tolalpa, L., Jiménez, D. J., de Lima Brossi, M. J., Salles, J. F. & van Elsas, J. D. Different inocula produce distinctive microbial consortia with similar lignocellulose degradation capacity. Appl. Microbiol. Biotechnol. https://doi.org/10.1007/s00253-016-7516-6 (2016).22.Lee, D.-J., Show, K.-Y. & Wang, A. Unconventional approaches to isolation and enrichment of functional microbial consortium – a review. Bioresour. Technol. 136, 697–706 (2013).CAS 
    Article 

    Google Scholar 
    23.Lazuka, A., Auer, L., O’Donohue, M. & Hernandez-Raquet, G. Anaerobic lignocellulolytic microbial consortium derived from termite gut: enrichment, lignocellulose degradation and community dynamics. Biotechnol. Biofuels 11, 284 (2018).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    24.Puentes-Téllez, P. E. & Falcao Salles, J. Construction of effective minimal active microbial consortia for lignocellulose degradation. Microb. Ecol. 76, 419–429 (2018).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    25.He, X., McLean, J. S., Guo, L., Lux, R. & Shi, W. The social structure of microbial community involved in colonization resistance. ISME J. 8, 564–574 (2014).Article 

    Google Scholar 
    26.Jung, J., Philippot, L. & Park, W. Metagenomic and functional analyses of the consequences of reduction of bacterial diversity on soil functions and bioremediation in diesel-contaminated microcosms. Sci. Rep. 6, 23012 (2016).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    27.Franklin, R. B. & Mills, A. L. Structural and functional responses of a sewage microbial community to dilution-induced reductions in diversity. Microb. Ecol. 52, 280–288 (2006).Article 

    Google Scholar 
    28.Kang, D. et al. Enrichment and characterization of an environmental microbial consortium displaying efficient keratinolytic activity. Bioresour. Technol. 270, 303–310 (2018).CAS 
    Article 

    Google Scholar 
    29.Goodnight, C. J. Evolution in metacommunities. Phil. Trans. R. Soc. B 366, 1401–1409 (2011).Article 

    Google Scholar 
    30.Swenson, W., Wilson, D. S. & Elias, R. Artificial ecosystem selection. Proc. Natl Acad. Sci. USA 97, 9110–9114 (2000).CAS 
    Article 

    Google Scholar 
    31.Jochum, M. D., McWilliams, K. L., Pierson, E. A. & Jo, Y.-K. Host-mediated microbiome engineering (HMME) of drought tolerance in the wheat rhizosphere. PLoS ONE 14, e0225933 (2019).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    32.Mueller, U. G. et al. Artificial microbiome-selection to engineer microbiomes that confer salt-tolerance to plants. Preprint at bioRxiv https://doi.org/10.1101/081521 (2016).33.Panke-Buisse, K., Poole, A. C., Goodrich, J. K., Ley, R. E. & Kao-Kniffin, J. Selection on soil microbiomes reveals reproducible impacts on plant function. ISME J. 9, 980–989 (2015).CAS 
    Article 

    Google Scholar 
    34.Panke-Buisse, K., Lee, S. & Kao-Kniffin, J. Cultivated sub-populations of soil microbiomes retain early flowering plant trait. Microb. Ecol. https://doi.org/10.1007/s00248-016-0846-1 (2016).35.Arora, J., Mars Brisbin, M. A. & Mikheyev, A. S. Effects of microbial evolution dominate those of experimental host-mediated indirect selection. PeerJ 8, e9350 (2020).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    36.Swenson, W., Arendt, J. & Wilson, D. S. Artificial selection of microbial ecosystems for 3-chloroaniline biodegradation. Environ. Microbiol. 2, 564–571 (2000).CAS 
    Article 

    Google Scholar 
    37.Wright, R. J., Gibson, M. I. & Christie-Oleza, J. A. Understanding microbial community dynamics to improve optimal microbiome selection. Microbiome 7, 85 (2019).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    38.Blouin, M., Karimi, B., Mathieu, J. & Lerch, T. Z. Levels and limits in artificial selection of communities. Ecol. Lett. 18, 1040–1048 (2015).Article 

    Google Scholar 
    39.Raynaud, T., Devers, M., Spor, A. & Blouin, M. Effect of the reproduction method in an artificial selection experiment at the community level. Front. Ecol. Evol. 7, 416 (2019).Article 

    Google Scholar 
    40.Chang, C.-Y., Osborne, M. L., Bajic, D. & Sanchez, A. Artificially selecting bacterial communities using propagule strategies. Evolution https://doi.org/10.1111/evo.14092 (2020).41.Arias-Sánchez, F. I., Vessman, B. & Mitri, S. Artificially selecting microbial communities: if we can breed dogs, why not microbiomes? PLoS Biol. 17, e3000356 (2019).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    42.Day, M. D., Beck, D. & Foster, J. A. Microbial communities as experimental units. BioScience 61, 398–406 (2011).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    43.Wade, M. J. Group selections among laboratory populations of Tribolium. Proc. Natl Acad. Sci. USA 73, 4604–4607 (1976).CAS 
    Article 

    Google Scholar 
    44.Wade, M. J. An experimental study of group selection. Evolution 31, 134–153 (1977).Article 

    Google Scholar 
    45.Wade, M. J. A critical review of the models of group selection. Q. Rev. Biol. 53, 101–114 (1978).Article 

    Google Scholar 
    46.Goodnight, C. J. Experimental studies of community evolution I: The response to selection at the community level. Evolution 44, 1614–1624 (1990).Article 

    Google Scholar 
    47.Guo, X. & Boedicker, J. High-order interactions between species strongly influence the activity of microbial communities. Biophys. J. 110, 143a (2016).Article 

    Google Scholar 
    48.Stein, R. R. et al. Computer-guided design of optimal microbial consortia for immune system modulation. eLife 7, e30916 (2018).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    49.Arnold, F. H. Innovation by evolution: bringing new chemistry to life (Nobel lecture). Angew. Chem. Int. Ed. 58, 14420–14426 (2019).CAS 
    Article 

    Google Scholar 
    50.Tracewell, C. A. & Arnold, F. H. Directed enzyme evolution: climbing fitness peaks one amino acid at a time. Curr. Opin. Chem. Biol. 13, 3–9 (2009).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    51.Williams, H. T. P. & Lenton, T. M. Artificial selection of simulated microbial ecosystems. Proc. Natl Acad. Sci. USA 104, 8918–8923 (2007).CAS 
    Article 

    Google Scholar 
    52.Williams, H. T. P. & Lenton, T. M. in Advances in Artificial Life ECAL 2007. Lecture Notes in Computer Science, vol. 4648 (eds Almeida e Costa, F. et al.) 93–102 (Springer, 2007).53.Doulcier, G., Lambert, A., De Monte, S. & Rainey, P. B. Eco-evolutionary dynamics of nested Darwinian populations and the emergence of community-level heredity. eLife 9, e53433 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    54.Xie, L., Yuan, A. E. & Shou, W. Simulations reveal challenges to artificial community selection and possible strategies for success. PLoS Biol. 17, e3000295 (2019).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    55.Wilson, D. S. Complex interactions in metacommunities, with implications for biodiversity and higher levels of selection. Ecology 73, 1984–2000 (1992).Article 

    Google Scholar 
    56.Marsland, R. III et al. Available energy fluxes drive a transition in the diversity, stability, and functional structure of microbial communities. PLoS Comput. Biol. 15, e1006793 (2019).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    57.Marsland, R., Cui, W., Goldford, J. & Mehta, P. The Community Simulator: a Python package for microbial ecology. PLoS ONE 15, e0230430 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    58.Marsland, R. III, Cui, W. & Mehta, P. A minimal model for microbial biodiversity can reproduce experimentally observed ecological patterns. Sci. Rep. 10, 3308 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    59.Advani, M., Bunin, G. & Mehta, P. Statistical physics of community ecology: a cavity solution to MacArthur’s consumer resource model. J. Stat. Mech. 2018, 033406 (2018).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    60.Goldford, J. E. et al. Emergent simplicity in microbial community assembly. Science 361, 469–474 (2018).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    61.Lu, N., Sanchez-Gorostiaga, A., Tikhonov, M. & Sanchez, A. Cohesiveness in microbial community coalescence. Preprint at bioRxiv https://doi.org/10.1101/282723 (2018).62.Faith, J. J., Ahern, P. P., Ridaura, V. K., Cheng, J. & Gordon, J. I. Identifying gut microbe-host phenotype relationships using combinatorial communities in gnotobiotic mice. Sci. Transl. Med. 6, 220ra11 (2014).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    63.Estrela, S. et al. Metabolic rules of microbial community assembly. Preprint at bioRxiv https://doi.org/10.1101/2020.03.09.984278 (2020).64.Friedman, J., Higgins, L. M. & Gore, J. Community structure follows simple assembly rules in microbial microcosms. Nat. Ecol. Evol. 1, 0109 (2017).Article 

    Google Scholar 
    65.Venturelli, O. S. et al. Deciphering microbial interactions in synthetic human gut microbiome communities. Mol. Syst. Biol. 14, e8157 (2018).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    66.Hall, B. G. Experimental evolution of a new enzymatic function. II. Evolution of multiple functions for ebg enzyme in E. coli. Genetics 89, 453–465 (1978).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    67.Smith, G. P. & Petrenko, V. A. Phage display. Chem. Rev. 97, 391–410 (1997).CAS 
    Article 

    Google Scholar 
    68.Bloom, J. D. & Arnold, F. H. In the light of directed evolution: pathways of adaptive protein evolution. Proc. Natl Acad. Sci. USA 106, 9995–10000 (2009).CAS 
    Article 

    Google Scholar 
    69.Romero, P. A., Krause, A. & Arnold, F. H. Navigating the protein fitness landscape with Gaussian processes. Proc. Natl Acad. Sci. USA 110, E193–E201 (2013).CAS 
    Article 

    Google Scholar 
    70.Ho, K.-L., Lee, D.-J., Su, A. & Chang, J.-S. Biohydrogen from cellulosic feedstock: dilution-to-stimulation approach. Int. J. Hydrog. Energy 37, 15582–15587 (2012).CAS 
    Article 

    Google Scholar 
    71.Shepherd, E. S., DeLoache, W. C., Pruss, K. M., Whitaker, W. R. & Sonnenburg, J. L. An exclusive metabolic niche enables strain engraftment in the gut microbiota. Nature 557, 434–438 (2018).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    72.Ting, S.-Y. et al. Targeted depletion of bacteria from mixed populations by programmable adhesion with antagonistic competitor cells. Cell Host Microbe https://doi.org/10.1016/j.chom.2020.05.006 (2020).73.Sheth, R. U., Cabral, V., Chen, S. P. & Wang, H. H. Manipulating bacterial communities by in situ microbiome engineering. Trends Genet. 32, 189–200 (2016).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    74.Lemon, K. P., Armitage, G. C., Relman, D. A. & Fischbach, M. A. Microbiota-targeted therapies: an ecological perspective. Sci. Transl. Med. 4, 137rv5 (2012).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    75.Harcombe, W. R. & Bull, J. J. Impact of phages on two-species bacterial communities. Appl. Environ. Microbiol. 71, 5254–5259 (2005).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    76.Chan, B. K. et al. Phage treatment of an aortic graft infected with Pseudomonas aeruginosa. Evol. Med. Public Health 2018, 60–66 (2018).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    77.Rillig, M. C., Tsang, A. & Roy, J. Microbial community coalescence for microbiome engineering. Front. Microbiol. 7, 1967 (2016).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    78.Sierocinski, P. et al. A single community dominates structure and function of a mixture of multiple methanogenic communities. Curr. Biol. 27, 3390–3395.e4 (2017).CAS 
    Article 

    Google Scholar 
    79.Tilman, D. The ecological consequences of changes in biodiversity: a search for general principles. Ecology 80, 1455–1474 (1999).
    Google Scholar 
    80.Shade, A. et al. Fundamentals of microbial community resistance and resilience. Front. Microbiol. 3, 417 (2012).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    81.Kang, D. et al. Construction of simplified microbial consortia to degrade recalcitrant materials based on enrichment and dilution-to-extinction cultures. Front. Microbiol. 10, 3010 (2019).Article 

    Google Scholar 
    82.Zanaroli, G. et al. Characterization of two diesel fuel degrading microbial consortia enriched from a non acclimated, complex source of microorganisms. Microb. Cell Factories 9, 10 (2010).Article 
    CAS 

    Google Scholar 
    83.Peter, H. et al. Function-specific response to depletion of microbial diversity. ISME J. 5, 351–361 (2011).CAS 
    Article 

    Google Scholar 
    84.Pacheco, A. R., Moel, M. & Segrè, D. Costless metabolic secretions as drivers of interspecies interactions in microbial ecosystems. Nat. Commun. 10, 103 (2019).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    85.West, S. A., Griffin, A. S., Gardner, A. & Diggle, S. P. Social evolution theory for microorganisms. Nat. Rev. Microbiol. 4, 597–607 (2006).CAS 
    Article 

    Google Scholar 
    86.Scheuerl, T. et al. Bacterial adaptation is constrained in complex communities. Nat. Commun. 11, 754 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    87.Lewontin, R. C. The units of selection. Annu. Rev. Ecol. Syst. 1, 1–18 (1970).Article 

    Google Scholar 
    88.Marsland, R., Cui, W., Goldford, J. & Mehta, P. The Community Simulator: a Python package for microbial ecology. PLoS ONE 15, e0230430 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    89.Shoemaker, W. R., Locey, K. J. & Lennon, J. T. A macroecological theory of microbial biodiversity. Nat. Ecol. Evol. 1, 0107 (2017).Article 

    Google Scholar 
    90.Degnan, P. H., Taga, M. E. & Goodman, A. L. Vitamin B12 as a modulator of gut microbial ecology. Cell Metab. 20, 769–778 (2014).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    91.Degnan, P. H., Barry, N. A., Mok, K. C., Taga, M. E. & Goodman, A. L. Human gut microbes use multiple transporters to distinguish vitamin B12 analogs and compete in the gut. Cell Host Microbe 15, 47–57 (2014).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar  More

  • in

    Reduced resilience of terrestrial ecosystems locally is not reflected on a global scale

    1.Smol, J. P. et al. Climate-driven regime shifts in the biological communities of arctic lakes. Proc. Natl Acad. Sci. USA 102, 4397–4402 (2005).CAS 
    Article 

    Google Scholar 
    2.Wernberg, T. et al. Climate-driven regime shift of a temperate marine ecosystem. Science 353, 169–172 (2016).CAS 
    Article 

    Google Scholar 
    3.Scheffer, M., Carpenter, S., Foley, J. A., Folke, C. & Walker, B. Catastrophic shifts in ecosystems. Nature 413, 591–596 (2001).CAS 
    Article 

    Google Scholar 
    4.Staver, A. C., Archibald, S. & Levin, S. A. The global extent and determinants of savanna and forest as alternative biome states. Science 334, 230–232 (2011).CAS 
    Article 

    Google Scholar 
    5.Su, H. et al. Long‐term empirical evidence, early warning signals and multiple drivers of regime shifts in a lake ecosystem. J. Ecol. https://doi.org/10.1111/1365-2745.13544 (2020).6.Barnosky, A. D. et al. Approaching a state shift in Earth’s biosphere. Nature 486, 52–58 (2012).CAS 
    Article 

    Google Scholar 
    7.Steffen, W. et al. Trajectories of the Earth system in the Anthropocene. Proc. Natl Acad. Sci. 115, 8252–8259 (2018).CAS 
    Article 

    Google Scholar 
    8.Holling, C. S. Resilience and stability of ecological systems. Ann. Rev. Ecol. Syst. 4, 1–23 (1973).Article 

    Google Scholar 
    9.Ratajczak, Z. et al. Abrupt change in ecological systems: inference and diagnosis. Trends Ecol. Evol. 33, 513–526 (2018).Article 

    Google Scholar 
    10.Pimm, S. L. The complexity and stability of ecosystems. Nature 307, 321–326 (1984).Article 

    Google Scholar 
    11.Holling, C. S. Engineering resilience versus ecological resilience. Eng. Ecol.Constraints 31, 32 (1996).
    Google Scholar 
    12.Li, X. et al. Temporal trade-off between gymnosperm resistance and resilience increases forest sensitivity to extreme drought. Nat. Ecol. Evol. 4, 1075–1083 (2020).Article 

    Google Scholar 
    13.Carpenter, S. R. & Brock, W. A. Rising variance: a leading indicator of ecological transition. Ecol. Lett. 9, 311–318 (2006).CAS 
    Article 

    Google Scholar 
    14.Dakos, V. et al. Slowing down as an early warning signal for abrupt climate change. Proc. Natl Acad. Sci. USA 105, 14308–14312 (2008).CAS 
    Article 

    Google Scholar 
    15.Guttal, V. & Jayaprakash, C. Changing skewness: an early warning signal of regime shifts in ecosystems. Ecol. Lett. 11, 450–460 (2008).Article 

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

    Google Scholar 
    17.Drake, J. M. & Griffen, B. D. Early warning signals of extinction in deteriorating environments. Nature 467, 456 (2010).CAS 
    Article 

    Google Scholar 
    18.Wang, R. et al. Flickering gives early warning signals of a critical transition to a eutrophic lake state. Nature 492, 419–422 (2012).Article 
    CAS 

    Google Scholar 
    19.Clements, C. F. & Ozgul, A. Including trait-based early warning signals helps predict population collapse. Nat. Commun. 7, 10984 (2016).CAS 
    Article 

    Google Scholar 
    20.Chevalier, M. & Grenouillet, G. Global assessment of early warning signs that temperature could undergo regime shifts. Sci. Rep. 8, 10058 (2018).Article 
    CAS 

    Google Scholar 
    21.Cole, L. E., Bhagwat, S. A. & Willis, K. J. Recovery and resilience of tropical forests after disturbance. Nat. Commun. 5, 3906 (2014).CAS 
    Article 

    Google Scholar 
    22.Willis, K. J., Jeffers, E. S. & Tovar, C. What makes a terrestrial ecosystem resilient? Science 359, 988–989 (2018).CAS 
    Article 

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

    Google Scholar 
    24.Seddon, A. W., Macias-Fauria, M., Long, P. R., Benz, D. & Willis, K. J. Sensitivity of global terrestrial ecosystems to climate variability. Nature 531, 229–232 (2016).CAS 
    Article 

    Google Scholar 
    25.Ehleringer, J. R., Cerling, T. E. & Helliker, B. R. C4 photosynthesis, atmospheric CO2, and climate. Oecologia 112, 285–299 (1997).Article 

    Google Scholar 
    26.Higgins, S. I. & Scheiter, S. Atmospheric CO2 forces abrupt vegetation shifts locally, but not globally. Nature 488, 209 (2012).CAS 
    Article 

    Google Scholar 
    27.Holmgren, M., Hirota, M., Van Nes, E. H. & Scheffer, M. Effects of interannual climate variability on tropical tree cover. Nat. Clim. Chang. 3, 755–758 (2013).Article 

    Google Scholar 
    28.Thornton, P. K., Ericksen, P. J., Herrero, M. & Challinor, A. J. Climate variability and vulnerability to climate change: a review. Glob. Change Biol. 20, 3313–3328 (2014).Article 

    Google Scholar 
    29.Ray, D. K., Gerber, J. S., MacDonald, G. K. & West, P. C. Climate variation explains a third of global crop yield variability. Nat. Commun. 6, 5989 (2015).CAS 
    Article 

    Google Scholar 
    30.Jha, S., Das, J. & Goyal, M. K. Assessment of risk and resilience of terrestrial ecosystem productivity under the influence of extreme climatic conditions over India. Sci. Rep. 9, 18923 (2019).CAS 
    Article 

    Google Scholar 
    31.Li, D., Wu, S., Liu, L., Zhang, Y. & Li, S. Vulnerability of the global terrestrial ecosystems to climate change. Glob. Change Biol. 24, 4095–4106 (2018).Article 

    Google Scholar 
    32.Gonzalez, P., Neilson, R. P., Lenihan, J. M. & Drapek, R. J. Global patterns in the vulnerability of ecosystems to vegetation shifts due to climate change. Glob. Ecol. Biogeogr. 19, 755–768 (2010).Article 

    Google Scholar 
    33.Wang, S. & Loreau, M. Ecosystem stability in space: α, β and γ variability. Ecol. Lett. 17, 891–901 (2014).Article 

    Google Scholar 
    34.Stenseth, N. C. et al. The effect of climatic forcing on population synchrony and genetic structuring of the Canadian lynx. Proc. Natl Acad. Sci. USA 101, 6056–6061 (2004).CAS 
    Article 

    Google Scholar 
    35.Koenig, W. D. & Liebhold, A. M. Temporally increasing spatial synchrony of North American temperature and bird populations. Nat. Clim. Chang. 6, 614–617 (2016).Article 

    Google Scholar 
    36.Sheppard, L. W., Bell, J. R., Harrington, R. & Reuman, D. C. Changes in large-scale climate alter spatial synchrony of aphid pests. Nat. Clim. Chang. 6, 610–613 (2016).Article 

    Google Scholar 
    37.Dakos, V., van Nes, E. H., Donangelo, R., Fort, H. & Scheffer, M. Spatial correlation as leading indicator of catastrophic shifts. Theor. Ecol. 3, 163–174 (2010).Article 

    Google Scholar 
    38.Paruelo, J. M., Epstein, H. E., Lauenroth, W. K. & Burke, I. C. ANPP estimates from NDVI for the central grassland region of the United States. Ecology 78, 953–958 (1997).Article 

    Google Scholar 
    39.Piao, S., Fang, J., Zhou, L., Tan, K. & Tao, S. Changes in biomass carbon stocks in China’s grasslands between 1982 and 1999. Global Biogeochem. Cycles 21, 2 (2007).
    Google Scholar 
    40.Maurer, G. E., Hallmark, A. J., Brown, R. F., Sala, O. E. & Collins, S. L. Sensitivity of primary production to precipitation across the United States. Ecol. Lett. 23, 527–536 (2020).Article 

    Google Scholar 
    41.Brown, J. H. & Kodric-Brown, A. Turnover rates in insular biogeography: effect of immigration on extinction. Ecology 58, 445–449 (1977).Article 

    Google Scholar 
    42.Earn, D. J., Levin, S. A. & Rohani, P. Coherence and conservation. Science 290, 1360–1364 (2000).CAS 
    Article 

    Google Scholar 
    43.Hodgson, D., McDonald, J. L. & Hosken, D. J. What do you mean,‘resilient’? Trends Ecol. Evol. 30, 503–506 (2015).Article 

    Google Scholar 
    44.Seidl, R. et al. Forest disturbances under climate change. Nat. Clim. Chang. 7, 395–402 (2017).Article 

    Google Scholar 
    45.Bernstein, L. et al. IPCC, 2007: Climate Change 2007: Synthesis Report. (IPCC, Geneva, 2008)46.Myers-Smith, I. H. et al. Shrub expansion in tundra ecosystems: dynamics, impacts and research priorities. Environ. Res. Lett. 6, 045509 (2011).Article 

    Google Scholar 
    47.Myers-Smith, I. H. et al. Climate sensitivity of shrub growth across the tundra biome. Nat. Clim. Chang. 5, 887–891 (2015).Article 

    Google Scholar 
    48.Thompson, I., Mackey, B., McNulty, S. & Mosseler, A. Forest resilience, biodiversity, and climate change. In Secretariat of the Convention on Biological Diversity, Montreal. Technical Series 43, 1–67 (2009).
    Google Scholar 
    49.Carpenter, S. R. et al. Early warnings of regime shifts: a whole-ecosystem experiment. Science 332, 1079–1082 (2011).CAS 
    Article 

    Google Scholar 
    50.Gsell, A. S. et al. Evaluating early-warning indicators of critical transitions in natural aquatic ecosystems. Proc. Natl Acad. Sci. USA 113, E8089–E8095 (2016).CAS 
    Article 

    Google Scholar 
    51.Clements, C. F., Blanchard, J. L., Nash, K. L., Hindell, M. A. & Ozgul, A. Body size shifts and early warning signals precede the historic collapse of whale stocks. Nat. Ecol. Evol. 1, 0188 (2017).Article 

    Google Scholar 
    52.Dakos, V., Carpenter, S. R., van Nes, E. H. & Scheffer, M. Resilience indicators: prospects and limitations for early warnings of regime shifts. Philos. Trans. R. Soc. B, Biol. Sci. 370, 20130263 (2015).Article 

    Google Scholar 
    53.Zemp, D. C. et al. Self-amplified Amazon forest loss due to vegetation-atmosphere feedbacks. Nat. Commun. 8, 14681 (2017).CAS 
    Article 

    Google Scholar 
    54.Staal, A. et al. Forest-rainfall cascades buffer against drought across the Amazon. Nat. Clim. Chang. 8, 539–543 (2018).Article 

    Google Scholar 
    55.Poorter, L. et al. Biomass resilience of Neotropical secondary forests. Nature 530, 211–214 (2016).CAS 
    Article 

    Google Scholar 
    56.Locosselli, G. M. et al. Global tree-ring analysis reveals rapid decrease in tropical tree longevity with temperature. Proc. Natl Acad. Sci. USA 117, 33358–33364 (2020).CAS 
    Article 

    Google Scholar 
    57.Ruiz-Pérez, G. & Vico, G. Effects of temperature and water availability on Northern European boreal forests. Front. For. Glob.Change 3, 34 (2020).Article 

    Google Scholar 
    58.Kitzberger, T., Aráoz, E., Gowda, J. H., Mermoz, M. & Morales, J. M. Decreases in fire spread probability with forest age promotes alternative community states, reduced resilience to climate variability and large fire regime shifts. Ecosystems 15, 97–112 (2012).Article 

    Google Scholar 
    59.Scheffer, M., Hirota, M., Holmgren, M., Van Nes, E. H. & Chapin, F. S. Thresholds for boreal biome transitions. Proc. Natl Acad. Sci. USA 109, 21384–21389 (2012).CAS 
    Article 

    Google Scholar 
    60.Newbold, T. et al. Climate and land-use change homogenise terrestrial biodiversity, with consequences for ecosystem functioning and human well-being. Emerg. Top. Life Sci. 3, 207–219 (2019).Article 

    Google Scholar 
    61.Senior, R. A., Hill, J. K., González del Pliego, P., Goode, L. K. & Edwards, D. P. A pantropical analysis of the impacts of forest degradation and conversion on local temperature. Ecol. Evol. 7, 7897–7908 (2017).Article 

    Google Scholar 
    62.Wang, S. et al. An invariability-area relationship sheds new light on the spatial scaling of ecological stability. Nat. Commun. 8, 1–8 (2017).Article 
    CAS 

    Google Scholar 
    63.Mehrabi, Z. & Ramankutty, N. Synchronized failure of global crop production. Nat. Ecol. Evol. 3, 780–786 (2019).Article 

    Google Scholar 
    64.Post, E. & Forchhammer, M. C. Spatial synchrony of local populations has increased in association with the recent Northern Hemisphere climate trend. Proc. Natl Acad. Sci. 101, 9286–9290 (2004).CAS 
    Article 

    Google Scholar 
    65.Ripa, J. Analysing the Moran effect and dispersal: their significance and interaction in synchronous population dynamics. Oikos 89, 175–187 (2000).Article 

    Google Scholar 
    66.Peterson, G., Allen, C. R. & Holling, C. S. Ecological resilience, biodiversity, and scale. Ecosystems 1, 6–18 (1998).Article 

    Google Scholar 
    67.Wang, S. & Loreau, M. Biodiversity and ecosystem stability across scales in metacommunities. Ecol. Lett. 19, 510–518 (2016).Article 

    Google Scholar 
    68.Dakos, V. et al. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PloS ONE 7, e41010 (2012).CAS 
    Article 

    Google Scholar 
    69.R core team. R: a language and environment for statistical computing. R Foundation for Statistical Computing https://www.R-project.org/ (2019).70.Bivand, R., Keitt, T. & Rowlingson, B. rgdal: bindings for the ‘Geospatial’ Data Abstraction Library. R package version 1.5-16 https://CRAN.R-project.org/package=rgdal (2020).71.Tucker, C. J. et al. An extended AVHRR 8‐km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data. Int. J. Remote Sens. 26, 4485–4498 (2005).Article 

    Google Scholar 
    72.Pinzon, J. E. & Tucker, C. J. A non-stationary 1981–2012 AVHRR NDVI3g time series. Remote Sens. 6, 6929–6960 (2014).Article 

    Google Scholar 
    73.Holben, B. N. Characteristics of maximum-value composite images from temporal AVHRR data. Int. J. Remote Sens. 7, 1417–1434 (1986).Article 

    Google Scholar 
    74.Piao, S. et al. Changes in vegetation net primary productivity from 1982 to 1999 in China. Global Biogeochem. Cycles 19, 2 (2005).Article 
    CAS 

    Google Scholar 
    75.Olson, D. M. et al. Terrestrial ecoregions of the world: a new map of life on Earth. BioScience 51, 933–938 (2001).Article 

    Google Scholar 
    76.Harris, I., Osborn, T. J., Jones, P. & Lister, D. Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. Sci. Data 7, 1–18. (2020).Article 

    Google Scholar 
    77.Mitchell, A. The ESRI Guide to GIS Analysis: Spatial Measurements and Statistics (Environmental System Research Institute Press, 2005).78.Fang, J., Piao, S., He, J. & Ma, W. Increasing terrestrial vegetation activity in China, 1982–1999. Sci. China C Life Sci. 47, 229–240 (2004).
    Google Scholar 
    79.Peng, S. et al. Recent change of vegetation growth trend in China. Environ. Res. Lett. 6, 044027 (2011).Article 

    Google Scholar 
    80.Thenkabail, P. S. & Lyon, J. G. Hyperspectral Remote Sensing of Vegetation (CRC press, 2016).81.Feng, Y. et al. Changes in the trends of vegetation net primary productivity in China between 1982 and 2015. Environ. Res. Lett. 14, 124009 (2019).Article 

    Google Scholar 
    82.He, H. et al. Altered trends in carbon uptake in China’s terrestrial ecosystems under the enhanced summer monsoon and warming hiatus. Natl Sci. Rev. 6, 505–514 (2019).CAS 
    Article 

    Google Scholar  More

  • in

    Community composition of microbial microcosms follows simple assembly rules at evolutionary timescales

    Strains and mediaThe set of 16 strains used in this experiment contains environmental isolates along with strains from the ATCC collection (Supplementary Table 1). The strains were chosen based on two criteria: a distinct colony morphology that would enable visual identification when plated on an NB agar plate; and ability to coexist for ~60 generations with at least two other strains in our collection.All cultures were grown in M9 minimal salts media containing 1X M9 salts, 2 mM MgSO4, 0.1 mM CaCl2, 1X trace metal solution (Teknova), supplemented with 3 mM galacturonic acid (Sigma), 6.1 mM Serine (Sigma), and 9.1 mM sodium acetate as carbon sources, which correspond to 16.67 mM carbon atoms for each compound and 50 mM overall. We chose a combination of carbon sources representing three chemical groups—a carbohydrate, an amino acid, and a carboxylic acid—in order to promote the survival and coexistence of a diverse set of species. The media was prepared on the day of each transfer. A carbon source mixture was prepared ahead at 10X, and was kept in aliquots at 4 °C for up to four weeks.Evolution experimentFrozen stocks of individual species were streaked out on nutrient agar Petri plates and grown at 28 °C. After 48 h single colonies were picked and inoculated into 15 ml falcon tubes containing 3 ml nutrient broth (5 g/L peptone BD difco, BD Bioscience; 3 g/L yeast extract BD difco, BD Bioscience), and were grown overnight at 28 °C shaken at 250 rpm. Initial mixtures were prepared by diluting each species separately to an OD of ({10}^{-2}) and mixing the normalized cultures at equal volumes. OD measurements were done using a Epoch2 microplate reader (BioTek) and were recorded using the Gen5 v3.09 software (BioTek). After mixing, the cocultures were aliquoted to replicates and further diluted to a final OD of ({10}^{-4}), at which the evolutionary experiment was initialized. The number of replicates for each community varied between 3 and 18 (Supplementary Data 1).Communities were grown in 96-well plates containing 200 µl M9 at 28°C and were shaken at 900 rpm. Every 48 h cultures were diluted by a factor of 1500 into fresh M9 media, and OD600 was measured. For this dilution factor, each cycle corresponds to ~10.5 generations. As 1 OD600 ~ of ({10}^{9}) C.F.U/ml, and communities reached ~ 0.5 OD600 and were grown in 200 µl and was diluted by 1500, ~({10}^{5}) cells were transferred each dilution. To avoid cross contaminations, cultures were grown in a checkerboard formation, meaning that each community was surrounded by wells containing media but no bacteria.At transfers 0, 2, 5, 7, 10, 14, 19, 30, and 38 community composition was measured by plating on nutrient agar plates (5 g/L peptone BD difco, BD Bioscience; 3 g/L yeast extract BD difco, BD Bioscience, 15 g/L agar Bacto, BD Bioscience) and counting colonies. For that, the cultures were diluted to an OD of (2.4* {10}^{-8})− (1* {10}^{-8}) and 100 µl of the diluted culture was plated on NB plates and spread using glass beads. Plates were incubated at 28 °C for 48 h and colonies were counted manually. The distribution of the number of colonies counted at each plate to infer community composition is found in Supplementary Fig. 11.We chose the communities based on a preliminary experiment that was conducted by the same protocol for six transfers. In this experiment, 114 of 171 possible pairs of a set of 19 strains (3 strains were not included in the evolution experiment) were cocultured. Pairs that had coexisted for the duration of this experiment, and were confidently distinguishable by colony morphology, and trios that are composed of these pairs, were used for the coevolutionary experiment. We started the evolutionary experiment with 51 pairs and 51 trios, and removed communities that did not coexist for the first ~70 from the final analysis. If a replicate was suspected to be contaminated it was also excluded from further analysis.Ecological experimentsWe supplemented the data of the evolutionary experiment with two ecological competition experiments with the same experimental condition. In order to assess whether communities typically reach an ecological equilibrium within ~50–70 generations (Supplementary Fig. 3), we cultured eight of the pairs that were used in the evolutionary experiment. This experiment was initiated in the same way as the evolutionary experiment, only that after the species’ starters were normalized they were inoculated at the varying initial fractions – 9:1, 5:5, 1:9. Because the normalization depended on optical density, there is a variation in the actual initial fractions between different pairs. Community composition was then measured on six transfers during this experiment: 0, 1, 2, 4, 5, and 6.In order to assess whether changes in composition are due to heritable changes in species’ phenotypes, we used strains that were re-isolated from 31 evolved pairs, and 13 pairs of ancestral strains (Supplementary Figs. 6, 7). Strains were replicated from glycerol stocks into the experimental media and grown for 24 h. The starters were normalized to initiate the competition assay at ({rm{OD}}={10}^{-4}) in fresh M9 media. Species were mixed at equal volume and were propagated for five cycles. community composition was measured at initial conditions, and at the end of the final cycle (5).Quantification of repeatabilityIn order to quantify the qualitative repeatability of different replicate communities we first identified which species was the maximally increasing member at each replicate, that is, which species had increased its abundance by the largest factor between generation 70 and 400. Then, we quantified the frequency of the replicates that had the same maximally increasing member for each community. This measure always produces a value between 1 and 1/n where n is the number of species in the community. We checked the distribution of the repeatability scores against the null hypothesis that the factor by which a species’ abundance increases during evolution is independent of the species or the community. For this, we shuffled the factor of change in relative abundance across all samples, for pairs and trios separately, and quantified the new repeatability scores of the shuffled data. Data of the null hypothesis were generated over 2000 times, and the p value was given by the probability to get a mean equal or above the real data mean.We used the average Euclidean distance of replicates from the median replicate in order to quantify the variability between replicate communities. In order to check whether the distribution of variabilities is similar to what can be expected of random communities, in which each species in the community is just as likely to have any relative abundance, we replaced the real fractions with fractions drawn from a uniform Dirichlet distribution with (underline{{boldsymbol{alpha }}}=underline{1}). We then checked the statistical difference between the two distributions using one-sided Mann–Whitney U test.Trio composition predictionsWe used the formerly established method for predicting the composition of trios from the composition of pairs that was developed by Abreu et al.14 In this approach the fraction of a species when grown in a multispecies community is predicted as the weighted geometric mean of the fraction of the species in all pairwise cultures. The accuracy of the predictions was measured as the Euclidean distance between the prediction and the mean composition of the observed trio, normalized to the largest possible distance between each two communities, (sqrt{n}), where n is the number of species.We used the factors by which species increased their abundance during coevolution in pairs (between generations ~70 and ~400) to predict which species would increase by the largest factor in trios. The maximally increasing member in a given community was assigned to be the one that was the maximally increasing member in the most replicates of that community. If the same species was the maximally increasing member in both pairs it was a member of, then this species was predicted to be the maximally increasing member of the trio. If in every pair a different species was the maximally increasing member, then we predicted that the maximally increasing member of the trio would be the one with the highest mean increase. Only two trios had such transient topology, where in each pair a different species increases, thus we are unable to determine the general utility of the latter approach.Re-isolationEach ~50 generations all communities were frozen at −80 °C with 50% glycerol in a 96-deep well plate. In order to re-isolate strains, stocks were inoculated to a 96-well plate containing the experimental media using a 96-pin replicator, and grown for 24 h at 28 °C. After growth, cultures were diluted by a factor of (2.4* {10}^{-8}) and 100 µl were spread on a nutrient agar plate using glass beads. Plates were kept at room temperature for at least two days and no longer than a week before re-isolations. 5-15 colonies of each strain were picked using a sterile toothpick, and pooled together into 200 µl M9. Re-isolated strains were incubated at 28 °C and shaken at 900 rpm for 24 h and kept in 50% glycerol stock at −80 °C until further use.Growth rates and carrying capacities of individually evolved strainsRe-isolated strains were replicated from glycerol stocks into the experimental media and grown for 24 h. The starters were normalized to initiate the growth assay at (OD={10}^{-4}) in fresh M9 media. The optical density was measured in two automated plate readers simultaneously, Epoch2 microplate reader (BioTek) and Synergy microplate reader (BioTek), and was recorded using Gen5 v3.09 software (BioTek). Plates were incubated at 28 °C with a 1 °C gradient to avoid condensation on the lid, and were shaken at 250 cpm. OD was measured every 10 min. Each strain was measured in four technical replicates, evenly distributed between the two plates, and 2–3 evolutionary replicates were measured for each species (replicates that evolved separately for the duration of the experiment). Growth rates were quantified as the number of divisions it takes a strain to grow from the initial OD of ({10}^{-4}) to an OD of (8* {10}^{-2}) (({log }_{2}frac{0.08}{{10}^{-4}})) divided by the time it took the strain to reach this OD. This measure gives the average doubling time during the initial growth and also accounts for the lag times of the strain. The growth rates of evolutionary replicates were averaged after averaging technical replicates.Carrying capacity was defined as the OD a monoculture reached at the end of each growth cycle of the evolutionary experiment averaged across replicates. These measurements were done in an Epoch2 microplate reader (BioTek). In order to reduce noise, the trajectories of OD measurements were smoothed for each well using moving mean with an averaging window of three.Carrying capacities of coevolved strainsRe-isolated strains were replicated from glycerol stocks into the experimental media and grown for 48 h in M9 media at 28 °C. Cultures were then diluted by 1500 into 3 technical replicates in fresh M9-media, and were given another 48 h to reach carrying capacity. The strains used in this experiment were isolated from 1-3 evolutionary replicates (Supplementary Data 2).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

  • in

    The epidemicity index of recurrent SARS-CoV-2 infections

    Data and data processingThe modeling tools described in the following sections are applied to the Italian COVID-19 epidemic at the scale of second-level administrative divisions, i.e., provinces and metropolitan cities (as of 2020, 107 spatial units). Official data about resident population at the provincial level are produced yearly by the Italian National Institute of Statistics (Istituto Nazionale di Statistica, ISTAT; data available at http://dati.istat.it/Index.aspx?QueryId=18460). The January 2019 update has been used to inform the spatial distribution of the population.The data to quantify nation-wide human mobility prior to the pandemic come from ISTAT (specifically, from the 2011 national census; data available online at https://www.istat.it/it/archivio/139381). Mobility fluxes, mostly reflecting commuting patterns related to work and study purposes, are provided at the scale of third-level administrative units (municipalities)53,54. These fluxes were upscaled to the provincial level following the administrative divisions of 2019, and used to evaluate the fraction pi of mobile people and the fraction qij of mobile people between i and all other administrative units j (see Supplementary Material in Gatto et al.7).Airport traffic data for year 2019, used to inform the simulation shown in Fig. 4c, d, are from the Italian Airports Association (Assaeroporti; data available at http://assaeroporti.com/statistiche_201912/). Note that airports have been assigned to the main Metropolitan Area they serve, rather than to the province where they are geographically located (e.g., Malpensa Airport has been assigned to the Metropolitan City of Milano, rather than to the neighboring Varese province, where it actually lies).Model parameters are taken from a paper by Bertuzzo et al.14, where they were inferred in a Bayesian framework on the basis of the official epidemiological bulletins released daily by Dipartimento della Protezione Civile55 (data available online at https://github.com/pcm-dpc/COVID-19) and the bulletins of Epicentro, at ISS51,56. The parameters estimated for the initial phase of the Italian COVID-19 epidemic14, during which SARS-CoV-2 was spreading unnoticed in the population, reflect a situation of unperturbed social mixing and human mobility, absent any effort devoted to disease control. This parameterization, in which all parameters (including the transmission rates) are spatially homogeneous, is reported in Table 2 and has been used to produce all the results presented in the main text, except for those of Fig. 6. In this case, to account for the containment measures put in place by the Italian authorities and their effects on transmission rates and mobility patterns during the first months of the pandemic, a time-varying parameterization14 for the period February 24 to May 1, 2020 has been used. In this parameterization, the transmission rates were allowed to take different values over different time windows, corresponding to the timing of the implementation of the main nation-wide restrictions, or lifting thereof. Specifically, the effect of the containment measures was parameterized by assuming that the transmission parameters had a sharp decrease after the containment measures announced at the end of February and the beginning of March, and that they were further reduced in the following weeks as the country was effectively entering full lockdown. As a by-product, these time-varying transmission rates can also at least partially account for seasonal effects on disease transmission. Due to the emerging nature of the pathogen, seasonality has not been given further consideration in this work; however, it may become a key component of future modeling efforts aimed at studying post-pandemic SARS-CoV-2 transmission dynamics3, i.e., if/when the pathogen establishes as endemic. Spatial connectivity too was modified with respect to the baseline scenario to reflect the disruption of mobility patterns induced by the pandemic and the associated containment measures14. Specifically, between-province mobility was progressively reduced as the epidemic unfolded according to estimates obtained through mobility data from mobile applications53,57.Spatially explicit SEPIAR with distributed controlsWe consider a set of n communities connected by human mobility fluxes. In each community, the human population is subdivided according to infection status into the epidemiological compartments of susceptible, exposed (latently infected), post-latent (incubating infectious, also termed pre-symptomatic7), symptomatic infectious, asymptomatic infectious (including paucisymptomatic), and recovered individuals. The present model utilizes previous work aimed to describe the first wave of COVID-19 infections7,14. In particular, it allows us to account for three widely adopted types of containment measures: reduction of local transmission (as a result of the use of personal protections, social distancing, and local mobility restriction), travel restriction, and isolation of infected individuals. To describe the effects of isolation, each infected compartment (exposed, post-latent, symptomatic and asymptomatic) is actually split into two, which allows keeping track of the abundances of infected individuals who are still in the community vs. those who are removed from it (i.e., either in isolation at a hospital, if symptomatic, or quarantined at home, if exposed, post-latent, or asymptomatic). The state variables of the model are summarized in Table 1. Supplementary Figure 1 recapitulates the structure of the model.COVID-19 transmission dynamics are thus described by the following set of ordinary differential equations:$${dot{S}}_{i} =mu ({N}_{i}-{S}_{i})-{lambda }_{i}{S}_{i}\ {dot{E}}_{i} ={lambda }_{i}{S}_{i}-(mu +{delta }^{E}+{chi }_{i}^{E}){E}_{i}\ {dot{P}}_{i} ={delta }^{E}{E}_{i}-(mu +{delta }^{P}+{chi }_{i}^{P}){P}_{i}\ {dot{I}}_{i} =sigma {delta }^{P}{P}_{i}-(mu +alpha +{gamma }^{I}+eta +{chi }_{i}^{I}){I}_{i}\ {dot{A}}_{i} =(1-sigma ){delta }^{P}{P}_{i}-(mu +{gamma }^{A}+{chi }_{i}^{A}){A}_{i}\ {dot{E}}_{i}^{{rm{q}}} ={chi }_{i}^{E}{E}_{i}-(mu +{delta }^{E}){E}_{i}^{{rm{q}}}\ {dot{P}}_{i}^{{rm{q}}} ={chi }_{i}^{P}{P}_{i}+{delta }^{E}{E}_{i}^{{rm{q}}}-(mu +{delta }^{P}){P}_{i}^{{rm{q}}}\ {dot{I}}_{i}^{{rm{h}}} =(eta +{chi }_{i}^{I}){I}_{i}+sigma {delta }^{P}{P}_{i}^{{rm{q}}}-(mu +alpha +{gamma }^{I}){I}_{i}^{{rm{h}}}\ {dot{A}}_{i}^{{rm{q}}} ={chi }_{i}^{A}{A}_{i}+(1-sigma ){delta }^{P}{P}_{i}^{{rm{q}}}-(mu +{gamma }^{A}){A}_{i}^{{rm{q}}}\ {dot{R}}_{i} ={gamma }^{I}({I}_{i}+{I}_{i}^{{rm{h}}})+{gamma }^{A}({A}_{i}+{A}_{i}^{{rm{q}}})-mu {R}_{i}.$$
    (3)
    Susceptible individuals are recruited into community i (i = 1…n) at a constant rate μNi, with μ and Ni being the average mortality rate of the population and the size of the community in the absence of disease, respectively, and die at rate μ. In this way, the equilibrium size of community i without disease amounts to Ni. Susceptible individuals get exposed to the pathogen at rate λi, corresponding to the force of infection for community i (detailed below), thus becoming latently infected (but not infectious yet). Exposed individuals die at rate μ and transition to the post-latent, infectious stage at rate δE. If containment measures including mass testing and preventive isolation of positive cases are in place, exposed individuals may be removed from the general population and quarantined at rate ({chi }_{i}^{E}). Post-latent individuals die at rate μ, progress to the next infectious classes at rate ηP, developing an infection that can be either symptomatic—with probability σ—or asymptomatic, including the case in which only mild symptoms are present—with probability 1 − σ, and may be tested and quarantined at rate ({chi }_{i}^{P}). Symptomatic infectious individuals die at rate μ + α, with α being an extra-mortality term associated with disease-related complications, recover from infection at rate γI, may spontaneously seek treatment at a hospital at rate η, and may be identified through mass screening and hospitalized at rate ({chi }_{i}^{I}). Asymptomatic individuals die at rate μ, recover at rate γA, and may be quarantined at rate ({chi }_{i}^{A}). Infected individuals who are either hospitalized or quarantined at home are subject to the same epidemiological dynamics as those who are still in the community, but are considered to be effectively removed from it, thus not contributing to disease transmission. Individuals who recover from the infection die at rate μ, and are assumed to have permanent immunity to reinfection. This last assumption is not fundamental, as loss of immunity can be easily included in the model. However, immunity to SARS-CoV-2 reinfection is reported to be relatively long-lasting (a few months at least), hence its loss cannot alter transmission dynamics over epidemic timescales14.The cornerstone of model (Eq. (3)) is the force of infection, λi, which in a spatially explicit setting must account not only for locally acquired infections but also for the role played by human mobility. We assume that, at the spatiotemporal scales of interest for our problem, human mobility mostly depicts daily commuting flows (also coherently with the data available for parameterization; see above) and does not actually entail a permanent relocation of individuals. We thus describe human mobility (and the associated social contacts possibly conducive to disease transmission) by means of instantaneous spatial-mixing matrices ({M}_{c,ij}^{X}) (with X ∈ {S, E, P, I, A, R}), i.e.,$${M}_{c,ij}^{X}=left{begin{array}{ll}{r}^{X}{p}_{i}{q}_{ij}(1-{xi }_{ij})hfill&,{text{if}},i,ne, jhfill\ (1-{p}_{i})+(1-{r}^{X}){p}_{i}+{r}^{X}{p}_{i}{q}_{ij}(1-{xi }_{ij})&,{text{if}},i=j,end{array}right.$$
    (4)
    where pi (0 ≤ pi ≤ 1 for all i’s) is the fraction of mobile people in community i, qij (0 ≤ qij ≤ 1 for all i’s and j’s) represents the fraction of people moving between i and j (including j = i, (mathop{sum }nolimits_{j = 1}^{n}{q}_{ij}=1) for all i’s), rX (0 ≤ rX ≤ 1 for all X’s) quantifies the fraction of contacts occurring while individuals in epidemiological compartment X are traveling, and ξij (0 ≤ ξij ≤ 1 for all i’s and j’s) represents the effects of travel restrictions that may be imposed between any two communities i and j as a part of the containment response. Therefore, the probability that residents from i have social contacts while being in j (independently of with whom) is assumed to be proportional to the fraction rX of the mobility-related contacts of the individuals in epidemiological compartment X, multiplied by the probability pi that people from i travel (independently of the destination) and the probability qij that the travel occurs between i and j, possibly reduced by a factor 1 − ξij accounting for travel restrictions. All other contacts contribute to mixing within the local community (i in this case). Note also that if ξij = 0 for all i’s and j’s, then ({M}_{c,ij}^{X}) reduces to ({M}_{ij}^{X}), i.e., to the mixing matrix in the absence of disease-containment measures. In this case, (mathop{sum }nolimits_{j = 1}^{n}{M}_{ij}^{X}=1) for all i’s and X’s. It is important to remark, though, that the epidemiologically relevant contacts between the residents of two different communities, say i and j, may not necessarily occur in either i or j; in fact, they could happen anywhere else, say in community k, between residents of i and j simultaneously traveling to k. On this basis, we define the force of infection as$${lambda }_{i}=mathop{sum }limits_{j=1}^{n}{M}_{c,ij}^{S}frac{(1-{epsilon }_{j})left({beta }_{j}^{P}mathop{sum }nolimits_{k = 1}^{n}{M}_{c,kj}^{P}{P}_{k}+{beta }_{j}^{I}mathop{sum }nolimits_{k = 1}^{n}{M}_{c,kj}^{I}{I}_{k}+{beta }_{j}^{A}mathop{sum }nolimits_{k = 1}^{n}{M}_{c,kj}^{A}{A}_{k}right)}{mathop{sum }nolimits_{k = 1}^{n}left({M}_{c,kj}^{S}{S}_{k}+{M}_{c,kj}^{E}{E}_{k}+{M}_{c,kj}^{P}{P}_{k}+{M}_{c,kj}^{I}{I}_{k}+{M}_{c,kj}^{A}{A}_{k}+{M}_{c,kj}^{R}{R}_{k}right)},$$
    (5)
    where the parameters ({beta }_{j}^{X}) (X ∈ {P, I, A}) are the community-dependent rates of disease transmission from the three infectious classes, ϵj (0 ≤ ϵj ≤ 1 for all j’s) represents the reduction of transmission induced by social distancing, the use of personal protective equipment, and local mobility restrictions if such containment measures are in fact in place, and the terms ({M}_{c,ij}^{X}) (with X ∈ {S, E, P, I, A, R}) describe the epidemiological effects of mobility between i and j in the presence of disease-containment measures. Note that transmission has been assumed to be frequency-dependent.The parameters μ, δX (X ∈ {E, P}), σ, α, η, γX (X ∈ {I, A}), and rX (X ∈ {S, E, P, I, A, R}) are assumed to be community-independent, for they pertain to population demography at the country scale or the clinical course of the disease. By contrast, the transmission rates ({beta }_{i}^{X}) (X ∈ {P, I, A}) and the control parameters, namely the isolation rates ({chi }_{i}^{X}) (X ∈ {E, P, I, A}), the reductions of transmission due to personal protection, social distancing, and local mobility restriction ϵi, and the travel restrictions ξij, are assumed to be possibly community-dependent, thereby reflecting spatial heterogeneities in disease transmission prior to the implementation of containment measures (({beta }_{i}^{X})), testing effort and/or strategy (({chi }_{i}^{X})), local transmission reduction (ϵi), and travel restriction (ξij).Derivation of the basic and control reproduction numbersClose to the DFE, a state in which all individuals are susceptible to the disease (Si = Ni, with Ni being the baseline population size of community i) and all the other epidemiological compartments are empty (({E}_{i}={P}_{i}={I}_{i}={A}_{i}={E}_{i}^{{rm{q}}}={P}_{i}^{{rm{q}}}={I}_{i}^{{rm{h}}}={A}_{i}^{{rm{q}}}={R}_{i}=0) for all i’s), the dynamics of model (Eq. (3)) is described by the linearized system (dot{{bf{x}}}={{bf{J}}}_{{bf{c}}}{bf{x}}), where ({bf{x}}={[{S}_{i},{E}_{i},{P}_{i},{I}_{i},{A}_{i},{E}_{i}^{{rm{q}}},{P}_{i}^{{rm{q}}},{I}_{i}^{{rm{h}}},{A}_{i}^{{rm{q}}},{R}_{i}]}^{T}) (where i = 1…n and the superscript T denotes matrix transposition) and Jc is the spatial Jacobian matrix$${{bf{J}}}_{{bf{c}}}=left[begin{array}{llllllllll}-mu {bf{I}}&{bf{0}}&-{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&-{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&-{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{delta }^{E}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&sigma {delta }^{P}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{{boldsymbol{chi }}}^{{bf{E}}}&{bf{0}}&{bf{0}}&{bf{0}}&-(mu +{delta }^{E}){bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{{boldsymbol{chi }}}^{{bf{P}}}&{bf{0}}&{bf{0}}&{delta }^{E}{bf{I}}&-(mu +{delta }^{P}){bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&eta {bf{I}}+{{boldsymbol{chi }}}^{{bf{I}}}&{bf{0}}&{bf{0}}&sigma {delta }^{P}{bf{I}}&-(mu +alpha +{gamma }^{I}){bf{I}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{{boldsymbol{chi }}}^{{bf{A}}}&{bf{0}}&(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-(mu +{gamma }^{A}){bf{I}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{gamma }^{I}{bf{I}}&{gamma }^{A}{bf{I}}&{bf{0}}&{bf{0}}&{gamma }^{I}{bf{I}}&{gamma }^{A}{bf{I}}&-mu {bf{I}}end{array}right],$$
    (6)
    where I and 0 are the identity and null matrices of size n, respectively, ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{X}}}) (X ∈ {E, P, I, A}) are diagonal matrices whose non-zero elements are (mu +{delta }^{E}+{chi }_{i}^{E}) (for ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}})), (mu +{delta }^{P}+{chi }_{i}^{P}) (for ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}})), (mu +alpha +eta +{gamma }^{I}+{chi }_{i}^{I}) (for ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})), and (mu +{gamma }^{A}+{chi }_{i}^{A}) (for ({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})), and the matrices ({{boldsymbol{theta }}}_{{bf{c}}}^{{bf{X}}}) (X ∈ {P, I, A}) are given by$${{boldsymbol{theta }}}_{{bf{c}}}^{{bf{X}}}={bf{N}}{{bf{M}}}_{{bf{c}}}^{{bf{S}}}({bf{I}}-{boldsymbol{epsilon }}){{boldsymbol{beta }}}^{{bf{X}}}{({{boldsymbol{Delta }}}_{{bf{c}}})}^{-1}{({{bf{M}}}_{{bf{c}}}^{{bf{X}}})}^{T},$$
    (7)
    where N is a diagonal matrix whose non-zero elements are the population sizes Ni, ({{bf{M}}}_{{bf{c}}}^{{bf{X}}}=[{M}_{c,ij}^{X}]) (X ∈ {S, P, I, A}) are sub-stochastic matrices representing the spatially explicit contact terms in the presence of containment measures, ϵ is a diagonal matrix whose non-zero entries are the transmission reductions ϵi, βX (X ∈ {P, I, A}) are diagonal matrices whose non-zero elements are the contact rates ({beta }_{i}^{X}), and Δc is a diagonal matrix whose non-zero entries are the elements of vector ({bf{u}}{bf{N}}{{bf{M}}}_{{bf{c}}}^{{bf{S}}}), with u being a unitary row vector of size n.Because of its block-triangular structure, it is immediate to see that Jc has 6n strictly negative eigenvalues, namely −μ, with multiplicity 2n, and −(μ + δE),−(μ + δP), −(μ + α + γI), and −(μ + γA), each with multiplicity n. Therefore, the asymptotic stability properties of the DFE of model (Eq. (3)), which determine whether long-term disease circulation in the presence of controls is possible, are linked to the eigenvalues of a reduced-order spatial Jacobian associated with the infection subsystem, i.e., the subset of state variables directly related to disease transmission, in this case {E1, …, En, P1, …, Pn, I1, …, In, A1, …, An}. Note that introducing waning immunity would not change the spectral properties of the Jacobian matrix evaluated at the DFE. The reduced-order Jacobian ({{bf{J}}}_{{bf{c}}}^{* }) thus reads$${{bf{J}}}_{{bf{c}}}^{* }=left[begin{array}{llll}-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}\ {delta }^{E}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}&{bf{0}}&{bf{0}}\ {bf{0}}&sigma {delta }^{P}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}}&{bf{0}}\ {bf{0}}&(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}}end{array}right].$$
    (8)
    The asymptotic stability properties of the DFE can be assessed through a NGM approach22,37. In fact, the spectral radius of the NGM provides an estimate of the so-called control reproduction number58, ({{mathcal{R}}}_{{rm{c}}}), which can be thought of as the average number of secondary infections produced by one infected individual in a completely susceptible population in the presence of disease-containment measures. Clearly, if ({{mathcal{R}}}_{{rm{c}}}, > , 1) the pathogen can invade the population in the long run, and endemic transmission will eventually be established despite the implementation of disease-containment measures. To evaluate ({{mathcal{R}}}_{{rm{c}}}) for model (Eq. (3)), the Jacobian of the infection subsystem can be decomposed into a spatial transmission matrix$${{bf{T}}}_{{bf{c}}}=left[begin{array}{llll}{bf{0}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}end{array}right],$$
    (9)
    and a transition matrix$${{boldsymbol{Sigma }}}_{{bf{c}}}=left[begin{array}{llll}-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}&{bf{0}}&{bf{0}}&{bf{0}}\ {delta }^{E}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}&{bf{0}}&{bf{0}}\ {bf{0}}&sigma {delta }^{P}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}}&{bf{0}}\ {bf{0}}&(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}}end{array}right],$$
    (10)
    so that Jc = Tc + Σc. The spatial NGM with large domain ({{bf{K}}}_{{bf{c}}}^{{bf{L}}}), including variables other than the states-at-infection59 (i.e., the exposed individuals Ei) thus reads$${{bf{K}}}_{{bf{c}}}^{{bf{L}}}=-{{bf{T}}}_{{bf{c}}}{({{mathbf{Sigma }}}_{{bf{c}}})}^{-1}=left[begin{array}{llll}{{bf{K}}}_{{bf{c}}}^{{bf{1}}}&{{bf{K}}}_{{bf{c}}}^{{bf{2}}}&{{bf{K}}}_{{bf{c}}}^{{bf{3}}}&{{bf{K}}}_{{bf{c}}}^{{bf{4}}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}end{array}right],$$
    (11)
    with$${{bf{K}}}_{{bf{c}}}^{{bf{1}}} ={delta }^{E}left[{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}+sigma {delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})}^{-1}+(1-sigma ){delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})}^{-1}right]{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}})}^{-1}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}})}^{-1}\ {{bf{K}}}_{{bf{c}}}^{{bf{2}}} =left[{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}+sigma {delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})}^{-1}+(1-sigma ){delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})}^{-1}right]{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}})}^{-1}\ {{bf{K}}}_{{bf{c}}}^{{bf{3}}} ={{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})}^{-1}\ {{bf{K}}}_{{bf{c}}}^{{bf{4}}} ={{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})}^{-1}.$$
    (12)
    Because of the peculiar block-triangular structure of ({{bf{K}}}_{{bf{c}}}^{{bf{L}}}), the spatial NGM with small domain (Kc, accounting only for Ei) is simply ({{bf{K}}}_{{bf{c}}}^{{bf{1}}}) (see again Diekmann et al.59). The control reproduction number can thus be found as the spectral radius of the NGM (with either large or small domain), i.e.,$${{mathcal{R}}}_{{rm{c}}}=rho ({{bf{K}}}_{{bf{c}}}^{{bf{L}}})=rho ({{bf{K}}}_{{bf{c}}})=rho ({{bf{G}}}_{{bf{c}}}^{{bf{P}}}+{{bf{G}}}_{{bf{c}}}^{{bf{I}}}+{{bf{G}}}_{{bf{c}}}^{{bf{A}}}),$$
    (13)
    where$${{bf{G}}}_{{bf{c}}}^{{bf{P}}} ={delta }^{E}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}})}^{-1}\ {{bf{G}}}_{{bf{c}}}^{{bf{I}}} =sigma {delta }^{E}{delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}})}^{-1}\ {{bf{G}}}_{{bf{c}}}^{{bf{A}}} =(1-sigma ){delta }^{E}{delta }^{P}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}{({{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}})}^{-1}$$
    (14)
    are three spatially explicit generation matrices describing the contributions of post-latent infectious people, infectious symptomatic people, and asymptomatic/paucisymptomatic infectious people to the next generation of infections in a neighborhood of the DFE in the presence of disease-containment measures.In the absence of controls, i.e., if the isolation rates ({chi }_{i}^{X}) (X ∈ {E, P, I, A}), the transmission reductions ϵi, and the travel restrictions ξij are equal to zero for all i’s and j’s, then the control reproduction number ({{mathcal{R}}}_{{rm{c}}}) reduces to the basic reproduction number ({{mathcal{R}}}_{0}), defined as the average number of secondary infections produced by one infected individual in a population that is completely susceptible to the disease and where no containment measures are in place. ({{mathcal{R}}}_{0}) can be evaluated as the spectral radius of matrix GP + GI + GA, where$${{bf{G}}}^{{bf{P}}} ={delta }^{E}{{boldsymbol{theta }}}^{{bf{P}}}{({{boldsymbol{phi }}}^{{bf{E}}}{{boldsymbol{phi }}}^{{bf{P}}})}^{-1}\ {{bf{G}}}^{{bf{I}}} =sigma {delta }^{E}{delta }^{P}{{boldsymbol{theta }}}^{{bf{I}}}{({{boldsymbol{phi }}}^{{bf{E}}}{{boldsymbol{phi }}}^{{bf{P}}}{{boldsymbol{phi }}}^{{bf{I}}})}^{-1}\ {{bf{G}}}^{{bf{A}}} =(1-sigma ){delta }^{E}{delta }^{P}{{boldsymbol{theta }}}^{{bf{A}}}{({{boldsymbol{phi }}}^{{bf{E}}}{{boldsymbol{phi }}}^{{bf{P}}}{{boldsymbol{phi }}}^{{bf{A}}})}^{-1}.$$
    (15)
    In the previous set of expressions, ϕX (X ∈ {E, P, I, A}) are diagonal matrices whose non-zero elements are μ + δE (for ϕE), μ + δP (for ϕP), μ + α + η + γI (for ϕI), and μ + γA (for ϕA), while matrices θX (X ∈ {P, I, A}) are given by ({bf{N}}{{bf{M}}}^{{bf{S}}}{{boldsymbol{beta }}}^{{bf{X}}}{({boldsymbol{Delta }})}^{-1}{({{bf{M}}}^{{bf{X}}})}^{T}), with ({{bf{M}}}^{{bf{X}}}=[{M}_{ij}^{X}]) (X ∈ {S, P, I, A}) and ({M}_{ij}^{X}={M}_{c,ij}^{X}) evaluated with ξij = 0 for all i’s and j’s, and Δ is a diagonal matrix whose non-zero entries are the elements of vector uNMS.Derivation of basic and control epidemicity indicesThe concept of epidemicity26 extends previous work24,25 where a reactivity index was defined and applied to study the transient dynamics of ecological systems characterized by steady-state behavior. To explain, in physical terms, the meaning of reactivity and of the Hermitian matrix used to derive it, consider a linear system dx/dt = Ax, where ({bf{x}}={({x}_{1},ldots ,{x}_{n})}^{T}) is the state vector and A is a n × n real state matrix. The system is subject to pulse perturbations x(0) = x0  > 0. Reactivity is defined as the gradient of the Euclidean norm (| | {bf{x}}| | =sqrt{{x}_{1}^{2}+cdots +{x}_{n}^{2}}=sqrt{{{bf{x}}}^{T}{bf{x}}}) of the state vector, evaluated for the fastest-growing initial perturbation, and corresponds to the spectral abscissa ({{{Lambda }}}_{max }^{{rm{Re}}}(cdot )) of the Hermitian part (A + AT)/2 of matrix A24. Following Mari et al.25, an asymptotically stable equilibrium is characterized by positive generalized reactivity if there exist small perturbations that can lead to a transient growth in the Euclidean norm of a suitable system output y = Wx, with matrix W describing a linear transformation of the system state.In epidemiological applications, W should include the variables of the infection subsystem26. Therefore, a suitable output transformation for the problem at hand is$${bf{W}}=left[begin{array}{llllllllll}{bf{0}}&{w}^{E}{bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{w}^{P}{bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{w}^{I}{bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}\ {bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{w}^{A}{bf{I}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}&{bf{0}}end{array}right],$$
    (16)
    where wE, wP, wI, wA are the weights assigned to the variables of the infection subsystem in the output ({bf{y}}=[{w}^{E}{E}_{1},ldots ,{w}^{E}{E}_{n},{w}^{P}{P}_{1},ldots ,{w}^{P}{P}_{n},{w}^{I}{I}_{1},ldots ,{w}^{I}{I}_{n},{w}^{A}{A}_{1},ldots ,{w}^{A}{A}_{n}]^{T}). Generalized reactivity for the DFE of system (Eq. (3)) is positive if the spectral abscissa of a suitable Hermitian matrix (either H0 or Hc, depending on whether the spread of disease is uncontrolled or some containment measures are in place) is also positive. In SEPIAR, the expressions of matrices H0 and Hc are far from trivial, as shown below, and the evaluation of spectral abscissae typically requires numerical techniques. Note also that, since recovered individuals are not accounted for in the system output, including waning immunity would not alter the epidemicity properties of the DFE.Let us consider the most general case of disease-containment measures being in place (which includes as a limit case also uncontrolled pathogen spread). If we note that (ker ({bf{W}})=ker ({bf{W}}{{bf{J}}}_{{bf{c}}})), with Jc being the Jacobian of SEPIAR at the DFE in the presence of controls, matrix Hc can be defined25,27 as the Hermitian part of WJc(W)+, i.e.,$${{bf{H}}}_{{bf{c}}}=H({bf{W}}{{bf{J}}}_{{bf{c}}}{({bf{W}})}^{+})=frac{1}{2}left{{bf{W}}{{bf{J}}}_{{bf{c}}}{({bf{W}})}^{+}+{[{({bf{W}})}^{+}]}^{T}{({{bf{J}}}_{{bf{c}}})}^{T}{({bf{W}})}^{T}right},$$
    (17)
    where (W)+ is the right pseudo-inverse (a generalization of the concept of inverse for non-square matrices) of W, and can be evaluated as$${({bf{W}})}^{+}={({bf{W}})}^{T}{[{bf{W}}{({bf{W}})}^{T}]}^{-1}.$$
    (18)
    Matrix$${{bf{H}}}_{{bf{c}}}=left[begin{array}{llll}-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{E}}}&frac{{w}^{P}}{2{w}^{E}}{delta }^{E}{bf{I}}+frac{{w}^{E}}{2{w}^{P}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&frac{{w}^{E}}{2{w}^{I}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&frac{{w}^{E}}{2{w}^{A}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}\ frac{{w}^{P}}{2{w}^{E}}{delta }^{E}{bf{I}}+frac{{w}^{E}}{2{w}^{P}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{P}}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{P}}}&frac{{w}^{I}}{2{w}^{P}}sigma {delta }^{P}{bf{I}}&frac{{w}^{A}}{2{w}^{P}}(1-sigma ){delta }^{P}{bf{I}}\ frac{{w}^{E}}{2{w}^{I}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{I}}}&frac{{w}^{I}}{2{w}^{P}}sigma {delta }^{P}{bf{I}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{I}}}&{bf{0}}\ frac{{w}^{E}}{2{w}^{A}}{{boldsymbol{theta }}}_{{bf{c}}}^{{bf{A}}}&frac{{w}^{A}}{2{w}^{P}}(1-sigma ){delta }^{P}{bf{I}}&{bf{0}}&-{{boldsymbol{phi }}}_{{bf{c}}}^{{bf{A}}}end{array}right]$$
    (19)
    is Hermitian, hence real and symmetric. Therefore all eigenvalues are real and the spectral abscissa ({e}_{{rm{c}}}={{{Lambda }}}_{max }^{{rm{Re}}}({{bf{H}}}_{{bf{c}}})) coincides with the largest eigenvalue, which corresponds to the fastest-growing perturbation in the system output. Thus, ec can be interpreted as a control epidemicity index: if ec  > 0, there must exist some small perturbations to the DFE that are temporarily amplified in the system output, thus generating a transient, subthreshold epidemic wave.Absent any containment measures, the control epidemicity index, ec, reduces to the basic epidemicity index, ({e}_{0}={{{Lambda }}}_{max }^{{rm{Re}}}({{bf{H}}}_{{bf{0}}})), where$${{bf{H}}}_{{bf{0}}}=H({bf{W}}{{bf{J}}}_{{bf{0}}}{({bf{W}})}^{+})=frac{1}{2}left{{bf{W}}{{bf{J}}}_{{bf{0}}}{({bf{W}})}^{+}+{[{({bf{W}})}^{+}]}^{T}{({{bf{J}}}_{{bf{0}}})}^{T}{({bf{W}})}^{T}right}$$
    (20)
    and the Jacobian matrix J0 can be obtained from Jc by setting equal to zero the isolation rates ({chi }_{i}^{X}) (X ∈ {E, P, I, A}), the transmission reductions ϵi, and the travel restrictions ξij for all i’s and j’s.The effective reproduction number and the effective epidemicity indexThe reproduction numbers and the epidemicity indices defined above can be rigorously applied only to characterize the spread of disease in a fully naïve population (Si = Ni ∀ i). As soon as the pathogen begins to circulate within the population, the state of the system gradually departs from the DFE. Under these circumstances, it is customary19,21 to define a time-dependent, effective reproduction number, ({mathcal{R}}(t)), to track the number of secondary infections caused by a single infectious individual in a population in which the pool of susceptible individuals is progressively depleted, and control measures are possibly in place58. Similarly, it is possible to define an effective epidemicity index, e(t), to evaluate the likelihood that transient epidemic waves may occur even if ({mathcal{R}}(t), More

  • in

    Spatial and temporal analysis of cumulative environmental effects of offshore wind farms in the North Sea basin

    The area of study (Fig. 6) was the Greater North Sea ecoregion, which includes the EEZs of six countries (England, Scotland, the Netherlands, Denmark, Norway and Germany). The Kattegat area, the English Channel, and the Belgium EEZ were omitted from the study area. The North Sea Marine Ecosystem is a large semi-closed continental sea situated on the continental shelf of North-western Europe, with a dominant physical division between the comparatively deep northern part (50–200 m, with the Norwegian Trench dropping to 700 m) and the shallower southern part (20–50 m)48. The North Sea is one of the most varied coastal regions in the world, which is characterised by, among others, rocky, fjord and mountainous shores as well as sandy beaches with dunes48. Apart from the marine seabirds feeding primarily in the coastal areas, under 5 km from the coast (e.g., terns, sea-ducks, grebes), the North Sea basin also hosts pelagic birds feeding further offshore, with some also diving for food (guillemot, razorbill, etc.). The North Sea basin is also a major habitat for four marine mammal species, of which the harbour porpoise and harbour seal are the most common. Moreover, fish ecology has been a widely studied topic, especially for commercial species, due to evidence of a decline in the fish stock, such as sprat, whiting, bib, and mackerel. Fish communities, and in particular the small pelagic fish group (such as European sprat, European pilchard), play also a key ecologic role, constituting the main pray for most piscivorous fishes, cetacean and seabirds49, Based on early surveys, the predominant species divided by the three North Sea fish communities are: saithe (43.6% in the shelf edge), haddock (42.4% in the central North Sea, 11.6% in the shelf edge), whiting (21.6% in the eastern North Sea, 13.9% central North Sea), and dab (21.8% in the eastern North Sea)34. More recent assessments of North Sea fish community are emphasizing the clear geographical distinction between the fish species living in the southern part of the North Sea, a shallow area with high primary production and pronounced seasonality, and northern part, a deeper area with lower primary production and lower seasonal variation in temperature and salinity. The southern North Sea fish community is represented by fish species such as lesser weever, while the northern North Sea fish community is represented by species such as saithe, with species like whitting, haddock representative for the North–West subdivision, and the European plaice having the highest abundance in the South–East community50. The future fish stock and spatial distribution is however uncertain due to impacts of climate change related factors (e.g., growing temperatures)49 and overexploitation.Figure 6Offshore wind farm prospects (existing/authorised/planned) in the North Sea basin.Full size imageThe most prominent human activities in the North Sea basin are fishing, coastal construction, maritime transport, oil and gas exploration and production, tourism, military, and OWF construction38. Within this list, the construction of OWFs has seen a rapid increase, aiming to reach a total cumulative installed capacity of 61.8–66.8 GW by 203051. As indicated in Fig. 6, the new designated/search/scoping areas for the location of future OWFs will significantly increase the current space reserved for the offshore production of renewable energy in the North Sea basin.Spatio-temporal database of OWF developments in the North Sea basinFor the input of the geo-spatial layers with the location of OWF areas we compiled a comprehensive spatial data repository in QGIS containing the shapefiles of analysed OWF, from 1999 to 2027 (last year of available official information on OWF development, Appendix D). The analysis was performed for the North Sea geographic area, referred here as the basin scale, taking into account the cumulative pressures from individual OWF projects (project scale). The main data sources for geospatial information for OWF, for the entire North Sea basin, are EMODnet (Human Activities data portal) and OSPAR, which were complemented by data on the country level, where needed; i.e. from Crown Estate Scotland (Energy infrastructure, Legal Agreements), Rijkswaterstaat for the Netherlands. From the available geo-spatial data for OWF, we selected the OWF in our area of study (Fig. 6) with the status of consent-authorised, authorised, pre-construction, under construction, or fully commissioned (operational). Therefore, planned OWF such as Vesterhavet Syd and Vesterhavet Nord, for which the start date of construction is still unknown, were not included in the analysis. Similarly, for the Horns Rev 3 OWF no geo-referenced spatial footprint was available in the open-access data sets, and therefore it was not included in the analysis.The collected OWF geospatial data was aggregated to create a geospatial database, for the studied period of 1999–2050, composed by the following attributes: code name, country, name, production capacity (MW), area (({mathrm{km}}^{2})), number of turbines, start operation (year), installation time, and status in the period 1999–2050 (construction, operation, decommissioning). The created geospatial dataset was additionally cross-checked for integrity with the information provided through the online platform 4coffshore.com.The lack of data regarding the construction time was complemented with the methodology proposed by Lacal-Arántegui et al.36. Based on this research, we calculated the time required for OWF construction phase related activities multiplying 1.06 days by the known production capacity (total MW) for each analysed OWF.The average time of operation is considered to be 20 years, probably profitably extendable to 25 years, as stated in a number of studies on the cycle of offshore wind farms52. For this case study, the operation time considered is 20 years (subject to change). Since there is little experience with the decommissioning of offshore wind farms (only a few OWFs have so far been decommissioned in the UK and Denmark), the decommissioning time is not yet clear. There are a number of parameters that influence the decommissioning time, which are: the number of turbines, the foundation type, the distance to port, etc. It is estimated that the time taken for decommissioning should be around 50–60% less than the installation time37. Our study considers the decommissioning time as 50% of the construction time.Time-aware cumulative effects assessmentIn this study, Tools4MSP53,54, a Python-based Free and Open Source Software (FOSS) for geospatial analysis in support of Maritime Spatial Planning and marine environmental management, was used for the assessment of the impacts of OWFs on the marine ecosystem, in the three development stages. We applied the Tools4MSP CEA module to the OWF of the North Sea basin for the period 1999–2050, taking into account the full life cycle of the OWF development, namely the construction, operation and decommissioning phases. The modified methodology from Menegon et al.31 and subsequent implementation55, proposes to calculate the CEA score for each cell of analysis as follows (Eqs. 1, 2):$$CEA=sum_{k=1}^{n}d({E}_{k}) sum_{j=1}^{m}{s}_{i,j} eff({P}_{j}{E}_{k})$$
    (1)
    where eff is the effect of pressure P over the environmental component E and is defined as follows:$$eff left({P}_{j}{E}_{k}right)=(sum_{i=1}^{l}{w}_{i,j} i({U}_{i},{M}_{i,j,k})){^{prime}}$$
    (2)
    whereas,

    ({U}_{i}) defines the human activity, namely the OWF activity in the study area

    ({E}_{k}) defines the environmental components of the study area described in the Table 1

    ({d(E}_{k})) defines intensity or presence/absence of the k-th environmental component

    ({P}_{j}) defines the pressures exerted by human activities dependent on the three different OWF development phases (Annex B)

    ({w}_{i,j}) refers to the specific pressure weight according to the OWF phase

    ({s(P}_{j}, {E}_{k})) is the sensitivity of the k-th environmental component to the j-th pressure

    ({i({U}_{i, }M({U}_{i, }P}_{j}, {E}_{k}))) is the distance model propagating j-th pressure caused by i-th activity over the k-th environmental component

    ({M(U}_{i}, {P}_{j})) is the 2D Gaussian kernel function used for convolution, which considers buffer distances at 1 km, 5 km, 10 km, 20 km, and 50 km56.

    Table 1 Primary sources for the environmental component data sets.Full size tableIn Eq. (3), the CEA 1999–2050 describes the modelling over the time frame 1999–2050, whereas ({CEA}_{t}) is the cumulative effect of year t within the timeframe 1999–2050:$${CEA}_{1999-2050}= sum_{t=1999}^{2050}{CEA}_{t}$$In this study, each final CEA score was normalised. To normalise the value of each initial CEA score obtained using the Eq. (1), we calculated its percentage of the sum of all CEA scores for all OWFs in the three development phases, period spanning the period 1999–2050 (({CEA}_{1999-2050})).Environmental componentsThe selection of the environmental components (receptors) impacted by the identified pressures is an essential part of the scoping phase for OWF location, as monitoring the status (distribution, abundance) of different identified species represents a relevant indicator for the ecosystem status. For the evaluation of the habitats and species that can be affected by the cumulative ecological effects of OWF, we adapted the methodology of Meissl et al.14. Therefore, we selected the environmental components based first on their: (1) ecological value, supported by legal documents identifying species protected by law or through various national and international agreements (e.g. EU Habitats Directive, Wild Mammals (Protection) Act (UK), see Table 1 in Appendix E), to which we added species with (2) commercial value, but also with a (3) broad geographic-scale habitat occurrence of the species in the studied area, based on previous studies35 and on 35 EIA studies for OWF in the North Sea basin.Among the five fish species selected, sprat and sandeel play key roles in the marine food web (small pelagic fish), as prey source for piscivorous fish, cetacean and birds. The ecological value of sandeel, sprat, whiting and saither is also highlighted through EU or national protection agreements such as Priority Marine Features—PMF or Scottish/UK Biodiversity list (see Appendix E, Table 2). The list is completed by haddock, one of the fish species with commercial importance, highly dominant in the Central North Sea. With regards to the spatial occurrence at the basin level, the fish species selected are representative for both of the two distinct North Sea communities50, the southern part of the North Sea (sprat), and the northern and north-west part (haddock, whiting, saithe).The three selected seabird species are of ecological importance for the marine ecosystem, as indicated through the European, national and international protection agreements, such as the EU Birds Directive Migratory Species or the IUCN Red List (see Appendix E, Table 1). While razorbill and guillemot have similar feeding and flying patterns (low flight, catch pray underwater), there is evidence of different behaviors towards OWFs, with relatively more avoidance from razorbill compared to guillemot. In relation to the spatial distribution of the three selected species, there is a clear distinction between razorbill, highly present in the coastal areas of west North Sea basin, guillemot, with a relatively even distribution across the marine basin, and fulmar, one of the 4 most common seabirds in the studied area, in particular in the central and N–E parts.In the marine mammals category we selected the harbor porpoise, indicated to be one of the most impacted species in this category57, with a high occurrence in the North Sea basin. Its ecological value is emphasized by its presence in European and international lists for habitat protection, such as EU Habitats Directive58, OSPAR List of Threatened and/or Declining Species59, the Agreement on the Conservation of Small Cetaceans in the Baltic and the North Seas (ASCOBANS)60. The harbor porpoise is the protected species in numerous Natura 2000 areas in the North Sea basin, such as the Spatial Area of Conservation Southern North Sea61 (British EEZ) or The Special area of Protection Kleverbank62 (Dutch EEZ).Among the selected fish species, sandeel had the highest occurrence in EIA studies of OWF developments (23 out of 35), while guillemot had the highest occurrence among seabird species (25 out of 35). With an occurrence of 26 out of the 35 analysed EIA document, the harbour porpoise is the most studied mammal in relation to the impact of OWF.As a result, we selected three EUNIS marine seabed habitat types (European Union Nature Information System)58 (Appendix E, Table 2), three seabird species, one mammal species and five fish species (Appendix E, Table 1). The list can be extended; however, for this exercise we considered it sufficient.The data sets used to represent the spatial distribution (presence/absence, intensity) of the environmental components in the studied area were obtained from multiple sources and were used in the Tools4MSP model either directly (EUNIS habitats, marine mammals, seabirds) or further processed using a predictive distribution model (fish species). In the case of EUNIS marine habitats, the data source was the online geo-portal EMODnet, through the Seabed Habitat service (Table 1), which provided GIS polygon layers for each habitat type and was further used to indicate presence/absence of a specific habitat.For the distribution of the selected mammal species, the harbour porpoise, we used the modelling results of Waggit et al.16, translated into maps for the prediction of densities (nr. animals/({mathrm{km}}^{2})). The mapping approach starts with collating data from available surveys, which are further standardised with regards to transect length, number of platform sides, and the effective strip width. Finally, the standardised data sets were used in a binomial and a Poisson model, in association with environmental conditions (Table 1), in order to deliver a homogenous cover of species distribution maps, on 10 km × 10 km spatial resolution grid16.For the distribution of the selected seabird species (razorbill, fulmar, guillemot), we used the results of the SEAPOP program (http://www.seapop.no/en/distribution-status/), through the open-source data portal (https://www2.nina.no/seapop/seapophtml/). The proposed methodology for creating the occurrence density prediction maps, on a 10 × 10 km spatial resolution grid, starts with the modelling of the presence/absence of birds using a binomial distribution and “logit link”. This was followed by the modelling of the number of birds using a Gamma distribution with a “log link” function, which also took into account geographically fixed explanatory variables (geographic position, water depth, and distance to coast).The predictive model for the spatial distribution of fish species biomass (haddock, sandeel, whiting, saithe, sprat) was developed using AI4Blue software, an open-source, python-based library for Artificial Intelligence based geospatial analysis of Blue Growth settings (AI4Blue, 2021)63. The model was based on two types of inputs: (1) the observation data on the presence of species and (2) data on the absence of species (absence data) for the period 2000–2019. Both data types were extracted by the ICES North Sea International Bottom Trawl Survey (NSI-IBTS, extracted survey year 2000–2019 including all available quarters) for commercial fish species, which was accessed on the online ICES-DATRAS database64. Data was extracted using two DATRAS web service Application Programming Interfaces (APIs): (1) the HHData, that returns detailed haul-based meta-data of the survey (e.g. haul position, sampling method etc.) and (2) the CPUEPerLengthPerHaulPerHour for the catch/unit of effort per length of sampled species.The presence data were represented by the catch/unit of effort (CPUE), expressed in kg of biomass of the specified species per one hour of hauling. The biomass was estimated by using the SAMLK (sex-maturity-age-length keys) dataset for ICES standard species. This approach is a viable alternative to presence-only data models, as it tackles the biased outcomes resulting from an non-uniform marine coverage of the data sets (mainly along the shipping routes)65. The absence data were estimated using the methodology presented by Coro et al.65, which detects absence location for the chosen species as the locations in which repeated surveys (with the selected species on the survey’s species target list) report information only on other species.Additionally, the predictive model automatically correlates the presence/absence data with environmental conditions (Appendix E, Table 3) data to more accurately estimate the likelihood of species presence in the North Sea basin. Intersecting a large number of surveys containing observation data on the presence of selected species can return the true absence data locations, which represent a valuable indicator for geographical areas with unsuitable habitat (see methodology by Coro et al.65). Those locations were estimated from abiotic and biotic parameters and differed to the sampling absences which were estimated from surveys without presence data65. The environmental conditions (Appendix E, Table 2) data were accessed through direct queries using the MOTU Client option from the Marine Copernicus database. In order to input the layers to the CEA calculation, the input layer for the biomass was transformed using log[x + 1] to avoid an over-dominance of extreme values and all datasets rescaled from 0 to 1 in order to allow direct comparison on a single, unit-less scale55.The rescaled special distribution of biomass for the selected species are presented in Appendix F (Fig. a–j).OWF pressures and relative weightsA systematic literature review was conducted to reach a first quantification of the OWF pressure weights (({w}_{i,j}),) in the construction, operation, and decommissioning phases (({U}_{i})). The OWF-related pressures specific to each of the phases of the OWF life cycle were based on the comprehensive analysis of all the existing Environmental Impact Assessment (EIA) methodologies used in the North Sea countries14. The review enabled the collection of 18 pressures that were subsequently compared and merged with the pressures established in the Marine Strategy Framework Directive, applied by the EU countries in the assessment of environmental impacts66. Figure 7 illustrates the impact chain linking the three OWF development phases with the exerted 18 pressures and the 12 selected environmental components impacted.Figure 7Impact chain defining OWF phases-pressure-environmental components analysed in the North Sea (the strength of the link between pressures and environmental components is proportional to the sensitivity scores. The order is descending from the pressures with highest impact, as well as from the environmental components most affected).Full size imageSensitivity in this research is defined as the likelihood of change when a pressure is applied to a receptor (environmental component) and is a function of the ability of the receptor to adapt, tolerate or resist change and its ability to recover from the impact67. The criteria for assessing the sensitivities of environmental components is based on MarLIN (Marine Life Information Network) detailed criteria (https://www.marlin.ac.uk/sensitivity/sensitivity_rationale).We validated the weights of pressures (({w}_{i,j}) from 0 to 5) and scores of environmental components sensitivities (({s(P}_{j}, {E}_{k})) from 0 to 5), as well as the distance of pressure propagation (≤1000 m to ≥ 25,000 m), through a series of 4 questionnaires for the marine mammals, seabirds, fish and seabed habitats. The compiled questionnaires were further validated through semi-interviews of 9 experts in the field of marine ecology, spatial planning, environmental impact assessment and offshore wind energy development. The expert-based questionnaires also included a confidence level for the proposed scores, which ranged between 0.2 (very low confidence: based on expert judgement; proxy assessment) and 1 (very high confidence: based on peer reviewed papers, report, assessment on the same receptor). The confidence level was used in determining the final scores for the pressure weights and species sensitivities. The final scores for weights and sensitivity scores were identified either by calculating the mean value (for cases where literature review scores and expert scores differed by  > 2 units) or selecting the higher value—precautionary principle (for cases where scores from different sources differed by  More