Ecology

The LCA approach includes the potential effects of depriving humans and ecosystems of water resources, as well as the specific potential effects of pollutants affecting water and thus the environment49. Water stress is commonly defined as the ratio of total freshwater consumption to the level of its hydrological availability. ISO 14046 presents a new concept, i.e., WF, which is associated with the LCA approach. The standard’s “water scarcity footprint” refers to the potential impacts associated with the quantitative aspect of water use50. Figure 2 shows the WF per cultivation area of conventional and organic carrot production. In general, there are significant differences in the total value of the WF in question. For conventional carrot production technology, it is 10.25 m3 ha−1, while for organic technology, it is only 1.96 m3 ha−1. In the case of conventional production, treatments using significant amounts of chemicals have the greatest impact on the WF, i.e., fertilization (mainly mineral) (WF = 6.85 m3 ha−1), and chemical plant protection (WF = 1.19 m3 ha−1). The analysis of WF in organic farming showed that its highest value (WF = 0.84 m3 ha−1) concerns the harvesting of carrots, while soil preparation ranks second (WF = 0.45 m3 ha−1). A slightly lower WF of 0.38 m3 ha−1 was recorded in the case of transporting the harvested carrots to the farm buildings. It can therefore be concluded that in organic farming, it is (diesel) fuel consumption that has the greatest impact on WF level.Figure 2Water footprint in conventional and organic carrot production (m3 ha−1).Full size imageIn terms of production volume, in conventional technology, the WF is 0.196 m3 t−1. On the other hand, in organic technology the value of WF is approx. four times lower and amounts to 0.049 m3 t−1 of harvested carrots. For comparison, the WF in tomato production is 160 m3 per 1 tonne of produce51. Such a high value results mainly from irrigation of the plants.In order to explain in detail the impact of individual agricultural treatments on the water deficit in carrot production, a detailed WF analysis was carried out for the treatments that demonstrated the highest values. In the case of conventional technology, it was fertilization (Fig. 3). Upon analyzing Fig. 3, it can be observed that the use of urea, and hence nitrogen, has the greatest impact on WF with regard to fertilization. Most nitrogen mineral fertilizers have a negative impact on the environment, causing ozone depletion in the stratosphere, groundwater pollution, global warming, and water eutrophication52,53. The largest water footprint associated with the use of mineral fertilizers in conventional cultivation is mainly due to a very energy-intensive fertilizer production process. Depending on the type of fertilizer and the technology used, the production process involves machines and equipment for cleaning, grinding, drying, sieving, extruding, granulating, packing, pumping, evaporation (crystallization) and transport. The vast majority of these treatments are powered by electricity. In contrast, conventional power production, regardless of the technology and fuel used (nuclear, natural gas, or coal), is characterized by very high water consumption. Mineral fertilizers are also a material whose consumed mass is relatively high compared to other production materials (seeds, pesticides, and diesel fuel). These two factors mentioned above have a decisive impact on the largest water footprint associated with the use of mineral fertilizers in conventional carrot cultivation. Processes requiring the use of machinery, i.e., fertilizer spreaders (1%) and the consumption of diesel fuel (0.1%) have the lowest impact on the level of fertilization-induced WF. Such a low impact of diesel fuel results mainly from its relatively low consumption during fertilization, most often using very efficient centrifugal spreaders. For comparison, WF related only to the use of carrot irrigation water is 20 m3 t−1 of harvested crops54.Figure 3Water footprint related to carrot fertilization in conventional production (m3 ha−1).Full size imageIn the case of organic technology, WF of harvesting was analyzed in detail (Fig. 4). Carrots were excavated with harvesters, which cut the aboveground parts, cleaned the roots and collected them in a hopper. Sometimes the excavation was preceded by mowing the carrot leaves with mowers. Carrot harvesters are machines that require farm tractors with high-power combustion engines, and the harvesting procedure itself is very time-consuming, hence such a large impact of fuel consumption on WF in carrot harvesting. Despite the above, the share of diesel consumption in the total value of WF related to carrot harvesting is only 11%. However, when comparing the WF related to fuel consumed during harvesting and during fertilization, it can be noticed that in the case of harvesting, WF is approx. 15 times higher.Figure 4Water footprint related to carrot harvest in organic production (m3 ha−1).Full size imageIn LCA, the potential effects of water pollution have traditionally been addressed in impact categories such as (eco) toxicity, acidification, and eutrophication42,43. In the WF analysis, the impact of water consumption is generally related to specific goals within a given conservation area, such as: Human Health, Ecosystems Quality and Resources43. The impact of water consumption on human health is expressed in DALY and is obtained by modeling the cause-effect chain of water scarcity (lack of irrigation water) leading to malnutrition. Ecosystem quality is assessed by modeling the cause-effect chain of freshwater consumption with the quality of the terrestrial ecosystem, based on the number of species disappearing each year (species * year). On the other hand, the impact of water consumption in the resources category is assessed by modeling the cause-effect chain of freshwater consumption in relation to the depletion of water resources, along with the cost ($) of extracting an additional cubic meter of water46. The data in Table 2 shows WF in conventional carrot production related to the three impact categories, and Fig. 5 shows its structure. The total impact of individual processes in the Human Health category is 1.15E−05 DALY, in the Ecosystem Quality category—1.53E−07 species * year, and in the Resources category—2.97 $ surplus. For comparison, WF in the above-mentioned impact areas per 1 ha of tomatoes is, respectively: Human Health—5.00E−03 DALY, Ecosystem Quality—2.50E−05 species * year55. When analyzing Fig. 5, it can be observed that in all impact categories, fertilization has the greatest environmental impact, the share of which in individual categories is at approx. 67.0–67.7%. Chemical plant protection ranks second, the impact of which in the three categories ranges from 11.9 to 12.6%. In addition to the treatments related to fertilizers and chemicals, treatments associated with high consumption of diesel fuel, i.e., soil preparation and harvest, have a significant impact on the value of individual categories in carrot production. This confirms the results of many studies, i.e. that the extraction, production and, above all, the use of diesel fuel bring significant damage to the environment56,57.Table 2 Environmental impact related to the use of water in conventional carrot production per area unit (ha).Full size tableFigure 5The structure of WF in individual impact categories in conventional carrot production.Full size imageBearing in mind that carrot yield range in conventional cultivation is 43–65 t ha−1, the total impact of individual processes per 100 tons of harvested carrots is as follows: Human Health: 2.17E−05 DALY, Ecosystem Quality: 2.88E−07 species * year and Resources: 5.57 $ surplus.Upon comparing the obtained results with the research presented in the literature and conducted with a similar methodology, it can be concluded that for the production of 100 tons of tomatoes, the total environmental footprint for the above-mentioned impact areas, including factors other than water, is respectively: Human Health: 2.7E−01 DALY, Ecosystems Quality: 1.45E−03 Species * year, Resources: 1.05E + 06 $58. On the other hand, the WF for green beans, per 100 tons of harvest was reported as follows: Human Health: from 2.00E−2 to 1.08E−1 DALY, Ecosystem Quality: from1.10E−3 to 1.80E−3 species * year, Resources: from 1.90E+2 to 1.40E+3 $ surplus59.Detailed WF results for the fertilization process in conventional carrot production are presented in Fig. 6. Among the individual factors shaping the environmental impact, what stands out is the consumption of urea, i.e. nitrogen (44.9–47.0% of the total impact in individual categories) and of phosphorus fertilizers, the impact of which is at 31.4–32.4%.Figure 6Structure of WF of the fertilization process in conventional carrot production as per individual impact categories.Full size imageIn endpoint analysis, the impact of water use is generally related to specific endpoints in a given conservation area: Human Health, Ecosystems Quality or Resources43. The data in Table 3 shows WF in organic carrot production as per the three impact categories, and Fig. 7 shows its structure. The total WF values in each category are as follows: in the Human Health category—2.11E−06 DALY, in the Ecosystem Quality category—3.00E−08 species * year and in the Resources category—0.56 $ surplus. The above results are over five times lower compared to the footprint in conventional production (Table 2), and therefore it can be concluded that organic production not only enables the production of healthy carrot, but also has a very positive impact on the broadly understood environment. Upon analyzing the data from Tables 2 and 3, it can be observed that the environmental impact of fertilization treatment in organic production is over thirty times lower compared to the impact of fertilization in conventional production. Moreover, the fact that no pesticides are used means that the impact of chemical plant protection treatments is 0. Upon analyzing Fig. 7, it can be observed that the largest share in the total value of WF in individual impact categories is that of carrot harvest, from 41.9% (Ecosystem Quality) to 43.1% (Resources). The reason for such a significant environmental impact of carrot harvesting technology is the use of complex harvesters, as explained in Fig. 8. When calculating WF, both direct water consumption in a technology is taken into account, as well as indirect, related e.g., to the production of agricultural equipment used in the technology. The complexity of the machinery, the type of materials it is made of and the type of technology used in its production determine the WF.Table 3 Environmental impact related to the use of water in organic carrot production per area unit (ha).Full size tableFigure 7The structure of WF in individual impact categories in organic carrot production.Full size imageFigure 8Structure of WF of the harvest process in organic carrot production as per individual impact categories.Full size imageThe methodology for calculating WF is very diverse and includes many methods. Moreover, the results of research on WF related to the production of vegetable species presented in the literature often differ in terms of the analyzed system boundaries, production technology, irrigation, etc. Therefore, the possibility of a broad discussion of the results of WF of conventional and organic carrot production is limited.The detailed structure of WF of the carrot harvesting process in organic farming is shown in Fig. 8. The use of machines, i.e. harvesters, has a decisive share (90.3–96.6%) in the total value of individual impact categories. The reason for this significant impact was explained above. Relatively high consumption of fuel during harvester operation contributes little to the water footprint structure. The use of diesel fuel has the highest impact on damage caused in the Ecosystem Quality category, with the lowest impact in the Human Health category. More

Rahmstorf, S. & Coumou, D. Increase of extreme events in a warming world. Proc. Natl Acad. Sci. USA 108, 17905–17909 (2011).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Somero, G. N. The physiology of climate change: how potentials for acclimatization and genetic adaptation will determine ‘winners’ and ‘losers’. J. Exp. Biol. 213, 912–920 (2010).Hoffmann, A. A. & Sgrò, C. M. Climate change and evolutionary adaptation. Nature 470, 479–485 (2011).ADS 
CAS 
PubMed 

Google Scholar 
Chevin, L.-M., Lande, R. & Mace, G. M. Adaptation, plasticity, and extinction in a changing environment: towards a predictive theory. PLoS Biol. 8, e1000357 (2010).PubMed 
PubMed Central 

Google Scholar 
Kawecki, T. J. & Ebert, D. Conceptual issues in local adaptation. Ecol. Lett. 7, 1225–1241 (2004).
Google Scholar 
Campbell-Staton, S. C. et al. Winter storms drive rapid phenotypic, regulatory, and genomic shifts in the green anole lizard. Science 357, 495–498 (2017).ADS 
CAS 
PubMed 

Google Scholar 
Barrett, R. D. H. et al. Linking a mutation to survival in wild mice. Science 363, 499–504 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Therkildsen, N. O. et al. Contrasting genomic shifts underlie parallel phenotypic evolution in response to fishing. Science 365, 487–490 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Brennan, R. S., Garrett, A. D., Huber, K. E., Hargarten, H. & Pespeni, M. H. Rare genetic variation and balanced polymorphisms are important for survival in global change conditions. Proc. R. Soc. B: Biol. Sci. 286, 20190943 (2019).CAS 

Google Scholar 
Stearns, S. C. The evolutionary significance of phenotypic plasticity. Bioscience 39, 436–445 (1989).Thompson, J. D. Phenotypic plasticity as a component of evolutionary change. Trends Ecol. Evol. 6, 246–249 (1991).CAS 
PubMed 

Google Scholar 
Kelly, M. Adaptation to climate change through genetic accommodation and assimilation of plastic phenotypes. Philos. Trans. R. Soc. Lond. B Biol. Sci. 374, 20180176 (2019).PubMed 
PubMed Central 

Google Scholar 
Chevin, L. M., Collins, S. & Lefèvre, F. Phenotypic plasticity and evolutionary demographic responses to climate change: taking theory out to the field. Funct. Ecol. https://doi.org/10.1111/j.1365-2435.2012.02043.x (2013).Hendry, A. P. Key questions on the role of phenotypic plasticity in eco-evolutionary dynamics. J. Hered. 107, 25–41 (2016).PubMed 

