Ecology

The logic of the model
Our spatially-resolved model was simulated in discrete grid boxes of a 100 × 100 array, which included four basic assumptions: (1) Initial individuals were assumed to secrete three public goods but may randomly mutate to lose any of those functions with a certain probability; (2) Secreting a public good created a corresponding metabolic burden, therefore in losing a function the individual would gain a benefit; (3) All public goods were essential for growth. The net growth rates of individuals were dependent on the local concentrations of public goods; (4) Substrate and public goods diffused between two grid boxes at rates proportional to the concentration gradient.
For the 1st assumption, we included three functions because it is the minimal unit and tersest design to simulate complex communities, allows for the emergence of three categories of interaction patterns, and a single cooperative LOF genotype might evolve from differential evolutionary paths (Fig. 1A, B). The genotypes were described by bit strings containing 1 and 0 which indicated the genotype could produce the corresponding public good or not, respectively. Eight genotypes could emerge during the simulations, which were the initial autonomous producer [1, 1, 1], three one-function loss genotypes (OFLGs, i.e., [1, 1, 0], [1, 0, 1], and [0, 1, 1]), three two-function loss genotypes (TFLGs, i.e., [1, 0, 0], [0, 1, 0], and [0, 0, 1]), and a nonproducing cheater [0, 0, 0] (Fig. 1A).
Fig. 1: Logic of the individual-based model.

A Possible genotypes and evolutionary relationships among them emerging from reductive evolution when starting with an autonomous genotype that performs three essential public functions. Note that in this three-function model, some genotypes, i.e., Two-function loss genotypes and cheaters, might evolve from different mother genotypes. B Interaction patterns that could possibly be established in the spatially structured communities. C Schematic of the individual-based simulations. A 100 × 100 array initialization with all autonomous phenotypic individuals (left) was conducted with a long-term stepwise iteration to investigate if diverse interaction patterns could form (right). At each time step, calculations were done from the level of individual grids (top) to whole lattice (bottom). Within each grid box, Monod equation modified by basic assumptions of the Black Queen Hypothesis was used to calculate the microbial growth, while minimum and maximum thresholds of biomass were defined to decide the division and death of individuals (top middle). Microbial individuals were allowed to randomly mutate to lose functions (top middle). Classical discretization of the diffusion equation gave local rules for updating the concentrations of public goods and nutrients in each box (middle). State changes at the individual level lead to the evolutionary dynamics of the communities, which may give rise to the formation of diverse interaction patterns (bottom).

Full size image

The 2nd and 3rd assumptions were developed from the basic mathematical assumption of the BQH [19], and defined individual growth by integrating the benefit and cost of function loss (Fig. 1C). To conceptualize the cost of performing a function, we supposed a parameter (α) which is the fraction of biomass used to produce a public good per unit time of an individual. In addition, we defined a second parameter (β) as the ratio of the amount of public goods required during each step to account for the produced public goods. Therefore the redundant fraction of public goods production was 1−βj, and lower βj reflected a higher amount of redundant public goods that could be gained from the producers by the LOF genotypes, resulting in decreased risk in association with function loss (see Supporting Information S1 for more details). During the model simulation, spatiotemporal dynamic variables, i.e., positioning of genotypes and the time points at which genotypes evolved, would be collected. We initiated the simulations by randomly distributing 100 ancestor cells [1, 1, 1] into the grid boxes and iterated for at least 1,500,000 time steps. During each time step, individuals grew, decayed, reproduced, and mutated according to the previously mentioned assumptions (Fig. 1C). We paid attention to whether stable communities with various interdependent patterns could be formed after a specified number of iterations, as well as recorded the spatiotemporal dynamics of the communities.
Diverse interdependent patterns emerged with high level of function cost and varied level of functional redundancy
For model simulations, the function cost (parameter α) and functional redundancy (parameter β) were assigned to 0.0001, 0.0005, 0.001, and 0.4, 0.6, 0.8, respectively. A total of 2891 independent simulations with 9 parameter sets displayed different community structures (Fig. 2A). When the function cost was assigned to a low level, i.e., 0.0001, the autonomous ancestor dominated the community. When function costs were assigned to higher levels, 0.0005 and 0.001, new genotypes evolved and later interacted to form three distinct types of interdependent patterns even within the same α and β combination, i.e., asymmetric functional complementation (AFC), complete functional division pattern, and one-way dependency, with the relative amounts of 1677/2891, 143/2891, and 48/2891, respectively. In addition, higher functional redundancies favored the loss of more functions, increasing complexity of the community structures.
Fig. 2: Reductive evolution shapes diverse interdependent patterns in microbial communities.

A The final (steady state) community structures across gradients of function cost (α) and functional redundancy (1-β). Results were summarized from at least 300 interdependent runs for each parameter set. Community structures were assessed after simulation for 170,000 iterations, where 98.9% (2891/2923) of runs reached steady state. According to the structures, replicates were clustered into several scenarios for each parameter set, which are shown separately in the area plots. Note that the values of β is the proportion of public goods that is required for growth, and thus 1-β reflects the level of function redundancy. B, C Six representative community dynamics on the spatial lattices were selected from one interdependent simulation with the given conditions (mut = 10−5, α = 0.001, β = 0.8), showing the evolution of three types of asymmetric functional complementary pairs (AFCPs) (B), three different paths for the evolution of pairs [0, 0, 1] & [1, 1, 0] (C). Left images indicate the distribution of different genotypes at different points in evolutionary time. Curve plot in the middle describes the community dynamics of the corresponding simulation. Schematics at right briefly summarize the spatiotemporal dynamics of each simulation: the arrays in (B) indicate one type of AFCP directly dominated the communities without competition from others; the boxes in (C) indicates the composition of ancestor or AFCP in the related time points, while the windows inside indicate the spatial coexistence of multiple AFCPs and the size of the windows represents the relative fraction of different AFCPs.

Full size image

