Philippa Kaur

Effects of seasonality on physiochemical properties of waterThe characteristics of the physical and chemical factors of the river are shown in Table 1. Variance analysis showed that the environmental conditions during the four seasons were significantly different (P  More

Coral collection and laboratory acclimationCorals were obtained from the nursery at Mote Marine Laboratory’s Elizabeth Moore International Center for Coral Reef Research & Restoration. Acropora cervicornis fragments (apical tips ~ 5–6 cm in length) were obtained from Mote’s in situ nursery and mounted to a plastic base with marine epoxy (Instant Ocean Holdfast). The A. cervicornis genotypes used in this study were Mote Marine Laboratory genotype numbers 7, 31, 50, 57, and 70 (n = 12 ramets from each of the 5 genets). The coral O. faveolata was obtained from Mote’s ex situ nursery as fragments (~ 3–4 cm2) attached to individual carbonate plugs. Six genotypes of O. faveolata were used in deoxygenation treatments from genotype numbers F3A, F8, F27, F61, F125, and F132 (n = 12 ramets from each of the 6 genets). The immediate environmental conditions at the ex situ coral nursery are not known in detail, however the corals were maintained in an open, flow-through system with high flow rates and continuous bubbling with air stones.Corals were transported to the wet lab facilities at the Smithsonian Marine Station (SMS) in Fort Pierce, FL and experiments were conducted between October and December 2019. Corals were acclimated to laboratory conditions for 4–6 weeks before exposure to deoxygenation treatments. During this period, corals were kept in an indoor, closed seawater system. Holding tanks contained ~ 570 L of seawater that was recirculated through a sump containing rigorous aeration (ambient air) and a heater/chiller that maintained the temperature at ~ 27 ºC (see Table S3 for all parameters). Water was pumped from the sump through a UV sterilizer (Coralife) and then back into the holding tank. Water flow and circulation was maintained by two aquarium pumps (AquaTop MaxFlow MCP-5), and yielded a full water exchange through the sump ~ 8.5 times per hour. Light was provided by LED aquarium lights (HQD) with a maximum photosynthetically active radiation (PAR) of ~ 300 µmol photon m−2 s−1. Full water changes were conducted weekly, and non-living surfaces on the coral fragments (e.g., plastic or carbonate bases) were cleaned at least once a week to reduce algal growth. Conditions within the acclimation tanks closely matched the ambient (i.e., normoxic) conditions at the collection locations and in the oxygen treatments (Table S3, Fig. S4).Field dissolved oxygen conditionsOxygen concentrations were monitored at Mote Marine Laboratory’s in situ nursery (24.56 latitude, − 81.40 longitude) at the depth of coral outplants (~ 5–6 m) to parameterize oxygen conditions used in laboratory experiments. A dissolved oxygen sensor (miniDot, PMEL) was calibrated according to manufacturer protocols and attached to the benthos adjacent to coral outplants. Dissolved oxygen and temperature were logged every ten minutes for the duration of deployment. Hourly and daily average DO concentrations are presented for one month encompassing the approximate time of coral retrievals (Fig. S4). The 6.25 mg L−1 treatment in the laboratory experiments simulates the average normoxic conditions at the site of collection (6.18 ± 0.10, SD), while the most severe deoxygenation treatment (1.00 mg L−1) simulated extreme oxygen depletion that has occurred during acute deoxygenation events on coral reefs21,22,52 (Fig. S4). The two intermediate deoxygenation treatments represented incremental increases by ~ 2 mg L−1 to capture coral responses over a full range of oxygen conditions (Fig. S4).Experimental designDeoxygenation experiments were conducted in the wet lab facilities at SMS using a mesocosm array consisting of 12 tanks (50 L, AquaLogic Systems), with each tank functioning as a closed system with fully independent temperature and oxygen control. Targeted oxygen levels for treatments were 1.00, 2.25, 4.25, and 6.25 mg L−1 DO, with three independent tank replicates for each oxygen level. Oxygen treatments are referred to by the targeted DO concentrations. Species were run in separate, sequential experiments, and genotypes within a species were co-mingled in treatment tanks with one replicate of each genotype per tank.Oxygen concentrations were maintained in each independent tank by bubbling seawater with nitrogen gas and ambient air, with gas injection controlled through a DO feedback and solenoid valves (Neptune Systems). Oxygen levels were monitored every 15 s in each tank by an OxyGuard DO probe connected to an aquarium controller (Neptune Systems, Apex Aquacontroller), which opened or closed the respective solenoid valves to add nitrogen or air in order to maintain the programmed treatment conditions. Oxygen probes were calibrated at the start of experiments, and then again after five days, following the manufacturer’s protocol.Temperature control was maintained in treatment aquaria through an independent heating/chilling loop in each tank and monitored by AquaLogic temperature probes (sensu53). The average temperature at the collection site in the month prior to collection was 29.01 ± 0.62 ℃ (SD) (Fig. S4), which represents seasonally warm temperatures of the Florida Keys54,55. To reduce the potential confounding effects of temperature on responses to deoxygenation, we acclimated corals to 27 ℃ during the holding period, and conducted experiments at ~ 27 ℃. This represents the average annual temperature at the collection site in the Florida Keys, is the temperature that O. faveolata fragments were maintained at in the ex situ nursery prior to collection (27.0 ± 0.62 ℃), and is commonly used as the ambient temperature for corals from this location54,56. Each tank contained an additional probe (Neptune Systems) that logged temperature every 10 s for the duration of the experiment. In addition to continuous monitoring of DO and temperature, discrete measurements were taken each day for DO, temperature, pH, and salinity with a handheld multi-parameter water quality meter (YSI ProDSS with optical DO sensor) (Tables 1, 2, Tables S1, S2). The YSI DO probe was calibrated at the start of every day following manufacturer protocols and used to validate DO measurements from OxyGuard probes and to adjust programmed values to maintain treatment targets when necessary.Individual tanks were supplied with a 7-color LED aquarium light (Aquaillumination, Hydra 52), programmed to simulate a diel light cycle over a 12:12 h photoperiod, and set to maximum irradiance of ~ 300 µmol photon m−2 s−1 (PAR). This maximum intensity is likely lower than what corals may experience in situ during peak midday irradiances, however, it matches the irradiance levels these corals were maintained at in the ex situ nursery (320 ± 182 SD, n = 45) and is sufficient to stimulate maximal photosynthesis without causing light stress45. Light levels were measured with a light meter (Licor, LI-1400) and an underwater spherical quantum sensor (LI-193SA) submerged at the center of each tank at the start of each experiment and again after five days. Oxygen, temperature, and light conditions were effectively maintained at or near targeted levels for the duration of each experiment (Fig. 1, Table 1), with minimal differences between replicate tanks within a treatment (Tables S1, S2). We did not simultaneously manipulate pH in treatment tanks, which resulted in differences in pH between treatments due to bubbling with nitrogen gas (Tables 1, 2).Coral conditionCorals were monitored daily during the experiments for changes in live tissue cover and photophysiology, and then destructively sampled at the end of the experiment for symbiont densities. Fragments were visually evaluated at each time point for evidence of bleaching, tissue loss, and mortality. Tissue loss was quantified as the percent of each fragment with exposed skeleton, resulting from the sloughing of tissue. Fragments were categorized as dead when all living tissue was lost.The A. cervicornis experiment ended after five days, at which point 50% of corals in the lowest DO treatment were dead or displayed partial mortality. We ended the O. faveolata experiment after 11 days, more than twice the duration of the A. cervicornis experiment (Fig. 1), at which point no observable signs of deterioration were detected in corals from any treatment (Fig. 2a).PAM fluorometryPAM fluorometry is a non-destructive method that directly measures chlorophyll fluorescence and the activity of photosystem II (PSII)57, and the quantum yield can be used as a proxy for coral stress36,58. For all PAM readings, one measurement was taken from the same position on each fragment with the probe held at a 90º angle ~ 0.5 cm from the coral surface (angle and distance were maintained with a Walz probe holder). To optimize initial fluorescence (F0) to between 300–500, the following PAM settings were used: gain = 2, damp = 2, saturation intensity = 8, saturation width = 0.8, and measuring light intensity = 10 for A. cervicornis and 6 for O. faveolata. PAM measurements were taken pre-dawn daily, at 0500–0600 h. The position on A. cervicornis fragments was shifted if tissue sloughing occurred so that measurements were taken on living tissue. After the final PAM measurements, corals were snap frozen in liquid nitrogen and stored at − 80 ºC for subsequent analyses.Symbiont densitySymbiont density analyses were conducted at the University of Florida in Gainesville, FL following standard protocols (sensu59). In brief, tissue was stripped from the coral skeleton with an air brush (Master Airbrush S68) and filtered seawater. The tissue slurry was then homogenized with a handheld electric tissue homogenizer (Tissue-Tearor) and subsampled for symbiont counts. Symbiont cells were counted with a hemocytometer, with six replicate counts per fragment. Replicate counts were averaged per fragment for all analyses.Surface area of A. cervicornis fragments was determined by wax dipping60 and the surface area of O. faveolata fragments was determined through image analysis in ImageJ. For A. cervicornis, surface area was then corrected by the percent tissue loss recorded for each fragment. Symbiont densities were normalized to fragment live tissue surface area and expressed as cells per cm−2.Statistical analysesAll analyses were conducted in R (v 4.0.2)61. The effects of fixed and random factors were evaluated with linear mixed effects models using the package lme462. Normality and homoscedasticity of variances of response variables were evaluated by visual inspection of residuals and Levene’s tests, respectively. All variables met model assumptions, and results are reported for full models that included all fixed and random effects.For maximum quantum yield and tissue loss, the factor corresponding to measurements repeated daily from the start to end of the experiment (i.e., Day 1, 2, 3, etc.) is referred to as “day”. Day, genotype, and treatment were treated as fixed factors in the analysis of maximum quantum yield and tissue loss, with individual included as a random factor to account for repeated measures, and tank as random factor to account for fragments of different genotypes within a tank (Tables 2, 3). For symbiont density, genotype and treatment were analyzed as fixed factors, and tank as a random factor (Table 4). The significance of fixed effects was evaluated with type II ANOVA tables using Satterthwaite’s method, and Tukey’s post-hoc tests were used where necessary to determine significant differences between levels of a factor using the package emmeans63. Data in main figures in which genotype was not a significant effect (i.e., all response variables) are presented as treatment averages (± SE), pooled across genotypes. More

1.de Souza Dias, V., Pereira da Luz, M., Medero, G. M. & Tarley Ferreira Nascimento, D. An overview of hydropower reservoirs in Brazil: Current situation, future perspectives and impacts of climate change. Water 10, 592 (2018).
Google Scholar 
2.Patias, J., Zuquette, L. V. & Rodrigues-Carvalho, J. A. Piezometric variations in the basaltic massif beneath the Itaipu hydroelectric plant (Brazil/Paraguay border): Right Buttress Dam. Bull. Eng. Geol. Environ. 74, 207–231 (2015).CAS 

Google Scholar 
3.Agostinho, A. A. Pesquisas, monitoramento e manejo da fauna aquática em empreendimentos hidrelétricos. In Seminário Sobre Fauna Aquática E O Setor Elétrico Brasileiro 38–59 (Brasil, 1994).4.Makrakis, S., Gomes, L. C., Makrakis, M. C., Fernandez, D. R. & Pavanelli, C. S. The Canal da Piracema at Itaipu Dam as a fish pass system. Neotrop. Ichthyol. 5, 185–195 (2007).
Google Scholar 
5.Dos Reis, R. B., Frota, A., Depra, G. D. C., Ota, R. R. & Da Graca, W. J. Freshwater fishes from Paraná State, Brazil: An annotated list, with comments on biogeographic patterns, threats, and future perspectives. Zootaxa 4868, 451–494 (2020).
Google Scholar 
6.Becker, R. A., Sales, N. G., Santos, G. M., Santos, G. B. & Carvalho, D. C. DNA barcoding and morphological identification of neotropical ichthyoplankton from the Upper Paraná and São Francisco. J. Fish Biol. 87, 159–168 (2015).CAS 
PubMed 

Google Scholar 
7.Milan, D. T. et al. New 12S metabarcoding primers for enhanced Neotropical freshwater fish biodiversity assessment. Sci. Rep. 10, 1–12 (2020).ADS 

Google Scholar 
8.Agostinho, A. A., Pelicice, F. M. & Gomes, L. C. Dams and the fish fauna of the Neotropical region: Impacts and management related to diversity and fisheries. Braz. J. Biol. 68, 1119–1132 (2008).CAS 
PubMed 

Google Scholar 
9.Bonar, S. A., Hubert, W. A. & Willis, D. W. Standard methods for sampling North American freshwater fishes. American Fisheries Society, Bethesda, (USA, 2009).
Google Scholar 
10.Shaw, J. L. A. et al. Comparison of environmental DNA metabarcoding and conventional fish survey methods in a river system. Biol. Conserv. 197, 131–138 (2016).
Google Scholar 
11.Reis, R. E. et al. Fish biodiversity and conservation in South America. J. Fish Biol. 89, 12–47 (2016).CAS 
PubMed 

Google Scholar 
12.Baumgartner, G. et al. Peixes do baixo rio Iguaçu. (Eduem, 2012).13.Taberlet, P., Bonin, A., Coissac, E. & Zinger, L. Environmental DNA: For Biodiversity Research and Monitoring (Oxford University Press, 2018).
Google Scholar 
14.Taberlet, P., Coissac, E., Pompanon, F., Christian, B. & Willerslev, E. Towards next-generation biodiversity assessment using DNA metabarcoding. Mol Ecol 33, 2045–2050 (2012).
Google Scholar 
15.Ritter, C. D. et al. The pitfalls of biodiversity proxies: Differences in richness patterns of birds, trees and understudied diversity across Amazonia. Sci. Rep. 9, 1–3 (2019).CAS 

Google Scholar 
16.Sales, N. G. et al. Space-time dynamics in monitoring neotropical fish communities using eDNA metabarcoding. Sci. Total Environ. 754, 142096 (2021).ADS 
CAS 
PubMed 

Google Scholar 
17.Zinger, L. et al. Body size determines soil community assembly in a tropical forest. Mol. Ecol. 28, 528–543 (2019).CAS 
PubMed 

Google Scholar 
18.Baird, D. J. & Hajibabaei, M. Biomonitoring 2.0: A new paradigm in ecosystem assessment made possible by next-generation DNA sequencing. Mol. Ecol. 21, 2039–2044 (2012).PubMed 

Google Scholar 
19.Zinger, L. et al. Advances and prospects of environmental DNA in neotropical rainforests. Adv. Ecol. Res. 62, 331–373 (2020).
Google Scholar 
20.Cilleros, K. et al. Unlocking biodiversity and conservation studies in high-diversity environments using environmental DNA (eDNA): A test with Guianese freshwater fishes. Mol. Ecol. Resour. 19, 27–46 (2019).CAS 
PubMed 

Google Scholar 
21.Sales, N. G., Wangensteen, O. S., Carvalho, D. C. & Mariani, S. Influence of preservation methods, sample medium and sampling time on eDNA recovery in a neotropical river. Environ. DNA 119–130. https://doi.org/10.1002/edn3.14 (2020).Article 

Google Scholar 
22.Blaxter, M. et al. Defining operational taxonomic units using DNA barcode data. Philos. Trans. R. Soc. B Biol. Sci. 360(1462), 1935–1943. https://doi.org/10.1098/rstb.2005.1725 (2005).CAS 
Article 

Google Scholar 
23.Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
24.Edgar, R. C. UNOISE2: Improved error-correction for Illumina 16S and ITS amplicon sequencing. BioRxiv 81257 (2016).
25.Muha, T. P., Rodriguez-Barreto, D., O’Rorke, R., Garcia de Leaniz, C. & Consuegra, S. Using eDNA metabarcoding to monitor changes in fish community composition after barrier removal. Front. Ecol. Evol. 9, 28 (2021).
Google Scholar 
26.Kitano, T., Umetsu, K., Tian, W. & Osawa, M. Two universal primer sets for species identification among vertebrates. Int. J. Legal Med. 121, 423–427 (2007).PubMed 

Google Scholar 
27.Stoeckle, M. Y., Soboleva, L. & Charlop-Powers, Z. Aquatic environmental DNA detects seasonal fish abundance and habitat preference in an urban estuary. PLoS One 12, e0175186 (2017).PubMed 
PubMed Central 

Google Scholar 
28.Bylemans, J. et al. An environmental DNA-based method for monitoring spawning activity: A case study, using the endangered Macquarie perch (Macquaria australasica). Methods Ecol. Evol. 8, 646–655 (2017).
Google Scholar 
29.De Souza, L. S., Godwin, J. C., Renshaw, M. A. & Larson, E. Environmental DNA (eDNA) detection probability is influenced by seasonal activity of organisms. PLoS One 11, e0165273 (2016).PubMed 
PubMed Central 

Google Scholar 
30.Ritter, C. D. et al. Locality or habitat? Exploring predictors of biodiversity in Amazonia. Ecography (Cop.) 42, 321–333 (2019).
Google Scholar 
31.CFMV-Resolução no 1000 de 11 de maio de 2012—Dispõe sobre procedimentos e métodos de eutanásia em animais e dá outras providências. (2012).32.Britski, H. A., de Silimon, K. Z. S. & Lopes, B. S. Peixes do Pantanal: manual de identificação, ampl. Brasília, DF, Embrapa Informação Tecnológica (2007).33.Ota, R. R., Deprá, G. de C., Graça, W. J. da & Pavanelli, C. S. Peixes da planície de inundação do alto rio Paraná e áreas adjacentes: revised, annotated and updated. Neotrop. Ichthyol. 16(2). https://www.scielo.br/j/ni/a/tScwvm8JLhKnbxKjtBQLPBx/abstract/?lang=en (2018).34.Neris, N., Villalba, F., Kamada, D. & Viré, S. Guía de peces del Paraguay/Guide of fishes of Paraguay. Zamphiropolos, (Paraguay, 2010).
Google Scholar 
35.Pie, M. R. et al. Development of a real-time PCR assay for the detection of the golden mussel (Limnoperna fortunei, Mytilidae) in environmental samples. An. Acad. Bras. Cienc. 89, 1041–1045 (2017).CAS 
PubMed 

Google Scholar 
36.Miya, M. et al. MiFish, a set of universal PCR primers for metabarcoding environmental DNA from fishes: Detection of more than 230 subtropical marine species. R. Soc. Open Sci. 2, 150088 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
37.Boeger, W. A. et al. Testing a molecular protocol to monitor the presence of golden mussel larvae (Limnoperna fortunei) in plankton samples. J. Plankton Res. 29, 1015–1019 (2007).CAS 

Google Scholar 
38.Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 17, 10–12 (2011).
Google Scholar 
39.Van Rossum, G. & Drake, F. L. Python 3 References Manual. Scotts Valley CA: CreateSpace. (2009).40.R Core Team. R: the R project for statistical computing. 2019. https://www.r-project.org/ (accessed 30 Mar 2020).41.Edgar, R. C. UPARSE: Highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996–998 (2013).CAS 
PubMed 

Google Scholar 
42.Edgar, R. C. & Flyvbjerg, H. Error filtering, pair assembly and error correction for next-generation sequencing reads. Bioinformatics 31, 3476–3482 (2015).CAS 
PubMed 

Google Scholar 
43.Camacho, C. et al. BLAST+: Architecture and applications. BMC Bioinform. 10, 421 (2009).
Google Scholar 
44.Team, Rs. RStudio: integrated development for R. RStudio, Inc., Boston, MA https://www.rstudio.com42, 84 (2015).45.Wickham, H. tidyverse: Easily Install and Load “Tidyverse” Packages (Version R package version 1.1. 1). (2017).46.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).MATH 

Google Scholar 
47.Tang, Y., Horikoshi, M. & Li, W. ggfortify: Unified interface to visualize statistical results of popular R packages. R J. 8, 474 (2016).
Google Scholar 
48.Auguie, B. & Antonov, A. gridExtra: Miscellaneous functions for “grid” graphics (Version 2.2. 1)[Computer software]. (2016).49.Kassambara, A. & Kassambara, M. A. Package ‘ggpubr’. (2020).50.Oksanen, J. et al. Vegan: Community ecology package. R package version 1.17-4. https://cran.r-project.org. Acesso em 23, 2010 (2010).51.McMurdie, P. J. & Holmes, S. Waste not, want not: Why rarefying microbiome data is inadmissible. PLoS Comput. Biol. 10, e1003531 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
52.Jost, L. Entropy and diversity. Oikos 113, 363–375 (2006).
Google Scholar 
53.Marcon, E., Herault, B. & Marcon, M. E. Package ‘entropart’. (2021).54.Mächler, E., Walser, J.-C. & Altermatt, F. Decision making and best practices for taxonomy-free eDNA metabarcoding in biomonitoring using Hill numbers. BioRxiv (2020).55.McMurdie, P. J. & Holmes, S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8, e61217 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
56.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 
57.León, A., Reyes, J., Burriel, V. & Valverde, F. Data quality problems when integrating genomic information. In International Conference on Conceptual Modeling 173–182 (Springer, 2016).58.Callahan, B. J., McMurdie, P. J. & Holmes, S. P. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 11, 2639–2643 (2017).PubMed 
PubMed Central 

Google Scholar 
59.Stahlhut, J. K. et al. DNA barcoding reveals diversity of hymenoptera and the dominance of parasitoids in a sub-arctic environment. BMC Ecol. 13, 2 (2013).PubMed 
PubMed Central 

Google Scholar 
60.Gillet, B. et al. Direct fishing and eDNA metabarcoding for biomonitoring during a 3-year survey significantly improves number of fish detected around a South East Asian reservoir. PLoS One 13, e0208592 (2018).PubMed 
PubMed Central 

Google Scholar 
61.Barrett, M. et al. Living planet report 2018: Aiming higher. WWF. Available at: https://www.globallandscapesforum.org/publication/living-planet-report-2018-aiming-higher/ (2018).62.Díaz, S. M. et al. The global assessment report on biodiversity and ecosystem services: Summary for policy makers. Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. 56, (2019).63.Dudgeon, D. Asian river fishes in the Anthropocene: Threats and conservation challenges in an era of rapid environmental change. J. Fish Biol. 79, 1487–1524 (2011).CAS 
PubMed 

Google Scholar 
64.Dudgeon, D. Multiple threats imperil freshwater biodiversity in the Anthropocene. Curr. Biol. 29, R960–R967 (2019).CAS 
PubMed 

Google Scholar 
65.He, F. et al. Disappearing giants: A review of threats to freshwater megafauna. Wiley Interdiscip. Rev. Water 4, e1208 (2017).
Google Scholar 
66.Agostinho, A. A., Thomaz, S. M. & Gomes, L. C. Threats for biodiversity in the floodplain of the Upper Paraná River: Effects of hydrological regulation by dams. (2018). Int. J. Ecohydrol. Hydrobiol Warsaw. 4(3), 267–280 (2004).
Google Scholar 
67.Santana, M. L., Carvalho, F. R. & Teresa, F. B. Broad and fine-scale threats on threatened Brazilian freshwater fish: Variability across hydrographic regions and taxonomic groups. Biota Neotrop. 21 (2). https://www.scielo.br/j/bn/a/YqFbWSy5vbfHy3QK9kNpdKp/?format=html&lang=en (2021).68.Matthews, W. J. Patterns in Freshwater Fish Ecology. (Springer Science & Business Media, 2012).69.de Oliveira Bueno, E., Alves, G. J. & Mello, C. R. Hydroelectricity water footprint in Parana hydrograph region, Brazil. Renew. Energy 162, 596–612 (2020).
Google Scholar 
70.Camacho Guerreiro, A. I., Amadio, S. A., Fabre, N. N. & da Silva Batista, V. Exploring the effect of strong hydrological droughts and floods on populational parameters of Semaprochilodus insignis (Actinopterygii: Prochilodontidae) from the Central Amazonia. Environ. Dev. Sustain. 23, 3338–3348 (2021).
Google Scholar 
71.Jespersen, H., Rasmussen, G. & Pedersen, S. Severity of summer drought as predictor for smolt recruitment in migratory brown trout (Salmo trutta). Ecol. Freshw. Fish 30, 115–124 (2021).
Google Scholar 
72.Pool, T. K., Grenouillet, G. & Villéger, S. Species contribute differently to the taxonomic, functional, and phylogenetic alpha and beta diversity of freshwater fish communities. Divers. Distrib. 20, 1235–1244 (2014).
Google Scholar 
73.de Oliveira, E. F., Goulart, E. & Minte-Vera, C. V. Fish diversity along spatial gradients in the Itaipu Reservoir, Paraná, Brazil. Braz. J. Biol. 64, 447–458 (2004).CAS 
PubMed 

Google Scholar 
74.Daga, V. S. et al. Homogenization dynamics of the fish assemblages in Neotropical reservoirs: Comparing the roles of introduced species and their vectors. Hydrobiologia 746, 327–347 (2015).
Google Scholar 
75.Vitule, J. R. S. Introdução de peixes em ecossistemas continentais brasileiros: revisão, comentários e sugestões de ações contra o inimigo quase invisível. Neotrop. Biol. Conserv. 4, 111–122 (2009).
Google Scholar 
76.Mariac, C. et al. Species‐level ichthyoplankton dynamics for 97 fishes in two major river basins of the Amazon using quantitative metabarcoding. Mol. Ecol. https://onlinelibrary.wiley.com/action/showCitFormats?doi=10.1111%2Fmec.15944 (2021).77.Jackman, J. M. et al. eDNA in a bottleneck: Obstacles to fish metabarcoding studies in megadiverse freshwater systems. Environ. DNA 3, 837–849 (2021).
Google Scholar 
78.Bessey, C. et al. Maximizing fish detection with eDNA metabarcoding. Environ. DNA 2, 493–504 (2020).
Google Scholar 
79.Evans, N. T. et al. Fish community assessment with eDNA metabarcoding: Effects of sampling design and bioinformatic filtering. Can. J. Fish. Aquat. Sci. 74, 1362–1374 (2017).CAS 

Google Scholar 
80.Prodan, A. et al. Comparing bioinformatic pipelines for microbial 16S rRNA amplicon sequencing. PLoS One 15, e0227434 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
81.Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 2460–2461 (2010).CAS 

Google Scholar 
82.Elbrecht, V. & Leese, F. Can DNA-based ecosystem assessments quantify species abundance? Testing primer bias and biomass—sequence relationships with an innovative metabarcoding protocol. PLoS ONE 10, e0130324 (2015).PubMed 
PubMed Central 

Google Scholar 
83.Pawluczyk, M. et al. Quantitative evaluation of bias in PCR amplification and next-generation sequencing derived from metabarcoding samples. Anal. Bioanal. Chem. 407, 1841–1848 (2015).CAS 
PubMed 

Google Scholar 
84.Holman, L. E., Chng, Y. & Rius, M. How does eDNA decay affect metabarcoding experiments? Environ. DNA https://onlinelibrary.wiley.com/action/showCitFormats?doi=10.1002%2Fedn3.201 (2021).85.Edgar, R. C. UNCROSS2: identification of cross-talk in 16S rRNA OTU tables. BioRxiv 400762 (2018).86.MacArthur, R. H. Geographical Ecology: Patterns in the Distribution of Species. (Princeton University Press, 1984).87.Leray, M. & Knowlton, N. Random sampling causes the low reproducibility of rare eukaryotic OTUs in Illumina COI metabarcoding. PeerJ 5, e3006 (2017).PubMed 
PubMed Central 

Google Scholar 
88.Team, Q. D. QGIS geographic information system. Open Source Geospatial Found. Proj. Versão 2, (2015). More

Microbial community samplingHippo gut microbiomeWe characterized the microbial communities in the hippo gut by collecting ten samples of fresh hippo feces adjacent to four hippo pools early in the morning (prior to desiccation by the sun) in September 2017. We collected feces from different pools and locations adjacent to the pool to include the feces of different individuals so we could estimate the similarity of the gut microbiome among individuals and across the landscape. The four hippo pools are sufficiently far apart that there was likely no intermixing of hippos among them.Each individual hippo feces sample was gently homogenized by hand and then the liquid was gently squeezed from the coarse particulate organic matter. A portion of the liquid (approximately 10 mL) was vacuum filtered through a Supor polysulfone membrane (0.2-µm pore size; Pall, Port Washington, NY, USA). After approximately 10 mL had filtered through and the filter was dry, 15 mL of RNALater was gently poured onto the filter and allowed to contact the collected biomass on the filter for 15 min before being removed by filtration. The filter was stored dry in a sterile petri dish and transferred to a refrigerator within several hours, then to a − 20 °C freezer for storage within several days.During the July 2016 survey of hippo pools, we collected an additional two samples of fresh hippo feces near a high-subsidy hippo pool and filtered approximately 10 mL of the liquid portion after homogenization as detailed above. The filter was then folded twice to preserve the biomass on the filter and stored in 14 mL of RNALater.Aquatic ecosystemWe characterized the microbial communities in the water column of hippo pools across a gradient of hippo subsidy (July 2016, N = 12 pools). We collected samples from the upstream, downstream, surface, and bottom of both pools containing hippos and pools that lacked hippos. Subsamples were also analyzed for biogeochemical variables (details provided below). We also collected water samples in four of the high-subsidy hippo pools every 2–3 days starting immediately after a flushing event until the next flushing event (August and September 2017, Supplementary Fig. S1)23. The number of hippos, discharge and volume for each pool are presented in Dutton et al (2020)23.We sampled the aquatic microbial community and biogeochemical variables along a longitudinal transect down both the Mara and Talek rivers (Supplementary Fig. S1, Supplementary Table S2). For the Mara River, we sampled an approximately 100-km transect along a gradient of hippo numbers (N = 10 locations, from 0 to ~ 4000 hippos). For the Talek River, we sampled an approximately 30-km transect to the confluence with the Mara (N = 8 locations, from 0 to 700 hippos). Mara River sites 9 and 10 are downstream of the confluence with the Talek River. Water samples were collected from each site in a well-mixed flowing section away from any hippo pools.Aquatic microbial samples were collected by filtering water samples through a Supor polysulfone filter (0.2-µm pore size; Pall, Port Washington, NY, USA) and then preserving the filter in RNALater Stabilization Solution (Ambion, Inc., Austin, TX, USA). In 2017, the filters were preserved with RNA Later and then frozen for analysis.Mesocosm experimentWe collected river water from the Mara River upstream of the distribution of hippos and placed it in 45 1-L bottles in a large water basin covered by a dark tarp to help regulate temperature and prevent algal production. Bottles were randomly assigned to the control, bacteria, and bacteria + virus treatments. We collected fresh hippo feces from multiple locations adjacent to the Mara River. After homogenization, half of the hippo feces was sterilized in a pressure cooker, which testing confirmed had similar sterilization results as an autoclave53 (see Supplementary Materials). Five grams of sterilized hippo feces was placed into each bottle to provide an organic matter substrate without viable bacteria or viruses. The unsterilized hippo feces was expressed, and the resulting liquid was filtered through 0.7-µm GF/F filters (0.7-µm pore size; Whatman, GE Healthcare Life Sciences, Pittsburgh, PA, USA) and 0.2-µm Supor filters to physically separate the bacteria (on the filter papers) from the viruses (in the filtrate). Half the filtrate was then sterilized with a UV light treatment (Supplementary Fig. S4). The UV light treatment did not significantly alter DOC quality (see Supplementary Materials).We prepared 15 bottles for each of three treatments—control, bacteria, and bacteria + virus—as follows: Control Unfiltered river water, 5 g wet weight sterilized hippo feces, and two blank Supor filters; Bacteria Unfiltered river water, 5 g wet weight sterilized hippo feces, two Supor filters containing bacteria, and 4 mL sterilized filtrate; Virus Unfiltered river water, 5 g wet weight sterilized hippo feces, two Supor filters containing bacteria, and 4 mL unsterilized filtrate containing viruses.We conducted the experiment for 27 days from September to October 2017. We terminated the experiment after 27 days because we were trying to replicate the microbial communities in hippo pools as best as we could and the hippo pools rarely go more than 1 month before they are flushed out by a flood25. Initial microbial samples of the river water, hippo feces bacteria and hippo fecal liquid filtrate were taken on day 0, and three replicate samples per treatment were destructively sampled on day 3, 9, 15, 21, and 27. During each time step, the microbial communities were sampled using the methods detailed above, and chemical analyses were done on the water samples as described below. We also measured chlorophyll a, dissolved oxygen, temperature, conductivity, total dissolved solids, turbidity, and pH with a Manta2 water quality sonde (Eureka Water Probes, Austin, TX, USA).Microbial community characterizationWe used 16S rRNA sequencing to characterize the active microbial communities. We extracted both DNA and RNA from our preserved samples, then used RNA to synthesize cDNA to represent the “active” microbial community and the total DNA in the sample to represent the “total” microbes present, including those that may not be actively replicating54. Due to the continual loading of hippo feces into pools and the long half-life of DNA, we would expect there to be significant quantities of microbial DNA derived from hippo feces within the pools. However, there would be less accumulation of RNA because of RNA’s shorter half-life. The active communities identified through this RNA-based approach are the ones that would potentially contribute to ecosystem function55 as indicated by the protein synthesis potential, although relationships between activity and rRNA concentrations in individual taxa within mixed communities can vary56. Nevertheless, this method provides an overall characterization of the microbial community’s potential activity.We used the Qiagen RNeasy Powerwater Kit (Qiagen, Hilden, Germany) to extract the DNA and RNA from the material on the filter using a slightly modified manufacturer’s protocol to allow for the extraction of both DNA and RNA. After extraction, we split the total extracted volume (100 µL per sample) into two groups. We treated one group with the DNase Max Kit (Qiagen, Hilden, Germany) to remove all DNA and serve as the RNA group of samples.We used the RNA group of samples to synthesize cDNA using the SuperScript III First Strand Synthesis Kit (Invitrogen, Carlsbad, CA, USA). DNA and cDNA were quantified using the PicoGreen dsDNA Assay Kit (Molecular Probes, Eugene, OR, USA) then normalized to 5 ng/µL. Amplicon library preparation was done using a dual-index paired-end approach57. We amplified the V4 region of the 16S rRNA gene using dual-index primers (F515/R805) and AccuPrime Pfx SuperMix (Invitrogen, Carlsbad, CA, USA) in duplicate for each sample using the manufacturer’s recommended thermocycling routine.Samples were then pooled, purified and normalized using the SequelPrep Normalization Plate Kit (Invitrogen, Carlsbad, CA, USA). Barcoded amplicon libraries were then sequenced at the Yale Center for Genome Analysis (New Haven, CT, USA) using an Illumina Miseq v2 reagent kit (Illumina, San Diego, CA, USA) to generate 2 × 250 base pair paired-end reads.Sampling took place in 2016 and 2017 and involved two separate sequencing runs. The first sequencing run included negative controls and a mock community (D6306, Zymo Research, Irvin, CA, USA). The second sequencing run included negative controls, a mock community (D6306), and a single E. coli strain. In both runs, the mock community and single E. coli strain were well reconstructed from the sequences, and there was minimal contamination in the negative controls, mock community and E. coli strain.From those two sequencing campaigns, we received over 2 million raw sequences from the first campaign and over 7 million raw sequences for the second campaign. For the microbial community analyses, only samples collected and sequenced during the same campaign are analyzed together to prevent preservation or sequencing biases. However, samples within the two separate campaigns were preserved and sequenced using identical methods with only a minor modification (mentioned above) to increase the preservation of genetic material.We de-multiplexed sequenced reads then removed barcodes, indexes, and primers using QIIME258. We used DADA2 with a standard workflow in R59 to infer exact sequence variants (ESV) for each sample60. We assigned taxonomy using a naïve Bayesian classifier and the SILVA training set v. 128 database61,62. We removed potential contamination in samples from both campaigns by using the statistical technique in the R package, decontam63. We used Phyloseq to characterize, ordinate, and compare microbial communities64 with their standard workflow59.Chemical analysesAll water samples collected in the field and in the experiment were analyzed for dissolved ferrous iron (Fe(II)), hydrogen sulfide (H2S), dissolved organic carbon (DOC), inorganic nutrients, major ions, dissolved gases, and biochemical oxygen demand following the standard methods provided in detail in Dutton et al (2020)23.Statistical analysesWe computed all statistical analyses in the R 4.1.1 statistical language in RStudio 2021.09.0 using α = 0.05 to determine significance65,66. Error bars in the figures represent standard deviation of the means. All data and R code for the statistics and data treatments are provided in the Mendeley Data Online Repository67.We used the Bray–Curtis dissimilarity matrix followed by ordination with NMDS to examine differences between individual hippo gut microbiomes; between low-, medium-, and high-subsidy hippo pools; and between a gradient of hippo pools and the environment. We used a CCA to test for the influence of biogeochemical drivers on microbial community composition using biogeochemical data that were previously published but collected concurrently with this study23. We constrained the CCA ordination by soluble reactive phosphorus, nitrate, methane, BOD, and sulfate, which were all previously shown to be important drivers in the variation between pools23. We used PERMANOVA and PERMDISP to test for significant differences between groups68.
We compared aquatic microbial communities from the bottom of high-subsidy hippo pools (N = 15), from hippo feces (N = 10, the hippo gut microbiome) and upstream of high-subsidy hippo pools (N = 15, free of hippo gut microbiome influence) using the Bray–Curtis dissimilarity matrix on the relative abundances for the active aquatic microbial communities collected from the different sample types followed by ordination with NMDS. 95% confidence ellipses were generated. We then determined the active taxa that were shared between the hippo gut microbiome (hippo feces) and the bottom of the high-subsidy hippo pools and not present in the upstream samples from high-subsidy hippo pools.We used LEfSe to calculate the differential abundance of microbial taxa between upstream (N = 14), downstream (N = 16), at the surface (N = 17) and at the bottom (N = 14) of hippo pools and calculated their effect size69. We then calculated the correlation of microbial taxa to the measured biogeochemistry using Pearson’s correlation coefficient with a false discovery rate corrected p-value in the microeco R package70.
We used SourceTracker to quantify the contribution of the hippo gut, upstream waters, or unknown sources to the active aquatic microbial communities in the bottom waters of three of the high-subsidy hippo pools between flushing flows71. We also used the Bray–Curtis dissimilarity matrix followed by ordination with NMDS to examine changes in the active aquatic microbial communities in one of the high subsidy hippo pools through time after flushing flows.For the experiment, we calculated the Bray–Curtis dissimilatory matrix followed by ordination with NMDS for the active aquatic microbial communities over time in each of the three experimental treatments. We used SourceTracker to determine the proportion of the active aquatic microbial community in each treatment that originated from the hippo gut, the river water, or unknown sources71. We analyzed the biogeochemical differences among experimental treatments by fitting a linear mixed effects model for each of the biogeochemical variables throughout the experiment with the nlme package in R72. We fit the model with the restricted maximum likelihood method and a continuous autoregressive temporal correlation structure with sample day as the repeated factor. Treatment and time were fixed effects and individual bottles were treated as random effects. We conducted a pairwise post-hoc test with an ANOVA and the emmeans package in R to estimate marginal means with a Tukey adjusted p-value for multiple comparisons73,74. More

MPs abundanceIn Table 1 MP abundance (mean value ± standard deviation) values are presented by shape, size range, color and polymer type categories for each sampling site. MP were found in all analyzed salt samples including pellets, fibers, fragments, films and lines (Fig. 3). MP total abundance values per site ranged from 74.7 to 136.7 particles kg−1 in the following order of increasing abundance: S3  black (17%)  > blue (15%)  > green and transparent (10% each)  > pink (6%)  > colorless (5%). In terms of size, most particles were in the category 500–1000 µm, except for S3 (1000–5000 µm) (Table 1). The distribution of MP particles based on size range was: 500–1000 µm (40%)  > 1000–5000 µm (34%)  > 250–500 µm (26%). For salts from the Atlantic and the Pacific Ocean, originating from Brazil, the United Kingdom, and the USA, Kim et al.12 reported a higher abundance of MP in size range 100–1000 µm while sizes in the range 100–5000 µm were reported for salt samples from the Black Sea. Seth and Shriwastay20 found that 80% of fibers found in salt samples from the Indian Sea were smaller than 2000 μm in length. MP size range differences among the various studies are suggested to depend on the degree of weathering for a given environment30, different climatic conditions such as wind, rain, temperature, salinity, and waves influencing size range composition. Also, for runoff, rivers, and atmospheric fallout transportation, smaller MP size ranges can be expected to be associated with a longer range from the initial plastic sources31,32,33. Nevertheless, more detailed information about MP polymer/color features within the size ranges are needed to achieve stronger conclusions about potential long/short-range sources.Figure 6Microplastics abundance (particles kg−1) by color in sea salt samples from stations S1 to S8.Full size imageFigure 7Microplastics abundance (particles kg−1) by size range in sea salt samples from stations S1 to S8.Full size imageMP polymer compositionFour types of polymer, namely polypropylene (PP), polystyrene (PS), polyethylene (PE), and polyethylene terephthalate (PET), were identified with FT-MIR-NIR (Supplementary Figure S1). These results are in accordance with those reported for salt samples in other studies worldwide (Table 1). These polymer types are widely used in daily life products, packaging, single-use plastics, and clothes, contributing to plastic pollution worldwide21. PET presented the highest contribution at all sampling sites, at ~ 48%, whereas PS was found to be least, at ~ 15% (Fig. 8, Table 1). Iñiguez et al.34 also reported PET predominance (83.3%) in Spanish table salt samples. PET predominance could be explained by its high density (1.30 g cm−3), making particles prone to sedimentation during the salt crystallization process19. PE (0.94 g cm−3), PP (0.90 g cm−3), and PS (1.05 cm−3) presented lower or similar densities to seawater (~ 1.02 g cm−3), making these more prone to flotation and possible loss due to wind during desiccation.Figure 8Microplastics abundances (particles kg−1) by polymer composition in sea salt samples from stations S1 to S8.Full size imageRisks assessmentDuring degradation, MP tends to emit monomers and different types of additives, these having the potential to cause harm to ecological systems and health18, 35. Results for the polymeric risks indices are presented in Fig. 9. According to polymer risk classification, all salts samples showed low risks, being similar to the entire study area. To date, none of the published studies have applied chemometric models in evaluating MP pollution in salts, posing difficulties when comparing our results. Information on the hazards of MP from ingestion to human health is still highly unclear. Other than exposure, the destiny and transit of ingested MP in the human body, including intestinal digestion and biliary discharge, have not been determined in previous research and remained largely unclear36. Some studies conducted impact assessments based on in vitro models37,38. However, whether the exposure concentrations used in such studies indicate the MP consumed and collected in humans is inconclusive. Previous studies found that toxicity, oxidative stress, and inflammation could result from MP exposure, including immune disruption and neurotoxicity effects, among others39. Therefore, an immediate effort is required to assess the health consequences of these MP when they reach the human body.Figure 9Polymeric risk indices for MP types in salts from stations S1 to S8.Full size image More

AffiliationsDepartment of Ocean Sciences, University of California, Santa Cruz, CA, USARachel A. FosterDepartment of Biogeochemistry, Max Planck Institute for Marine Microbiology, Bremen, GermanyRachel A. Foster, Daniela Tienken, Sten Littmann & Marcel M. M. KuypersDepartment of Ecology, Environment and Plant Sciences, Stockholm University, Stockholm, SwedenRachel A. FosterDepartment of Geosciences, Swedish Museum of Natural History, Stockholm, SwedenMartin J. WhitehouseDepartment of Oceanography, University of Hawai’i at Mānoa, Honolulu, HI, USAAngelicque E. WhiteAuthorsRachel A. FosterDaniela TienkenSten LittmannMartin J. WhitehouseMarcel M. M. KuypersAngelicque E. WhiteCorresponding authorCorrespondence to
Rachel A. Foster. More

Well-mixed populationAs shown in the Methods section, in a well-mixed population, the model can be described in terms of the replicator-mutator dynamics. I begin, by a case where the quality of the two resources are similar, r1 = r2 = r, and plot the frequency (solid blue), the average payoffs from the game (dashed red), and the amplitude of fluctuations (dotted blue) for different strategies in Fig. 1a–d. Here, the replicator dynamic is solved starting from a uniform initial condition in which all the strategies’ initial frequency is equal. The results of simulations in finite populations are in good agreement with the replicator dynamics results (see Supplementary Note 2 and Supplementary Figs. 1, 2, 3, and 4 for comparison to simulations). Throughout this manuscript I fix c = 1.Fig. 1: The frequency, game payoffs, and the amplitude of fluctuations for different strategies.The frequency (solid blue), game payoffs (dashed red), and amplitude of fluctuations (dotted blue) of costly cooperators (a), costly defectors (b), non-costly cooperators (c), and non-costly defectors (d), as a function of the enhancement factor, r. As r increases, above a first threshold (r* = 1 + cg), cooperation in the costly institution evolves, and above a second threshold (approximately r = 2) cooperation in both the costly and the free institutions evolves. For medium r, the system shows periodic fluctuations. Parameter values: g = 5, nu = 10−3, π0 = 2, cg = 0.398. The replicator dynamic, derived in the Methods Section, is solved for 9000 time steps, and the time averages are taken over the last 2000 time steps.Full size imageFor small enhancement factors r, the dynamics settle in a fixed point where only defectors in the free institution survive. The advantage of the costly institution becomes apparent as r increases beyond r* = 1 + cg, such that the maximum possible payoff of the costly institution, which is achieved when nobody defects and is equal to r − 1 − cg becomes positive. As shown in the Methods, a focal defector in a group composed of ({n}_{C}^{1}) costly cooperators and ({n}_{D}^{1}) other costly defectors (and ({n}_{C}^{2}+{n}_{D}^{2}) individuals who prefer the free resource), obtains a payoff of ({n}_{C}^{1}r/({n}_{D}^{1}+{n}_{C}^{1}+1)-{c}_{g}). This payoff becomes negative for a small enough value of ({n}_{C}^{1}) (or a large enough value of ({n}_{D}^{1})). Since groups with a small number of costly cooperators are drawn with a high probability when ({rho }_{C}^{1}) is small (these probabilities can be derived in terms of multinational coefficients, see the Methods), the average payoff of a costly defector remains negative in this regime (for instance, given at the transition ({rho }_{C}^{1}approx nu), a mutant costly defector finds itself in a group with no costly cooperator with probability ({(1-{rho }_{C}^{1})}^{g-1}approx {(1-nu )}^{g-1}), which is close to 1 for low mutation rates, and pays a pure cost of −cg). On the other hand, a costly cooperator’s payoff is equal to (({n}_{C}^{1}+1)r/({n}_{D}^{1}+{n}_{C}^{1}+1)-{c}_{g}-1), which for small enough ({n}_{D}^{1}) becomes positive. As in this region, the frequency of costly defectors, ({rho }_{D}^{1}), is small, such group compositions occur with a high probability (at the transition, ({rho }_{D}^{1}approx nu), and thus the probability that a costly cooperator joins a group with no costly defectors is ({(1-{rho }_{D}^{1})}^{g-1}approx {(1-nu )}^{g-1}), which is close to 1 for low mutation rates). Consequently, the average payoff of costly cooperators from the game becomes positive, and thus, larger than the dominant non-costly defectors’ payoff (who receive a payoff of zero). Consequently, the frequency of costly cooperators rapidly increases at r*. However, due to the rapid increase in the frequency of costly cooperators at r*, the probability of formation of such mixed groups increases, and costly defectors start to appear in the system. Further increasing r in this region, the frequency of costly cooperators and costly defectors increases at the expense of non-costly defectors.As r increases above a second threshold, cyclic fluctuations set in, and the dynamics settle in a periodic orbit. An example of this periodic orbit is presented in Fig. 2a, b. Interestingly, the average payoff of costly cooperators, costly defectors, and non-costly defectors in this region remains close to zero despite the evolution of cooperation. Although individuals constantly update their strategy to overcome others, no strategy wins in the evolution. Instead, individuals engage in a winnerless red queen dynamic. The game payoffs of costly cooperators and costly defectors fluctuate around zero (which is equal to the game payoff of non-costly defectors). The dynamics of the system in this regime resembles the frequency-dependent selection in the host-parasite evolution, coined the red queen dynamic based on the fact that no matter how much they run, all end up in the same place53,54. On this basis, I call this periodic orbit the red queen periodic orbit.Fig. 2: Red queen and black queen orbits.The frequency of different strategies (a) and the game payoffs (b) in the red queen, and the black queen (c, d) periodic orbits. In the red queen orbit, cooperators in the costly institution survive. However, the payoff of the surviving strategies fluctuates around zero, and none dominate others. In contrast, cooperators in both institutions evolve in the black queen orbit, and cooperators of each type suppress defection in their opposite institution. Consequently, the payoff of all the strategies starts to deviate from zero. Parameter values: g = 5, nu = 10−3, π0 = 2, and cg = 0.398. In (a, b) r = 1.7, and in (c, d) r = 2.2.Full size imageThe existence of a costly institution can facilitate the evolution of cooperation in its competing free institution too. As the amplitude of fluctuations increases, episodes where most of the individuals prefer the costly institution occur. During these episodes, ({rho }_{D}^{2}) drops to a small value. Consequently, the probability that a mutant non-costly cooperator finds itself in a group devoid of non-costly defectors ({(1-{rho }_{D}^{2})}^{g-1}), increases. In such groups, non-costly cooperators receive a payoff of r−1, which is larger than the payoff of all the other strategies and outcompete other strategies. At this point, a second periodic orbit emerges in which cooperation in both the costly and free institutions evolves. The evolution of cooperation in the free institution can, in turn, have a positive impact on cooperation in the costly institution. This is the case because above the point where cooperation in the free institution evolves, the frequency of individuals who prefer the free institution starts to increase by increasing r. This effect decreases the frequency of those who prefer the costly institution and its effective size. This decreases the mixing probability between costly cooperators and costly defectors and increases the costly cooperators’ payoffs. Consequently, a functional complementation between cooperators with different game preferences emerges, which is reminiscent of a black queen dynamics in which different types crucially depend on each other for performing vital functions36,55. While vulnerable to defectors in their own institution, cooperators complement each other by beating defectors in their opposite institution. Synergistically thus, they can suppress defection in the population and alternately dominate the population (see Fig. 2c, d). At this stage, the game payoff of all the strategies starts to increase beyond zero. I call this periodic orbit the black queen orbit.The picture depicted above is the typical behavior of the model for large enough values of the cost. To see this, in Fig. 3a, I plot the phase diagram of the model in the cg − r plane. Here, the frequency of cooperators in the population, ({rho }_{C}={rho }_{C}^{1}+{rho }_{C}^{2}), is color plotted as well (see Supplementary Fig. 1 for the frequency of different strategies). Red dashed lines show the boundary of the region where the system settles in a periodic orbit. For high costs, as r increases, the system shows a series of successive cross-overs from a defective fixed point to the red queen periodic orbit, black queen periodic orbit, and finally a cooperative fixed point. On the other hand, for small costs, the system possesses a bistable region where both the red queen and black queen periodic orbits (or a partially cooperative fixed point and black queen periodic orbit to the left of the red dashed line in the bistable region) are stable, and the system shows a discontinuous transition between these two orbits. Orange circles show the lower boundary of the bistable region, below which the black queen orbit is unstable. Its upper boundary, above which the red queen orbit becomes unstable, is plotted by red squares, in Fig. 3. The transition between the two periodic orbits becomes a continuous transition at a single critical point (see Supplementary Fig. 4).Fig. 3: Evolution of cooperation.a Time average total frequency of cooperators, ({rho }_{C}={rho }_{C}^{1}+{rho }_{C}^{2}) in the r − cg plane is color plotted. The dynamics can settle in fixed point (FP) (small and large enhancement factors), or two different periodic orbits, red queen periodic orbit (RPO) where cooperation only in the costly institution evolve and black queen periodic orbit (BPO) where cooperation in both institutions evolve. b Time average difference between the probability that an individual in the costly institution is a cooperator from the probability that an individual in the free institution is a cooperator, (gamma ={rho }_{C}^{1}/({rho }_{C}^{1}+{rho }_{D}^{1})-{rho }_{C}^{2}/({rho }_{C}^{2}+{rho }_{D}^{2})). Individuals are more likely to be cooperators in a costly institution. c The time average total frequency of cooperators in the r − cg plane under pure selection dynamic (ν = 0). Red queen and black queen periodic orbit can occur for, respectively, small and large enhancement factors. In other regions, the dynamics settle in a fixed point where either non-costly defectors (small enhancement factors), costly cooperators (inside the region marked with dashed black line), or non-costly cooperators survive. Parameter values: g = 5, and π0 = 2. In (a, b) ν = 10−3, and in (c) ν = 0. In (a, b) the replicator dynamic is solved for 8000 time steps, and the time average is taken over the last 2000 steps. In (c) the replicator dynamic is solved for 200,000 time steps, and the time average is taken over the last 150,000 time steps.Full size imageExamination of the overall cooperation in the population shows that an entrance cost has a contrasting effect on population cooperation for large and small enhancement factors. An entrance cost keeps free-riders away from a costly institution. This fact makes the relative frequency of cooperators to defectors higher in the costly institution than that in the free institution. To see that defectors are less likely to join the costly institution, I plot the difference between the probabilities that an individual in the costly institution is a cooperator and the probability that an individual in the free institution is a cooperator, (gamma ={rho }_{C}^{1}/({rho }_{C}^{1}+{rho }_{D}^{1})-{rho }_{C}^{2}/({rho }_{C}^{2}+{rho }_{D}^{2})) in Fig. 3b, where it can be seen it is always positive. Intuitively, as a costly defector’s payoff in a group with ({n}_{C}^{1}) cooperators and ({n}_{D}^{1}) other defectors is equal to (r{n}_{C}^{1}/({n}_{C}^{1}+{n}_{D}^{1})-{c}_{g}), a costly defector can reach a positive payoff only when ({n}_{C}^{1}) is large. Otherwise, costly defectors are better off hedging the risk of obtaining a negative payoff by joining the free institution, where their payoff is necessarily non-negative. Consequently, the expected number of cooperators in the costly institution, ({rho }_{C}^{1}g), sets a bound for the frequency of costly defectors. This fact increases a costly institution’s profitability, especially for small enhancement factors, and positively impacts cooperation in the population. On the other hand, for high enhancement factors, a large entrance cost is detrimental to cooperation. This is because, although the frequency of defectors in the costly institution remains close to zero, fewer individuals are willing to choose a costly institution with a high cost. This increases the effective size of the free institution and the mixing between cooperators and defectors in the free institution. Since defectors can better exploit cooperators in well-mixed groups, the increased mixing between cooperators and defectors in the free institution hinders cooperation.As shown in the Supplementary Note 3, while the phenomenology of the model remains the same for lower mutation rates, lower mutation rates increase the size of the region where the dynamics settle in a periodic orbit (see Supplementary Fig. 5). Regarding the dependence of the dynamics on the mutation rate, an interesting case is the zero mutation rate, where selection is the sole driver of the dynamics. The time average cooperation for zero mutation rate, starting from an initial condition where all the strategies are equal, is plotted in Fig. 3c (See Supplementary Fig. 6 for the frequency of different strategies). Both the red queen (for small enhancement factors) and the black queen (for large enhancement factors) periodic orbits are observed in this case. However, for zero mutation rate, both the amplitude and period of fluctuations increase: The fluctuating dynamics go through periods where one of the surviving strategies reaches a frequency close to 1 only to be later replaced by another strategy (see Supplementary Fig. 7). The dynamics can also settle in different fixed points. For cg = 0, depending on the enhancement factor, either cooperators or defectors in both institutions survive in equal densities. For nonzero cg, however, only one of the strategies survives. For small enhancement factors, non-costly defectors dominate the population. For larger enhancement factors, either costly cooperators (the region marked with a dashed black line) or non-costly cooperators dominate the population.In the Supplementary Note 2, I consider a case where the two institutions have different productivities, i.e., different enhancement factors, and show that similar phases are at work in this case (see the Supplementary Figs. 2 and 3). For instance, I show that a large entrance cost destabilizes full defection, removes the system’s bistability, and ensures the evolution of cooperation starting from all the initial compositions of the population. In addition, I study the continuous replicator dynamics and show similar phenomenology is at work in this case (see Supplementary Notes 1.4 and 4, and Supplementary Figs. 8 and 9).Finally, I note that a similar phenomenology is at work in a context where instead of a costly and a cost-free institution, two costly institutions interact. To see this, assume institution 1 has a cost cg and institution 2 has a cost ({c}_{g}^{0}). Without loss of generality, assume ({c}_{g} , > , {c}_{g}^{0}). Writing ({c}_{g}=({c}_{g}-{c}_{g}^{0})+{c}_{g}^{0}), it is easy to see that it is possible to absorb ({c}_{g}^{0}) in the base payoff b (as all the individual pay a cost ({c}_{g}^{0}) irrespective of their institution choice). Thus, the model is equivalent to a context where resource 2 has zero cost, resource 1 has a cost of ({c}_{g}-{c}_{g}^{0}), and all the individuals receive a shifted base payoff of (b-{c}_{g}^{0}) (see Supplementary Note 5 and Supplementary Fig. 10).Structured populationIn contrast to the well-mixed population, the model shows no bistability in a structured population, and the fate of the dynamics is independent of the initial condition. To see why this is the case, I note that in a well-mixed population, a situation where all the individuals are defectors, and randomly prefer one of the two institutions, is the worst case for the evolution of cooperation, as in this case, mutant cooperators are in a disadvantage in both institutions. However, in a structured population, starting from such an initial condition, blocks of defectors, most of whom prefer the same institution, form due to spatial fluctuations. A mutant cooperator who prefers the minority institution in these blocks obtains a high payoff and proliferates. This removes the bistability of the dynamics in a structured population.To study the model’s behavior in a structured population, I perform simulations starting from an initial condition in which all the individuals are defectors and prefer one of the two institutions at random. The model shows similar behavior in a structured population to that in a well-mixed population. This can be seen in Fig. 4a–d, where the densities of different strategies are color plotted in the cg − r plane (see Supplementary Note 6 and Supplementary Figs. 11 and 12 for further details). As was the case in a well-mixed population, cooperation does not evolve for too small values of r. As r increases beyond a threshold, cooperation does evolve in the costly institution but not in the free institution. In this region, for a fixed enhancement factor, an optimal cost, approximately equal to cg = r − 1, optimizes the cooperation in the population. On the other hand, cooperation in both the costly and the free institutions evolves for large enhancement factors. In this region, increasing the cost can slightly increases defection in the free institution and have a detrimental effect on the evolution of cooperation, but not as much as it does in a well-mixed population.Fig. 4: The frequency of different strategies in the cg − r plane in a structured population.The time average frequencies of costly cooperators (a), costly defectors (b), non-costly cooperators (c), and non-costly defectors (d) in the cg − r plane are color plotted. The system shows a red queen dynamic in which cooperators only in the costly institution survive in large numbers (for smaller enhancement factors), or a black queen dynamic, where cooperators in both institutions survive and help each other to suppress defection (for larger enhancement factors). Parameter values: g = 5, nu = 10−3, and π0 = 2. The population resides on a 200 × 200 first nearest neighbor square lattice with von Neumann connectivity and periodic boundaries. The simulation is performed for 5000 time steps starting from an initial condition in which all the individuals are defectors and prefer one of the two institutions at random. The time average is taken over the last 2000 steps.Full size imageInstead of periodic orbits observed in the well-mixed population, on a spatial structure the model’s dynamic is governed by the cyclic dominance of different strategies through spatiotemporal fluctuations manifested by traveling waves. In Fig. 5, I present snapshots of the population’s stationary state in different phases. In this figure, I consider a model in which individuals reproduce with a probability proportional to the exponential of their payoff, π, times a selection parameter, β, (exp (beta pi )) (see the Supplementary Note 1.3), with β = 5. The situation in the model where individuals reproduce with a probability proportional to their payoff is similar. In Fig. 5a, I have set r1 = r2 = 1.7, and cg = 0.6. This phase corresponds to the red queen periodic orbit in the well-mixed population case. Here, the majority of the population are non-costly defectors. Costly cooperators experience an advantage over the former and can proliferate in the sea of non-costly defectors. However, costly cooperators are vulnerable to both costly defectors and non-costly cooperators. The former can only survive in small bands around costly cooperators, as they rapidly get replaced by non-costly defectors once they eliminate costly cooperators. This phenomenon shows that spatial competition between defectors with differing institution preferences can positively impact the evolution of cooperation. Non-costly cooperators, in turn, can survive by forming compact domains where they reap the benefit of cooperation among themselves. However, as the effect of network reciprocity is too small to promote cooperation in this region, non-costly cooperators get eliminated by non-costly defectors once costly cooperators are out of the picture. Consequently, the system’s dynamic is governed by traveling waves of costly cooperators followed by small trails of costly defectors and non-costly cooperators in a sea of non-costly defectors (see the Supplementary Video, SV.156, and Supplementary Note 7 for an illustration of the dynamics in this regime).Fig. 5: Snapshots of the population in the stationary state for different parameter values.Different strategies are color codded (legend). In (a,) r1 = 1.7, r2 = 1.7, in (b,) r1 = 3.5 and r2 = 3.5, and in c, r1 = 3, and r2 = 1.8. In all the cases cg = 0.6. For small enhancement factors (a), the red queen dynamics in which cooperators only in the costly institute survive in large numbers occur. For larger enhancement factors (b), the black queen dynamics in which cooperators in both institutions survive and help each other suppress defection occur. By increasing the enhancement factors (c), non-costly cooperators dominate. However, a small frequency of costly cooperators survives and purge the population from defectors by moving along the bands of non-costly defectors. Here, individuals reproduce with a probability proportional to the exponential of their payoff with a selection parameter equal to β = 5. The population resides on a 400 × 400 square lattice with von Neumann connectivity and periodic boundaries. Parameter values: g = 5 and ν = 10−3.Full size imageFigure 5b shows a snapshot of the population for r1 = r2 = 2.2. This phase corresponds to the black queen periodic orbit in a well-mixed population. In this phase, cooperators in both the costly and free institutions evolve. Cooperators are vulnerable to defectors in their institution and lose their territory to defectors of similar type. Defectors are in turn vulnerable to cooperators in their opposite institution and are replaced by them. Consequently, the dynamic of the model is governed by traveling waves of cooperators, chased by defectors of similar type, who are in turn extincted by cooperators of the opposite type. Thus, while cooperators of different types on their own either can not survive (in the case of non-costly cooperators) or are doomed to a winnerless competition with defectors (in the case of costly cooperators), they complement each other to efficiently suppress defection in the population (see the Supplementary Video, SV.256, for an illustration of the dynamics in this regime).Another manifestation of functional complementation between cooperators of different types can be seen in the regime of large enhancement factors. An example of this situation is plotted in Fig. 5c. Here, r1 = r2 = 3.5 and cg = 0.6. In this region, non-costly cooperators dominate the population. However, non-costly defectors can survive in small bands in the sea of non-costly cooperators. While at a disadvantage in the sea of non-costly cooperators, costly cooperators beat non-costly defectors. Consequently, small blocks of costly cooperators are formed within the bands of non-costly defectors. These blocks of costly cooperators move along the bands of non-costly defectors and purge the population from non-costly defectors. In this way, although costly cooperators exist only in small frequency, they play a constructive role in helping non-costly cooperators to dominate the population.In summary, the analysis of the spatial patterns reveals that competition or synergistic relation between individuals with different institution preferences plays an essential role in the evolution of cooperation in the system. Defectors with different institution preferences always appear as competitors who compete over space. By eliminating each other, they play a surprisingly constructive role in the evolution of cooperation. Cooperators, on the other hand, while having direct competition over scarce sites, can also act synergistically and help the evolution of cooperation in their opposite institution since they can eliminate defectors in their opposite institution. In this way, by purging defectors with an opposite game preference, cooperators help fellow cooperators with an opposite game preference. Consequently, cooperators with different game preferences can engage in a mutualistic relation to efficiently suppress defection in the population.Finally, as shown in the Supplementary Note 6, the spatial model shows similar phases in the case where the two public resources have heterogeneous profitability, that is, when r1 ≠ r2 (see the Supplementary Fig. 12). More

Data for all analysed radionuclides are presented in the “Supplementary Material”. Cryoconite samples were collected on Nalli Glacier (Supplementary Fig. S1) on Sept. 25, 2017 (samples 1701–1714) and on Sept. 10, 2018 (samples 1801–1814) at 28 spots (Fig. 2, Supplementary Table S1). Gamma spectrometric analysis of samples showed the presence of anthropogenic radionuclides 137Cs, 241Am, and 207Bi. All quoted radioactivity values were recalculated for the sampling date, except those for 241Am since the concentration of the parent 241Pu isotope is unknown. However, for this isotope, the correction for decay is negligible. The activity of 137Cs reached 8093 (± 69) Bq kg−1 of dry weight, that of 241Am reached 58.3 (± 2.3) Bq kg−1 and that of 207Bi reached 6.3 (± 0.6) Bq kg−1. The natural radionuclides 210Pb and 7Be were also present in all samples. The activity of 210Pb varied in the range of 1280–9750 Bq kg−1. In addition, in the investigated samples, a significant amount of short-lived cosmogenic radionuclide 7Be was found, whose specific activity reached 2418 (± 76) Bq kg−1 (Fig. 3, Supplementary Table S2). To evaluate the contribution of atmospheric components to the total 210Pb activity, 226Ra activity was determined and found to be 17–27 Bq kg−1 (Supplementary Table S2). Based on the 210Pb/226Ra ratio, we conclude that more than 98% of 210Pb was of atmospheric provenance.Figure 2Location of sampling points on Nalli Glacier. A—137Cs activity zone  95%) of corresponding rocks and numerous outcrops likely promoted entrapment of these elements into explosion clouds, and their subsequent precipitation with radionuclides. This feature of the geological structure of the area explains the extremely high enrichment of surface waters in elements such as Zn, Pb, Sr, Ni, As, Cr, Co, Se, Te, Cd, W, Cu, Sb, and Sn; for many of them, the excess reaches 10-fold with respect to the Clrake values51. This hypothesis is supported by obvious correlations between the concentrations of Bi, Ag, Sn, Sb, Pb, Cd, W, and Cu and the activity of anthropogenic radionuclides 137Cs, 241Am and 207Bi. This relationship is obviously related to the simultaneous release of elements and radionuclides from the contaminated ice layer and their entrapment in cryoconite holes. More