Google Scholar 
Calosi, P., De Wit, P., Thor, P. & Dupont, S. Will life find a way? Evolution of marine species under global change. Evol. Appl. 9, 1035–1042 (2016).PubMed 
PubMed Central 

Google Scholar 
Fox, R. J., Donelson, J. M., Schunter, C., Ravasi, T. & Gaitán-Espitia, J. D. Beyond buying time: the role of plasticity in phenotypic adaptation to rapid environmental change. Philos. Trans. R. Soc. Lond. B Biol. Sci. 374, 20180174 (2019).PubMed 
PubMed Central 

Google Scholar 
Lande, R. Adaptation to an extraordinary environment by evolution of phenotypic plasticity and genetic assimilation. J. Evol. Biol. 22, 1435–1446 (2009).PubMed 

Google Scholar 
Murren, C. J. et al. Constraints on the evolution of phenotypic plasticity: limits and costs of phenotype and plasticity. Heredity 115, 293–301 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Posavi, M., Gulisija, D., Munro, J. B., Silva, J. C. & Lee, C. E. Rapid evolution of genome-wide gene expression and plasticity during saline to freshwater invasions by the copepod Eurytemora affinis species complex. Mol. Ecol. 29, 4835–4856 (2020).CAS 
PubMed 

Google Scholar 
Ghalambor, C. K. et al. Non-adaptive plasticity potentiates rapid adaptive evolution of gene expression in nature. Nature 525, 372–375 (2015).ADS 
CAS 
PubMed 

Google Scholar 
Kelly, M. W., Pankey, M. S., DeBiasse, M. B. & Plachetzki, D. C. Adaptation to heat stress reduces phenotypic and transcriptional plasticity in a marine copepod. Funct. Ecol. 31, 398–406 (2017).
Google Scholar 
Sikkink, K. L., Reynolds, R. M., Ituarte, C. M., Cresko, W. A. & Phillips, P. C. Rapid evolution of phenotypic plasticity and shifting thresholds of genetic assimilation in the nematode Caenorhabditis remanei. G3 4, 1103–1112 (2014).PubMed 
PubMed Central 

Google Scholar 
Brennan, R. S., Galvez, F. & Whitehead, A. Reciprocal osmotic challenges reveal mechanisms of divergence in phenotypic plasticity in the killifish Fundulus heteroclitus. J. Exp. Biol. 218, 1212–1222 (2015).PubMed 

Google Scholar 
Kelly, M. W., Pankey, M. S. & DeBiasse, M. B. Adaptation to heat stress reduces phenotypic and transcriptional plasticity in a marine copepod. Funct. Ecol. https://doi.org/10.1111/1365-2435.12725 (2017).Waddington, C. H. Genetic assimilation of an acquired character. Evolution 7, 118–126 (1953).
Google Scholar 
Schlötterer, C., Kofler, R., Versace, E., Tobler, R. & Franssen, S. U. Combining experimental evolution with next-generation sequencing: a powerful tool to study adaptation from standing genetic variation. Heredity 114, 431–440 (2015).PubMed 

Google Scholar 
Munday, P. L., Warner, R. R., Monro, K., Pandolfi, J. M. & Marshall, D. J. Predicting evolutionary responses to climate change in the sea. Ecol. Lett. 16, 1488–1500 (2013).PubMed 

Google Scholar 
Huang, Y. & Agrawal, A. F. Experimental evolution of gene expression and plasticity in alternative selective regimes. PLoS Genet. 12, e1006336 (2016).PubMed 
PubMed Central 

Google Scholar 
Mallard, F., Nolte, V. & Schlötterer, C. The evolution of phenotypic plasticity in response to temperature stress. Genome Biol. Evol. 12, 2429–2440 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Schaum, C. E. & Collins, S. Plasticity predicts evolution in a marine alga. Proc. Biol. Sci. 281, 20141486 (2014).Kelly, S. A., Czech, P. P., Wight, J. T., Blank, K. M. & Garland, T. Jr Experimental evolution and phenotypic plasticity of hindlimb bones in high-activity house mice. J. Morphol. 267, 360–374 (2006).PubMed 

Google Scholar 
Garland, T. Jr & Kelly, S. A. Phenotypic plasticity and experimental evolution. J. Exp. Biol. 209, 2344–2361 (2006).PubMed 

Google Scholar 
Gibbin, E. M., Massamba N’Siala, G., Chakravarti, L. J., Jarrold, M. D. & Calosi, P. The evolution of phenotypic plasticity under global change. Sci. Rep. 7, 17253 (2017).ADS 
PubMed 
PubMed Central 

Google Scholar 
McCairns, R. J. S. & Bernatchez, L. Adaptive divergence between freshwater and marine sticklebacks: insights into the role of phenotypic plasticity from an integrated analysis of candidate gene expression. Evolution 64, 1029–1047 (2010).CAS 
PubMed 

Google Scholar 
Whitehead, A. The evolutionary radiation of diverse osmotolerant physiologies in killifish (Fundulus sp.). Evolution 64, 2070–2085 (2010).PubMed 

Google Scholar 
Lind, M. I. & Johansson, F. The degree of adaptive phenotypic plasticity is correlated with the spatial environmental heterogeneity experienced by island populations of Rana temporaria. J. Evol. Biol. 20, 1288–1297 (2007).CAS 
PubMed 

Google Scholar 
Lázaro-Nogal, A. et al. Environmental heterogeneity leads to higher plasticity in dry-edge populations of a semi-arid Chilean shrub: insights into climate change responses. J. Ecol. 103, 338–350 (2015).
Google Scholar 
Gianoli, E. Plasticity of traits and correlations in two populations of Convolvulus arvensis (Convolvulaceae) differing in environmental heterogeneity. Int. J. Plant Sci. 165, 825–832 (2004).
Google Scholar 
Fischer, E. K., Song, Y., Hughes, K. A., Zhou, W. & Hoke, K. L. Nonparallel transcriptional divergence during parallel adaptation. Mol. Ecol. 30, 1516–1530 (2021).PubMed 

Google Scholar 
Gunter, H. M., Schneider, R. F., Karner, I., Sturmbauer, C. & Meyer, A. Molecular investigation of genetic assimilation during the rapid adaptive radiations of East African cichlid fishes. Mol. Ecol. 26, 6634–6653 (2017).CAS 
PubMed 

Google Scholar 
Bitter, M. C. et al. Fluctuating selection and global change: a synthesis and review on disentangling the roles of climate amplitude, predictability and novelty. Proc. Biol. Sci. 288, 20210727 (2021).CAS 
PubMed 

Google Scholar 
Skliris, N. et al. Salinity changes in the World Ocean since 1950 in relation to changing surface freshwater fluxes. Clim. Dyn. 43, 709–736 (2014).
Google Scholar 
Collins, M. et al. Long-term climate change: projections, commitments and irreversibility. in Climate Change 2013-The Physical Science Basis: Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. 1029–1136 (Cambridge University Press, 2013).Sunday, J. M. et al. Evolution in an acidifying ocean. Trends Ecol. Evol. 29, 117–125 (2014).PubMed 

Google Scholar 
Reusch, T. B. H. & Boyd, P. W. Experimental evolution meets marine phytoplankton. Evolution 67, 1849–1859 (2013).PubMed 

Google Scholar 
Palumbi, S. R., Evans, T. G., Pespeni, M. H. & Somero, G. N. Present and future adaptation of marine species assemblages. Oceanography https://doi.org/10.5670/oceanog.2019.314 (2019).Helmuth, B. et al. Long-term, high frequency in situ measurements of intertidal mussel bed temperatures using biomimetic sensors. Sci. Data 3, 160087 (2016).MathSciNet 
PubMed 
PubMed Central 

Google Scholar 
Feely, R. A., Sabine, C. L., Hernandez-Ayon, J. M., Ianson, D. & Hales, B. Evidence for upwelling of corrosive ‘acidified’ water onto the continental shelf. Science 320, 1490–1492 (2008).ADS 
CAS 
PubMed 

Google Scholar 
Cowen, R. K. & Sponaugle, S. Larval dispersal and marine population connectivity. Ann. Rev. Mar. Sci. 1, 443–466 (2009).PubMed 

Google Scholar 
Huys, R. & Boxshall, G. A. Copepod Evolution. (marinespecies.org, 1991).Langer, J. A. F. et al. Acclimation and adaptation of the coastal calanoid copepod Acartia tonsa to ocean acidification: a long-term laboratory investigation. Mar. Ecol. Prog. Ser. 619, 35–51 (2019).ADS 
CAS 

Google Scholar 
Dam, H. G. Evolutionary adaptation of marine zooplankton to global change. Ann. Rev. Mar. Sci. 5, 349–370 (2013).PubMed 

Google Scholar 
De Wit, P., Dupont, S. & Thor, P. Selection on oxidative phosphorylation and ribosomal structure as a multigenerational response to ocean acidification in the common copepod Pseudocalanus acuspes. Evol. Appl. 9, 1112–1123 (2016).PubMed 

Google Scholar 
Thor, P. & Dupont, S. Transgenerational effects alleviate severe fecundity loss during ocean acidification in a ubiquitous planktonic copepod. Glob. Chang. Biol. 21, 2261–2271 (2015).ADS 
PubMed 

Google Scholar 
Donelson, J. M. et al. Understanding interactions between plasticity, adaptation and range shifts in response to marine environmental change. Philos. Trans. R. Soc. Lond. B Biol. Sci. 374, 20180186 (2019).PubMed 
PubMed Central 

Google Scholar 
Gibbin, E. M. et al. Can multi-generational exposure to ocean warming and acidification lead to the adaptation of life history and physiology in a marine metazoan? J. Exp. Biol. 220, 551–563 (2017).PubMed 

Google Scholar 
Mauchline, J. The Biology of Calanoid Copepods (Academic Press, 1998).Steinberg, D. K. & Landry, M. R. Zooplankton and the ocean carbon cycle. Ann. Rev. Mar. Sci. 9, 413–444 (2017).PubMed 

Google Scholar 
Gobler, C. J. & Baumann, H. Hypoxia and acidification in ocean ecosystems: coupled dynamics and effects on marine life. Biol. Lett. 12, 20150976 (2016).Rice, E., Dam, H. G. & Stewart, G. Impact of climate change on estuarine zooplankton: surface water warming in Long Island Sound is associated with changes in copepod size and community structure. Estuaries Coasts 38, 13–23 (2015).
Google Scholar 
IPCC. Climate Change 2014: Mitigation of Climate Change. Contribution of Working Group III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Vol. 1454 (IPCC, 2014).Caldeira, K. & Wickett, M. E. Oceanography: anthropogenic carbon and ocean pH. Nature 425, 365 (2003).ADS 
CAS 
PubMed 

Google Scholar 
Dam, H. G. et al. Rapid, but limited, zooplankton adaptation to simultaneous warming and acidification. Nat. Clim. Chang. 11, 780–786 (2021).ADS 

Google Scholar 
Behrenfeld, M. J. et al. Climate-driven trends in contemporary ocean productivity. Nature 444, 752–755 (2006).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Barghi, N., Hermisson, J. & Schlötterer, C. Polygenic adaptation: a unifying framework to understand positive selection. Nat. Rev. Genet. 21, 769–781 (2020).CAS 
PubMed 

Google Scholar 
Láruson, Á. J., Yeaman, S. & Lotterhos, K. E. The importance of genetic redundancy in evolution. Trends Ecol. Evol. 35, 809–822 (2020).PubMed 

Google Scholar 
Tobler, R. et al. Massive habitat-specific genomic response in D. melanogaster populations during experimental evolution in hot and cold environments. Mol. Biol. Evol. 31, 364–375 (2014).CAS 
PubMed 

Google Scholar 
Belhadj Slimen, I. et al. Reactive oxygen species, heat stress and oxidative-induced mitochondrial damage. A review. Int. J. Hyperth. 30, 513–523 (2014).CAS 

Google Scholar 
Downs, C. A. & Heckathorn, S. A. The mitochondrial small heat-shock protein protects NADH:ubiquinone oxidoreductase of the electron transport chain during heat stress in plants. FEBS Lett. 430, 246–250 (1998).CAS 
PubMed 