Among the three possible kinds of interactions, the AFC pattern was the most widespread, which was the combination of a two-function-loss genotype (TFLG) and its complementary one-function-loss genotype (OFLG). For example, [0, 0, 1], which produced a single essential public good, depended on its functional complement one-function-loss partner [1, 1, 0], for the other two public goods. Specifically, three types of the asymmetric functional complementary pairs (AFCPs), that is, [0, 0, 1] coupled with [1, 1, 0], [0, 1, 0] coupled with [1, 0, 1], and [1, 0, 0] coupled with [0, 1, 1], colonized most of the grid with a similar frequency of emergence. Interestingly, under the condition of high level of cost, the emergence of AFC patterns was accompanied by some nonproducing cheaters, whose relative abundance rose with the increase in functional redundancy (Fig. 2A top row). The addition of cheaters significantly reduced the total biomass of the communities, suggesting that high functional redundancy favors the evolution of cheaters which may decrease the community productivity. In addition, function loss happened more easily with high function cost. As the function cost parameter α increased from 0.0005 to 0.001, relative abundance of TFLGs increased approximately from 55 to 70% (Fig. 2A).
Besides the AFC patterns, two additional types of interdependent patterns evolved at a relatively lower frequency. The complete functional division pattern, that is, coexistence of [0, 0, 1], [0, 1, 0], and [1, 0, 0], only evolved when both factors were at high levels (α = 0.001, β = 0.4) with a frequency of approximately 45% (143 of 319 simulations, Fig. 2A, top right), which described a scenario with high benefit and low cost of function loss, favoring the loss of more functions and consequently more likely to maintain the evolution of TFLPs. Another form of interactions that emerged was one-way dependency, where one partner performs all functions and other none (i.e., coexistence of [1, 1, 1] and [0, 0, 0]). This form emerged at a low frequency (48 out of all 2891 simulations shown in Fig. 2A), but evolved with a higher probability under the condition of a mid-level function cost and low level of functional redundancy (α = 0.0005, β = 0.6, Fig. 2A, middle left), where the extinction of [1, 1, 1] was ~2.5 times slower than in other scenarios (Supplementary Fig. 1), leading to a higher potential for the spatial proximity between [1, 1, 1] and [0, 0, 0] during evolution.
Taken together, these phenomena demonstrated that the mutualistic exchange of complementary functions happened only when function cost was high. The emergence of different interdependent interaction patterns was related to the function cost and function redundancy, especially for the complete functional division and one-way dependency pattern, which only emerged within a limited parameter range. However, even for a given combination of α and β, it still remained possible for the evolution of distinct interaction patterns, suggesting that stochastic processes may play a role.
Same interdependent patterns might evolve via different modes
Because the evolution of three kinds of AFCPs were the most common scenarios in our simulations, we then focused on the role of stochastic processes, i.e., the key random events, in deciding the winning complementary pair among the three similar but different AFCPs. As a first step, we traced the variation in the spatiotemporal dynamics, trying to cluster the numerous evolutionary dynamics into limited modes and divide the complex evolutionary courses into several stages. These simplifications would facilitate the search for key random events.
Therefore, we analyzed the dynamics of 296 simulations with a typical parameter set (α = 0.001, β = 0.8), because under this condition, only the three types of AFCPs evolved, with a similar frequency of emergence (Fig. 2A, Top left), in order to avoid interference from the other interaction patterns. As described above, any of the three types of AFCPs could potentially take over the final community under this condition (Fig. 2B; Supplementary video 1–3). Using the emergence of AFCP [0, 0, 1] & [1, 1, 0] as an example, three categories of dynamic modes could give rise to its final domination. (1) After pair [0, 0, 1] & [1, 1, 0] emerged and formed a spatial aggregation, it rapidly expanded and took over the entire grid (Fig. 2C, first line; Supplementary video 3). (2) In addition to the pair [0, 0, 1] & [1, 1, 0], spatial aggregations of another AFCP also emerged (e.g., pair [0, 1, 0] & [1, 0, 1] in Fig. 2C, second line and Supplementary video 4). In this scenario, a special spatial pattern was established in a short period after the evolution of both AFCPs e, where pairs of two complementary members exhibited strong spatial mixing, while the two different AFCPs were totally segregated. Community succession was then governed by spatial competition between the two AFCPs. If pair [0, 0, 1] & [1, 1, 0] won the competition, it would dominate the final community. (3) Spatial aggregations of all three AFCPs emerged, and then pair [0, 0, 1] & [1, 1, 0] dominated the community after outcompeting the other two AFCPs (Fig. 2C, third line; Supplementary video 5). The clustering of these three possible modes of AFC patterns was also shown by the temporal dynamics of the α-diversity across different parameter sets (Supplementary Fig. 2), where the evolution modes of the AFC patterns were clearly clustered into three possible categories, suggesting that this clustering is independent of the determined factors α and β.
In sum, the succession of interdependent patterns could be divided into two stages: (1) the emergence of spatial aggregations composed of two interdependent members with strong connections; (2) spatial competition among different aggregations drive the community to evolve to the final state, composed of only one type of interdependent interactions. Of course, if only one type of AFCP emerged, the spatial competition stage would be unnecessary during succession.
Evolutionary random events play important roles in deciding the dominant AFCP in equilibrium communities
The presence of two evolutionary stages lead us to hypothesize that the random events affecting ecological outcomes should arise from two aspects. First, in the initial evolutionary stage, the emergence of interdependent spatial aggregations should be related to the order in which new genotypes emerge. Second, the outcome of the spatial competition should be also influenced by the initial positioning of the new genotypes.
The fact that each TFLP had two possible evolutionary paths (e.g., [1, 0, 0] could inherit its function from [1, 1, 0] or [1, 0, 1]), suggested that the effects of the random order of emergence for different genotypes were highly correlated with the evolutionary lineage. Therefore, to investigate the effects of this, we analyzed the evolutionary lineage of emergence, colonization, and loss of every genotype within the 296 simulations with the typical parameter set (α = 0.001, β = 0.8). In total, there were 24 evolutionary branches leading to the evolution of the three forms of AFC patterns (8 for each, Fig. 3). Among all these branches, we summarized two key random events (Fig. 3, red and blue boxes).
Fig. 3: The evolutionary trajectories of 296 independent simulations with the typical parameter set (mut = 10−5, α = 0.001, β = 0.8).

We analyzed the evolutionary trajectories of every interdependent run and clustered them into 24 types of branches (top, see Methods). The area plot shows the final community structures and the frequencies of each branch (bottom). Blue dashed box shows the evolutionary diversification into four scenarios after the first key event occurs, while the red dashed box indicates the 24 different evolutionary trajectories that diverged after the second key event occurs. Solid boxes with colored circles represent the genotypic composition of communities at different evolutionary time points. Red arrows indicate the branches where one type of asymmetric functional complementary pair (AFCP) directly dominated the communities without competition with other AFCPs, while the blue arrows indicate the branches where one type of AFCP took over the entire space after competitions with other AFCPs. Dashed boxes at the figure labels (right) indicate different AFCPs.

Full size image