Google Scholar 
Harada, A. E., Healy, T. M. & Burton, R. S. Variation in thermal tolerance and its relationship to mitochondrial function across populations of Tigriopus californicus. Front. Physiol. 10, 213 (2019).PubMed 
PubMed Central 

Google Scholar 
Chung, D. J. & Schulte, P. M. Mitochondria and the thermal limits of ectotherms. J. Exp. Biol. 223 (2020).Mathew, A. N. U. & Morimoto, R. I. Role of the heat-shock response in the life and death of proteins. Ann. N. Y. Acad. Sci. 851, 99–111 (1998).ADS 
CAS 
PubMed 

Google Scholar 
Evans, T. G., Pespeni, M. H., Hofmann, G. E., Palumbi, S. R. & Sanford, E. Transcriptomic responses to seawater acidification among sea urchin populations inhabiting a natural pH mosaic. Mol. Ecol. 26, 2257–2275 (2017).CAS 
PubMed 

Google Scholar 
Bailey, A. et al. Regulation of gene expression is associated with tolerance of the Arctic copepod Calanus glacialis to CO2-acidified sea water. Ecol. Evol. 7, 7145–7160 (2017).PubMed 
PubMed Central 

Google Scholar 
Tenaillon, O. et al. The molecular diversity of adaptive convergence. Science 335, 457–461 (2012).ADS 
CAS 
PubMed 

Google Scholar 
Anjum, R. & Blenis, J. The RSK family of kinases: emerging roles in cellular signalling. Nat. Rev. Mol. Cell Biol. 9, 747–758 (2008).CAS 
PubMed 

Google Scholar 
Marshall, D. J. Transgenerational plasticity in the sea: context-dependent maternal effects across the life history. Ecology 89, 418–427 (2008).PubMed 

Google Scholar 
Vehmaa, A., Brutemark, A. & Engström-Öst, J. Maternal effects may act as an adaptation mechanism for copepods facing pH and temperature changes. PLoS ONE 7, e48538 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Skinner, M. K. What is an epigenetic transgenerational phenotype? F3 or F2. Reprod. Toxicol. 25, 2–6 (2008).CAS 
PubMed 

Google Scholar 
Sasaki, M. C. & Dam, H. G. Integrating patterns of thermal tolerance and phenotypic plasticity with population genetics to improve understanding of vulnerability to warming in a widespread copepod. Glob. Chang. Biol. 25, 4147–4164 (2019).ADS 
PubMed 

Google Scholar 
Sasaki, M. C. & Dam, H. G. Genetic differentiation underlies seasonal variation in thermal tolerance, body size, and plasticity in a short‐lived copepod. Ecol. Evol. 90, 193 (2020).
Google Scholar 
Ho, W.-C., Li, D., Zhu, Q. & Zhang, J. Phenotypic plasticity as a long-term memory easing readaptations to ancestral environments. Sci. Adv. 6, eaba3388 (2020).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Caswell, H. Matrix population models. Encyclopedia of Environmetrics 3, https://doi.org/10.1002/9781118445112.stat07481 (2006).Huey, R. B., Wakefield, T., Crill, W. D. & Gilchrist, G. W. Within- and between-generation effects of temperature on early fecundity of Drosophila melanogaster. Heredity 74, 216–223 (1995). Pt 2.PubMed 

Google Scholar 
Zwaan, B., Bijlsma, R. & Hoekstra, R. F. Direct selection on life span in Drosophila melanogaster. Evolution 49, 649–659 (1995).PubMed 

Google Scholar 
Reznick, D. A., Bryga, H. & Endler, J. A. Experimentally induced life-history evolution in a natural population. Nature 346, 357–359 (1990).ADS 

Google Scholar 
Jerison, E. R., Nguyen Ba, A. N., Desai, M. M. & Kryazhimskiy, S. Chance and necessity in the pleiotropic consequences of adaptation for budding yeast. Nat. Ecol. Evol. 4, 601–611 (2020).PubMed 
PubMed Central 

Google Scholar 
Zhong, S., Khodursky, A., Dykhuizen, D. E. & Dean, A. M. Evolutionary genomics of ecological specialization. Proc. Natl Acad. Sci. USA 101, 11719–11724 (2004).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
MacLean, R. C., Bell, G. & Rainey, P. B. The evolution of a pleiotropic fitness tradeoff in Pseudomonas fluorescens. Proc. Natl Acad. Sci. USA 101, 8072–8077 (2004).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Bettencourt, B. R., Feder, M. E. & Cavicchi, S. Experimental evolution of HSP70 expression and thermotolerance in Drosophila melanogaster. Evolution 53, 484–492 (1999).CAS 
PubMed 

Google Scholar 
Schaum, C.-E., Buckling, A., Smirnoff, N., Studholme, D. J. & Yvon-Durocher, G. Environmental fluctuations accelerate molecular evolution of thermal tolerance in a marine diatom. Nat. Commun. 9, 1719 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
Orr, H. A. Adaptation and the cost of complexity. Evolution 54, 13–20 (2000).CAS 
PubMed 

Google Scholar 
Chen, P. & Zhang, J. Antagonistic pleiotropy conceals molecular adaptations in changing environments. Nat. Ecol. Evol. 4, 461–469 (2020).PubMed 
PubMed Central 

Google Scholar 
Cox, P. M., Betts, R. A., Jones, C. D., Spall, S. A. & Totterdell, I. J. Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model. Nature 408, 184–187 (2000).ADS 
CAS 
PubMed 

Google Scholar 
Mayor, D. J., Sommer, U., Cook, K. B. & Viant, M. R. The metabolic response of marine copepods to environmental warming and ocean acidification in the absence of food. Sci. Rep. 5, 13690 (2015).ADS 
PubMed 
PubMed Central 

Google Scholar 
Pedersen, S. A. et al. Multigenerational exposure to ocean acidification during food limitation reveals consequences for copepod scope for growth and vital rates. Environ. Sci. Technol. 48, 12275–12284 (2014).ADS 
CAS 
PubMed 

Google Scholar 
Bono, L. M., Smith, L. B. Jr, Pfennig, D. W. & Burch, C. L. The emergence of performance trade-offs during local adaptation: insights from experimental evolution. Mol. Ecol. 26, 1720–1733 (2017).PubMed 

Google Scholar 
Masel, J., King, O. D. & Maughan, H. The loss of adaptive plasticity during long periods of environmental stasis. Am. Nat. 169, 38–46 (2007).PubMed 

Google Scholar 
Bay, R. A. et al. Genomic signals of selection predict climate-driven population declines in a migratory bird. Science 359, 83–86 (2018).ADS 
CAS 
PubMed 

Google Scholar 
Bay, R. A. et al. Predicting responses to contemporary environmental change using evolutionary response architectures. Am. Nat. 189, 463–473 (2017).PubMed 

Google Scholar 
Bush, A. et al. Incorporating evolutionary adaptation in species distribution modelling reduces projected vulnerability to climate change. Ecol. Lett. 19, 1468–1478 (2016).PubMed 

Google Scholar 
Valladares, F. et al. The effects of phenotypic plasticity and local adaptation on forecasts of species range shifts under climate change. Ecol. Lett. 17, 1351–1364 (2014).PubMed 

Google Scholar 
Feinberg, L. R. & Dam, H. G. Effects of diet on dimensions, density and sinking rates of fecal pellets of the copepod Acartia tonsa. Mar. Ecol. Prog. Ser. 175, 87–96 (1998).ADS 

Google Scholar 
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
Jørgensen, T. S. et al. The genome and mRNA transcriptome of the cosmopolitan calanoid copepod Acartia tonsa Dana improve the understanding of copepod genome size evolution. Genome Biol. Evol. https://doi.org/10.1093/gbe/evz067 (2019).Haas, B. J. et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc. 8, 1494–1512 (2013).CAS 

Google Scholar 
Davidson, N. M., Hawkins, A. D. K. & Oshlack, A. SuperTranscripts: a data driven reference for analysis and visualisation of transcriptomes. Genome Biol. 18, 148 (2017).PubMed 
PubMed Central 

Google Scholar 
Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Preprint at https://arxiv.org/abs/1303.3997 (2013).Faust, G. G. & Hall, I. M. SAMBLASTER: fast duplicate marking and structural variant read extraction. Bioinformatics 30, 2503–2505 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
Koboldt, D. C. et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Soneson, C., Love, M. I. & Robinson, M. D. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res. 4, 1521 (2016).Patro, R., Duggal, G., Love, M. I., Irizarry, R. A. & Kingsford, C. Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 14, 417–419 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Core Team, 2019).Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).PubMed 
PubMed Central 

Google Scholar 
Kenkel, C. D. & Matz, M. V. Gene expression plasticity as a mechanism of coral adaptation to a variable environment. Nat. Ecol. Evol. 1, 14 (2016).PubMed 

Google Scholar 
Campbell-Staton, S. C., Velotta, J. P. & Winchell, K. M. Selection on adaptive and maladaptive gene expression plasticity during thermal adaptation to urban heat islands. Nat. Commun. 12, 6195 (2021).ADS 
PubMed 
PubMed Central 

Google Scholar 
Jombart, T. & Ahmed, I. adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinformatics https://doi.org/10.1093/bioinformatics/btr521 (2011).Hadfield, J. D. MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R Package. J. Stat. Softw. 33, 1–22 (2010).
Google Scholar 
Orozco-terWengel, P. et al. Adaptation of Drosophila to a novel laboratory environment reveals temporally heterogeneous trajectories of selected alleles. Mol. Ecol. 21, 4931–4941 (2012).PubMed 
PubMed Central 

Google Scholar 
Kofler, R. et al. PoPoolation: a toolbox for population genetic analysis of next generation sequencing data from pooled individuals. PLoS ONE 6, e15925 (2011).ADS 
CAS 

Google Scholar 
Wright, R. M., Aglyamova, G. V., Meyer, E. & Matz, M. V. Gene expression associated with white syndromes in a reef building coral, Acropora hyacinthus. BMC Genomics 16, 371 (2015).PubMed 
PubMed Central 

Google Scholar 
Therneau, T. M. & Grambsch, P. M. Modeling Survival Data: Extending the Cox Model (Springer, 2013).Therneau, T. A Package for Survival Analysis in S. version 2.38. (Mayo Foundation, 2015).Kassambara, A., Kosinski, M., Biecek, P. & Fabian, S. Package ‘survminer’. Drawing Survival Curves using ‘ggplot2’. (R package version 0. 3. 1.) (2017).Houde, S. E. L. & Roman, M. R. Effects of food quality on the functional ingestion response of the copepod Acartia tonsa. Mar. Ecol. Prog. Ser. 40, 69–77 (1987).ADS 

Google Scholar 
Brennan, R. S. et al. Code repository for ‘Loss of transcriptional plasticity but sustained adaptive capacity after adaptation to global change conditions in a marine copepod’. Zenodo https://doi.org/10.5281/zenodo.5840148 (2022). More