The first event occurred after two types of OFLGs emerged. After this evolutionary time point, all three public functions were included in OFLGs. With the benefit of the function loss, these two OFLGs would expand and gradually outcompete the autonomous genotype [1, 1, 1]. Thus, the first key event was whether all three OFLGs could emerge before the autonomous genotype entirely disappeared (Fig. 3, blue box). If not, the third type of AFCP would never evolve; if so, all three types of AFCPs would still have a chance to dominate the final community. In the 296 simulations, the frequencies of these two scenarios were nearly same, that is, 147 simulations were clustered to the former, while 149 simulations were clustered to the latter. The 147 simulations, where the third type of AFCP never evolved, could be then divided into three categories with similar frequencies, where two of the three OFLGs occupied the whole space and excluded the ancestral population.
The second key evolutionary event was the emergence of TFLGs (Fig. 3, red box). After the two or three types of OFLGs successfully colonized, whose functional complementary TFLGs first to emerge in the next evolutionary time would lead to the prior formation of the spatial aggregation of the AFCP. It is obvious that if no other AFCP aggregations formed later, this AFCP would dominate the final community (Fig. 3, red arrow indicated branches). Alternatively, if other AFCP aggregations formed during the expansion process, the spatial competition between different AFCPs would decide the dominant AFCP in the equilibrium communities (Fig. 3, blue arrow indicated branches). In our analysis, the chance of only one AFCP evolving reached 64.7% (198 of the 296 simulations). If only two OFLGs evolved after the first event, the frequency of only one AFCP evolving reached 79.6% (121 of the 152 simulations). In contrast, if three OFLGs evolved after the first event, there could be a relative higher possibility of two or three AFCPs evolving (47.4%), meaning that spatial competition could then be an important process.
What decided the winner of the competition? We observed that after the segregated interdependent spatial pattern was newly established, the relative region sizes occupied by different AFCPs were the key to determining the winner (Fig. 2C, the second and third lines; Supplementary video 4 and 5). We analyzed the time gaps between the emergence of the two AFCPs in the second categories of succession modes and the size of the regions they occupied (Fig. 4A). The result indicated a significantly positive correlation between the length of the time gaps and the region size the prior AFCP occupied (t-test, p  1 indicates pair [0, 0, 1] & [1, 1, 0] is more spatially associated than pair [0, 1, 0] & [1, 0, 1]. Applying these definitions, the simulation results where the advantage of prior space occupancy was not significant (left side of blue line in Fig. 4C, 33 replicates) were selected for analysis, and we found a significantly positive correlation between the relative PAD at the beginning of spatial self-organization and the ‘region size advantage’ (Fig. 5A; p  1 means the prior emerged AFCPs are more spatially associated than the second AFCPs. Red dots indicate the first to emerge AFCP won the competition in the corresponding replicate, while the green dots indicate the second to emerge AFCP won the competition. B Two typical examples of simulations initialized with premixing the two types of AFCPs, [0, 0, 1] & [1, 1, 0] and [0, 1, 0] & [1, 0, 1], which represent scenarios when initial PAI001:010 1, respectively. C The significant positive correlation between the winning frequency of pair [0, 0, 1] & [1, 1, 0] and the initial value of PAI001:010. When initial PAI001:010  > 1, final communities were more likely to be dominated by pair [0, 0, 1] & [1, 1, 0], oppositely, pair [0, 1, 0] & [1, 0, 1] were more favorable when PAI001:010  More

1.
Viver T, Orellana LH, Díaz S, Urdiain M, Ramos‐Barbero MD, González‐Pastor JE, et al. Predominance of deterministic microbial community dynamics in salterns exposed to different light intensities. Environ Microbiol. 2019;21:4300–15.
CAS  PubMed  Article  Google Scholar 
2.
Tettelin H, Masignani V, Cieslewicz MJ, Donati C, Medini D, Ward NL, et al. Genome analysis of multiple pathogenic isolates of Streptococcus agalactiae: implications for the microbial “pan-genome”. Proc Natl Acad Sci USA. 2005;102:13950–5.
CAS  PubMed  Article  Google Scholar 

3.
Amann R, Rosselló-Móra R. After all, only millions? MBio. 2016;7:e00999–16.
PubMed  PubMed Central  Google Scholar 

4.
Konstantinidis KT, Ramette A, Tiedje JM. The bacterial species definition in the genomic era. Philos Trnas R Soc Lond B Biol Sci. 2006;361:1929–40.
Article  Google Scholar 

5.
Shapiro BJ, Polz MF. Ordering microbial diversity into ecologicaly and genetically cohesive units. Trends Microbiol. 2014;22:235–47.
CAS  PubMed  PubMed Central  Article  Google Scholar 

6.
McInerney JO, McNally A, O’Connell MJ. Why prokaryotes have pangenomes. Nat Microbiol. 2017;2:17040.
CAS  PubMed  Article  Google Scholar 

7.
Andreani NA, Hesse E, Vos M. Prokaryote genome fluidity is dependent on effective population size. ISMEJ. 2017;11:1719–21.
CAS  Article  Google Scholar 

8.
Cohan FM. What are bacterial species? Annu Rev Microbiol. 2002;56:457–87.
CAS  PubMed  Article  Google Scholar 

9.
Lan R, Reeves PR. When does a clone deserve a name? A perspective on bacterial species based on population genetics. Trends Microbiol. 2001;9:419–24.
CAS  PubMed  Article  Google Scholar 

10.
Fraser C, Alm EJ, Polz MF, Spratt BG, Hanage WP. The bacterial species challenge: making sense of genetic and ecological diversity. Science. 2009;323:741–46.
CAS  PubMed  Article  Google Scholar 

11.
Vázquez DP, Simberloff D. Ecological specialization and susceptibility to disturbance: conjectures and refutations. Am Nat. 2002;159:606–23.
PubMed  Article  Google Scholar 

12.
Prosser JI, Bohannan BJM, Curtis TP, Ellis RJ, Firestone MK, Freckleton RP, et al. The role of ecological theory in microbial ecology. Nat Rev Microbiol. 2007;5:384–92.
CAS  PubMed  Article  Google Scholar 

13.
Shade A, Peter H, Allison SD, Baho DL, Berga M, Bürgmann H, et al. Fundamentals of microbial community resistance and resilience. Front Microbiol. 2012;3:417.
PubMed  PubMed Central  Article  Google Scholar 

14.
Rodriguez-R LM, Overholt WA, Hagan C, Huettel M, Kostka JE, Konstantinidis KT. Microbial community successional patterns in beach sands impacted by the Deepwater Horizon oil spill. ISME J. 2015;9:1928–40.
PubMed  PubMed Central  Article  Google Scholar 

15.
Petraitis PS, Latham RE, Niesenbaum RA. The maintenance of species diversity by disturbance. Q Rev Biol. 1989;64:393–418.
Article  Google Scholar 

16.
Narasingarao P, Podell S, Ugalde JA, Brochier-Armanet C, Emerson JB, Brocks JJ, et al. De novo metagenomic assembly reveals abundant novel major lineage of Archaea in hypersaline microbial communities. ISME J. 2012;6:81–93.
CAS  PubMed  Article  Google Scholar 

17.
Antón J, Rosselló-Móra R, Rodriguez-Valera F, Amann R. Extremely halophilic bacteria in crystallizer ponds from solar salterns. Appl Environ Microbiol. 2000;66:3052–57.
PubMed  PubMed Central  Article  Google Scholar 

18.
Gomariz M, Martínez-García M, Santos F, Rodriguez F, Capella-Gutiérrez S, Gabaldón T, et al. From community approaches to single-cell genomics: the Discovery of ubiquitous hyperhalophilic Bacteroidetes generalists. ISME J. 2015;9:1–16.
Article  CAS  Google Scholar 

19.
Mora-Ruiz MR, Font-Verdera F, Díaz-Gil C, Urdiain M, Rodríguez-Valdecantos G, González G, et al. Moderate halophilic bacteria colonizing the phylloplane of halophytes of the subfamily Salicornioideae (Amaranthaceae). Syst Appl Microbiol. 2015;38:406–16.
CAS  Article  Google Scholar 

20.
Antón J, Lucio M, Peña A, Cifuentes A, Brito-Echeverría J, Moritz, F, et al. High metabolomic microdiversity within co-occurring isolates of the extremely halophilic bacterium Salinibacter ruber. PLoS ONE. 2013;8:e64701.
PubMed  PubMed Central  Article  CAS  Google Scholar 

21.
Conrad EC, Viver T, Hatt JK, Rosselló-Móra R, Konstantinidis KT. Unrestricted but ecologically-important gene-content diversity within a natural sequence-discrete population as revealed by sequencing of 112 isolates. 2020. In review.

22.
Cuadros-Orellana S, Martin-Cuadrado AB, Legault B, D’Auria G, Zhaxybayeva O, Papke RT, et al. Genomic plasticity in prokaryotes: the case of the square haloarchaeon. ISME J. 2007;1:235–45.
CAS  PubMed  Article  Google Scholar 

23.
Konopka A, Lindemann S, Fredrickson J. Dynamics in microbial communities: unraveling mechanisms to identify principles. ISME J. 2015;9:1488–95.
PubMed  Article  Google Scholar 

24.
Millán MM, Estrela MJ, Miró J. Rainfall components: variability and spatial distribution in a Mediterranean Area (Valencia Region). J Clim. 2005;18:2682–705.
Article  Google Scholar 

25.
Santos F, Moreno-Paz M, Meseguer I, López C, Rosselló-Móra R, Parro V, et al. Metatranscriptomic analysis of extremely halophilic viral communities. ISME J. 2011;5:1621–33.
CAS  PubMed  PubMed Central  Article  Google Scholar 

26.
Begon M, Townsend CR, Harper JL, editors. Ecology: from individuals to ecosystems. 4th ed. Malten, MA, USA: Blackwell Publishing Ltd; 2006.

27.
Rodriguez-R LM, Konstantinidis KT. Nonpareil: a redundancy-based approach to assess the level of coverage in metagenomics datasets. Bioinform. 2014;30:629–35.
CAS  Article  Google Scholar 

28.
Ondov BD, Treangen TJ, Melsted P, Mallonee AB, Bergman NH, Koren S, et al. Mash: fast genome and metagenome distance estimation using MinHash. Genome Biol. 2016;17:132.
PubMed  PubMed Central  Article  CAS  Google Scholar 

29.
Oksanen J, Kindt R, Legendre P.O’Hara B. Vegan: community ecology package. Com Ecol Pack. 2007;10:631–37.

30.
Pruesse E, Quast C, Knittel K, Fuchs BM, Ludwig W, Peplies J, et al. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 2007;35:7188–96.
CAS  PubMed  PubMed Central  Article  Google Scholar 

31.
Ludwig W, Strunk O, Westram R, Richter L, Meier H, Yadhukumar, et al. ARB; a software environment for sequence data. Nucleic Acids Res. 2004;32:1363–71.
CAS  PubMed  PubMed Central  Article  Google Scholar 

32.
Wu YW, Tang YH, Tringe SG, Simmons BA, Singer SW. MaxBin: an automated binning method to recover individual genomes from metagenomes using an expectation-maximization algorithm. Microbiome. 2014;2:26.
CAS  PubMed  PubMed Central  Article  Google Scholar 

33.
Rodriguez-R LM, Gunturu S, Harvey WT, Rosselló-Móra R, Tiedje JM, Cole JR, et al. The Microbial Genomes Atlas (MiGA) webserver: taxonomic and gene diversity analysis of Archaea and Bacteria at the whole genome level. Nucleic Acids Res. 2018;43:W282–8.
Article  CAS  Google Scholar 

34.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.
CAS  PubMed  Article  Google Scholar 

35.
Eren AM, Esen ÖC, Quince C, Vineis JH, Morrison HG, Sogin ML, et al. Anvi’o: an advanced analysis and visualization platform for ‘omics data. PeerJ. 2015;3:e1319.
PubMed  PubMed Central  Article  Google Scholar 

36.
Hyatt D, Chen GL, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinform. 2010;11:119.
Article  CAS  Google Scholar 

37.
UniProt Consortium. Uniprot: a hub for protein information. Nucleic Acids Res. 2015;43:D204–12.
Article  CAS  Google Scholar 

38.
Rodriguez-R LM, Konstantinidis KT. The enveomics collection: a toolbox for specialized analyses of microbial genomes and metagenomes. PeerJ Prepr. 2016;4:e1900v1.
Google Scholar 

39.
Caro-Quintero A, Konstantinidis KT. Bacterial species may exist, metagenomics reveal. Environ Microbiol. 2012;14:347–55.
CAS  PubMed  Article  Google Scholar 

40.
Goris J, Konstantinidis KT, Klappenbach JA, Coenye T, Vandamme P, Tiedje JM. DNA-DNA hybridization values and their relationship to whole-genome sequence similarities. Int J Syst Evol Microbiol. 2007;57:81–91.
CAS  PubMed  Article  Google Scholar 

41.
Viver T, Orellana LH, Hatt JK, Urdiain M, Díaz S, Richter M, et al. The low diverse gastric microbiome of the jellyfish Cotylorhiza tuberculata is dominated by four novel taxa. Environ Microbiol. 2017;19:3039–58.
PubMed  Article  Google Scholar 

42.
Haynes WM, Lide DR, Bruno TJ, editors. CRC handbook of chemistry and physics, 94th ed. London, UK: CRC Press; 2013. p. 4–89.

43.
Griebler C, Lueders T. Microbial biodiversity in groundwater ecosystems. Freshw Biol. 2009;54:649–77.
Article  Google Scholar 

44.
Konstantinidis KT, Tiedje JM. Towards a genome-based taxonomy for prokaryotes. J Bacteriol. 2005;187:6258–64.
CAS  PubMed  PubMed Central  Article  Google Scholar 

45.
Oren A. Microbial life at high salt concentrations: phylogenetic and metabolic diversity. Saline Syst. 2008;4:2.
PubMed  PubMed Central  Article  CAS  Google Scholar 

46.
Pedrós-Alió C. Marine microbial diversity: can it be determined? Trends Micribiol. 2006;14:257–63.
Article  CAS  Google Scholar 

47.
Azua-Bustos A, Fairén AG, González-Silva C, Ascaso C, Carrizo D, Fernández-Martínez MÁ, et al. Unprecedented rains decimate surface microbial communities in the hyperarid core of the Atacama Desert. Sci Rep. 2018;8:16706.
CAS  PubMed  PubMed Central  Article  Google Scholar 

48.
Allison SD, Martiny JBH. Resistance, resilience, and redundancy in microbial communities. Proc Natl Acad Sci U S A. 2008;105:11512–19.
CAS  PubMed  PubMed Central  Article  Google Scholar 

49.
Uritskiy G, Getsin S, Munn A, Gomez-Silva B, Davila A, Glass B, et al. Halophilic microbial community compositional shift after a rare rainfall in the Atacama Desert. ISME J. 2019;13:2737–49.
CAS  PubMed  PubMed Central  Article  Google Scholar 

50.
Ghai R, Pašić L, Fernández AB, Martin-Cuadrado AB, Mizuno CM, McMahon KD, et al. New abundant microbial groups in aquatic hypersaline environments. Sci Rep Nat 2011;1:135.
Article  CAS  Google Scholar 

51.
Burns DG, Janssen PH, Itoh T, Minegishi H, Usami R, Kamekura M, et al. Natronomonas moolapensis sp. nov., non-alkaliphilic isolates recovered from a solar saltern crystallizer pond, and emended description of the genus Natronomonas. Int J Syst Evol Microbiol. 2010;60:1173–76.
CAS  PubMed  Article  Google Scholar 

52.
López-Pérez M, Ghai R, Leon MJ, Rodríguez-Olmos Á, Copa-Patiño JL, Soliveri J, et al. Genomes of “Spiribacter”, a streamlined, successful halophilic bacterium. BMC Genom. 2013;14:787.
Article  CAS  Google Scholar 

53.
Martin-Cuadrado AB, Pašić L, Rodriguez-Valera F. Diversity of the cell-wall associated genomic island of the archaeon Haloquadratum walsbyi. BMC Genom. 2015;16:603.
Article  CAS  Google Scholar 

54.
Mirete S, Mora-Ruiz MF, Lamprecht-Grandío M, de Figueras CG, Rosselló-Móra R, González-Pastor J. Salt resistance genes revealed by functional metagenomics from brines and moderate-salinity rhizosphere within a hypersaline environment. Front Microbiol. 2015;6:1121.
PubMed  PubMed Central  Article  Google Scholar 

55.
Cray JA, Bell AN, Bhaganna P, Mswaka AY, Timson DJ, Hallsworth JE. The biology of habitat dominance; can microbes behave as weeds? Microb Biotechnol. 2013;6:453–92.
PubMed  PubMed Central  Article  Google Scholar  More

Study area
This project was conducted in the downstream segment of the Karnali River basin of Nepal (Fig. 1), which is the largest of Nepal’s three major river systems and is characterized by the steep terrain of the Himalayan Mountains. The highest runoff occurs during the monsoon season (e.g., June–October), and the lowest occurs during the winter season (e.g., December–May). Below the Siwalik Mountain range (a physiographic zone, Fig. 1), a vast network of small tributaries combines to form a single narrow channel of the Karnali River with well-defined banks. Originating from the Tibetan Plateau, the Karnali River is the largest tributary to the Ganges River in India, which harbours the most significant density of GRDs in the world. The lower Karnali River basin provides the furthest upstream range for GRDs, critically endangered gharials (Gavialis gangeticus), smooth Indian otters (Lutrogale perspicillata), and 36 native fish species49. The GRD population size in the Karnali River has declined from 26 to six individuals50. Such a sharp decline in the GRD population is due to the effects of habitat degradation, mainly from water-based development projects (i.e., water diversion,51). Concurrently, several upstream development projects are proposed, under construction, or completed [e.g., planned: the Karnali Chisapani multipurpose dam, 10,800 megawatt (MW); under construction: the upper Karnali hydropower project, 900 MW, and Bhari Babai diversion project; completed: Rani Jamara Kulariya irrigation intakes] and further threaten downstream aquatic life. All projects adopt traditional preconstruction environmental impact assessments procedure to define flow proportions (generally 10–20% of natural regimes) anecdotally and unscientifically. Thus, traditional flow proportions might be inadequate to sustain native aquatic biodiversity. Our study focused on the lower catchment area of the Karnali River basin, which is downstream from all megaprojects. All measurement protocols, including dolphin observation methods, were carried out in accordance with the Department of National Parks and Wildlife Conservation, Government of Nepal, guidelines and regulations. Habitat measurement protocols, including dolphin observation methods, were approved by the Department of National Parks and Wildlife Conservation, Government of Nepal (No 1129; 12 December 2016).
Available habitat assessment
Reduced water levels during the low-water season (e.g., December–May) escalate threats to aquatic biota by limiting physical habitat availability. Here, habitat refers to the hydro-physical habitat, which is defined by the flow and depth interactions at a particular geomorphic condition over space and time. Therefore, habitat availability (i.e. the area accessible to species) is assumed to be the greatest bottleneck, critically limiting species reproduction and survival42,51. We measured the available habitats in the low-water season when suitable habitat is critically limited (i.e., December–May in 2018/2019), excluding the monsoon season (June–November). Further, to capture dynamic flow variation within the dry (i.e., low) water season, we selected three temporal periods—March (mid dry season), May (late dry season), and December (early dry season)—based on 39 years of flow records available from the Department of Hydrology, Government of Nepal. Assessing available habitat includes habitat mapping and bank and instream surveys. We divided the study area into three segments [upper segment (S1): length = 11 km, average width = 218 m; middle segment (S2): length = 29 km, average width = 121 m; and lower segment (S3): length = 10 km, average width = 198 m; Fig. 1] based on uniform flow and channel geomorphology mapped along the selected stretch of the river. The three segments vary hydrologically and structurally. S1 consists of river channels with natural flows without any infrastructural diversion. Because of water diversion operations (e.g., Rani Jamuna irrigation intake and several traditional agricultural irrigation channels) and distributaries, the natural flow volume in S2 was low compared to that in S1. S3 benefited slightly from distributaries and received more water than S2.
Within each segment, the study reach (the linear segment where cross-sections are established) was established in such a way that the length of each reach was at least higher than the mean width (so the number varies among segments) of the respective segment. We also tried to maintain relatively similar flow at the top and bottom of the reach. Within each reach, random cross-sections were established to capture the hydraulic properties based on flow variation. As the flow variability of the stream increased, the number of cross-sections increased, and each section was kept at least 300 m apart from the other sections. Therefore, the number of cross-sections was based on the flow variation within a reach instead of the length of the reach. Bank and instream surveys started in an upstream direction, wherein directional readings of the cross-sections were noted. For the bank and instream measurements, pin heights were established at either side of the cross-section using GPS and a permanent reference marker for repeated flow measurements. Water surface elevations were estimated using a total station (an optical instrument for land surveying; Leica 772737 Builder 503) across the pin heights for each cross-section required for hydrological simulation. A new benchmark was established for each effort to measure the water surface elevation at each cross-section. The total number of cross-sections examined for the available habitats was 177 (March = 60, May = 47, and December = 60). The hydraulic parameters (see habitat characterization section below) at each cross-section were measured using a RiverSurveyor S5 acoustic water current profile reader [Sontek, Acoustic Doppler Profiler (ADP S5)], which records hydro-parameters continuously at a cell size between 0.02 to 0.5 m offering complete underwater available hydro-physical profile.
Occupied habitat assessment
We conducted a GRD population survey to capture occupied (selectivity) habitat characteristics (n = 97) at three temporal scales (previously mentioned) using the developed approach24. Within each temporal scale, we conducted three replications to capture the temporal and spatial variability in the characteristics of the occupied habitat. When we first detected dolphins, we observed surfacing behaviours for at least five minutes before establishing a cross-section. The habitats that were used for at least five minutes were considered occupied habitats, and then cross-sections were established to measure habitat characteristics using the ADP. If the dolphins disappeared after the location of the first sighting in less than five minutes, we excluded those habitats from our analysis. The Dolphin observation (only observation done) protocols were approved and permitted by the Department of National Parks and Wildlife Conservation, Government of Nepal.
Data analysis
Data preparation and software
The ADP S5 hydraulic data were imported into Excel databases (Microsoft v. 2010) to format for System for Environmental Flow Analysis (SEFA, version 1.5; Aquatic Habitat Analysts Inc.) software. All the hydraulic properties [depth (m), velocity (m/s), wetted perimeter-WP (m), width (m), cross-sectional area-CSA (m2), Froude number, and discharge (m3/s)], suitability, and flow regime determination were calculated using SEFA software and analysed at the cross-section and segment levels. The average flow of each segment was used as a base flow while running the habitat simulation model for the respective segment. We found critical flows ( 417 m3/s, excess flow with a negative contribution to habitat suitability) in May. Therefore, the habitat retention hydraulic simulation model was performed only with excess flow (for May) using 39 years of 90% exceedance flow (the flow that is equaled or exceeded 90% of the time).
Habitat characterization
The cross-sectional hydro-physical parameters—width, flow, depth, velocity, wetted perimeter, cross-sectional area, and habitat (types)—were reported spatially and temporally. The habitat type (e.g., pool, run, and riffle) was classified based on the Froude number (Fr), where Froude is an index of hydraulic turbulence (the ratio of velocity by the acceleration of gravity). Points with Froude numbers exceeding 0.41 were considered riffles, points with Froude numbers less than 0.18 were considered pools, and intermediate values were classified as run habitats. The proportion of run, riffle, and pool habitats within each study reach was calculated from the Froude numbers. The GRD’s seasonal hydro-physical habitats were characterized using basic descriptive statistics (mean and 95% CI). The variation in these hydraulic parameters among seasons, habitat types, and segments was examined by an analysis of variance (ANOVA), and post hoc pairwise comparisons were performed using Tukey’s honestly significant difference (HSD) test. A two-way ANOVA test was used to investigate any interactive effects of season and habitat on hydraulic variations. The level of significance was set at p  one and zero to those categories for which w ≤ one. By assigning one and zero to each group, we developed an HSC to calculate the area weighted suitability (AWS) at each measured point. Hydraulic habitat suitability is expressed as AWS in terms of usable area in metres of width or square metres per metre of reach (m2/m).
To obtain the AWS value for the reach, we multiplied the combined suitability index (CSI, which is the product of the suitability of depth and velocity at a point) and the proportion of the reach area represented by that point. Using a 39-year average base flow of 536.11 m3/s (90% exceedance flow) in May, we predicted the fluctuation (decrease by 10%) in the currently available maximum AWS (i.e., 22.718 m2/m, AWS of May) in the range of flows from 200–900 m3/s. We simulated the AWS in this particular range because this range represents the 39-year low and maximum values of the 90% exceedance flow for the low-water season (November–May). Covering this variation over a broader scale increases the applicability of our ecological thresholds across time. Using the same base flow and range, we also estimated the minimum flows that retain various standards (%) of habitat protection. Further, we also determined the minimum flow that provides the maximum AWS for the low-water season.
Ecological thresholds using flow-ecology relationships
As water depth and velocity are the result of instream habitat features, such as pools, riffles, and runs, we only incorporated depth and velocity when estimating the hydraulic habitat suitability. Additionally, GRD habitat selection is strongly guided by the depth and velocity of a river section24,51. Generalized linear models (GLMs) using logit functions were used to examine the relationship between GRD presence and hydraulic properties (depth and velocity). Four different GLMs (depth, velocity, depth*velocity, and depth + velocity) were developed, and the Akaike information criterion (AIC) was used to select the best models. The additive effect of depth and velocity on the GRD presence was found in the model with the best performance; therefore, we further used a generalized additive model (GAM) to capture the possible non-linear influence of depth and velocity on GRD presence. Because of the possibility of both linear and non-linear relationships11, we again used a GAM to capture the functional relationships between ecology (AWS) and flow. The degree of smoothness for all the GAMs identified by the iterative approach (up to 25 smoothing factors were checked) and the selected smoothing parameter (i.e., 20 for all the GAM models) that yielded a significant covariate (at the 0.005 level of significance) explained the maximum deviance and adjusted R2. Both the GLM and GAM models were fitted using the lm and mgcv packages in R Studio. More

1.
Belden, L. K. et al. Panamanian frog species host unique skin bacterial communities. Front. Microbiol. 6, 1171 (2015).
PubMed  PubMed Central  Article  Google Scholar 
2.
Jani, A. J. & Briggs, C. J. Host and aquatic environment shape the amphibian skin microbiome but effects on downstream resistance to the pathogen Batrachochytrium dendrobatidis are variable. Front. Microbiol. 9, 487 (2018).
PubMed  PubMed Central  Article  Google Scholar 

3.
McFall-Ngai, M. et al. Animals in a bacterial world, a new imperative for the life sciences. Proc. Natl. Acad. Sci. USA 110, 3229–3236 (2013).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

4.
Woodhams, D. C., Bletz, M., Kueneman, J. & McKenzie, V. Managing amphibian disease with skin microbiota. Trends Microbiol. 24, 161–164 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
Flechas, S. V. et al. Current and predicted distribution of the pathogenic fungus Batrachochytrium dendrobatidis in Colombia, a hotspot of amphibian biodiversity. Biotropica 49, 685–694 (2017).
Article  Google Scholar 

6.
Rollins-Smith, L. A. & Conlon, J. M. Antimicrobial peptide defenses against chytridiomycosis, an emerging infectious disease of amphibian populations. Dev. Comp. Immunol. 29, 589–598 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

7.
Jiménez, R. R. & Sommer, S. The amphibian microbiome: natural range of variation, pathogenic dysbiosis, and role in conservation. Biodiv. Conserv. 26, 763–786 (2017).
Article  Google Scholar 

8.
Bletz, M. C. et al. Host ecology rather than host phylogeny drives amphibian skin microbial community structure in the biodiversity hotspot of Madagascar. Front. Microbiol. 8, 1530 (2017).
PubMed  PubMed Central  Article  Google Scholar 

9.
Romano-Bertrand, S., Licznar-Fajardo, P., Parer, S. & Jumas-Bilak, E. Impact de l’environnement sur les microbiotes: focus sur l’hospitalisation et les microbiotes cutanés et chirurgicaux. Revue Francophone des Laboratoires 469, 75–82 (2015).
Article  Google Scholar 

10.
Woodhams, D. C. et al. Host-associated microbiomes are predicted by immune system complexity and climate. Genome Biol. 21, 23 (2020).
PubMed  PubMed Central  Article  Google Scholar 

11.
Cheng, Y. et al. The Tasmanian devil microbiome—implications for conservation and management. Microbiome 3, 1–11 (2015).
Article  Google Scholar 

12.
Lemieux-Labonté, V., Tromas, N., Shapiro, B. J. & Lapointe, F. J. Environment and host species shape the skin microbiome of captive neotropical bats. PeerJ 4, e2430 (2016).

13.
Grice, E. A. & Segre, J. A. The skin microbiome. Nat. Rev. Microbiol. 9, 244–253 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

14.
Daskin, J. H., Bell, S. C., Schwarzkopf, L. & Alford, R. A. Cool temperatures reduce antifungal activity of symbiotic bacteria of threatened amphibians–implications for disease management and patterns of decline. PLoS ONE 9, e100378 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

15.
Duellman, W. E. & Trueb, L. Integumentary, Sensory, and Visceral Systems. Biology of Amphibians (McGraw-Hill, New York, 1986).
Google Scholar 

16.
Bataille, A. et al. Susceptibility of amphibians to chytridiomycosis is associated with MHC class II conformation. Proc. R. Soc. B Biol. Sci. 282, 20143127 (2015).
Article  Google Scholar 

17.
Longo, A. V., Savage, A. E., Hewson, I. & Zamudio, K. R. Seasonal and ontogenetic variation of skin microbial communities and relationships to natural disease dynamics in declining amphibians. R. Soc. Open Sci. 2, 140377 (2015).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

18.
Kueneman, J. G. et al. The amphibian skin-associated microbiome across species, space and life history stages. Mol. Ecol. 23, 1238–1250 (2014).
PubMed  Article  Google Scholar 

19.
Kueneman, J. G. et al. Community richness of amphibian skin bacteria correlates with bioclimate at the global scale. Nat. Ecol. Evol. 3, 381–389 (2019).
PubMed  Article  Google Scholar 

20.
Rollins-Smith, L. A., Ramsey, J. P., Pask, J. D., Reinert, L. K. & Woodhams, D. C. Amphibian immune defenses against chytridiomycosis: impacts of changing environments. Integr. Comp. Biol. 51, 552–562 (2011).

21.
Sanchez, E. et al. Cutaneous bacterial communities of a poisonous salamander: a perspective from life stages, body parts and environmental conditions. Microb. Ecol. 73, 455–465 (2017).
CAS  PubMed  Article  Google Scholar 

22.
Antwis, R. E. et al. Ex situ diet influences the bacterial community associated with the skin of red-eyed tree frogs (Agalychnis callidryas). PLoS ONE 9, e85563 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

23.
Jani, A. J. & Briggs, C. J. The pathogen Batrachochytrium dendrobatidis disturbs the frog skin microbiome during a natural epidemic and experimental infection. Proc. Natl. Acad. Sci. USA 111, E5049–E5058 (2014).
ADS  CAS  PubMed  Article  Google Scholar 

24.
Medina, D. et al. Variation in metabolite profiles of amphibian skin bacterial communities across elevations in the Neotropics. Microb. Ecol. 74, 227–238 (2017).
CAS  PubMed  Article  Google Scholar 

25.
Loudon, A. H. et al. Vertebrate hosts as islands: dynamics of selection, immigration, loss, persistence, and potential function of bacteria on salamander skin. Front. Microbiol. 7, 333 (2016).
PubMed  PubMed Central  Article  Google Scholar 

26.
Bates, K. A. et al. Amphibian chytridiomycosis outbreak dynamics are linked with host skin bacterial community structure. Nat. Commun. 9, 1–11 (2018).
ADS  CAS  Article  Google Scholar 

27.
Woodhams, D. C. et al. Interacting symbionts and immunity in the amphibian skin mucosome predict disease risk and probiotic effectiveness. PLoS ONE 9, e96375 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

28.
Costa, S., Lopes, I., Proença, D. N., Ribeiro, R. & Morais, P. V. Diversity of cutaneous microbiome of Pelophylax perezi populations inhabiting different environments. Sci. Total Environ. 572, 995–1004 (2016).
ADS  CAS  PubMed  Article  Google Scholar 

29.
Sabino-Pinto, J. et al. Composition of the cutaneous bacterial community in Japanese amphibians: effects of captivity, host species, and body region. Microb. Ecol. 72, 460–469 (2016).
PubMed  Article  Google Scholar 

30.
Kueneman, J. G. et al. Inhibitory bacteria reduce fungi on early life stages of endangered Colorado boreal toads (Anaxyrus boreas). ISME J. 10, 934–944 (2016).
PubMed  Article  PubMed Central  Google Scholar 

31.
Becker, C. G., Longo, A. V., Haddad, C. F. B. & Zamudio, K. R. Land cover and forest connectivity alter the interactions among host, pathogen and skin microbiome. Proc. R. Soc. B Biol. Sci. 284, 20170582 (2017).
Article  Google Scholar 

32.
Belasen, A. M., Bletz, M. C., Leite, D. D. S., Toledo, L. F. & James, T. Y. Long-term habitat fragmentation is associated with reduced MHC IIB diversity and increased infections in amphibian hosts. Front. Ecol. Evol. 6, 236 (2019).
Article  Google Scholar 

33.
Greenspan, S. E. et al. Arthropod–bacteria interactions influence assembly of aquatic host microbiome and pathogen defense. Proc. R. Soc. B 286, 20190924 (2019).
CAS  PubMed  Article  Google Scholar 

34.
Becker, C. G. et al. Low-load pathogen spillover predicts shifts in skin microbiome and survival of a terrestrial-breeding amphibian. Proc. Roy. Soc. B 286, 20191114 (2019).
Article  Google Scholar 

35.
Christian, K., Weitzman, C., Rose, A., Kaestli, M. & Gibb, K. Ecological patterns in the skin microbiota of frogs from tropical Australia. Ecol. Evol. 8, 10510–10519 (2018).
PubMed  PubMed Central  Article  Google Scholar 

36.
McKenzie, V. J., Bowers, R. M., Fierer, N., Knight, R. & Lauber, C. L. Co-habiting amphibian species harbor unique skin bacterial communities in wild populations. ISME J. 6, 588–596 (2012).
CAS  PubMed  Article  Google Scholar 

37.
Walke, J. B. et al. Amphibian skin may select for rare environmental microbes. ISME J. 8, 2207–2217 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

38.
Varela, B. J., Lesbarrères, D., Ibáñez, R. & Green, D. M. Environmental and host effects on skin bacterial community composition in Panamanian frogs. Front. Microbiol. 9, 298 (2018).
PubMed  PubMed Central  Article  Google Scholar 

39.
Loudon, A. H. et al. Microbial community dynamics and effect of environmental microbial reservoirs on red-backed salamanders (Plethodon cinereus). ISME J. 8, 830–840 (2014).
CAS  PubMed  Article  Google Scholar 

40.
Toledo, L. F. & Batista, R. F. Integrative study of Brazilian anurans: geographic distribution, size, environment, taxonomy, and conservation. Biotropica 44, 785–792 (2012).
Article  Google Scholar 

41.
Haddad, C. F. B. et al. Guide to the Amphibians of the Antic Forest: Diversity and Biology (Anolisbooks, São Paulo, 2013).
Google Scholar 

42.
Toledo, L. F., Becker, C. G., Haddad, C. F. & Zamudio, K. R. Rarity as an indicator of endangerment in Neotropical frogs. Biol. Conserv. 179, 54–62 (2014).
Article  Google Scholar 

43.
Sasso, T. et al. Environmental DNA characterization of amphibian communities in the Brazilian Atlantic forest: potential application for conservation of a rich and threatened fauna. Biol. Conserv. 215, 225–232 (2017).
Article  Google Scholar 

44.
Ribeiro, M. C., Metzger, J. P., Martensen, A. C., Ponzoni, F. J. & Hirota, M. M. The Brazilian Atlantic Forest: how much is left, and how is the remaining forest distributed? Implications for conservation. Biol. Conserv. 142, 1141–1153 (2009).
Article  Google Scholar 

45.
Ledru, M. P., Montade, V., Blanchard, G. & Hély, C. Long-term spatial changes in the distribution of the Brazilian Atlantic Forest. Biotropica 48, 159–169 (2016).
Article  Google Scholar 

46.
Joly, C. A., Metzger, J. P. & Tabarelli, M. Experiences from the Brazilian Atlantic F orest: ecological findings and conservation initiatives. New Phytol. 204, 459–473 (2014).
PubMed  Article  PubMed Central  Google Scholar 

47.
Myers, N., Mittermeier, R. A., Mittermeier, C. G., Da Fonseca, G. A. & Kent, J. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).
ADS  CAS  Article  Google Scholar 

48.
Eterovick, P. C. et al. Amphibian declines in Brazil: an overview 1. Biotropica 37, 166–179 (2005).
Article  Google Scholar 

49.
Becker, C. G., Fonseca, C. R., Haddad, C. F. B., Batista, R. F. & Prado, P. I. Habitat split and the global decline of amphibians. Science 318, 1775–1777 (2007).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

50.
Both, C. et al. Widespread occurrence of the American bullfrog, Lithobates catesbeianus (Shaw, 1802) (Anura: Ranidae), Brazil. S. Am. J. Herpetol. 6, 127–134 (2011).
Article  Google Scholar 

51.
Carvalho, T., Becker, C. G. & Toledo, L. F. Historical amphibian declines and extinctions in Brazil linked to chytridiomycosis. Proc. R. Soc. B Biol. Sci. 284, 20162254 (2017).
Article  Google Scholar 

52.
Haddad, C. F., Toledo, L. F. & Prado, C. P. Anfíbios da Mata Atlântica: guia dos anfíbios anuros da Mata Atlântica. Editora Neotropica (2008).

53.
Carnaval, A. C. O. Q., Toledo, L. F., Haddad, C. F. B. & Britto, F. B. (2005). Chytrid fungus infects high-altitude stream-dwelling Hylodes magalhaesi (Leptodactylidae) in the Brazilian Atlantic rainforest. Froglog 70, 3–4 (2005).

54.
Carnaval, A. C. O., Puschendorf, R., Peixoto, O. L., Verdade, V. K. & Rodrigues, M. T. Amphibian chytrid fungus broadly distributed in the Brazilian Atlantic Rain Forest. EcoHealth 3, 41–48 (2006).
Article  Google Scholar 

55.
Toledo, L. F., Britto, F. B., Araújo, O. G., Giasson, L. M. & Haddad, C. F. The occurrence of Batrachochytrium dendrobatidis in Brazil and the inclusion of 17 new cases of infection. S. Am. J. Herpetol. 1, 185–191 (2006).
Article  Google Scholar 

56.
Toledo, L. F., Haddad, C. F. B., Carnaval, A. C. O. Q. & Britto, F. B. A Brazilian anuran (Hylodes magalhaesi: Leptodactylidae) infected by Batrachochytrium dendrobatidis: a conservation concern. Amphib. Reptile Conserv. 4, 17–21 (2006).
Google Scholar 

57.
de Oliveira Ramalho, A. C., De Paula, C. D., Catao-Dias, J. L., Vilarinho, B. & Brandao, R. A. First record of Batrachochytrium dendrobatidis in two endemic Cerrado hylids, Bokermannohyla pseudopseudis and Bokermannohyla sapiranga, with comments on chytridiomycosis spreading in Brazil. North West. J. Zool. 9, 145–150 (2013).
Google Scholar 

58.
Rodriguez, D., Becker, C. G., Pupin, N. C., Haddad, C. F. B. & Zamudio, K. R. Long-term endemism of two highly divergent lineages of the amphibian-killing fungus in the Atlantic F orest of B razil. Mol. Ecol. 23, 774–787 (2014).
CAS  PubMed  Article  Google Scholar 

59.
Preuss, J. F., Lambertini, C., da Silva Leite, D., Toledo, L. F. & Lucas, E. M. Batrachochytrium dendrobatidis in near threatened and endangered amphibians in the southern Brazilian Atlantic Forest. North West. J. Zool 11, 360–362 (2015).
Google Scholar 

60.
Preuss, J. F., Lambertini, C., Leite, D. D. S., Toledo, L. F. & Lucas, E. M. Crossing the threshold: an amphibian assemblage highly infected with Batrachochytrium dendrobatidis in the southern Brazilian Atlantic forest. Stud. Neotrop. Fauna E 51, 68–77 (2016).
Article  Google Scholar 

61.
Valencia-Aguilar, A., Toledo, L. F., Vital, M. V. & Mott, T. Seasonality, environmental factors, and host behavior linked to disease risk in stream-dwelling tadpoles. Herpetologica 72, 98–106 (2016).
Article  Google Scholar 

62.
Becker, C. G., Rodriguez, D., Lambertini, C., Toledo, L. F. & Haddad, C. F. Historical dynamics of Batrachochytrium dendrobatidis in Amazonia. Ecography 39, 954–960 (2016).
Article  Google Scholar 

63.
Skerratt, L. F. et al. Spread of chytridiomycosis has caused the rapid global decline and extinction of frogs. EcoHealth 4, 125–134 (2007).
Article  Google Scholar 

64.
Fisher, M. C. et al. Emerging fungal threats to animal, plant and ecosystem health. Nature 484, 186–194 (2012).
ADS  CAS  PubMed  Article  Google Scholar 

65.
Bataille, A., Lee-Cruz, L., Tripathi, B. & Waldman, B. Skin bacterial community reorganization following metamorphosis of the fire-bellied toad (Bombina orientalis). Microb. Ecol. 75, 505–514 (2018).
CAS  PubMed  Article  Google Scholar 

66.
Scheele, B. C. et al. Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity. Science 363, 1459–1463 (2019).
ADS  CAS  PubMed  Article  Google Scholar 

67.
Scheele, B. C. et al. Response to Comment on “Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity”. Science 367, eaay2905 (2020).
CAS  PubMed  Article  Google Scholar 

68.
Woodhams, D. C. et al. Antifungal isolates database of amphibian skin-associated bacteria and function against emerging fungal pathogens: ecological archives E096–059. Ecology 96, 595 (2015).
Article  Google Scholar 

69.
Harris, R. N. et al. Skin microbes on frogs prevent morbidity and mortality caused by a lethal skin fungus. ISME J. 3, 818–824 (2009).
CAS  PubMed  Article  PubMed Central  Google Scholar 

70.
Muletz-Wolz, C. R. et al. Inhibition of fungal pathogens across genotypes and temperatures by amphibian skin bacteria. Front. Microbiol. 8, 1551 (2017).
PubMed  PubMed Central  Article  Google Scholar 

71.
Niederle, M. V. et al. Skin-associated lactic acid bacteria from North American bullfrogs as potential control agents of Batrachochytrium dendrobatidis. PLoS ONE 14, e0223020 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

72.
Becker, M. H. et al. Composition of symbiotic bacteria predicts survival in Panamanian golden frogs infected with a lethal fungus. Proc. R. Soc. B Biol. Sci. 282, 20142881 (2015).
Article  CAS  Google Scholar 

73.
Rebollar, E. A. et al. Skin bacterial diversity of Panamanian frogs is associated with host susceptibility and presence of Batrachochytrium dendrobatidis. ISME J. 10, 1682–1695 (2016).
CAS  PubMed  PubMed Central  Article  Google Scholar 

74.
Longo, A. V. & Zamudio, K. R. Environmental fluctuations and host skin bacteria shift survival advantage between frogs and their fungal pathogen. ISME 11, 349–361 (2017).
Article  Google Scholar 

75.
Longo, A. V. & Zamudio, K. R. Temperature variation, bacterial diversity and fungal infection dynamics in the amphibian skin. Mol. Ecol. 26, 4787–4797 (2017).
PubMed  Article  PubMed Central  Google Scholar 

76.
Lambertini, C. et al. Biotic and abiotic determinants of Batrachochytrium dendrobatidis infections in amphibians of the Brazilian Atlantic Forest. Fung. Ecol. 49, 100995 (2021).
Article  Google Scholar 

77.
Becker, C. G. et al. Variation in phenotype and virulence among enzootic and panzootic amphibian chytrid lineages. Fung. Ecol. 26, 45–50 (2017).
Article  Google Scholar 

78.
Fierer, N. & Jackson, R. B. The diversity and biogeography of soil bacterial communities. Proc. Natl. Acad. Sci. USA 103, 626–631 (2006).
ADS  CAS  PubMed  Article  Google Scholar 

79.
Lozupone, C. A. & Knight, R. Global patterns in bacterial diversity. Proc. Natl. Acad. Sci. USA 104, 11436–11440 (2007).
ADS  CAS  PubMed  Article  Google Scholar 

80.
Bahram, M. et al. Structure and function of the global topsoil microbiome. Nature 560, 233–237 (2018).
ADS  CAS  PubMed  Article  Google Scholar 

81.
Boyle, A. H. D. et al. Diagnostic assays and sampling protocols for the detection of Batrachochytrium dendrobatidis. Dis. Aquat. Organ. 73, 175–192 (2007).
PubMed  Article  Google Scholar 

82.
Boyle, D. G., Boyle, D. B., Olsen, V., Morgan, J. A. T. & Hyatt, A. D. Rapid quantitative detection of chytridiomycosis (Batrachochytrium dendrobatidis) in amphibian samples using real-time Taqman PCR assay. Dis. Aquat. Organ. 60, 141–148 (2004).
CAS  PubMed  Article  PubMed Central  Google Scholar 

83.
Lambertini, C., Rodriguez, D., Brito, F. B., Leite, D. S. & Toledo, L. F. Diagnóstico do fungo Quitrídio: Batrachochytrium dendrobatidis. Herpetol. Bras. 2, 12–17 (2013).
Google Scholar 

84.
Kozich, J. J., Westcott, S. L., Baxter, N. T., Highlander, S. K. & Schloss, P. D. Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Appl. Envrion. Microbiol. 79, 5112–5120 (2013).
CAS  Article  Google Scholar 

85.
Bletz, M. C. et al. Amphibian gut microbiota shifts differentially in community structure but converges on habitat-specific predicted functions. Nat. Commun. 7, 13699 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

86.
Kriger, K. M. & Hero, J. M. The chytrid fungus Batrachochytrium dendrobatidis is non-randomly distributed across amphibian breeding habitats. Divers. Distrib. 13, 781–788 (2007).
Article  Google Scholar 

87.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. J. R. Meteorol. Soc. 25, 1965–78 (2005).
Article  Google Scholar 

88.
Kwon, S., Park, S., Lee, B. & Yoon, S. In-depth analysis of interrelation between quality scores and real errors in illumina reads. In 2013 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 635–638 (2013).

89.
Rideout, J. R. et al. Subsampled open-reference clustering creates consistent, comprehensive OTU definitions and scales to billions of sequences. PeerJ 2, e545 (2014).
PubMed  PubMed Central  Article  Google Scholar 

90.
Amir, A. et al. Deblur rapidly resolves single-nucleotide community sequence patterns. MSystems 2, e00191–16 (2017).

91.
Bokulich, N. A. et al. Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing. Nat. Methods 10, 57–59 (2013).
CAS  PubMed  Article  Google Scholar 

92.
Caporaso, J. G. et al. PyNAST: a flexible tool for aligning sequences to a template alignment. Bioinformatics 26, 266–267 (2010).
CAS  PubMed  Article  Google Scholar 

93.
Team, R. C. R: A Language and Environment for Statistical Computing. (2013).

94.
Wickham, H. Ggplot: Using the Grammar of Graphics with R. (2009)

95.
Calcagno, V., Calcagno, M. V., Java, S. & Suggests, M. A. S. S. Package ‘glmulti’ (2020).

96.
Bates, K. A. et al. Captivity and infection by the fungal pathogen Batrachochytrium salamandrivorans perturb the amphibian skin microbiome. Front. Microbiol. 10, 1834 (2019).
PubMed  PubMed Central  Article  Google Scholar 

97.
Lin, D., Foster, D. P. & Ungar, L. H. VIF regression: a fast regression algorithm for large data. J. Am. Stat. Assoc. 106, 232–247 (2011).
MathSciNet  CAS  MATH  Article  Google Scholar 

98.
Kruskal, J. B. Nonmetric multidimensional scaling: a numerical method. Psychometrika 29, 115–129 (1964).
MathSciNet  MATH  Article  Google Scholar 

99.
Bray, J. R. & Curtis, J. T. An ordination of the upland forest communities of southern Wisconsin. Ecol. Monogr. 27, 326–349 (1957).
Article  Google Scholar 

100.
Oksanen, J., Blanchet, F. G., Kindt, R., Legendre, P. & Minchin, P. R. In 2012: Vegan: Community Ecology Package. R Package Version 2.0-5 (eds. Hara, O. et al.) (2014).

101.
De Caceres, M., Jansen, F. & De Caceres, M. M. Indicspecies: relationship between species and groups of sites. R package Version 1, (2016).

102.
Longo, A. V., Burrowes, P. A. & Zamudio, K. R. Genomic studies of disease-outcome in host–pathogen dynamics. Am. Zool. 54, 427–438 (2014).
Google Scholar 

103.
Assis, A. B. D., Barreto, C. C. & Navas, C. A. Skin microbiota in frogs from the Brazilian Atlantic forest: species, forest type, and potential against pathogens. PLoS ONE 12, e0179628 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

104.
Estrada, A. et al. Skin bacterial communities of neotropical treefrogs vary with local environmental conditions at the time of sampling. PeerJ 7, e7044 (2019).
PubMed  PubMed Central  Article  Google Scholar 

105.
Muletz-Wolz, C. R., Fleischer, R. C. & Lips, K. R. Fungal disease and temperature alter skin microbiome structure in an experimental salamander system. Mol. Ecol. 28, 2917–2931 (2019).
CAS  PubMed  Google Scholar 

106.
Xu, L. L. et al. Changes in the community structure of the symbiotic microbes of wild amphibians from the eastern edge of the Tibetan Plateau. Microbiol. Open 9, e1004 (2020).
Article  Google Scholar 

107.
Puschendorf, R. et al. Environmental refuge from disease-driven amphibian extinction. Conserv. Biol. 25, 956–964 (2011).
PubMed  Article  Google Scholar 

108.
Whitfield, S. M., Kerby, J., Gentry, L. R. & Donnelly, M. A. Temporal variation in infection prevalence by the amphibian chytrid fungus in three species of frogs at La Selva, Costa Rica. Biotropica 44, 779–784 (2012).
Article  Google Scholar 

109.
Ruggeri, J. et al. Seasonal variation in population abundance and chytrid infection in stream-dwelling frogs of the Brazilian Atlantic forest. PLoS ONE 10, e0130554 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

110.
Longo, A. V., Burrowes, P. A. & Joglar, R. L. Seasonality of Batrachochytrium dendrobatidis infection in direct-developing frogs suggests a mechanism for persistence. Dis. Aquat. Organ. 92, 253–260 (2010).
PubMed  Article  PubMed Central  Google Scholar 

111.
Ellison, S., Knapp, R. A., Sparagon, W., Swei, A. & Vredenburg, V. T. Reduced skin bacterial diversity correlates with increased pathogen infection intensity in an endangered amphibian host. Mol. Ecol. 28, 127–140 (2019).
PubMed  Article  PubMed Central  Google Scholar 

112.
Piovia-Scott, J. et al. Greater species richness of bacterial skin symbionts better suppresses the amphibian fungal pathogen Batrachochytrium dendrobatidis. Microb. Ecol. 74, 217–226 (2017).
PubMed  Article  PubMed Central  Google Scholar 

113.
Flechas, S. V. et al. Surviving chytridiomycosis: differential anti-Batrachochytrium dendrobatidis activity in bacterial isolates from three lowland species of Atelopus. PLoS ONE 7, e44832 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

114.
Muletz, C. R., Myers, J. M., Domangue, R. J., Herrick, J. B. & Harris, R. N. Soil bioaugmentation with amphibian cutaneous bacteria protects amphibian hosts from infection by Batrachochytrium dendrobatidis. Biol. Conserv. 152, 119–126 (2012).
Article  Google Scholar 

115.
Woodhams, D. C., Ramsey, J. P. & Rollins-Smith, L. A. Effects of cold temperature on antimicrobial peptide synthesis and release in northern leopard frogs, Rana pipiens. Integr. Comp. Biol. 45, 1099–1099 (2005).
Google Scholar 

116.
Bovo, R. P. et al. Physiological responses of Brazilian amphibians to an enzootic infection of the chytrid fungus Batrachochytrium dendrobatidis. Dis. Aquat. Organ. 117, 245–252 (2016).
PubMed  Article  PubMed Central  Google Scholar 

117.
Familiar López, M., Rebollar, E. A., Harris, R. N., Vredenburg, V. T. & Hero, J. M. Temporal variation of the skin bacterial community and Batrachochytrium dendrobatidis infection in the terrestrial cryptic frog Philoria loveridgei. Front. Microbiol. 8, 2535 (2017).
PubMed  PubMed Central  Article  Google Scholar  More