Quantification of causal networksWe first compared the relative strengths of causal links across systems (Supplementary Fig. S3). Phytoplankton species richness was the major controlling factor for phytoplankton biomass (significant in 16 of 19 sites, Fig. 2a) in these diverse aquatic systems, consistent with experimental studies17. However, the averaged linkage strength for this effect was not significantly different from that of NO3 (i.e., BD → EF vs. NO3 → EF; permutation test P = 0.501), highlighting that nitrogen availability was equally important in affecting phytoplankton biomass in natural systems.Fig. 2: Relative strengths of various modules.Standardized linkage strengths of causal variables affecting (a) phytoplankton biomass and (b) species richness (here, BD) and loop weights for various types of (c) pairwise feedbacks and (d) triangular feedbacks. All statistics were calculated from the 19 independent sites (n = 19) and depicted as joint violins and box plots to present the empirical distribution that labels the maxima and minima at the top and bottom of the violins, respectively, and shows 25, 50, and 75% quantiles in the boxes with whiskers presenting at most 1.5 * interquartile range. The two numbers within the parentheses (S; R1) above each violin plot report the number of significant results in CCM (S; labeled blue) and the number of systems in which a particular module had the greatest strength (i.e., rank 1; R1; labeled red). Source data are provided as a Source Data file.Full size imageIn the opposite direction, phytoplankton biomass was a significant driver of phytoplankton species richness in most ecosystems (15 of 19 sites, Fig. 2b). However, NO3 more often had a stronger effect, appearing as the most important driver in 11 of 19 sites compared to phytoplankton biomass (4 of 19 sites) (Fig. 2b). Although the difference in effect strength was not significant (permutation test, P = 0.162), these results implicated nitrogen availability as an essential determinant affecting both phytoplankton diversity and biomass. As a sensitivity test, we also examined the effects of Shannon diversity. The results suggest that the importance of nutrients is robust to the use of other diversity indexes (e.g., Shannon diversity in Supplementary Fig. S4), although the causal effects from phytoplankton biomass became relatively more important compared to biomass effects on species richness (Fig. 2b). Based on these findings, we inferred that processes influencing nutrients (e.g., external loadings and internal cycling38) need to be considered when investigating aquatic biodiversity. Changes in those processes (e.g., climatic39 or anthropogenic40 driven nutrient changes) may indeed substantially impact phytoplankton biodiversity, and subsequent ecosystem functioning.The importance of NO3 uncovered in our analyses might not be a counter-intuitive result, as many systems analyzed in this study were P-rich. For instance, the average phosphate concentration was 57.5 and 41.7 μgP/L for Lake Mendota (Me) and Lake Monona (Mo) (Supplementary Table S1), respectively. In addition, there were also high total phosphorus (TP) concentrations in shallow lake systems, e.g., average TP was 106.1, 112.5, and 126.4 μgP/L in Lake Inba (Ib), Lake Kasumigaura (Ks), and Müggelsee (Mu), respectively. Phosphorus was not always a limiting factor in eutrophic and mesotrophic systems, e.g., Lake Kasumigaura41 and Lake Geneva (Gv)42. In addition, nitrogen was deficient and limited cyanobacteria bloom in Müggelsee (Mu)43. Nonetheless, we cannot exclude the possibility of colimitation44 in N and P and the possibility that P availability also depends on N45, which warrants further investigation.Apart from nutrients and temperature, the causal effects of other important drivers on phytoplankton biomass and diversity were also examined, though not in all 19 systems due to data limitation. The causal effects of physical environmental factors, such as irradiance and water column stability, were presented in Supplementary Fig. S5; the results indicated that the quantified causal strengths on average were not as strong as the effects of diversity and nutrients. Moreover, the effects of consumers (e.g., zooplankton), which have been suggested as important drivers affecting species diversity of phytoplankton communities46, were also examined. Based on our analysis of zooplankton, the causal effects of herbivorous crustaceans on phytoplankton biomass and diversity were significant in most of the analyzed systems. However, these effects were on average not as strong as the effects of phytoplankton diversity and nutrients, respectively (Supplementary Fig. S6). Nonetheless, these findings were not generalized to all 19 systems due to a lack of complete datasets as shown in Supplementary Table S3, and thus warrant more detailed investigation in future studies.In addition to individual causal effects, we investigated feedbacks across systems. Pairwise feedbacks (e.g., BD ↔ EF and NO3 ↔ EF) were common (Fig. 2c). However, the averaged linkage strength was often stronger in one direction when involving BD (Fig. 3). Specifically, the average strength of BD → EF was stronger than for the opposite direction of EF → BD (permutation test P = 0.015); BD → EF was stronger than EF → BD in 14 of the 19 systems (Fig. 3). In addition, biodiversity effects on nutrients (BD → NO3 and BD → PO4) were also stronger than their reversed effects (NO3 → BD and PO4 → BD) in 12 and 13 systems, respectively. In comparison, the interactions between nutrients and productivity were more symmetrical: nutrient effects on biomass (NO3 → EF and PO4 → EF) were stronger than biomass effects on nutrients (EF → NO3 and EF → PO4) in only 9 and 8 of 19 systems, respectively. These results supported the previous findings8 that biodiversity effects more often operate at short-term scales, which makes effects more observable in our monthly-scale analyses than feedback effects on diversity, which are expected to occur on a more prolonged timescale, e.g., through slowly changing nutrient cycling31 or decomposition47. Nevertheless, the timescale dependence of causal interactions in ecosystem networks is a topic that needs further study.Fig. 3: Directional bias in pairwise feedbacks.The difference in standardized linkage strengths between the two directions was computed for each pairwise feedback and depicted as joint violin and box plots. All statistics were calculated from the 19 independent sites (n = 19) and depicted as joint violins and box plots to present the empirical distribution that labels the maxima and minima at the top and bottom of the violins, respectively, and shows 25, 50, and 75% quantiles in the boxes with whiskers presenting at most 1.5 * interquartile range. The number above the plot indicates the number of systems with a positive difference in linkage strength. For example, BD → EF was stronger than its feedback, EF → BD, in 14 of the systems. In general, the strength of diversity effects (BD → EF, BD → NO3, BD → PO4) was usually stronger than feedback effects (EF → BD, NO3 → BD, PO4 → BD). Source data are provided as a Source Data file.Full size imageSubsequently, we quantified the strengths of pairwise feedbacks as the geometric mean of the linkage strengths in each direction, following a previous study9 (see more details in Methods). Among these feedbacks (Fig. 2c and Supplementary Fig. S7), BD ↔ NO3 had the highest median and average strength (0.78 and 0.68, respectively) across systems. However, strengths of BD ↔ NO3 were highly variable among systems (large interquartile range in Fig. 2c), and thus were only significant in 11 of 19 systems, compared to BD ↔ EF (15 of 19 systems). These findings reinforced the importance of nutrients as key determinants for aquatic biodiversity and implied that nutrient effects are context-dependent. In other words, BD ↔ NO3 was less common than BD ↔ EF across systems, despite its stronger average strength. The prevalence of BD ↔ EF indicated a need for more long-term experiments and process-based/theoretical modeling accounting for bidirectional interactions between diversity and biomass16, because bidirectional interactions and feedbacks may challenge our simple predictions for ecosystem dynamics, based on knowledge of unidirectional interactions30.Quantification of the causal network also allowed us to analyze triangular feedbacks. Within the conceptual framework of Fig. 1b, there are four kinds of triangular feedbacks involving biodiversity, ecosystem functioning, and either nitrate or phosphate (Type I: BD → EF → NO3 and BD → EF → PO4; Type II: EF → BD → NO3 and EF → BD → PO4). There was at least one significant triangular feedback in 14 of 19 sites (Fig. 2d). More specifically, NO3-associated feedbacks (Type I-N and Type II-N) were usually stronger than PO4-associated feedbacks (Type I-P and Type II-P) (Fig. 2d), although the difference in strength among the four types of feedbacks was not significant (Fig. 2d; Kruskal–Wallis test, P = 0.59). The dominance of NO3-associated feedbacks in our study was attributed to many of the sites being marine and eutrophic lakes, which are likely to be N-limited due to an imbalance in external loadings48 or strong denitrification49. Among both NO3- and PO4-associated feedbacks, there were no significant differences in strength between Type I and Type II feedbacks (Supplementary Fig. S7), suggesting that biodiversity can directly influence biomass (Type I), as well as through a pathway that involves endogenous nutrient variables (Type II) and eventually feeds back on itself.Causal networks under environmental contextsOur empirical analyses revealed state dependency of the causal links and feedbacks among biodiversity, biomass, and environmental factors in natural systems; that is, their strengths were highly dependent on the state of other variables. Based on a cross-system comparison (Methods), strengths of individual links (e.g., BD → EF), pairwise feedbacks (e.g., BD ↔ EF), and triangular feedbacks (e.g., BD → EF → NO3 → BD) varied systematically, depending on environmental characteristics (Fig. 4 and Supplementary Fig. S8). Ecosystems with higher species diversity (long-term average species richness) and lower average PO4 concentrations had stronger BD → EF links (Fig. 4a; correlation coefficient r = 0.600 and −0.513; P = 0.007 and 0.025 for species diversity and PO4, respectively). These results were further confirmed by stepwise regression, indicating that the ecosystems characterized by higher diversity, lower average temperature, and oligotrophic conditions had stronger BD → EF (best-fit regression model: BD → EF strength = 0.663 + 0.171*BD − 0.139*T − 0.096*PO4; F3, 15 = 9.958 and P  More

Partensky F, Garczarek L. Prochlorococcus: advantages and limits of minimalism. Ann Rev Mar Sci. 2010;2:305–31.PubMed 

Google Scholar 
Flombaum P, Gallegos JL, Gordillo RA, Rincón J, Zabala LL, Jiao N, et al. Present and future global distributions of the marine cyanobacteria Prochlorococcus and Synechococcus. Proc Natl Acad Sci. 2013;110:9824–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
Follows MJ, Dutkiewicz S, Grant S, Chisholm SW. Emergent biogeography of microbial communities in a model ocean. Science. 2007;315:1843–6.CAS 
PubMed 

Google Scholar 
Biller SJ, Berube PM, Lindell D, Chisholm SW. Prochlorococcus: the structure and function of collective diversity. Nat Rev Microbiol. 2014;13:13–27.PubMed 

Google Scholar 
Farrant GK, Doré H, Cornejo-Castillo FM, Partensky F, Ratin M, Ostrowski M, et al. Delineating ecologically significant taxonomic units from global patterns of marine picocyanobacteria. Proc Natl Acad Sci. 2016;113:E3365–74.CAS 
PubMed 
PubMed Central 

Google Scholar 
Kashtan N, Roggensack SE, Rodrigue S, Thompson JW, Biller SJ, Coe A, et al. Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus. Science. 2014;344:416–20.CAS 
PubMed 

Google Scholar 
Tagliabue A, Bowie AR, Boyd PW, Buck KN, Johnson KS, Saito MA. The integral role of iron in ocean biogeochemistry. Nature. 2017;543:51–9.CAS 
PubMed 

Google Scholar 
Gledhill M, van den Berg CMG. Determination of complexation of iron(III) with natural organic complexing ligands in seawater using cathodic stripping voltammetry. Mar Chem. 1994;47:41–54.CAS 

Google Scholar 
Rue EL, Bruland KW. The role of organic complexation on ambient iron chemistry in the equatorial Pacific Ocean and the response of a mesoscale iron addition experiment. Limnol Oceanogr 1997;42:901–10.Hassler CS, van den Berg CMG, Boyd PW. Toward a regional classification to provide a more inclusive examination of the ocean biogeochemistry of iron-binding ligands. Front Mar Sci. 2017;4:1–19.
Google Scholar 
Gledhill M, Buck KN. The organic complexation of iron in the marine environment: a review. Front Microbiol. 2012;3:1–17.
Google Scholar 
Bundy RM, Boiteau RM, McLean C, Turk-Kubo KA, McIlvin MR, Saito MA, et al. Distinct siderophores contribute to Iron cycling in the mesopelagic at station ALOHA. Front. 2018;5:1–15.
Google Scholar 
Boiteau RM, Mende DR, Hawco NJ, McIlvin MR, Fitzsimmons JN, Saito MA, et al. Siderophore-based microbial adaptations to iron scarcity across the eastern Pacific Ocean. Proc Natl Acad Sci. 2016;113:14237–42.Shaked Y, Lis H. Disassembling iron availability to phytoplankton. Front Microbiol. 2012;3:123.PubMed 
PubMed Central 

Google Scholar 
Morrissey J, Bowler C. Iron utilization in marine cyanobacteria and eukaryotic algae. Front Microbiol. 2012;3:43.PubMed 
PubMed Central 

Google Scholar 
Hogle SL, Cameron Thrash J, Dupont CL, Barbeau KA. Trace metal acquisition by marine heterotrophic bacterioplankton with contrasting trophic strategies. Appl Environ Microbiol. 2016;82:1613–24.CAS 
PubMed 
PubMed Central 

Google Scholar 
Hopkinson B, Barbeau K. Iron transporters in marine prokaryotic genomes and metagenomes. Environ Microbiol. 2012;14:114–28.CAS 
PubMed 

Google Scholar 
Webb EA, Moffett JW, Waterbury JB. Iron stress in open-ocean Cyanobacteria (Synechococcus, Trichodesmium, and Crocosphaera spp.): Identification of the IdiA protein. Appl Environ Microbiol. 2001;67:5444–52.CAS 
PubMed 
PubMed Central 

Google Scholar 
Hopkinson BM, Morel FMM. The role of siderophores in iron acquisition by photosynthetic marine microorganisms. Biometals. 2009;22:659–69.CAS 
PubMed 

Google Scholar 
Malmstrom RR, Rodrigue S, Huang KH, Kelly L, Kern SE, Thompson A, et al. Ecology of uncultured Prochlorococcus clades revealed through single-cell genomics and biogeographic analysis. ISME J. 2013;7:184–98.CAS 
PubMed 

Google Scholar 
Ustick LJ, Larkin AA, Garcia CA, Garcia NS, Brock ML, Lee JA, et al. Metagenomic analysis reveals global-scale patterns of ocean nutrient limitation. Science. 2021;372:287–91.CAS 
PubMed 

Google Scholar 
Garcia CA, Hagstrom GI, Larkin AA, Ustick LJ, Levin SA, Lomas MW, et al. Linking regional shifts in microbial genome adaptation with surface ocean biogeochemistry. Philos Trans R Soc Lond B Biol Sci. 2020;375:20190254.CAS 
PubMed 
PubMed Central 

Google Scholar 
Hogle SL, Dupont CL, Hopkinson BM, King AL, Buck KN, Roe KL, et al. Pervasive iron limitation at subsurface chlorophyll maxima of the California Current. Proc Natl Acad Sci. 2018;115:13300–5.CAS 
PubMed 
PubMed Central 

Google Scholar 
Hawco NJ, Fu F, Yang N, Hutchins DA, John SG. Independent iron and light limitation in a low-light-adapted Prochlorococcus from the deep chlorophyll maximum. ISME J. 2020;15:359–62.PubMed 
PubMed Central 

Google Scholar 
Biller SJ, Berube PM, Dooley K, Williams M, Satinsky BM, Hackl T, et al. Marine microbial metagenomes sampled across space and time. Sci Data. 2018;5:180176.CAS 
PubMed 
PubMed Central 

Google Scholar 
Pesant S, Not F, Picheral M, Kandels-Lewis S, Le Bescot N, Gorsky G, et al. Open science resources for the discovery and analysis of Tara Oceans data. Sci Data. 2015;2:150023.CAS 
PubMed 
PubMed Central 

Google Scholar 
Berube PM, Biller SJ, Hackl T, Hogle SL, Satinsky BM, Becker JW, et al. Single cell genomes of Prochlorococcus, Synechococcus, and sympatric microbes from diverse marine environments. Sci Data. 2018;5:180154.CAS 
PubMed 
PubMed Central 

Google Scholar 
Biller SJ, Berube PM, Berta-Thompson JW, Kelly L, Roggensack SE, Awad L, et al. Genomes of diverse isolates of the marine cyanobacterium Prochlorococcus. Sci Data. 2014;1:140034.CAS 
PubMed 
PubMed Central 

Google Scholar 
Karsenti E, Acinas SG, Bork P, Bowler C, De Vargas C, Raes J, et al. A holistic approach to marine eco-systems biology. PLoS Biol. 2011;9:e1001177.CAS 
PubMed 
PubMed Central 

Google Scholar 
Schlitzer R, Anderson RF, Dodas EM, Lohan M, Geibert W, Tagliabue A, et al. The GEOTRACES intermediate data product 2017. Chem Geol. 2018;493:210–23.CAS 

Google Scholar 
Salt LA, van Heuven SMAC, Claus ME, Jones EM, de Baar HJW. Rapid acidification of mode and intermediate waters in the southwestern Atlantic Ocean. Biogeosciences. 2015;12:1387–401.
Google Scholar 
Rijkenberg MJA, Middag R, Laan P, Gerringa LJA, van Aken HM, Schoemann V, et al. The distribution of dissolved iron in the West Atlantic Ocean. PLoS One. 2014;9:e101323.PubMed 
PubMed Central 

Google Scholar 
Wyatt NJ, Milne A, Woodward EMS, Rees AP, Browning TJ, Bouman HA, et al. Biogeochemical cycling of dissolved zinc along the GEOTRACES South Atlantic transect GA10 at 40°S: Dissolved zinc in the Atlantic at 40o S. Glob Biogeochem Cycles. 2014;28:44–56.CAS 

Google Scholar 
Ashkezari MD, Hagen NR, Denholtz M, Neang A, Burns TC, Morales RL, et al. Simons collaborative marine atlas project (Simons CMAP): an open‐source portal to share, visualize, and analyze ocean data. Limnol Oceanogr Methods. 2021;19:488–96.
Google Scholar 
Acker M, Hogle SL, Berube PM, Hackl T, Stepanauskas R, Chisholm SW, et al. Phosphonate production by marine microbes: exploring new sources and potential function. Proc Natl Acad Sci. 2022;119:e2113386119.Becker JW, Hogle SL, Rosendo K, Chisholm SW. Co-culture and biogeography of Prochlorococcus and SAR11. ISME J. 2019;13:1506–19.CAS 
PubMed 
PubMed Central 

Google Scholar 
Malmstrom RR, Coe A, Kettler GC, Martiny AC, Frias-Lopez J, Zinser ER, et al. Temporal dynamics of Prochlorococcus ecotypes in the Atlantic and Pacific oceans. ISME J. 2010;4:1252–64.PubMed 

Google Scholar 
Fitzsimmons JN, Hayes CT, Al-Subiai SN, Zhang R, Morton PL, Weisend RE, et al. Daily to decadal variability of size-fractionated iron and iron-binding ligands at the Hawaii Ocean Time-series Station ALOHA. Geochim Cosmochim Acta. 2015;171:303–24.CAS 

Google Scholar 
Buck KN, Sohst B, Sedwick PN. The organic complexation of dissolved iron along the U.S. GEOTRACES (GA03) North Atlantic Section. Deep Sea Res Part 2. 2015;116:152–65.CAS 

Google Scholar 
Mawji E, Gledhill M, Milton JA, Tarran GA, Ussher S, Thompson A, et al. Hydroxamate siderophores: occurrence and importance in the Atlantic Ocean. Environ Sci Technol. 2008;42:8675–80.CAS 
PubMed 

Google Scholar 
Buck KN, Bruland KW. The physicochemical speciation of dissolved iron in the Bering Sea, Alaska. Limnol Oceanogr. 2007;52:1800–8.CAS 

Google Scholar 
Boiteau RM, Fitzsimmons JN, Repeta DJ, Boyle EA. Detection of iron ligands in seawater and marine cyanobacteria cultures by high-performance liquid chromatography-inductively coupled plasma-mass spectrometry. Anal Chem. 2013;85:4357–62.CAS 
PubMed 

Google Scholar 
Huerta-Cepas J, Szklarczyk D, Forslund K, Cook H, Heller D, Walter MC, et al. eggNOG 4.5: a hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences. Nucleic Acids Res. 2016;44:D286–93.CAS 
PubMed 

Google Scholar 
Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, von Mering C, et al. Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper. Mol Biol Evol. 2017;34:2115–22.CAS 
PubMed 
PubMed Central 

Google Scholar 
Blin K, Shaw S, Steinke K, Villebro R, Ziemert N, Lee SY, et al. antiSMASH 5.0: updates to the secondary metabolite genome mining pipeline. Nucleic Acids Res. 2019;47:W81–W87.CAS 
PubMed 
PubMed Central 

Google Scholar 
Menzel P, Ng KL, Krogh A. Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat Commun. 2016;7:11257.CAS 
PubMed 
PubMed Central 

Google Scholar 
Wright M, Ziegler A. ranger: A fast implementation of random forests for high dimensional data in C++ and R. J Stat Softw, Artic. 2017;77:1–17.
Google Scholar 
Jolliffe IT. Principal Component Analysis. 1986. Springer-Verlag.Kursa M, Rudnicki W. Feature selection with the Boruta package. J Stat Softw, Artic. 2010;36:1–13.
Google Scholar 
Milligan AJ, Harrison PJ. Effects of non-steady-state iron limitation on nitrogen assimilatory enzymes in the marine diatom Thalassiosira weissflogii (Bacillariophyceae). J Phycol. 2000;36:78–86.CAS 

Google Scholar 
Lomas MW, Lipschultz F. Forming the primary nitrite maximum: Nitrifiers or phytoplankton? Limnol Oceanogr. 2006;51:2453–67.CAS 

Google Scholar 
Ahlgren NA, Belisle BS, Lee MD. Genomic mosaicism underlies the adaptation of marine Synechococcus ecotypes to distinct oceanic iron niches. Environ Microbiol. 2019;22:1801–15.PubMed 

Google Scholar 
Martiny AC, Coleman ML, Chisholm SW. Phosphate acquisition genes in Prochlorococcus ecotypes: evidence for genome-wide adaptation. Proc Natl Acad Sci. 2006;103:12552–7.CAS 
PubMed 
PubMed Central 

Google Scholar 
Coleman ML, Chisholm SW. Ecosystem-specific selection pressures revealed through comparative population genomics. Proc Natl Acad Sci. 2010;107:18634–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
Berube PM, Rasmussen A, Braakman R, Stepanauskas R, Chisholm SW. Emergence of trait variability through the lens of nitrogen assimilation in Prochlorococcus. eLife. 2019;8:e41043.PubMed 
PubMed Central 

Google Scholar 
Olgun N, Duggen S, Croot PL, Delmelle P, Dietze H, Schacht U, et al. Surface ocean iron fertilization: The role of airborne volcanic ash from subduction zone and hot spot volcanoes and related iron fluxes into the Pacific Ocean. Global Biogeochem Cycles. 2011;25:GB4001.Manck LE, Park J, Tully BJ, Poire AM, Bundy RM, Dupont CL, et al. Petrobactin, a siderophore produced by Alteromonas, mediates community iron acquisition in the global ocean. ISME J. 2021;16:358–69.Mackey KRM, Post AF, McIlvin MR, Saito MA. Physiological and proteomic characterization of light adaptations in marine Synechococcus. Environ Microbiol. 2017;19:2348–65.CAS 
PubMed 

Google Scholar 
Raven JA. Predictions of Mn and Fe use efficiencies of phototrophic growth as a function of light availability for growth and of C assimilation pathway. N. Phytol. 1990;116:1–18.CAS 

Google Scholar 
Sunda WG, Huntsman SA. Interrelated influence of iron, light and cell size on marine phytoplankton growth. Nature. 1997;390:389–92.CAS 

Google Scholar 
Hawco NJ, Barone B, Church MJ, Babcock-Adams L, Repeta DJ, Wear EK, et al. Iron depletion in the deep chlorophyll maximum: Mesoscale eddies as natural iron fertilization experiments. Global Biogeochem Cycles. 2021;35:e2021GB007112.Barbeau K, Rue EL, Bruland KW, Butler A. Photochemical cycling of iron in the surface ocean mediated by microbial iron(III)-binding ligands. Nature. 2001;413:409–13.CAS 
PubMed 

Google Scholar 
van den Engh GJ, Doggett JK, Thompson AW, Doblin MA, Gimpel CNG, Karl DM. Dynamics of Prochlorococcus and Synechococcus at Station ALOHA revealed through flow cytometry and high-resolution vertical sampling. Front Mar Sci. 2017;4:1–14.
Google Scholar 
Karl DM, Church MJ. Microbial oceanography and the Hawaii Ocean Time-series programme. Nat Rev Microbiol. 2014;12:699–713.CAS 
PubMed 

Google Scholar 
Neuer S, Davenport R, Freudenthal T, Wefer G, Llinás O, Rueda M-J, et al. Differences in the biological carbon pump at three subtropical ocean sites. Geophys Res Lett 2002;29:32–1.Giovannoni SJ, Vergin KL. Seasonality in ocean microbial communities. Science. 2012;335:671–6.CAS 
PubMed 

Google Scholar 
Jickells T, Moore CM. The importance of atmospheric deposition for ocean productivity. Annu Rev Ecol Evol Syst. 2015;46:481–501.
Google Scholar 
Rii YM, Karl DM, Church MJ. Temporal and vertical variability in picophytoplankton primary productivity in the North Pacific Subtropical Gyre. Mar Ecol Prog Ser. 2016;562:1–18.CAS 

Google Scholar 
Boyle EA, Bergquist BA, Kayser RA, Mahowald N. Iron, manganese, and lead at Hawaii Ocean Time-series station ALOHA: Temporal variability and an intermediate water hydrothermal plume. Geochim Cosmochim Acta. 2005;69:933–52.CAS 

Google Scholar 
Caprara S, Buck KN, Gerringa LJA, Rijkenberg MJA, Monticelli D. A compilation of iron speciation data for open oceanic waters. Front Mar Sci. 2016;3:1–7.
Google Scholar 
Gledhill M, Gerringa LJA. The effect of metal concentration on the parameters derived from complexometric titrations of trace elements in seawater: A model study. Front Mar Sci. 2017;4:1–15.
Google Scholar 
Jickells TD, Baker AR, Chance R. Atmospheric transport of trace elements and nutrients to the oceans. Philos Trans R Soc Lond A. 2016;374:20150286.
Google Scholar 
Sedwick PN, Church TM, Bowie AR, Marsay CM, Ussher SJ, Achilles KM, et al. Iron in the Sargasso Sea (Bermuda Atlantic Time-series Study region) during summer: Eolian imprint, spatiotemporal variability, and ecological implications. Glob Biogeochem Cycles. 2005;19:GB4006.
Google Scholar 
Ohnemus DC, Rauschenberg S, Cutter GA, Fitzsimmons JN, Sherrell RM, Twining BS. Elevated trace metal content of prokaryotic communities associated with marine oxygen deficient zones. Limnol Oceanogr. 2016;62:3–25.Browning TJ, Achterberg EP, Rapp I, Engel A, Bertrand EM, Tagliabue A, et al. Nutrient co-limitation at the boundary of an oceanic gyre. Nature. 2017;551:242–6.Louropoulou E, Gledhill M, Achterberg EP, Browning TJ, Honey DJ, Schmitz RA, et al. Heme b distributions through the Atlantic Ocean: evidence for ‘anemic’ phytoplankton populations. Sci Rep. 2020;10:4551.CAS 
PubMed 
PubMed Central 

Google Scholar 
Mark Moore C. Diagnosing oceanic nutrient deficiency. Philos Trans R Soc Lond A. 2016;374:20150290.
Google Scholar 
Moore CM, Mills MM, Arrigo KR, Berman-Frank I, Bopp L, Boyd PW, et al. Processes and patterns of oceanic nutrient limitation. Nat Geosci. 2013;6:701–10.CAS 

Google Scholar 
Moore JK, Doney SC, Lindsay K. Upper ocean ecosystem dynamics and iron cycling in a global three-dimensional model. Global Biogeochem Cycles 2004;18.Rafter PA, Sigman DM, Mackey KRM. Recycled iron fuels new production in the eastern equatorial Pacific Ocean. Nat Commun. 2017;8:1100.PubMed 
PubMed Central 

Google Scholar 
Boyd PW, Ellwood MJ, Tagliabue A, Twining BS. Biotic and abiotic retention, recycling and remineralization of metals in the ocean. Nat Geosci. 2017;10:167–73.CAS 

Google Scholar 
Ferreira D, Cessi P, Coxall HK, de Boer A, Dijkstra HA, Drijfhout SS, et al. Atlantic-Pacific asymmetry in deep water formation. Annu Rev Earth Planet Sci. 2018;46:327–52.CAS 

Google Scholar 
Rusch DB, Martiny AC, Dupont CL, Halpern AL, Venter JC. Characterization of Prochlorococcus clades from iron-depleted oceanic regions. Proc Natl Acad Sci. 2010;107:16184–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
Dutkiewicz S, Ward BA, Monteiro F, Follows MJ. Interconnection of nitrogen fixers and iron in the Pacific Ocean: Theory and numerical simulations. Glob Biogeochem Cycles. 2012;26:GB1012.
Google Scholar 
Resing JA, Sedwick PN, German CR, Jenkins WJ, Moffett JW, Sohst BM, et al. Basin-scale transport of hydrothermal dissolved metals across the South Pacific Ocean. Nature. 2015;523:200–3.CAS 
PubMed 

Google Scholar 
Sela I, Wolf YI, Koonin EV. Theory of prokaryotic genome evolution. Proc Natl Acad Sci. 2016;113:11399–407.CAS 
PubMed 
PubMed Central 

Google Scholar 
Zakem EJ, Al-Haj A, Church MJ, van Dijken GL, Dutkiewicz S, Foster SQ, et al. Ecological control of nitrite in the upper ocean. Nat Commun. 2018;9:1206.PubMed 
PubMed Central 

Google Scholar 
Lauderdale JM, Braakman R, Forget G, Dutkiewicz S, Follows MJ. Microbial feedbacks optimize ocean iron availability. Proc Natl Acad Sci. 2020;117:4842–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
Hogle SL, Barbeau KA, Gledhill M. Heme in the marine environment: from cells to the iron cycle. Metallomics. 2014;6:1107–20.CAS 
PubMed 

Google Scholar 
Hogle SL, Brahamsha B, Barbeau KA. Direct heme uptake by phytoplankton-associated Roseobacter bacteria. mSystems. 2017;2:e00124–16.CAS 
PubMed 
PubMed Central 

Google Scholar 
Coale KH, Fitzwater SE, Michael Gordon R, Johnson KS, Barber RT. Control of community growth and export production by upwelled iron in the equatorial Pacific Ocean. Nature. 1996;379:621–4.CAS 

Google Scholar 
Bressac M, Guieu C, Ellwood MJ, Tagliabue A, Wagener T, Laurenceau EC, et al. Resupply of mesopelagic dissolved iron controlled by particulate iron composition. Nat Geosci. 2019;12:995–1000.CAS 

Google Scholar 
Basu S, Gledhill M, de Beer D, Prabhu Matondkar SG, Shaked Y. Colonies of marine cyanobacteria Trichodesmium interact with associated bacteria to acquire iron from dust. Commun Biol. 2019;2:284.PubMed 
PubMed Central 

Google Scholar  More

Study site and contextAmbatovy is a very large nickel, cobalt and ammonium sulphate mine in central-eastern Madagascar owned by a consortium of international mining companies50. It represents the largest ever foreign investment in the country24 (US$8 billion by 201650) and a substantial source of fiscal income49. In 2018, the company contributed ~US$50 million in taxes, tariffs, royalties and other payments49 and employed >9,000 people (93% of whom were Malagasy)51. Commercial production began in January 201424 (Supplementary Fig. 1). As key components in batteries, supplies of nickel and cobalt are critical to the green energy transition and demand for these metals is predicted to increase notably in future52.The mining concession covers an area of 7,700 ha located in the eastern rainforests of Madagascar (Fig. 1) which have very high levels of biodiversity and endemism53,54. After avoidance and minimization measures were applied (Supplementary Methods) the mine was predicted to clear or substantially degrade 2,064 ha of high-quality natural forest at the mine footprint and upper pipeline24. Any impacts on plantations or secondary habitat are not included in this estimate. Losses at the impact site were not discounted in relation to a background rate of decline, meaning that the company took responsibility for the full area of forest lost25. Independent verification by our team (by measuring the size of the mine footprint on Google Earth) confirms the extent of forest loss at the mine footprint (Supplementary Fig. 2). Clearance of the footprint accounts for most of the forest loss associated with the mine as losses associated with the pipeline are small54.Ambatovy aims to generate biodiversity gains to offset the mine-induced losses by slowing deforestation driven by shifting agriculture elsewhere26. To this end the company designated four sites, totalling 28,740 ha, to be protected as biodiversity offsets; Ankerana, Corridor Forestier Analamay-Mantadia (CFAM), the Conservation Zone and Torotorofotsy54 (Fig. 1). The offsets are considered like-for-like30 and were selected on the basis of similarity to the impact site in terms of forest structure and type, geology, climate and altitude24. The large combined area of the offsets relative to the impacted area was designed to allow flexibility, account for uncertainty and incorporate as many of the affected biodiversity components as possible24. Ankerana is the flagship offset, selected on the basis of its size, connectivity to the Corridor Ankeniheny-Zahamena (CAZ) forest corridor and the presence of ultramafic outcrops thought to support the same rare type of azonal forest lost at the mine site54. Extensive surveys conducted within Ankerana to establish biological similarity concluded the offset to be of higher conservation significance than the forests of the mine site due to the presence of rare lowland tropical forest24.The Conservation Zone is directly managed by the company, given its location within the concession area, while the other offsets are managed in partnership with local and international NGOs24,25. Ambatovy funds the management of Ankerana by Conservation International and local NGO partners (although before 2015 Ankerana was directly managed by Ambatovy via a Memorandum of Understanding with Conservation International24), supports BirdLife partner Asity with the management of Torotorofotsy and a number of local NGOs including Voary Voakajy25 are involved in CFAM26. The company is also working to secure formal, legal protection for CFAM26 as part of a proposed Torotorofotsy–CFAM complex new protected area (although progress on this has stalled).Overview of methodsTo estimate the impact of the offsets on deforestation and determine whether this has prevented enough deforestation to offset forest loss at the mine site, we combined several complementary methods for robust impact evaluation. First, we used statistical matching to match a sample of pixels from each biodiversity offset to pixels from the wider forested landscape with similar exposure to drivers of deforestation. Then we used a site-based difference-in-differences regression for each matched offset–control sample and a fixed-effects panel regression on the pooled data, to estimate the effect of protection. We systematically explored how arbitrary modelling choices (including the statistical distance measure used in matching, caliper size, ratio of control to treated units, matching with or without replacement and which, if any, additional covariates were included) affected our inference, exploring the robustness of our results to 116 alternative model specifications.MatchingThe former province of Toamasina was selected as the geographic area from which control pixels were sampled as it encompasses forests of the same type as the concession area with varying degrees of intactness and accessibility. The four biodiversity offsets are located within this province (Fig. 1).The unit of analysis is a 30 × 30 m2 pixel that was forested in the baseline year 200045,55. It is important that the scale of analysis aligns with the scale at which the drivers of deforestation (in this case, small-scale shifting agriculture) operate56. The median agricultural plot size (from 564 measured plots) in the study region is ~36 × 36 m2 (ref. 57). We took a subsample of pixels to reduce computational effort while maintaining the capacity for robust statistical inference58,59. We used a grid-based sampling strategy ensuring a minimum distance between sample units to reduce spatial autocorrelation60 and equal coverage of the study area58. A 150 × 150 m2 resolution grid, aligned to the other 30-m resolution data layers (Fig. 1c), was overlaid on the province and the 30 × 30 m2 pixel at the centre of each grid square was extracted to produce a subsample of pixels that are 120 m away from their nearest neighbour. The 120 m is larger than the minimum distance between units used in another matching study in Madagascar (68 m; ref. 59) but smaller than that used in other studies (200 m; ref. 61) and so strikes an appropriate balance between the avoidance of spatial autocorrelation and maximizing the possible sample cells.Protected areas in the study area managed by Madagascar National Parks were excluded from our control sample as they are actively managed and therefore do not represent counterfactual outcomes for the biodiversity offsets in the absence of protection (Fig. 1). However, control pixels were sampled from within the CAZ new protected area as legal protection was only granted in 2015 and resources for management are limited and thinly spread62. Additionally, Ankerana and parts of CFAM overlap with the CAZ and would have experienced the same management, and likely trajectory, as the rest of the CAZ, had they not been designated biodiversity offsets. Areas within 10 km of an offset boundary were excluded from the control sample to reduce the chance of leakage (where pressures are displaced rather than avoided) biasing results17,29. The 10 km was selected as it is a commonly used buffer zone within the literature17,58.To test for leakage effects, we used Veronoi polygons to partition the buffer area for CFAM, the Conservation Zone and Torotorofotsy (which overlap) into three individual buffer areas according to the nearest offset centroid and took a subsample of pixels from each (Fig. 1). Areas that overlapped with the established protected areas of Mantadia National Park and Analamazotra Special Reserve were excluded from the buffer zones.The outcome variable is the annual deforestation rate sourced from the Global Forest Change (GFC) dataset34. Following Vieilledent et al.45 these data were restricted to only include pixels classed as forest in a forest cover map of Madagascar for the year 200045,55, reducing the probability of false positives (whereby tree loss is identified in pixels that were not forested). The resulting tree loss raster was snapped to the forest cover 2000 layer to align cells, resulting in a maximum spatial error of 15 m. The GFC product34 has been shown to perform reasonably well at detecting deforestation in humid tropical forests63. In the north-eastern rainforests of Madagascar, Burivalova et al.39 found that GFC data performed comparably to a local classification of very high resolution satellite imagery at detecting forest clearance for shifting agriculture (although it was not effective at detecting forest degradation from selective logging). As clearance for shifting agriculture is considered the principal agent of deforestation in the study area22 and the forests of the study area are tropical humid ( >75% canopy cover), the GFC data are an appropriate tool for quantifying forest loss. Although recent evidence suggests that GFC data may have temporal biases64, this phenomenon affects our control and treated samples equally and so is unlikely to impact our results.The choice of covariates is extremely important in matching analyses. They must include, or proxy, all important factors influencing selection to treatment and the outcome of interest so that the matched control sample is sufficiently similar to the treated sample in these characteristics to constitute a plausible counterfactual, otherwise the resulting estimates may not be valid33. On the basis of the literature and a local theory of change we selected five covariates that we believe capture or proxy for the aspects of accessibility, demand and agricultural suitability that drive deforestation in the study area22,59,65,66. These are slope, elevation, distance to main road, distance to forest edge and distance to deforestation (Supplementary Methods). These five essential covariates comprise the main matching specification and form the core set used in all alternative specifications that we tested in the robustness checks. We also defined five additional variables (annual precipitation, distance to river, distance to cart track, distance to settlement and population density) and tested the effect of including these in the robustness checks. The additional covariates were so defined because they were of poorer data quality (population density and distance to settlement), correlated with an essential variable (annual precipitation and population density) or simply considered less influential (distance to river and distance to cart track; Supplementary Methods).Statistical matching was conducted in R statistics using the MatchIt package v.4.1 (ref. 67). To improve efficiency and produce closer matches we cleaned the data before matching to remove control units with values outside the calipers of the treated sample in any of the essential covariates (see Supplementary Methods for caliper definition). Following the recommendations of Schleicher et al.68 we tested several matching specifications and selected the one that maximized the trade-off between the number of treated units matched and the closeness of matches as the main specification (Supplementary Table 7). This was 1:1 nearest-neighbour matching without replacement, using Mahalanobis distance and a caliper of 1 s.d. This specification produced acceptable matches (within 1 s.d. of the Mahalanobis distance) for all treated units within all offsets. The maximum postmatching standardized difference in mean covariate values between treated and control samples was 0.05, well below the threshold of 0.25 considered to constitute an acceptable match69. This indicates that, on average, treated and control units were very well matched across all covariates.Matching was run separately for each offset. The resulting matched datasets were aggregated by treated status (offset or control) and year to produce a matrix of the count of pixels that were deforested each year (2001–2019) in the offset and the matched control sample. Converting the outcome variable to a continuous measure of deforestation avoids the problem of attrition associated with binary measures of deforestation and is better suited to the framework of the subsequent regressions70.Robustness checksStatistical matching requires various choices to be made68, many of which are essentially arbitrary. There therefore exists a range of possible alternative specifications that are all a priori valid (although some may be better suited to the data and study objectives69) but which could influence the results20,28. We tested the robustness of our results to 116 different matching model specifications (Fig. 4). First, we tested the robustness of the estimates to the use of three alternative matching distance measures (Mahalanobis, standard propensity score matching using generalized linear model regressions with a logit distribution and propensity score matching using RandomForest), three different calipers (0.25, 0.5 and 1 s.d.), different ratios of control to treated units (one, five and ten nearest neighbours) and matching with/without replacement. Holding the choice of covariates constant (using only the essential covariates), the combination of these led to the estimation of 54 different models. Second, we tested the robustness of results to the inclusion of the five additional covariates. Holding the choice of distance measure and model parameters constant, we constructed 31 models comprising all possible combinations of additional covariates with the core set of essential covariates. Finally, we explore the robustness of results for 31 randomly selected combinations of distance measure, model parameters and additional covariates. All 116 specifications are a priori valid, assuming that the covariates capture or proxy for all important factors influencing outcomes, but may fail to satisfy the parallel trends condition or produce matches for insufficient numbers of treated observations ( More

Case-study selectionThe case studies are chosen to include different cities around the globe that are commonly used for urban studies50. From the 18 selected case-study cities, 12 cities are part of the C40 cities initiative, which aims to lead the way in urban sustainability. The selection ranges from large cities strongly following a grid street plan (for example, Mexico City and Tokyo) to smaller cities that are not dominated by a grid-like urban layout (for example, Zürich). This selection is not exhaustive, and comparatively more European cities are analysed. Barcelona is considered to test the methodology for the city where the superblock concept originates. The presented analysis can, however, be applied to further cities and is merely constrained by data availability.Street network data processingThe case-study extent for each city is 25 km2, for which the street network is downloaded with the help of Overpass API51 and represented as graphs with edges and nodes. For the graph-based algorithms and geometric processing steps, the python packages networkX52 and shapely53 are used. The downloaded raw data are processed in several steps to obtain more accurate street length estimations. First, network nodes that are within a 15 m distance of each other are clustered to obtain a simplified network. This street network abstraction particularly reduces the complexity of street intersections. Second, closed detached rings (loops) and isolated subgraphs are removed from the street network as well as long ( >300 m) tunnels and streets intersecting buildings. Further network cleaning is conceivable, depending on the case-study context, such as removing elevated streets. Third, very small network edges not forming part of longer streets are removed as they typically represent driveways (Supplementary Information section 2). The street network was not re-designed by extension. For example, the street network connectivity was not increased by adding new streets, which could potentially help to design more superblocks or miniblocks.Before searching for potential superblocks or miniblocks on the street network, the street network nodes are classified into higher- and lower-level nodes. Whenever a node forms part of an intersection (degree ≥3), the node is classified as a higher-level node. A lower-level node forms part of a street (edge) between two higher-level nodes. Typically, the street between two higher-level nodes consists of multiple lower-level nodes and edges as the streets are typically not perfectly straight. The lower-level nodes are considered for calculating distances on the street network and creating the blocks. However, when checking for the network topology criteria to identify superblocks and miniblocks, only the higher-level street network is considered. This differentiation is necessary as otherwise, superblocks and miniblocks often would not be identified on the street network due to the lower-level nodes (Supplementary Information section 1).Street hierarchyThe complexities surrounding superblock implementation are higher for streets that are essential to urban mobility. Contributors to OpenStreetMap can classify streets and assign different attribute labels to define a street network hierarchy. On the basis of the provided street hierarchy, the assumption is made that streets labelled as primary, secondary and trunk streets are not suitable for superblock design. Similarly, streets that form part of a trolleybus or tram route are excluded as potential superblock locations. For this analysis, all footways and private streets are ignored. Streets categorized as pedestrian or living streets are considered to be already not centred on car-based mobility and are ignored as the focus here is on transforming streets that are currently focused around car-based mobility.Calculation of population density and building coverageTo go beyond relying on street network characteristics for detecting superblock design, density values are calculated for the entire street network. For each network node, average density values are calculated considering a radius of 100 m and then averaged per edge by averaging the density values from the starting and end nodes. Alternatively, the population density could be calculated by buffering the street network edges. The population density calculation is based on population data provided by the Center for International Earth Science Information Network and Meta54, which are population estimates based on satellite data and census data at a ~30 m resolution. With the help of OpenStreetMap building footprints, the building footprint coverage is calculated first for every node on the street network considering the same radius of 100 m. Second, the average building coverage per edge is calculated by averaging the value from the start and end node of each edge.Detecting superblocks and miniblocks on the street networkThe street network forms the basis for locating potential superblock and miniblock candidates. Before searching potential candidates, the street network is characterized by street type, population density and building coverage as outlined in the previous two paragraphs for narrowing down the search. Then, the street network is first cleaned for cul-de-sac street elements, and the degree of all street network nodes is calculated. Second, all nodes with degree ≥3 are filtered as they could potentially be part of interior street intersections of superblocks or miniblocks. From this reduced filtered street network, all network cycles are identified as they could potentially be an interior street loop of a superblock. Isolated nodes with degree ≥3 not forming part of a network cycle or interior street loop of a superblock are further evaluated as a potential miniblock. Third, to find the exterior streets, all neighbouring nodes of the considered node(s) are identified, and a shortest-path algorithm is applied on the network where the considered nodes are removed to connect all the neighbouring nodes. If no path is found, the street network typology prevents the design of a miniblock or superblock. For miniblocks, all neighbours of the single node (single interior street intersection) are selected. Similarly, for superblocks, all neighbouring nodes not forming part of the interior street loop are connected along the shortest-path route to form the superblock polygon. If no path is found, or the path crosses a bridge, the node is not further considered. Finally, the geometric properties of the identified potential superblocks or miniblocks are checked to determine whether all the boundary conditions are fulfilled. In case a potential superblock does not fulfil the boundary conditions, the conditions for fulfilling a miniblock are tested. Table 1 lists all conditions, and the geometric scenario calculation is outlined in the next section.Geometric scenario calculationThe geometry of the superblock is calculated on the basis of the assumed 400 m length of the Barcelona superblock (lBCN) consisting of three individual blocks (block side length (l_{{mathrm{block}}} = frac{{l_{{mathrm{BCN}}}}}{3})). When identifying superblocks (S) or miniblocks (M), the following boundary conditions for the exterior street length (lext), interior street length (lint) and length of the interior street loop or ring (r) (for superblocks only) need to be fulfilled:$$l_{{mathrm{ext}},{mathrm{max}}}^{mathrm{M}} = left( {8 times l_{{mathrm{block}}} times z times left( {1 + f} right)} right)$$
(1)
$$l_{{mathrm{ext}},{mathrm{min}}}^{mathrm{M}} = left( {frac{{8 times l_{{mathrm{block}}} times z}}{{1 + f}}} right)$$
(2)
$$l_{{mathrm{int}},{mathrm{max}}}^{mathrm{M}} = left( {4 times l_{{mathrm{block}}} times z times left( {1 + f} right)} right)$$
(3)
$$l_{{mathrm{int}},{mathrm{min}}}^{mathrm{M}} = left( {frac{{3 times l_{{mathrm{block}}} times z}}{{1 + f}}} right)$$
(4)
$$l_{{mathrm{ext}},{mathrm{max}}}^{mathrm{S}} = left( {12 times l_{{mathrm{block}}} times {{{{z}}}} times left( {1 + f} right)} right)$$
(5)
$$l_{{mathrm{ext}},{mathrm{min}}}^{mathrm{S}} = left( {frac{{12 times l_{{mathrm{block}}} times {{{{z}}}}}}{{1 + f}}} right)$$
(6)
$$r_{{mathrm{int}},{mathrm{max}}}^{mathrm{S}} = left( {4 times l_{{mathrm{block}}} times {{{{z}}}} times left( {1 + f} right)} right)$$
(7)
$$r_{{mathrm{int}},{mathrm{min}}}^{mathrm{S}} = left( {frac{{4 times l_{{mathrm{block}}} times {{{{z}}}}}}{{1 + f}}} right)$$
(8)
where f denotes the deviation factor and z incorporates an uncertainty of ±20% of the Barcelona superblock and takes a value of 0.8 (min) or 1.2 (max). Example values for G0 and G1 are provided in Table 1.Street NDIThe Edmonds–Karp algorithm is applied to assess the importance of each network edge concerning traffic flow across the street network. This approach is inspired by ref. 55, where artificial ‘super’ sinks/sources are introduced. Super sinks/sources enable extending a network by adding a single node that feeds all the sources or drains all the sinks56. The following steps are performed to assess the street network edge criticality. In a first step, four supernodes are projected in each direction to the network and placed outside the street network extent. Then, 100 helper nodes are equally distributed along a straight axis on each side of the network extent. An auxiliary linking edge is established between each helper node and the supernode. An additional linking edge is next added between each helper node and the closest node on the street network. A visualization and a more detailed explanation of this step are provided in the Supplementary Information section 3. In a second step, the Edmonds–Karp algorithm is run consecutively for each direction between opposing super sinks and super sources. The resulting flows of each simulation run are summed and averaged for each edge (i) to obtain the average flow ((f_i^{{mathrm{avg}}})) per edge. The average edge flow is then normalized ((f_i^{{mathrm{norm}}})) with the calculated maximum average flow value of the network:$$f_i^{{mathrm{norm}}} = frac{{f_i^{{mathrm{avg}}}}}{{mathop {{max }}limits_j (f_j^{{mathrm{avg}}})}}$$
(9)
To consider local as well as regional network impacts, these outlined steps are performed on a raster with a resolution of 2.5 km as well as for the entire case-study area. The calculations on the raster provide local flows ((f_i^{{mathrm{norm}}})) as outlined in the preceding. As the overall analysed city extent is 25 km2, local flows are calculated across four regional cells across the city. In addition to these local flows, the same calculation is performed once for the entire case-study area using the 5 × 5 km2 as the input for the flow calculations in equation (9). This calculation using the entire street network reflects flows at a higher geographical level and is termed (F_i^{{mathrm{norm}}}). In a third step, the two calculations are equally weighted and combined into a single indicator (NDI) to calculate the relative importance of each edge concerning traffic flow:$${mathrm{NDI}}_i = frac{{f_i^{{mathrm{norm}}} + F_i^{{mathrm{norm}}}}}{2}$$
(10)
Edges with high NDI are edges that are critical to the street network; edges with low NDI have a low disruption potential as they do not form part of a critical network element and alternative paths exist for rerouting traffic. Calculating the NDI provides an approximate indication of the street disruption of a network towards urban mobility. The NDI is further used to derive classes indicating the street network disruption (low, middle, high). As identical geographical extents have been selected across our case studies, the obtained NDI values are comparable. More

Genetic diversity and population structureThe 614 bp length of mtCO1 sequences was successfully amplified and sequenced from 75 individuals of E. tetradactylum and 89 individuals of E. rhadinum from different sites. Based on the CO1 analysis, we detected 5 and 16 haplotypes, respectively, from E. tetradactylum and E. rhadinum (Table 1). Only one haplotype was inter-specifically shared in E. tetradactylum populations, as showed in the TCS haplotype networks (Fig. 2a). A total of 77 polymorphic sites was identified in E. rhadinum but 3 polymorphic sites in E. tetradactylum. Among these sites, a total of 3 and 11 parsimoniously informative sites was detected in E. tetradactylum and E. rhadinum, respectively. In E. tetradactylum, the number of CO1 haplotypes was 2 in ZS and 3 in PA and ZJ. The haplotype diversity was also much higher in ZJ (0.211) and PA (0.197) than ZS (0.105). In E. rhadinum, CO1 haplotypes varied from 3 (JH) to 8 (ZZ). The haplotype diversity was the highest in ZZ (0.663). The populations of ZJ and ZZ showed the statistically negative Tajima’s D value, which could signify the demographic expansion. The MDA revealed similar results (Fig. S3).Table 1 Genetic polymorphisms and neutrality tests of Eleutheronema tetradactylum and Eleutheronema rhadinum inferred from CO1 and 16s rRNA.Full size tableFigure 2The unrooted TCS haplotype networks were constructed based on the haplotypes of CO1 (a) and 16s rRNA (b) of Eleutheronema tetradactylum (left) and Eleutheronema rhadinum (right). Haplotype frequency was related to the size of the circle. Different colors within the nodes refer to various sampling sites.Full size imageThe mitochondrial 16s rRNA (574 bp in length) was also successfully sequenced from 75 and 89 individuals of E. tetradactylum and E. rhadinum (Table 1), which yielded 5 and 6 haplotypes, respectively (Fig. 2b). No haplotype was interspecifically shared of 16s rRNA both in E. tetradactylum and E. rhadinum. A total of 4 and 14 polymorphic sites of E. tetradactylum and E. rhadinum were identified, respectively, of which 3 and 4 were parsimoniously informative sites. Table 1 shows that only four haplotypes with 0.200 haplotype diversity were identified in E. tetradactylum from PA. In E. rhadinum, relatively high haplotype diversity (H = 0.481) and nucleotide diversity (π = 0.00170) were found in populations SA. Overall, the populations from Thailand showed higher genetic diversity than the China population both for E. tetradactylum and E. rhadinum.The TCS network was constructed to identify the phylogenetic relationships in E. tetradactylum and E. rhadinum between China and Thailand populations, as shown in Fig. 2. In E. tetradactylum, 5 haplotypes were closely related to a small number of mutation steps, and the Hap_1 was likely the most primitive haplotype, which evolved into others. In E. rhadinum, 16 haplotypes were distributed between the two branches, including China and Thailand branches. Only the Hap_7 was shared in ZJ and SA of the Thailand branch. One (hap_1) in E. tetradactylum and two (Hap_2 and Hap_8) in E. rhadinum were used as the central radiation distribution for most haplotypes. Other haplotypes were formed by a small number of mutations of these haplotypes. As shown in the TCS network of 16s rRNA haplotypes, the Hap_1 in E. tetradactylum and Hap_4 in E. rhadinum were the most primitive haplotype, which showed central radiation distributions. Also, in E. rhadinum, the haplotypes of China and Thailand populations were divided into two branches; only Hap_2 was shared in ZJ and SA.The level of population genetic differentiation between China and Thailand populations was also evaluated (Table S3). In E. tetradactylum, the average Fixation index (Fst) between PA and the other two sites was 0.81344 in ZS (p  More

Scat sample collection and preparationAt known harbour seal haulout sites individual scat samples were collected using a standardized protocol (Fig. 1). Disposable wooden tongue depressors were used to transfer deposited scats into 500 ml single-use jars or zip-style bags lined with 126 µm nylon mesh paint strainers18. Samples were either preserved immediately in the field by adding 300 ml 95% ethanol to the collection jar, or were taken to the lab and frozen at −20 °C within 6 hours of collection19. Later, samples were thawed and filled with ethanol before being manually homogenized with a disposable wooden depressor inside the paint strainer to separate the scat matrix material from hard prey remains (e.g. bones, cephalopod beaks). The paint strainer containing prey hard parts was then removed from the jar leaving behind the ethanol preserved scat matrix for genetic analysis20. The paint strainer containing prey hard parts was refrozen for subsequent parallel morphological prey ID.Fig. 1The 52 harbour seal scat collection sites in the Salish Sea represented in this dataset.Full size imageMolecular laboratory processingScat matrix samples were subsampled (approximately 20 mg), centrifuged and dried to remove ethanol prior to DNA extraction. DNA was extracted from scat with the QIAGEN QIAamp DNA Stool Mini Kit according to the manufacturer’s protocols. For additional details on the extraction process see Deagle et al.21 and Thomas et al.20.The metabarcoding marker we used to quantify fish and cephalopod proportions was a 16S mDNA fragment (~260 bp) previously described in Deagle et al.15 for pinniped scat analysis. We used the combined Chord/Ceph primer sets: Chord_16S_F (GATCGAGAAGACCCTRTGGAGCT), Chord_16S_R (GGATTGCGCTGTTATCCCT), and Ceph_16S_F (GACGAGAAGACCCTAWTGAGCT), Ceph_16S_R (AAATTACGCTGTTATCCCT). This multiplex PCR reaction is designed to amplify both chordate and cephalopod prey species DNA. A blocking oligonucleotide was included in the all 16S PCRs to limit amplification of seal DNA22. The oligonucleotide (32 bp: ATGGAGCTTTAATTAACTAACTCAACAGAGCA-C3) matches harbour seal sequence (GenBank Accession AM181032) and was modified with a C3 spacer so it is non-extendable during PCR22.A secondary metabarcoding marker was used in a separate PCR reaction to quantity the salmon portion of seal diet, because the primary 16S marker was unable to reliably differentiate between coho and steelhead DNA sequences. This marker was a COI “minibarcode” specifically for salmonids within the standard COI barcoding region: Sal_COI_F (CTCTATTTAGTATTTGGTGCCTGAG), Sal_COI_R (GAGTCAGAAGCTTATGTTRTTTATTCG). The COI amplicons were sequenced alongside 16S such that the overall salmonid fraction of the diet was quantified by 16S, and the salmon species proportions within that fraction were quantified by COI.To take full advantage of sequencing throughput, we used a two-stage labeling scheme to identify individual samples that involved both PCR primer tags and labeled MiSeq adapter sequences. The open source software package EDITTAG was used to create 96 primer sets each with a unique 10 bp primer tag and an edit distance of 5; meaning that to mistake one sample’s sequences for another, 5 insertions, substitutions or deletions would have to occur23.All PCR amplifications were performed in 20 μl volumes using the Multiplex PCR Kit (QIAGEN). Reactions contained 10 μl (0.5 X) master mix, 0.25 μM of each primer, 2.5 μM blocking oligonucleotide and 2 μl template DNA. Thermal cycling conditions were: 95 °C for 15 min followed by 34 cycles of: 94 °C for 30 s, 57 °C for 90 s, and 72 °C for 60 s.Amplicons from 96 individually labeled samples were pooled by running all samples on 1.5% agarose gels, and the luminosity of each sample’s PCR product was quantified using Image Studio Lite (Version 3.1). To combine all samples in roughly equal proportion (normalization), we calculated the fraction of each sample’s PCR product added to the pool based on the luminosity value relative to the brightest band. After 2013, amplicon normalization was performed using SequalPrep™ Normalization Plate Kits, 96-well.Sequencing libraries were prepared from pools of 96 samples using an Illumina TruSeq DNA sample prep kit which ligated uniquely labeled adapter sequences to each pool. Libraries were then pooled and DNA sequencing was performed on Illumina MiSeq using the MiSeq Reagent Kit v2 (300 cycle) for SE 300 bp reads. Samples were sequenced on multiple different runs as part of the larger study; however, typically between 4 and 6 libraries (each a pool of 96 individually identifiable samples) were sequenced on a single MiSeq run.BioinformaticsTo assign DNA sequences to a fish or cephalopod species, we created a custom BLAST reference database of 16S sequences by an iterative process. First, using a list of the fish species of Puget Sound, we searched Genbank for the 16S sequence fragment of all fishes known to occur in the region (71 fish families 230 species)24,25. Reference sequences for each prey species were included in the database if the entire fragment was available, and preference was given to sequences of voucher specimens. When the database was first generated (November, 2012) Genbank contained 16S sequences for 192 of the 230 fish species in the region, and the remaining 38 species were mostly uncommon species unlikely to occur in seal diets. Following a similar procedure, we added to this database sequences for all of the regional cephalopods for which 16S data were available (7 squid species, 2 octopus species). A separate reference database was generated for the COI salmon marker containing Genbank sequences for the nine salmonid species known to occur regionally: Oncorhynchus gorbuscha (Pink Salmon), Oncorhynchus keta (Chum Salmon), Oncorhynchus kisutch (Coho Salmon), Oncorhynchus mykiss (Steelhead), Oncorhynchus nerka (Sockeye Salmon), Oncorhynchus tshawytscha (Chinook Salmon), Oncorhynchus clarkii (Cutthroat Trout), Salmo salar (Atlantic Salmon), Salvelinus malma (Dolly Varden)24.To determine if some species in the database cannot be distinguished from each other at 16S (i.e. have identical sequences in the reference database) a distance matrix was performed on the complete database using the DistanceMatrix function in the R package DECIPHER26. Species with identical sequences were identified as having a distance of “0.00”. In some cases, one haplotype for a species was identical to another species but other haplotypes were not. When two species’ sequences were identical, we ultimately reported both species in the prey_ID field.Sequences were automatically sorted (MiSeq post processing) by amplicon pool using the indexed TruSeqTM adapter sequences. FASTQ sequence files for each library were imported into MacQIIME (version 1.9.1-20150604) for demultiplexing and sequence assignment to species27. For a sequence to be assigned to a sample, it had to match the full forward and reverse primer sequences and match the 10 bp primer tag for that sample (allowing for up to 2 mismatches in either primers or tag sequence).Next, we clustered the DNA sequences that were assigned to scat or tissue samples with USEARCH (similarity threshold = 0.99; minimum cluster size = 3; de novo chimera detection), and entered a representative sequence from each cluster into a GenBank nucleotide BLAST search28,29. If the top matching species for any cluster was not included in the existing database (or the sequence differed indicating haplotype variation), we put the top matching entry in the reference database. We repeated this procedure with every new batch of sequence data to minimize the potential for incorrect species assignment or prey species exclusion. This process was conducted for both the 16S and COI reference databases with each new batch of samples.For all DNA sequences successfully assigned to a sample, a BLAST search was performed against our custom 16S or COI reference databases. A sequence was assigned to a species based on the best match in the database (threshold BLASTN e-value  More