Ecology

Persian Gulf water and sediment samples along the oil pollution continuumWater and sediment samples were collected along the circulation current of the Persian Gulf from Hormuz Island [HW (SAMN12878178) and HS (SAMN12878113)), Asaluyeh area (AW (SAMN12878179) and AS (SAMN12878114)), and Khark Island (KhW (SAMN12878180) and KhS (SAMN12878115)] (Fig. 1). Physicochemical characteristics and Ionic content of the collected samples are presented in the Supplementary Table S1. The GC-FID analyses showed high TPH and polyaromatic hydrocarbon (PAH) concentrations in the Khark sediment (KhS) (Supplementary Table S2). The GC-SimDis analysis showed that C25–C38 HCs were dominant in the KhS (~ 60%), followed by  > C40 HCs (~ 14%) (Supplementary Fig. S1). Chrysene, fluoranthene, naphthalene, benzo(a)anthracene and phenanthrene were respectively the most abundant PAHs in KhS. This pollution could originate from oil spillage due to Island airstrikes during the imposed war (1980–1988), sub-sea pipeline failures, and discharge of oily wastewater or ballast water of oil tankers (ongoing for ~ 50 years)12. The TPH of other water and sediment samples was below the detection limit of our method ( 2%) representatives were enriched in KhS (Fig. 2B). The co-presence of the orders Methanosarcinals, Alteromanadales, and Thermotogae (Petrotogales) in the KhS hints at potential oil reservoir seepage around the sampling site since these taxa are expected to be present in oil reservoirs38. The main HC pollutants in the Asalouyeh are low molecular weight aromatic compounds that mainly influence the prokaryotic population in the water column and rarely precipitate into sediments hence the similarity of AS to HS microbial composition as they both experience low pollution rates.Apart from oil-degrading Proteobacteria (e.g. Alteromonadales, Rhodobacterales, and Oceanospirillales), a diversity of sulfur/ammonia-oxidizing chemolithoautotrophic Proteobacteria were present in these sediments although at lower abundances e.g., (Acidithiobacillales (KhS 1.8%), Chromatiales (HS 1.5, AS 1.1, KhS 0.85%), Ectothiorhodospirales (HS 3.75, AS 2.3, KhS 1.7%), Halothiobacillales (KhS 2.6%), Thiotrichales (HS 1.5, AS 1.1, KhS 0.3%), Thiohalorhabdales (HS 0.7, AS 1.2, KhS 0.5%), Thiomicrospirales (KhS 1.5%)) (Fig. 2B).Sulfate-reducing bacteria (SRB) in HS comprised up to 16.2% of the community (Desulfobacterales, NB1-j, Myxococcales, Syntrophobacterales, and Thermodesulfovibrionia). Similar groups along with Desulfarculales, comprised the SRB functional guild of the AS (~ 18.9%). In comparison, Desulfuromonadales and Desulfobacterales were the SRB representatives in KhS with a total abundance of only ~ 3.3%. The lower phylogenetic diversity and community contribution of SRBs in KhS hint at the potential susceptibility of some SRBs to oil pollution or that HC degraders might outcompete them (e.g., Deferribacterales). Additionally, KhS was gravel-sized sediment (particles ≥ 4 mm diameter), whereas HS and AS samples were silt and sand-sized sediments39. The higher oxygen penetration in gravel particles of KhS hampers anaerobic metabolism of sulfate/nitrate-reducing bacteria hence their lower relative abundance in this sample (Fig. 2B).Whereas in water samples, sulfur/ammonia-oxidizing chemolithoautotrophs such as Thiomicrospirales and sulfate/nitrate-reducing bacteria such as Desulfobacterales, NB1-j, Deferribacterales, Anaerolineales, Nitrosococcales, Nitrosopumilales, and Pirellulales were present in very small quantities (lower than 0.5% in each sample).Chronic exposure to oil pollution shapes similar prokaryotic communities as oil spill eventsWe analyzed the prokaryotic community composition of 41 oil-polluted marine water metagenomes (different depths in the water column) from Norway (Trondheimsfjord), Deepwater Horizon (Gulf of Mexico), the northern part of the Gulf of Mexico (dead zone) and Coal Oil Point of Santa Barbara; together with 65 oil exposed marine sediment metagenomes (beach sand, surface sediments and deep-sea sediments) originating from DWH Sediment (Barataria Bay), Municipal Pensacola Beach (USA) and a hydrothermal vent in Guaymas Basin (Gulf of California) in comparison with the PG water and sediment samples (in total 112 datasets) (Supplementary Table S3). This extensive analysis allowed us to get a comparative overview of the impact of chronic oil pollution on the prokaryotic community composition.Hydrocarbonoclastic bacteria affiliated to Oceanospirillales, Cellvibrionales (Porticoccaceae family), and Alteromonadales40 comprised a significant proportion of the prokaryotic community in samples with higher aliphatic compounds pollution e.g. DWHW.BD3 (sampled six days after the incubation of unpolluted water with Macondo oil), DWHW.he1, and DWHW.he2 (oil-polluted water samples incubated with hexadecane), DWHW.BM1, DWHW.BM2, DWHW.OV1 and DWHW.OV2 (sampled immediately after the oil spill in the Gulf of Mexico) (Fig. 3). Samples treated with Macondo oil, hexadecane, naphthalene, phenanthrene, and those taken immediately after the oil spill in the Gulf of Mexico had a significantly lower proportion of SAR11 due to the dominance of bloom formers and potential susceptibility of SAR11 to oil pollutants (Fig. 3).Figure 3The abundance of unassembled 16S rDNA reads from unassembled metagenomes of different oil-polluted water samples (41). Row names are microbial taxa at the order level. For taxa with lower frequency, the higher taxonomic level is shown (47 taxa in total). The right-hand dendrogram represents the clustering of rows based on the Pearson correlation. Columns are the name of water samples. Samples are clustered based on Pearson correlation and the color scale on the top left represents the row Z-score. Figure was plotted using “circlize” and “ComplexHeatmap” packages in R.Full size imageFlavobacteriales and Rhodobacterales were present in relatively high abundance in almost all oil-polluted water samples except for those with recent pollution. Samples named NTW5, NTW6, NTW11, NTW12, which were incubated with MC252 oil for 32–64 days, represented similar prokaryotic composition dominating taxa that are reportedly involved in degrading recalcitrant compounds like PAHs in the middle-to-late stages of the oil degradation process (Alteromonadales, Cellvibrionales, Flavobacteriales, and Rhodobacterales). Whereas at the earlier contamination stages, samples represented a different community composition with a higher relative abundance of Oceanospirillales (e.g., NTW8, NTW9, NTW15, NTW16, and NTW17 sampled after 0–8 days incubation) (Fig. 3).The non-metric multidimensional analysis of the prokaryotic community of 106 oil-polluted water and sediment samples, together with the PG samples, is represented in Fig. 4. Water and sediment samples expectedly represented distinct community compositions. The AW sample was placed near samples treated with phenanthrene and naphthalene in the NMDS plot showing the impact of aromatic compounds on its microbial community. The KhW sample was located near NTW13 in the plot, both of which had experienced recent oil pollution.Figure 4Non-metric multidimensional scaling (NMDS) of the Persian Gulf water and sediment metagenomes along with oil-polluted marine water and sediment metagenomes based on Bray–Curtis dissimilarity of the abundance of 16S rDNA reads in unassembled metagenomes at the order level. Samples with different geographical locations are shown in different colors. PG water and sediment samples are shown in red. Water and sediment samples are displayed by triangle and square shapes, respectively. Figure was plotted using “vegan” library in R.Full size imageThe orders Oceanospirillales, Alteromonadales, and Pseudomonadales were present in relatively high abundances in all oil-polluted water samples except for HW (PG input water) and samples collected from the northern Gulf of Mexico dead zone (GOMDZ) (Fig. 3). Persian Gulf was located in the proximity of the developing oxygen minimum zone (OMZ) of the Arabian Sea that is slowly expanding towards the Gulf of Oman41. Potential water exchange with OMZ areas could be the cause of higher similarity to the GOMDZ microbial community42.Our results suggest that water samples with similar contaminants and exposure time to oil pollution enrich for similar phylogenetic diversity in their prokaryotic communities (Fig. 3). Marine prokaryotes represent vertical stratification with discrete community composition across the depth profile. According to our analyses, the prokaryotic communities of the oil-polluted areas are consistently dominated by similar taxa regardless of sampling depth or geographical location. We speculate that the high nutrient input due to crude oil intrusion into the water presumably disturbs this stratification and HC degrading microorganisms are recruited to the polluted sites where their populations flourish.The inherent heterogeneity of the sediment prokaryotic communities is retained even after exposure to oil pollution, reflected in their higher alpha diversity (Supplementary Fig. S3). However, similar taxa dominate the community in response to oil pollution (Fig. 5).Figure 5The abundance of unassembled 16S rDNA reads from unassembled metagenomes of different oil-polluted sediment samples (65). Row names are microbial taxa at the order level. For taxa with lower frequency, the higher taxonomic level is shown (77 taxa in total). The right-hand dendrogram represents the clustering of rows based on the Pearson correlation. Columns are the name of sediment samples. Samples are clustered based on Pearson correlation and the color scale on the top left represents the row Z-score. Figure was plotted using “circlize” and “ComplexHeatmap” packages in R.Full size imageIn sediment samples, Deltaproteobacteria had the highest abundance, followed by Gammaproteobacteria representatives. Ectothiorhodospirales, Rhizobiales, Desulfobacterales, Myxococcales, and Betaproteobacteriales representatives were present in almost all samples at relatively high quantities (Fig. 5). Sulfate/nitrate-reducing bacteria were major HC degraders in sediment, showing substrate specificity for anaerobic HC degradation43. Desulfobacterales and Myxococcales were ubiquitous sulfate-reducers, present in almost all oil-polluted sediment samples44. Sulfate-reducing Deltaproteobacteria play a key role in anaerobic PAH degradation, especially in sediments containing recalcitrant HC types45. Members of Rhizobiales are involved in nitrogen fixation, which accelerates the HC removal process in the sediment samples46, and therefore their abundance increase in response to oil pollution (Fig. 5).Prokaryotes involved in nitrogen/sulfur cycling of sediments are defined by factors such as trace element composition, temperature, pressure, and more importantly, depth and oxygen availability. In oil-polluted sediment samples, the simultaneous reduction of available oxygen with an accumulation of recalcitrant HCs along the depth profile complicates the organic matter removal. However, anaerobic sulfate-reducing HC degrading bacteria will cope with this complexity47. Prokaryotic communities of HS and AS samples represented similar phylogenetic diversity (Figs. 4, 5). Their prokaryotic community involved in the nitrogen and sulfur cycling resembles the community of DWHS samples. The KhS sample had a similar prokaryotic community to deeper sediment samples collected from 30 to 40 cm depth (USFS3, USFS11, and USFS12) which could be due to our sampling method using a grab sampling device.Our results show that the polluted sediments’ sampling depth (surface or subsurface) defines the dominant microbial populations. Hydrocarbon degrading microbes had the ubiquitous distribution in almost all oil-polluted water and sediment samples including Oceanospirillales, Cellvibrionales, Alteromonadales, Flavobacteriales, Pseudomonadales, and Rhodobacterales. Mentioned orders along with Ectothiorhodospirales, Rhizobiales, Desulfobacterales, Myxococcales, and Betaproteobacteriales and also representatives of Deltaproteobacteria phylum dominated in sediment samples. However, their order of frequency varies depending on the type of oil pollution present at the sampling location and the exposure time.Genome-resolved metabolic analysis of the Persian Gulf’s prokaryotic community along the pollution continuumA total of 82 metagenome-assembled genomes (MAGs) were reconstructed from six sequenced metagenomes of the PG (completeness ≥ 40% and contamination ≤ 5%). Amongst them, eight MAGs belonged to domain Archaea and 74 to domain bacteria. According to GTDB-tk assigned taxonomy (release89) (https://data.gtdb.ecogenomic.org/releases/release89/), reconstructed MAGs were affiliated to Gammaproteobacteria (36.6%), Alphaproteobacteria (12.2%), Flavobacteriaceae (9.7%), Thermoplasmatota (5%) together with some representatives of other phyla (MAG stats in Supplementary Table S4).A collection of reported enzymes involved in the degradation of different aromatic and aliphatic HCs under both aerobic and anaerobic conditions was surveyed in the annotated MAGs of this study43,48,49,50. The KEGG orthologous accession numbers (KOs) of genes involved in HC degradation were collected, and the distribution of KEGG orthologues detected at least in one MAG (n = 76 genes) is represented in Fig. 6.Figure 6Hydrocarbon degrading enzymes present in recovered MAGs from the PG water and sediment metagenomes. Row names represent the taxonomy of recovered MAGs and their completeness is provided as a bar plot on the right side. The color indicates the MAG origin. The size of dots indicates the presence or absence of each enzyme in each recovered MAG. Columns indicate the type of hydrocarbon and in the parenthesis is the name of the enzyme hydrolyzing this compound followed by its corresponding KEGG orthologous accession number. Figure was plotted using “reshape2” and “ggplot2” packages in R.Full size imageA combination of different enzymes runs the oil degradation process. Mono- or dioxygenases are the main enzymes triggering the HC degradation process under aerobic conditions. Under anaerobic conditions, degradation is mainly started by the addition of fumarate or in some cases, by carboxylation of the substrate. Therefore, bacteria containing these genes will potentially initiate the degradation process that will be continued by other heterotrophs. Enzymes such as decarboxylase, hydroxylase, dehydrogenase, hydratase, and isomerases act on the products of initiating enzymes mentioned above through a series of oxidation/reduction reactions.Various microorganisms cooperate to cleave HCs into simpler compounds that could enter common metabolic pathways. Mono- or dioxygenases which are involved in the degradation of alkane (alkane 1-monooxygenase, alkB/alkM), cyclododecane (cyclododecanone monooxygenase, cddA), Biphenyl (Biphenyl 2, 3-dioxygenase subunit alpha/beta, bphA1/A2, Biphenyl-2, 3-diol 1, 2-dioxygenase, bphC), phenol (phenol 2-monooxygenase, pheA), toluene (benzene 1, 2-dioxygenase subunit alpha/beta todC1/C2, hydroxylase component of toluene-4-monooxygenase, todE), xylene (toluate/benzoate 1,2-dioxygenase subunit alpha/beta/electron transport component, xylX/Y/Z, hydroxylase component of xylene monooxygenase, xylM) and naphthalene/phenanthrene (catechol 1,2 dioxygenase, catA, a shared enzyme between naphthalene/phenanthrene /phenol degradation) were detected in recovered MAGs of the PG.The key enzymes including Alkylsuccinate synthase (I)/(II) (assA1/A2), benzylsuccinate synthase (BssA)/benzoyl-CoA reductase (BcrA), ethylbenzene dehydrogenase (EbdA), and 6-oxo-cyclohex-1-ene-carbonyl-CoA hydrolase (BamA) that are responsible for initiating the degradation of alkane, toluene, ethylbenzene and benzoate exclusively under anaerobic conditions were not detected in reconstructed MAGs of this study. Consequently, recovered MAGs of this study are not initiating anaerobic degradation via known pathways while they have the necessary genes to continue the degradation process started by other microorganisms.The MAG KhS_63 affiliated to Immundisolibacter contained various types of mono- or dioxygenases and had the potential to degrade a diverse range of HCs such as alkane, cyclododecane, toluene, and xylene (Fig. 6). Members of this genus have been reported to degrade high molecular weight PAHs51.Lutimaribacter representatives have been isolated from seawater and reported to be capable of degrading cyclohexylacetate52. We also detected enzymes responsible for alkane, cycloalkane (even monooxygenase enzymes), and naphthalene degradation under aerobic conditions and alkane, ethylbenzene, toluene, and naphthalene degradation under anaerobic conditions in KhS_39 affiliated to this genus (Fig. 6).The MAGs KhS_15 and KhS_26 affiliated to Roseovarius had the enzymes for degrading alkane (alkane monooxygenase, aldehyde dehydrogenase), cycloalkane, naphthalene, and phenanthrene under aerobic and toluene and naphthalene under anaerobic condition. PAHs degradation has been reported for other representatives of this taxa as well53.The MAGs KhS_11 (a representative of Rhodobacteraceae) and KhS_53 (Marinobacter) had alkB/alkM, KhS_27 (GCA-2701845), KhS_29 (UBA5862) and KhS_40 (from Porticoccaceae family) had cddA, KhS_13 and KhS_21 (UBA5335) and KhS_38 (Oleibacter) had both alkB/alkM and xylM genes. They were among microbes that were initiating the degradation of alkane, cycloalkane and xylene compounds. Other MAGs recovered from Khark sediment were involved in the continuation of the degradation pathway. For example, KhS_1 was affiliated to the genus Halomonas and had different enzymes to degrade intermediate compounds. Halomonas representatives have been frequently isolated from oil-polluted environments54. The phylum Krumholzibacteria has been first introduced in 2019 and reported to contain heterotrophic nitrite reducers55. Two MAGs, KhS_5 and KhS_10, were affiliated to this phylum and contained enzymes involved in the anaerobic degradation of toluene, phenol, and naphthalene (Fig. 6).The MAGs KhS_12 and KhW_31 affiliated to the genus Flexistipes, in Deferribacterales order, were reconstructed from both KhW and KhS samples. Deferribacterales are reported to be present in the medium to high-temperature oil reservoirs with HC degradation activity and also in high-temperature oil-degrading consortia56. The type strain of this species was isolated from environments with a minimum salinity of 3% and a temperature of 45–50 °C57. The presence of this genus in KhS could be due to natural oil seepage from the seabed as PG reservoirs mainly have medium to high temperature and high salinity. Enzymes involved in the degradation of alkane, phenol, toluene and naphthalene under anaerobic conditions were present in MAGs KhS_12 and KhW_31.As mentioned earlier, Flavobacteriales are potent marine indigenous HC degraders that bloom in response to oil pollution58. Flavobacteriales affiliated MAGs (KhW_2, KhW_3, AW_21, and AW_33) were recovered from KhW and AW and mostly contained enzymes that participate in the degradation of aromatic compounds under anaerobic conditions. KhW_2 and KhW_3 also had both alkB/M (alkane monooxygenase) and xylM enzyme, which initiates the alkane and xylene bioremediation in Khark water. Among other recovered MAGs from KhW sample, KhW_18 (UBA724), KhW_24 (clade SAR86), KhW_43 (UBA3478) had alkB/M, and xylM, KhW_24 (clade SAR86) had alkB/M and cddA, and KhW_28 (from Rhodobacteraceae family) had alkB/M and pheA genes in their genome to initiate the degradation process (Fig. 6).Marinobacter (KhW_15) was another MAG reconstructed from KhW sample. This genus is one of the main cultivable genera that play a crucial role in the bioremediation of a wide range of oil derivatives in polluted marine ecosystems54.Marine Group II (MGII) and Poseidonia representatives of Thermoplasmatota that have been reported to be nitrate-reducing Archaea59, were recovered from AW sample (AW_40, AW_45) and contained several enzymes contributing in alkane (alkane monooxygenase, aldehyde dehydrogenase) and naphthalene/phenanthrene/phenol/xylene degradation (decarboxylase) under aerobic conditions. The HC degradation potential of representatives of this phylum has been previously reported60.In the Asalouyeh water sample, MAGs AW_25 (UBA4421) and AW_38 (UBA8337) had cddA, AW_21 (UBA8444) had catA, AW_11 (Poseidonia) and AW_17 (from Rhizobiales order) had both alkB/M and xylM, and AW_4 (UBA8337) had catA and pheA genes and had potential to trigger the breakdown of their corresponding oil derivatives.Other recovered genomes had the potential to metabolize the product of initiating enzymes. For instance, AW_23 contained enzymes involved in the degradation of naphthalene, phenol and cyclododecane and was affiliated to the genus Alteromonas (Fig. 6).Three recovered MAGs of HW affiliated to Pseudomonadales (HW_23), Poseidoniales (HW_24), and Flavobacteriales (HW_30) contained some initiating enzymes to degrade cyclododecane/biphenyl/toluene, alkane/xylene, and alkane/xylene/naphthalene/phenanthrene, respectively. A representative of Heimdallarchaeia that are mainly recovered from sediment samples was reconstructed from the Hormuz water sample (HW_28). It had a completeness of 81% and contained enzymes involved in anaerobic degradation of alkanes. This archaeon could potentially be an input from the neighboring OMZ as this phyla include representatives adopted to microoxic niches61. Containing genes with the potential to initiate the oil derivative degradation in the input water with no oil exposure reiterates the intrinsic ability of marine microbiota for HC degradations and oil bioremediation.While 16S rRNA provides an overview of the community, MAGs provide the possibility to inspect the metabolic capability of the microbiota. We decided to provide both in this manuscript as we believe they are complementary. Having the full picture provided by the combination of these analyses allows for a better understanding of the community structure and their metabolic capabilities. This is even more evident for sediment samples as they are highly diverse, and reconstructing MAGs from sediment metagenomes is still a bottlenecks. In this case, we rely more on the 16S rRNA to provide an overall view of the community composition.This said, we see similar taxonomic distribution in the MAGs and 16S rRNA e.g., the prevalence of Flavobacteriales and Rhodobacterales in KhW and KhS, Synechococcales, and Desulfobacteriales and Flavobacteriales in HW, HS and AW samples, respectively.Additionally, some rare microbiota representatives were recovered among reconstructed MAGs. For example, the Immundisolibacterales showed an abundance of only 0.8% in the KhS sample based on 16S rRNA but the recovered KhS_63 MAG was affiliated to this taxon. Notably, this MAG contained many genes involved in hydrocarbon degradation having the highest potential in hydrocarbon degradation. More

Model-guided procedure guides the exploration of butyrate production landscapesWe aimed to explore the butyrate production landscape as a function of community composition to decipher microbial interactions shaping butyrate production. Exploring the butyrate production functional landscape is a major challenge because the number of sub-communities increases exponentially with the number of species43. To investigate the landscape, we developed a modeling framework to guide the iterative design of informative experiments (Fig. 1a, b). Microbial interactions can impact growth or metabolite production by influencing the availability of ecological niches or facilitating metabolite degradation. To capture these two types of interactions, we implemented a two-stage modeling framework to determine the contributions of microbial interactions to species growth and community assembly or metabolite production. In the first stage, a dynamic ecological model, referred to as the generalized Lotka–Volterra model (gLV), predicts community assembly. The second stage predicts metabolite production as a function of the resulting community composition (Fig. 1b). The gLV model is a set of coupled ordinary differential equations that capture the temporal change in species abundances due to monospecies growth parameters and inter-species growth interactions (see the “Methods” section)16. To estimate parameters for the gLV model, we use Bayesian parameter inference techniques to determine the uncertainty in our parameters based on biological and technical variability in the experimental data44.Fig. 1: Iterative modeling framework to predict microbial community assembly and metabolic function.a Two-stage modeling framework for predicting community assembly and function. The generalized Lotka–Volterra model (gLV) represents community dynamics. A Bayesian Inference approach was used to determine parameter uncertainties due to biological and technical variability. A linear regression model with interactions maps assembled community composition to metabolite concentration. Combining these two models enables prediction of a probability distribution of metabolite concentration from initial species abundances. b Design–Test–Learn cycle for model development. First, we use our model to explore the design space of possible experiments (i.e. different initial conditions of species presence/absence) and design communities that span a desired range of metabolite concentrations. Next, we use high-throughput experiments to measure species abundance and metabolite concentration. Finally, we evaluate the model’s predictive capability and infer an updated set of parameters based on the new experimental measurements. c Phylogenetic tree of the synthetic human gut microbiome composed of 25 highly prevalent and diverse species. Branch color indicates phylum and underlined species denote butyrate producers.Full size imageOur metabolite production model consists of a linear regression model with interaction terms mapping community composition (i.e. abundance of each species) at a specific time point to the concentration of an output metabolite at that time. This model was based on a phenomenological model of metabolite production used in bioprocess engineering expanded to microbial communities (see the “Methods” section). In the regression model, the first-order terms capture the monospecies production per unit biomass and the interaction terms represent the impact of inter-species interactions on metabolite production per unit biomass (i.e. deviations from constant metabolite production per unit biomass19). To estimate parameters for the regression model, we use Lasso regression to identify the most impactful interactions. Altogether, the composite gLV and regression model predicts the probability distribution of the metabolite concentration given an initial condition of species abundances (Fig. 1b, see the “Methods” section).In metabolic and protein engineering, a design–test–learn cycle (DTL) has been used to design biomolecules45 or metabolic pathways46 with properties that satisfy desired performance specifications. We hypothesized that this engineering-inspired approach could be used to explore community design spaces and understand the composition–function mapping for butyrate production. Each cycle consisted of: (1) a design phase wherein we used our model informed by experimental observations to simulate a vast number of potential community compositions to identify sub-communities that satisfied biological objectives (i.e. desired butyrate concentrations), (2) a test phase wherein the selected sub-communities were assembled and species abundance and butyrate concentration were measured, and (3) a learn phase wherein patterns in our experimental data were used to estimate model parameters and to extract information about the key microbial interactions influencing community assembly and butyrate production.Two-stage model enables efficient exploration of low richness community design spaceTo develop a system of microbes representing major metabolic functions in the gut, we selected 25 prevalent bacterial species from all major phyla in the human gut microbiome47 (Fig. 1c, Supplementary Data 1). This community contained five butyrate-producing Firmicutes which have been shown to play important roles in human health and protection from diseases (Fig. 1c, Supplementary Data 1). Due to the lack of a defined medium that universally supports the growth of gut microbes, most in vitro studies use undefined media, making it difficult to interrogate the effects of unknown components on community behaviors48. To maximize our knowledge of the substrates available to the communities, we developed a chemically defined medium to grow the synthetic communities (see the “Methods” section).Based on previous studies using pairwise communities to predict higher richness community behaviors16,18,49, we hypothesized that training our model on single and pairwise community measurements would provide an informative starting point for mapping composition–function relationships determining butyrate production. To do so, we first measured time-resolved growth of single species and observed a wide variety of growth dynamics within each phylum, including disparate growth rates and carrying capacities (Supplementary Fig. 1). We assembled each pairwise community containing at least one butyrate producer (the focal species of our system50) and measured species abundance and the concentrations of organic acid fermentation products (including butyrate, lactate, succinate, and acetate) after 48 h. The pairwise consortia displayed a broad range of butyrate concentrations of 0–50 mM (Fig. 2a).Fig. 2: Exploring the predicted butyrate production of 3–5 member communities with a model trained on 1–2 species communities.a Categorical scatter plot of butyrate production in 1–2 species and 24–25 species communities. Solid datapoints indicate the mean of the biological replicates which are represented by transparent datapoints connected to the mean with transparent lines. The colors indicate which butyrate producer was present in the community with green indicating the presence of multiple butyrate producers. DP− and AC− indicate the 24-species communities lacking Desulfovibrio piger (DP) and Anaerostipes caccae (AC), respectively. b Predicted medians (black line) and 60 percent confidence intervals (gray bars) of butyrate concentration for all 3–5 member communities containing at least one butyrate producer (46,591 community predictions). Colored lines indicate median and 60 percent confidence interval of butyrate production of communities chosen for the experimental design with the color indicating the number of species in the community (156 communities). Subplots separate groups of communities based on the identities of the combination of butyrate producers specified. c Scatter plot of measured butyrate versus predicted butyrate for 3–5 species communities. Colors indicate which butyrate producer was present in the community as in a. Biological replicates (n = 1–5, depending on the community, exact values in source data) are indicated by transparent squares connected to the corresponding mean, which is represented by the large data point. Prediction error bars (x-axis) indicate the 60% confidence interval of the predicted butyrate distribution as in b, with the center being the median prediction. Dashed line indicates the linear regression between the mean measured butyrate and the median predicted butyrate. Indicated statistics are for Pearson correlation (two-sided). Source data are available in the Source Data file.Full size imageSingle-species deletion communities have been used to investigate the contributions of individual species to a community function13,16. Therefore, we characterized the full 25-species community and each single-species deletion sub-community (i.e. 24-member consortia). In stark contrast to the pairwise communities, the 24- and 25-species communities exhibited similar low butyrate production (~2–22 mM Butyrate). The absence of only two species Desulfovibrio piger (DP) (~22 mM Butyrate) and Anaerostipes caccae (AC) (~2 mM Butyrate) resulted in a significant increase or decrease in butyrate concentration compared to the remaining 24-member and 25-member communities (Fig. 2a, Supplementary Fig. 2a). In addition, the concentrations of all measured organic acids spanned a much smaller range in the 24 and 25-member communities than the single and pairwise consortia (Supplementary Fig. 2b). These results suggest that high richness communities may trend towards a similar low butyrate-producing state that is difficult to change by the deletion of most single species and motivates a model-guided design strategy for exploring how community richness shapes butyrate production.To determine whether individual and pairwise communities could predict community composition and butyrate production of low richness communities (i.e. 3–5 species), we estimated the parameters of our model based on experimental measurements. Our initial model was informed only by pairwise communities that contained at least one butyrate producer (Supplementary Data 2, M1) and was thus naïve to all interactions between non-butyrate producers. We assumed that the unobserved growth interactions could be predicted based on trends in measured interactions across phylogenetic relatedness (see the “Methods” section)16. However, the resulting model was unable to predict butyrate production in the 24-and 25-member communities (Supplementary Fig. 3), which we attributed to missing information about non-butyrate producer interactions in our training data. Thus, we used our model to explore a low richness design space of 3–5 species communities based on the assumption that pairwise interactions would be more observable in low than high richness (i.e. >10 species) communities to identify an improved parameter estimate for non-butyrate producer interactions.We used our initial M1 model to predict the probability distributions of butyrate production for all 3–5 species communities containing at least one butyrate producer (46,591 communities). The predicted butyrate production varied substantially based on the combination of butyrate producers present in each community (Fig. 2b). In addition, we observed variations in the shapes of the probability distributions based on how the uncertainty in growth prediction propagated through the regression model. For instance, the butyrate concentration in the AC, Roseburia intestinalis (RI) pairwise community was lower than the AC monoculture, even though RI was low abundance, resulting in a high magnitude negative parameter in the regression model for a production interaction between AC and RI (Supplementary Data 3). Due to the uncertainty in the growth parameters, the model predicted that RI would grow substantially in a subset of the 3–5 member simulations containing both AC and RI. The variability in predicted RI growth combined with the high magnitude negative interaction parameter between AC and RI resulted in distributions where the median butyrate concentration was high (i.e. for simulations where RI did not grow substantially), and the 60 percent confidence interval extended to 0 mM butyrate (i.e. when RI grew substantially) (Fig. 2b). In sum, these results demonstrate that the shape of the predicted probability distributions can provide information about the uncertainty in species growth based on experimental observations.Based on the simulations, we selected 156 communities that spanned a broad range of predicted butyrate concentrations across the butyrate producer groups to evaluate experimentally (Fig. 2b). The model prediction exhibited good agreement with the rank order of butyrate production (Spearman rho = 0.84, p = 9*10−43) (Fig. 2c) and species abundance (Spearman rho = 0.76, p = 3*10−122) (Supplementary Fig. 4a–d), demonstrating that our initial model could predict a wide range of butyrate production in low richness communities.Composition–function landscape predicts contributions of growth and production interactionsEncouraged by our model’s predictive ability, we sought to explore composition–function relationships in higher richness communities (i.e. >10 species) using a model with updated parameters based on measurements of the 3–5 member communities (Supplementary Data 2, M2). Since the human gut microbiome exhibits functional redundancy in butyrate pathways51, we first used model M2 to simulate the assembly of all communities containing all five butyrate producers (5-butyrate producer or 5BP, 1,048,576 total) to map the composition–function landscape for butyrate production (Fig. 3a). In addition, we simulated the assembly of all communities containing the four butyrate producers excluding AC (4-butyrate producer or 4BP, 1,048,576 total) to understand how the composition–function landscape changes in the absence of the most productive butyrate producer (Fig. 3b). The majority of 5BP communities were predicted to have higher butyrate concentration than any of the 4BP communities (Fig. 3a, b), consistent with the substantial decrease in butyrate in the AC deletion community observed previously (Fig. 2a).Fig. 3: Community composition–function landscapes reveal key role of production interactions on A. caccae and negative impact of D. piger on butyrate production.a Scatter plot of predicted total butyrate producer abundance versus predicted butyrate concentration for all possible communities in which all five butyrate producers are present (1,048,576 communities). Histograms indicate the butyrate concentration distribution across the given axis. Communities are colored according to the presence (red) or absence (blue) of D. piger (DP). Blue and red dashed lines indicate the linear regression of communities with (red, y = −2.3x + 25.8, r = −0.34) or without (blue, y = 3.1x + 28.0, r = 0.76) DP. The white star indicates the full 25-member community and black star indicates the community of five butyrate producers alone. Large data points indicate communities chosen for experimental validation. Black triangles indicate leave-one-out communities, black circles indicate designed communities, and gray squares indicate random communities, with open/closed symbols indicate the absence/presence of DP. b Scatter plot of predicted total butyrate producer abundance versus predicted butyrate concentration for all possible communities in which all four butyrate producers excluding AC are present (1,048,576 communities). Histograms indicate the butyrate concentration distribution. Gray dashed line indicates the mean predicted butyrate concentration across all communities. Black dashed line indicates the linear regression of all communities (y = 13.2x−0.1, r = 0.64). The white star indicates the full 24-member community and the black star indicates the four butyrate producers alone. Large data points indicate communities chosen for experimental validation. c, d Scatter plots of experimental measurements of total butyrate producer abundance versus butyrate concentration for communities with (c) and without (d) AC. Data point shapes correspond to the legends in (a) and (b) and represent the mean of biological replicates, which are shown as small datapoints connected to the corresponding mean with lines (n = 2 except for 5 BP, 24 and 25 species communities where n = 5–8). Dashed line in (d) indicates the linear regression (y = 8.9x−0.5). e Comparison of butyrate concentration in communities from a and b with and without DP for both the designed and random experimental sets for the 5BP communities and designed set for the 4BP communities. Data point shapes correspond to the legends in a and b and represent the mean of biological replicates, which are shown as small datapoints connected to the corresponding mean with lines. Box and whisker plots represent the median (center line), quartiles (box), and range (whiskers) of the mean butyrate concentration for each community, excluding outliers (points outside 1.5 times the interquartile range). Indicated p-values are from a Mann–Whitney U test (5BP Designed: n = 28 for DP+ and n = 54 for DP−; 5BP Random: n = 27 for DP+ and n = 55 for DP−; 4BP: n = 42 for DP+ and n = 42 for DP−). f Butyrate concentration per unit biomass as a function of sulfide concentration after 24 h of growth. Butyrate concentration per biomass was normalized to the no sulfide condition. Circles indicate the mean of biological replicates, with individual replicates shown as transparent squares (n = 4). Inset: Butyrate concentration per biomass (mM OD600−1) for AC with and without the addition of 1.6 mM sulfide across time (n = 3). Endpoint sulfide concentrations were higher in the data shown in the inset than in the main figure (Supplementary Fig. 6). Source data are available in the Source Data file.Full size imageThe relationship between butyrate producer abundance and butyrate can provide insight into the contributions of growth and production interactions in the presence and absence of AC (Fig. 3a, b). If butyrate producer abundance correlates with butyrate, then growth interactions drive butyrate production, whereas the contributions of production interactions would reduce the strength of this correlation. The 5BP communities were predicted to have a large contribution of production interactions as evidenced by a weak correlation between butyrate concentration and butyrate producer abundance (Spearman rho = 0.17, p  More

1.Sala, O. E. et al. Global biodiversity scenarios for the year 2100. Science 287, 1770–1774 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
2.Dudgeon, D. et al. Freshwater biodiversity: Importance, threats, status and conservation challenges. Biol. Rev. 81, 163–182 (2006).PubMed 
Article 

Google Scholar 
3.Strayer, D. L. & Dudgeon, D. Freshwater biodiversity conservation: Recent progress and future challenges. J. N. Am. Benthol. Soc. 29, 344–358 (2010).Article 

Google Scholar 
4.Invasive Species Specialist Group IUCN guidelines for the prevention of biodiversity loss caused by alien invasive species. https://portals.iucn.org/library/node/12673 (2000).5.Clavero, M. & García-Berthou, E. Invasive species are a leading cause of animal extinction. Trends Ecol. Evol. 20, 110 (2005).PubMed 
Article 

Google Scholar 
6.Hassan, R., Scholes, R. J. & Ash, N. Ecosystems and human well-being: Current state and trends: Findings of the Condition and Trends working group (Millennium Ecosystem Assessment Series) (Island Press, 2005).7.Vitousek, P. M., D’Antonio, C. M., Loope, L. L., Rejmánek, M. & Westbrooks, R. Introduced species: A significant component of human-caused global change. N. Z. J. Ecol. 21, 1–16 (1997).
Google Scholar 
8.Mack, R. N. et al. Biotic invasions: Causes, epidemiology, global consequences, and control. Ecol. Appl. 10, 689–710 (2000).Article 

Google Scholar 
9.Hulme, P. E. Beyond control: Wider implications for the management of biological invasions. J. Appl. Ecol. 43, 835–847 (2006).Article 

Google Scholar 
10.Vander Zanden, M. J., Hansen, G. J. A., Higgins, S. N. & Kornis, M. S. A pound of prevention, plus a pound of cure: Early detection and eradication of invasive species in the Laurentian Great Lakes. J. Great Lakes Res. 36, 199–205 (2010).Article 

Google Scholar 
11.Myers, J. H., Simberloff, D., Kuris, A. M. & Carey, J. R. Eradication revisited: Dealing with exotic species. Trends Ecol. Evol. 15, 316–320 (2000).CAS 
PubMed 
Article 

Google Scholar 
12.Mehta, S. V., Haight, R. G., Homans, F. R., Polasky, S. & Venette, R. C. Optimal detection and control strategies for invasive species management. Ecol. Econ. 61, 237–245 (2007).Article 

Google Scholar 
13.McDonald, L. L. Sampling rare populations. In Sampling Rare or Elusive Species (ed. Thompson, W. L.) 11–42 (Island Press, 2004).
Google Scholar 
14.Harvey, C. T., Qureshi, S. A. & MacIsaac, H. J. Detection of a colonizing, aquatic, non-indigenous species. Divers. Distrib. 15, 429–437 (2009).Article 

Google Scholar 
15.Ficetola, G. F., Miaud, C., Pompanon, F. & Taberlet, P. Species detection using environmental DNA from water samples. Biol. Lett. 4, 423–425 (2008).PubMed 
PubMed Central 
Article 

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

Google Scholar 
17.Valentini, A., Pompanon, F. & Taberlet, P. DNA barcoding for ecologists. Trends Ecol. Evol. 24, 110–117 (2009).PubMed 
Article 

Google Scholar 
18.Thomsen, P. F. & Willerslev, E. Environmental DNA—An emerging tool in conservation for monitoring past and present biodiversity. Biol. Conserv. 183, 4–18 (2015).Article 

Google Scholar 
19.Rees, H. C., Maddison, B. C., Middleditch, D. J., Patmore, J. R. M. & Gough, K. C. The detection of aquatic animal species using environmental DNA—A review of eDNA as a survey tool in ecology. J. Appl. Ecol. 51, 1450–1459 (2014).CAS 
Article 

Google Scholar 
20.Brys, R. et al. Monitoring of spatio-temporal occupancy patterns of fish and amphibian species in a lentic aquatic system using environmental DNA. Mol. Ecol. https://doi.org/10.1111/mec.15742 (2021).Article 

Google Scholar 
21.Smart, A. S. et al. Assessing the cost-efficiency of environmental DNA sampling. Methods Ecol. Evol. 7, 1291–1298 (2016).Article 

Google Scholar 
22.Wilcox, T. M. et al. Understanding environmental DNA detection probabilities: A case study using a stream-dwelling char Salvelinus fontinalis. Biol. Conserv. 194, 209–216 (2016).Article 

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

Google Scholar 
24.Bohmann, K. et al. Environmental DNA for wildlife biology and biodiversity monitoring. Trends Ecol. Evol. 29, 358–367 (2014).PubMed 
Article 

Google Scholar 
25.Furlan, E. M., Gleeson, D., Hardy, C. M. & Duncan, R. P. A framework for estimating the sensitivity of eDNA surveys. Mol. Ecol. Resour. 16, 641–654 (2016).CAS 
PubMed 
Article 

Google Scholar 
26.Cristescu, M. E. & Hebert, P. D. N. Uses and misuses of environmental DNA in biodiversity science and conservation. Annu. Rev. Ecol. Evol. Syst. 49, 209–230 (2018).Article 

Google Scholar 
27.Sepulveda, A. J., Nelson, N. M., Jerde, C. L. & Luikart, G. Are environmental DNA methods ready for aquatic invasive species management?. Trends Ecol. Evol. 35, 668–678 (2020).PubMed 
Article 

Google Scholar 
28.Wilcox, T. M. et al. Robust detection of rare species using environmental DNA: The importance of primer specificity. PLoS One 8, e59520. https://doi.org/10.1371/journal.pone.0059520 (2013).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
29.Freeland, J. The importance of molecular markers and primer design when characterizing biodiversity from environmental DNA (eDNA). Genome 60, 358–374 (2016).PubMed 
Article 
CAS 

Google Scholar 
30.Goldberg, C. S. et al. Critical considerations for the application of environmental DNA methods to detect aquatic species. Methods Ecol. Evol. 7, 1299–1307 (2016).Article 

Google Scholar 
31.Veldhoen, N. et al. Implementation of novel design features for qPCR-based eDNA assessment. PLoS One 11, e0164907. https://doi.org/10.1371/journal.pone.0164907 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
32.Lin, M., Zhang, S. & Yao, M. Effective detection of environmental DNA from the invasive American bullfrog. Biol. Invasions 21, 2255–2268 (2019).Article 

Google Scholar 
33.Thalinger, B. et al. A validation scale to determine the readiness of environmental DNA assays for routine species monitoring. Environ. DNA https://doi.org/10.1002/edn3.189 (2021).Article 

Google Scholar 
34.Yates, M. C., Fraser, D. J. & Derry, A. M. Meta-analysis supports further refinement of eDNA for monitoring aquatic species-specific abundance in nature. Environ. DNA 1, 5–13 (2019).Article 

Google Scholar 
35.Hindson, B. J. et al. High-throughput droplet digital PCR system for absolute quantitation of DNA copy number. Anal. Chem. 83, 8604–8610 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Nathan, L. M., Simmons, M., Wegleitner, B. J., Jerde, C. L. & Mahon, A. R. Quantifying environmental DNA signals for aquatic invasive species across multiple detection platforms. Environ. Sci. Technol. 48, 12800–12806 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
37.Doi, H. et al. Use of droplet digital PCR for estimation of fish abundance and biomass in environmental DNA surveys. PLoS One 10, e0122763. https://doi.org/10.1371/journal.pone.0122763 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
38.Brys, R. et al. Reliable eDNA detection and quantification of the European weather loach (Misgurnus fossilis). J. Fish Biol. https://doi.org/10.1111/jfb.14315 (2020).Article 
PubMed 

Google Scholar 
39.Lacoursière-Roussel, A., Côté, G., Leclerc, V. & Bernatchez, L. Quantifying relative fish abundance with eDNA: A promising tool for fisheries management. J. Appl. Ecol. 53, 1148–1157 (2016).Article 
CAS 

Google Scholar 
40.Doi, H. et al. Environmental DNA analysis for estimating the abundance and biomass of stream fish. Freshw. Biol. 62, 30–39 (2017).CAS 
Article 

Google Scholar 
41.Buxton, A. S., Groombridge, J. J., Zakaria, N. B. & Griffiths, R. A. Seasonal variation in environmental DNA in relation to population size and environmental factors. Sci. Rep. 7, 46294. https://doi.org/10.1038/srep46294 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
42.Takahara, T., Iwai, N., Yasumiba, K. & Takeshi, I. Comparison of the detection of 3 endangered frog species by eDNA and acoustic surveys across 3 seasons. Freshw. Sci. 39, 18–27 (2020).Article 

Google Scholar 
43.Kats, L. B. & Ferrer, R. P. Alien predators and amphibian declines: Review of two decades of science and the transition to conservation. Divers. Distrib. 9, 99–110 (2003).Article 

Google Scholar 
44.Martel, A. et al. The novel ‘Candidatus Amphibiichlamydia ranarum’ is highly prevalent in invasive exotic bullfrogs (Lithobates catesbeianus). Environ. Microbiol. Rep. 5, 105–108 (2012).PubMed 
Article 

Google Scholar 
45.Blaustein, A. R. et al. Effects of invasive larval bullfrogs (Rana catesbeiana) on disease transmission, growth and survival in the larvae of native amphibians. Biol. Invasions 22, 1771–1784 (2020).Article 

Google Scholar 
46.Lowe, S., Browne, M., Boudjelas, S. & De Poorter, M. 100 of the world’s worst invasive alien species. A selection from the Global Invasive Species Database. Published by The Invasive Species Specialist Group (ISSG) a specialist group of the Species Survival Commission (SSC) of the World Conservation Union (IUCN), First published as special lift-out in Aliens 12 (2000).47.Adams, M. J. & Pearl, C. A. Problems and opportunities managing invasive bullfrogs: Is there any hope? In Biological Invaders in Waters: Profiles, Distribution and Threats (ed. Gherardi, F.) 679–693 (Springer, Paris, 2007).
Google Scholar 
48.Louette, G., Devisscher, S. & Adriaens, T. Combating adult invasive American bullfrog Lithobates catesbeianus. Eur. J. Wildl. Res. 60, 703–706 (2014).Article 

Google Scholar 
49.Kamoroff, C. et al. Effective removal of the American bullfrog (Lithobates catesbeianus) on a landscape level: Long term monitoring and removal efforts in Yosemite Valley, Yosemite National Park. Biol Invasions 22, 617–626 (2020).Article 

Google Scholar 
50.Jooris, R. Palmt de stierkikker uit Noord-Amerika ook Vlaanderen in?. Natuur. Focus 1, 13–15 (2001).
Google Scholar 
51.Adriaens, T., Devisscher, S. & Louette, G. Risk analysis of American bullfrog, Lithobates catesbeianus. Risk analysis report of non-native organisms in Belgium. Rapporten van het Instituut voor Natuur- en Bosonderzoek 41. https://doi.org/10.13140/2.1.2431.5688 (2013).52.Descamps, S. & De Vocht, A. Movements and habitat use of the invasive species Lithobates catesbeianus in the valley of the Grote Nete (Belgium). Belg. J. Zool. 146, 90–100 (2016).
Google Scholar 
53.Strickler, K. M., Fremier, A. K. & Goldberg, C. S. Quantifying effects of UV-B, temperature, and pH on eDNA degradation in aquatic microcosms. Biol. Conserv. 183, 85–92 (2015).Article 

Google Scholar 
54.Lefever, S., Pattyn, F., Hellemans, J. & Vandesompele, J. Single-nucleotide polymorphisms and other mismatches reduce performance of quantitative PCR assays. Clin. Chem. 59, 1470–1480 (2013).CAS 
PubMed 
Article 

Google Scholar 
55.Erligh, H. A., Gelfand, D. & Sninsky, J. J. Recent advances in the polymerase chain reaction. Science 252, 1643–1651 (1991).ADS 
Article 

Google Scholar 
56.Schrader, C., Schielke, A., Ellerbroek, L. & Johne, R. PCR inhibitors—Occurrence, properties and removal. J. Appl. Microbiol. 113, 1014–1026 (2012).CAS 
PubMed 
Article 

Google Scholar 
57.Lievens, A., Jacchia, S., Kagkli, D., Savini, C. & Querci, M. Measuring digital PCR quality: Performance parameters and their optimization. PLoS One 11, e0153317. https://doi.org/10.1371/journal.pone.0153317 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
58.Pecoraro, S. et al. Overview and recommendations for the application of digital PCR. EUR 29673 EN, Publications Office of the European Union. https://doi.org/10.2760/192883 (2019).59.Harper, L. R. et al. Prospects and challenges of environmental DNA (eDNA) monitoring in freshwater ponds. Hydrobiologia 826, 25–41 (2019).CAS 
Article 

Google Scholar 
60.Doi, H. et al. Droplet digital PCR outperforms real-time PCR in the detection of environmental DNA from an invasive fish species. Environ. Sci. Technol. 49, 5601–5608 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar 
61.Wells, K. D. (ed.) The Ecology and Behavior of Amphibians (The University of Chicago Press, 2007).
Google Scholar 
62.Willis, Y. L., Moyle, D. I. & Baskett, T. S. Emergence, breeding, hibernation, movements and transformation of the bullfrog, Rana catesbeiana Missouri. Copeia 1, 30–41 (1956).Article 

Google Scholar 
63.Maruyama, A., Nakamura, K., Yamanaka, H., Kondoh, M. & Minamoto, T. The release rate of environmental DNA from juvenile and adult fish. PLoS One 9, e114639. https://doi.org/10.1371/journal.pone.0114639 (2014).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
64.Barnes, M. A. et al. Environmental conditions influence eDNA persistence in aquatic systems. Environ. Sci. Technol. 48, 1819–1827 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
65.Lance, R. F. et al. Experimental observations on the decay of environmental DNA from bighead and silver carps. Manag. Biol. Invasions 8, 343–359 (2017).Article 

Google Scholar 
66.Hoorfar, J. Practical considerations in design of internal amplification controls for diagnostic PCR assays. J. Clin. Microbiol. 42, 1863–1868 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
67.Devisscher, S. et al. Beheer van de stierkikker in Vlaanderen en Nederland. Rapporten van het Instituut voor Natuur- en Bosonderzoek 52. https://www.researchgate.net/publication/235789235 (2012).68.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).Article 

Google Scholar 
69.Dunn, N., Priestley, V., Herraiz, A., Arnold, R. & Savolainen, V. Behavior and season affect crayfish detection and density inference using environmental DNA. Ecol. Evol. 7, 7777–7785 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
70.Bury, R. B. & Whelan, J. A. Ecology and management of the bullfrog. U.S. Fish and Wildlife Service 155 (1984).71.Gahl, M. K., Calhoun, A. J. K. & Graves, R. Facultative use of seasonal pools by American bullfrogs (Rana catesbeiana). Wetlands 29, 697–703 (2009).Article 

Google Scholar 
72.Biek, R., Funk, C., Maxell, B. A. & Mills, L. S. What is missing in amphibian decline research: Insights from ecological sensitivity analysis. Conserv. Biol. 16, 728–734 (2002).Article 

Google Scholar 
73.Govindarajulu, P., Altwegg, R. & Anholt, B. R. Matrix model investigation of invasive species control: Bullfrogs on Vancouver Island. Ecol. Appl. 15, 2161–2170 (2005).Article 

Google Scholar 
74.Carim, K. J. et al. Environmental DNA sampling informs fish eradication efforts: Case studies and lessons learned. N. Am. J. Fish. 40, 488–508 (2020).Article 

Google Scholar 
75.Riaz, T. et al. ecoPrimers: Inference of new DNA barcode markers from whole genome sequence analysis. Nucleic Acids Res. 39, e145. https://doi.org/10.1093/nar/gkr732 (2011).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
76.Moyer, G. R., Díaz-Ferguson, E., Hill, J. E. & Shea, C. Assessing environmental DNA detection in controlled lentic systems. PLoS One 9, e103767. https://doi.org/10.1371/journal.pone.0103767 (2014).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
77.Turner, C. R. et al. Particle size distribution and optimal capture of aqueous macrobial eDNA. Methods Ecol. Evol. 7, 676–684 (2014).Article 

Google Scholar 
78.U.S. Fish and Wildlife Service. Quality assurance project plan: eDNA monitoring of bighead and silver carps. https://www.fws.gov/midwest/fisheries/eDNA/documents/QAPP.pdf (2017).79.Spens, J. et al. Comparison of capture and storage methods for aqueous macrobial eDNA using an optimized extraction protocol: Advantage of enclosed filter. Methods Ecol. Evol. 8, 635–645 (2017).Article 

Google Scholar 
80.RStudio Team (2020) RStudio: Integrated Development Environment for R. RStudio, PBC. http://www.rstudio.com/. More

COVID-19 vs. MoMoThe COVID-19 official deaths and MoMo ED time series overlaped for the period from 3 March 2020 to 1 January 2021 for Spain and its 19 regions (Fig. 2). In general, there was good agreement between both datasets, meaning that most of MoMo ED were related to COVID-19 deaths. During the first wave, the most important differences were observed in Spain, Madrid, Cataluña, Castilla-La Mancha, and Castilla y León. Before 22 June in Spain, MoMo ED showed 15,445 accumulated deaths more than the official COVID-19 deaths, which is beyond the error band. That difference comes basically from the four regions with the largest numbers of deaths (Madrid, Cataluña, Castilla-La Mancha, and Castilla y León). Table 1 shows the accumulated values before 22 June, which were used to estimate the CFR for Spain and its 19 regions according to the third phase of the National Seroprevalence Study4,5. For all regions, the CFR estimated from MoMo ED was larger than the CFR estimated from COVID-19 deaths. In particular, Asturias, Canarias, and Murcia were twice as large. Ceuta and Melilla dramatically increased their CFR from MoMo ED, although that may be biased due to their small populations and numbers of deaths.Similarly, the same variables for the period from 23 June to 29 November 2020 are reported in Table 2. In Spain, MoMo ED showed 6173 accumulated deaths more than the official COVID-19 deaths. This difference is a third of the difference observed prior to 22 June; because this is within the error band, there was a significant improvement in the detection of COVID-19 deaths in this period. Figure 2 also shows a general agreement between MoMo ED and official COVID-19 deaths time series after the first wave, with the exception of late July and early August. These differences were due to two heat waves that were responsible for at least 25% of the MoMo ED16.Infections estimated from COVID-19 deathsTo illustrate the delay between official daily infections data and REMEDID estimated daily infections, we applied REMEDID from COVID-19 deaths assuming CFR = 100%. Figure 3 shows the current IO21 and the infections associated with COVID-19 deaths for the first wave. The latter in Spain reached a maximum on 13 March 2020 (Table 3), the day before the national government decreed a state of emergency and national lockdown. Thus, the adopted measures had an immediate effect, which was observed in the official data IO21 7 days later (20 March). This delay is similar to the incubation period (mean 5.78 days2), which could be explained because official infections were reported when symptoms appeared. This delay reached 16 days when we compared with earlier version IO20 (not shown), which highlights the usefulness of the methodology to reinterpret official data from very early stages of the pandemic. On the other hand, the maximum number of deaths was reached on 1 April, which was 19 days after the inferred infection maximum, bringing this delay close to the 20 days expected between infection and death (Figs. 1, 3).Figure 3Official COVID-19 infections and deaths, and estimated infections with case fatality ratio (CFR) of 100% in Spain during the first wave. Left y-axis: COVID-19 daily infections IO21 (blue curve). Right y-axis: COVID-19 deaths (orange curve) and its REMEDID-estimated infections with CFR = 100% (red curve). All curves are for Spain. Thin blue and orange curves are daily data, and thick curves are smoothed by 14-day running mean. Arrows show delays between the maximum of inferred infections and maxima from COVID-19 deaths (orange arrows) and COVID-19 infections (blue arrows). Solid arrows are expected delays, dotted while arrows are observed delays.Full size imageTable 3 Date of first infection for REMEDID estimated daily infections from COVID-19 deaths (IRO) and from MoMo Excess Deaths (ED) (IRM), and for official COVID-19 daily infections released on June 2020 (IO20) and on February 2021 (IO21).Full size tableWe applied REMEDID to the official COVID-19 deaths with the corresponding estimated CFRs (see “Data” section) to obtain the time series of estimated daily infections, hereafter referred to as IRO. Figure 4 shows IRO and the accumulated infections for Spain and its 19 regions. Note that in Spain, IRO are amplified versions of inferred infections in Fig. 3. In Spain, the first infection, according to IRO,, is on 8 January 2020 (Table 3), 43 days before the first infection was officially reported on 20 February 2020 according to IO20. By contrast, IO21 places the first infection on 1 January 2020. Spain reached the maximum number of IRO on 13 March, a day before the state of emergency and lockdown were enforced (Table 4). On 14 March, IO20 = 1832, and IO21 = 7478; however, IRO = 63,727 (CI 95% 60,050–67,403), 35 and 9 times IO20 and IO21, respectively (Table 5). This implies that on that day, IO20 and IO21 only reported 2.9% (CI 95% 2.7–3.1%) and 11.7% (CI 95% 11.1–12.5%) of new infections, respectively. Although detection of infections clearly improved from IO20 to IO21, almost 90% of the infections are still not documented in the peak of the first wave. The situation is similar for the accumulated infections before 22 June 2020, as reported by the National Seroprevalence Study4,5.Figure 4Daily and accumulated infections for official COVID-19 daily infections (IO21), and daily infections estimated from COVID-19 deaths (IRO). Lines are daily infections and refer to the y-axis on the right; bars are accumulated infections and refer to the y-axis on the left. Red lines and cyan bars are official COVID-19 data; orange lines and blue bars are inferred infections with case fatality ratio (CFR) in Table 1. Thin orange lines correspond to the CFR confidence interval.Full size imageTable 4 Date of the most prominent relative maxima, for Spain and the 19 regions, of the REMEDID estimated daily infections from COVID-19 deaths (IRO) and from MoMo excess deaths (ED) (IRM), and official COVID-19 daily infections (IO21).Full size tableTable 5 REMEDID estimated daily infections from COVID-19 deaths (IRO) and from MoMo excess deaths (ED) (IRM) on 14 March, and for official COVID-19 daily infections released on June 2020 (IO20) and released on February 2021 (IO21).Full size tableIn almost all regions, IO20 showed a delay of 1 month or more between the first infection and IRO (Table 3). No delay in IO21 occurred in Islas Baleares, Castilla-León, and Galicia, while in three regions (Cataluña, Madrid, and La Rioja), the first case occurred earlier than the first case of IRO. However, 6 regions had delays of 15 days and other 6 regions had delays of 1 month. According to IRO, all regions except Ceuta and Melilla had some infections in January, but in IO21 only 6 regions had infections in that month. In all scenarios, the first infections were in Madrid and Cataluña.During the first wave, according to IRO most of the regions had maximum daily infections around 14 March. In Madrid, the maximum was reached on 11 March, coinciding with the educational centres closing and an official warning by the regional government (Table 4). Asturias was the last region to reach peak infections (25–26 March). The maximum percentage of documented cases (12.6%, CI 95% 9.2–18.4%) occurred in Asturias on 14 March, but in the other regions, only between 1.2 and 8% of the infections were documented (Table 5).Figure 4 shows how the IO21 and IRO curves of Spain and the 19 different regions fluctuated following the same pattern until the middle of June 2020, but thereafter, they showed different patterns. This reflects the fact that the Spanish government had decreed the control measures for the whole nation until June, but thereafter, each regional government implemented its own control measures. For example, some regions (e.g., Aragón, Islas Baleares, Cantabria, Comunidad Valenciana, Extremadura, Galicia, Murcia, País Vasco, and La Rioja) had two peaks, but others had only one. An apparent maximum on 22 June in Islas Baleares is an artifact produced by the interpolation for transition from the two CFRs. Although beyond the scope of this work, it would be very interesting to investigate the effects of the different control measures implemented on the corresponding IRO for the 19 regions.The Spanish COVID-19s wave reached a maximum of daily infections on 22 October from IRO and on 26 October from IO21. The delay of 4 days is similar to the mean incubation period (5.78 days2). The estimated number of new infections is still larger than the documented cases, but the shapes of the two curves are more similar in the second wave than in the first wave (Fig. S1). The same is true for the 19 regions, most of which had the largest peak around 22–26 October, with the exceptions of Canarias and Madrid, which reached maxima in late August and early September, respectively.Infections from MoMo excess deathsAssuming that MoMo ED accounts for both recorded and non-recorded COVID-19 deaths, negative deaths are meaningless, and they were set to zero. Then, the associated daily infections can be estimated, as in “Infections estimated from COVID-19 deaths” section, with a CFR of 100% from MoMo ED for Spain (Fig. 5). Note two main differences between this time series and that estimated from official COVID-19 deaths: (1) MoMo data present an error band that was inherited by the estimated infections; (2) MoMo ED estimated infections reached a maximum of 1443 (CI 99% 1329–1547), doubling the 776 inferred daily infections from official COVID-19 deaths in Fig. 3. This is because maximum MoMo ED was 1,584 (CI 99% 1468–1686) and maximum COVID-19 official deaths was 828, both estimated from the 14-day running mean time series. The maximum of inferred infections was reached on 13 March, just one day prior to the state of emergency and lockdown. The expected and observed delays with respect to official infections and MoMo ED were similar to those observed for estimated infections from official COVID-19 deaths. Error bounds of the estimated infections in Fig. 5 were computed from the MoMo ED error bounds. However, it should be highlighted that the combination of the error bounds from MoMo ED and the estimated CFRs might lead to unrealistic error estimates. To avoid this, the error estimates in Fig. 6 were estimated from the MoMo ED time series (no error bounds) and the error bounds of the estimated CFRs.Figure 5Official COVID-19 infections, MoMo Excess Deaths (ED), and estimated infections with case fatality ratio (CFR) of 100% in Spain during the first wave. Left y-axis: COVID-19 daily infections IO21 (blue curve). Right y-axis: MoMo ED (orange curve) and its REMEDID-estimated infections with CFR = 100% (red curve). All curves are for Spain. Thin blue and orange curves are daily data, and thick curves are smoothed by 14-days running mean. Dashed curves represent the error estimate of MoMo ED (orange) and inferred infections (red). Arrows show delays between the maximum of inferred infections and maxima from MoMo ED (orange arrows) and COVID-19 infections (blue arrows). Solid arrows are expected delays, dotted while arrows are observed delays.Full size imageFigure 6Daily and accumulated infections for official COVID-19 daily infections (IO21), and daily infections estimated from MoMo Excess Deaths (ED) (IRM). Lines are daily infections and refer to the y-axis on the right; and bars are accumulated infections and refer to the y-axis on the left. Red lines and cyan bars are for official COVID-19 data; and orange lines and blue bars are for inferred infections with case fatality ratio (CFR) in Table 1. Thin orange lines represent the error estimate of inferred infections.Full size imageThe REMEDID was applied to the MoMo ED with the corresponding CFRs (see “Data” section) to obtain the estimated daily infections, which will be referred hereafter as IRM. The IRM were calculated for Spain and its 19 regions and are depicted in Fig. 6, as well as the accumulated IRM. In Spain, the first infection shown by IRM happened on 9 January, with an error estimate from 9 to 10 January, 41 to 42 days before the first documented infection of IO20 on 20 February 2020 (Table 3). The maximum IRM was 77,855 (CI 95% 73,364–82,347) reached on 13 March. On 14 March, IRM showed 14,128 infections more than IRO (Table 5). Notice that the CFR used with MoMo ED data was larger than the one used with official COVID-19 deaths data, which makes this difference even more remarkable, because the larger the CFR the lower the estimated infections. Therefore, if the true CFRs, which are unknown, were used in both cases, IRM would double IRO on 14 March, as happened when a CFR of 100% was used (Figs. 3, 5). Notice that with the CFRs used, the IRM and IRO resulted in the same accumulated infections on 22 June and 29 November, matching the results of the seroprevalence study. Nevertheless, IRM showed 42 times more cases than IO20 and 10 times more than IO21 on 14 March, detection of official cases of only 2.4% (2.2–2.5%) and 9.6% (9.1–10.2%), respectively.Table 3 shows the estimated date of first infection for Spain and by region. Note that the first cases of IRM in Spain were on 9 January and in Aragón, Canarias, and Navarra on 8 January, which is possible because significant excess deaths in a region may not become significant for the whole country. In general, the maxima of daily infections were closer to those on 14 March in IRM than in IRO. During the first wave, all regions showed a single maximum, except for Ceuta, Melilla, and Murcia, which showed two maxima (Fig. 6). In general, the IRM time series in all regions were similar during that period. The official data clearly under-detected infections during the first wave. On 14 March, IRM were comparable to IRO, overlapping CI in all regions, but not in Spain as a whole (Table 5). During the second wave, there was improved detection of cases with differences among regions. More

Study site and speciesThe sub-tropical forest analysed in this study is located in China, Yunnan, Xishuangbanna, Mengla county, Village Quingyanzhai (#94) N 21.802068, E 101.380214, geodetic datum WGS84 (Fig. 6). The site is characterized by a rocky outcrop rising 30–50 m over surrounding rubber plantations, harbouring about 20 ha of relict dry tropical forest. The outcrop sides are steep and mainly covered with bamboo. The top area is colonized by shrubs and 10–15 m high trees (a few trees on the slopes are much higher). The most conspicuous species is Quercus yiwuensis Y.C. Hsu & H.W. Jen. In March 2017 we selected four individual trees of Q. yiwuensis, and an equal number of Pistacia weinmannifolia Franch, and Beilschmiedia percoriacea C.K. Allen that were also numerous on the site. Plants were identified by the authors in the field and labelled. Botanical specimens were deposited in the School of Pharmaceutical Science and Technology, Tianjin University, Tianjin, China.Figure 6Map of the study site with approximate position of analysed trees (aerial perspective from Google Maps 2018). GPS positions were obtained less than 1 m from the tree trunks. The distance from N3 to B1 is approximately 60 m.Full size imageQ. yiwuensis was designated P-tree (P1, P2, P3, P4) because we consistently found the orchid Panisea uniflora Lindl. growing on this phorophyte species. P. weinmannifolia was designated B-tree (B1, B2, B3, B4) because it harboured the orchid species Bulbophyllum odoratissimum (Sw) Lindl. (Supplementary Fig. S2 a-d). On P. weinmannifolia trees no P. uniflora was observed, while B. odoratissimum was never found on Q. yiwuensis trees. Both tree species were richly colonized by several other orchid species. Beilschmiedia percoriacea trees were designated neutral tree, N-tree (N1, N2, N3, N4), because neither of the two target orchid species grew on them. The latter tree species carried several lichens and a single fern species (Lepisorus sp.), but only in one instance was observed to carry an orchid epiphyte (Coelogyne sp.).GPS positions of investigated trees were obtained less than 1 m from the trunk (Fig. 6). Accuracy is about 3 m. Accuracy in altitude readings is about 100 m. Distance between degrees of latitude is 111 km. At N 21.78978 the distance between degrees of longitude is 103 km, which means that the last digit in the 5-digit decimal degrees corresponds to 1.11 m in latitude and 1.03 m in longitude.The trees were labelled with different colours as follows:

P-trees (carrying P. uniflora and other epiphytes, but not B. odoratissimum), identified as Q. yiwuensis, with red labels (P1 N 21.79880, E 101.37909, 1073; P2 N 21.79882, E 101.37923, 1072; P3 N 21.79878, E 101.37947, 1074; P4 N 21.79878, E 101.37904, 1073).

B-trees (carrying B. odoratissimum and other epiphytes, but not P. uniflora), identified as P. weinmannifolia, with blue labels (B1 N 21.79878, E 101.37950, 1074; B2 N 21.79880, E 101.37938, 1078; B3 N 21.79881, E 101.37931, 1083; B4 N 21.79884, E 101.37923, 1076).

N-trees (carrying epiphytes, but neither B. odoratissimum nor P. uniflora), identified as B. percoriacea, with yellow labels (N1 N 21.79873, E 101.37905, 1064; N2 N 21.79868, E 101.37908, 1072; N3 N 21.79883, E 101.37893, 1071; N4 N 21.79879, E 101.37895, 1071).

The point of access to the outcrop top area was located at the Western edge (N 21.79880, E 101.37827, 1058, Fig. 6).

SamplingFor each of the twelve selected trees, breast height circumference (BH = 130 cm above ground) was measured. Approximate total height was determined by Nikon Laser Forestry Pro or estimated if sighting lines were interfered by other vegetation.The lowermost individual of the target orchid species was recorded in relation to BH. Bark samples were collected, and bark features recorded at BH, by target orchid, and 50 cm above target orchid or BH, whichever was highest point. In N-trees, where there were no target orchids, sampling was thus at BH, BH + 50 cm, and BH + 100 cm.Sampling on each tree involved approximately 12 cm2 bark cut out with a sterile knife and rubber gloves to prevent cross-contamination, for pH-analysis, metabarcoding, fungal isolation and chemical analysis. Besides, 3 bark cores were taken by trephor sampler (16 mm, 2 mm diam., Costruzioni Meccaniche Carabin Carlo) for water holding measurement.Roots of target orchids were sampled, from three adult individual plants on each P- and B-tree. No permissions were necessary to collect plant samples, using a protocol that avoided plant damages. All plants were left in the exact location where they were found in the sampling site, after collecting the small portions of bark and root material for the study. All experiments including the collection of plant material in this study are in compliance with relevant institutional, national, and international guidelines and legislation.All fresh material collected from the sampling site was first kept in cool boxes, brought to the laboratory, and processed within three days.Fungal isolation from barkFor each sample, half of the bark material and orchid roots were kept at − 80 °C for subsequent metabarcoding analysis. The rest of bark (about 2 g for each sample) was immediately processed for fungal isolation. The large bark portions were ground into powder using a sterile mortar and pestle; 5 ml were reserved for pH measurement, while the rest was suspended in a final volume 50 ml sterile water solution in a sterile centrifuge tube. The tube was shaken with Vortex vibration meter thoroughly and solution aliquots were spread homogenously onto isolation medium plates. For each bark sample, aliquots of 500, 300, 200, and 100 μl, were spread per triplicate to one plate each of PDA (Potato Dextrose Agar) medium, containing ampicillin (50 μg/mL) and streptomycin (50 μg/ml) to inhibit bacterial growth49,50. A diluted solution was also made by mixing 1 ml of the original solution with 9 ml sterile water and plated. Petri dishes were incubated at room temperature (23–25 °C) in the dark for up to 2 months to allow the development of slow-growing mycelia. Fast growing fungal strains started to grow after about two days. Colonies showing different morphology and appearance were transferred to fresh plates to obtain pure cultures. In the following days, other slower growing mycelia were available in the Petri dishes and were also regularly picked up and isolated onto new PDA plates every 2 days. All isolated fungal strains were stored at 4 °C for further analysis. All strains were deposited in the LP Culture Collection (personal culture collection held in the laboratory of Prof. Lorenzo Pecoraro), at the School of Pharmaceutical Science and Technology, Tianjin University, Tianjin, China.Molecular and morphological analysis of bark culturable fungiThe identification of isolated fungal colonies was performed using DNA sequencing combined with microscopy. Total genomic DNA from isolated fungi was extracted following the cetyltrimethyl ammonium bromide (CTAB) method modified from Doyle and Doyle51. Fungal ITS regions were PCR-amplified using the primer pair ITS1F/ITS452 following the procedure described in Pecoraro et al.37 for PCR reaction, thermal cycling, and purification of PCR products. Controls with no DNA were included in every amplification experiment in order to test for the presence of laboratory contamination from reagents and reaction buffers. Purified DNA amplicons were sequenced with the same primer pair used for amplification. DNA sequencing was performed at the GENEWIZ Company, Tianjin, China.Sequences were edited to remove vector sequences and to ensure correct orientation and assembled using Sequencher 4.1 for MacOsX (Genes Codes, Ann Arbor, MI). Sequence analysis was conducted with BLAST searches against the National Center for Biotechnology Information (NCBI) sequence database (GenBank; http://www.ncbi.nlm.nih. gov/BLAST/index.html) to determine the closest sequence matches that enabled taxonomic identification. DNA sequences were deposited in GenBank (Accession Nos. MW603206 – MW603451). Fungal morphological characters (hyphae, pseudohyphae, conidiophores, conidia, poroconidia, arthroconidia, etc.) were examined using a Nikon ECLIPSE Ci microscope for the identification of isolates following the standard taxonomic keys53,54,55,56,57.Assessment of bark and orchid associated fungal community using Illumina sequencingBark and orchid root samples were pulverized in a sterile mortar, and genomic DNA was extracted using the FastDNA® Spin Kit as described by the manufacturer (MP Biomedicals, Solon, OH, USA)58,59. In total, this resulted in 60 DNA samples, including 36 from bark (3 sampling points for each tree × 12 trees) and 24 from orchid roots (3 orchid individuals sampled on each P- and B-tree × 8 trees; the 4 individual N-trees were not used for orchid sampling because they did not carry the study orchid species). Subsequently, amplicon libraries were created using two primer combinations targeting the internal transcribed spacer 2 (ITS-2): ITS7F and ITS4R60 was used as universal fungal primer pair to target nearly all fungal species, while ITS361 and ITS4OF62 was used to more specifically target orchid mycorrhizal fungi. Previous research has shown that most universal fungal primers have multiple mismatches to many species of the orchid-associating basidiomycetes, in particular in Tulasnellaceae family46,58,63. Since the goal of the present work was to analyse the total fungal community in the orchid-phorophyte environment (bark and orchid roots), as well as more specifically detect the orchid mycorrhizal fungi in the studied samples, it was necessary to combine two different primer pairs to characterise the whole investigated fungal diversity47,64,65,66. Polymerase chain reaction (PCR) amplification was performed in 50 μl reaction volume, containing 38 μl steril distilled water, 5 μl 10 × buffer (100 mM Tris–HCl pH 8.3, 500 mM KCl, 11 mM MgCl2, 0.1% gelatin), 1 μl of dNTP mixture of 10 mM concentration, 0.25 μM of each primer, 1.5 U of RED TaqTM DNA polymerase (Sigma) and approximately 10 μg of extracted genomic DNA. PCR conditions were as follows: 1 cycle of 95 °C for 5 min initial denaturation before thermocycling, 30 cycles of 94 °C for 40 s denaturation, 45 s annealing at various temperatures following Taylor and McCormick62, 72 °C for 40 s elongation, followed by 1 cycle of 72 °C for 7 min extension. To minimize PCR bias, three PCRs were pooled for each sample. The resulting PCR products were electrophoresed in 1% agarose gel with ethidium bromide and purified with the QIAEX II Gel Extraction Kit (QIAGEN). Amplicon libraries were generated using the NEB Next Ultra DNA Library Prep Kit for Illumina (New England Biolabs, USA) following the manufacturer’s instructions to add index codes. Samples were sequenced using the Illumina MiSeq PE 250 sequencing platform (Illumina Inc., San Diego, CA) at Shanghai Majorbio Bio‐Pharm Technology Co., Ltd. (Shanghai, China).Bioinformatics of fungal sequencesSequences originated from the total (ITS7F and ITS4R primers) and orchid-associated (ITS3 and ITS4OF primers) fungi datasets were processed separately. Raw reads were merged with a minimum overlap of 30 nucleotides, and the primer sequences were trimmed off. Subsequently, reads were filtered by discarding all sequences with expected error  > 1. The quality-filtered reads were denoised using the UNOISE3 algorithm67 to create zero-radius operational taxonomic units (zOTUs), with chimera removal. All the steps were performed using USEARCH v.1168. Raw sequences have been deposited in the Sequences Read Archive (SRA) of NCBI as BioProject ID PRJNA702612. The fungal zOTUs were assigned to taxonomic groups using the Blast algorithm by querying against the UNITE + INSD fungal ITS database (version 7.2, released on 10 October 2017)69 using the sintax algorithm with 0.8 cutoff70. The zOTUs originated with the orchid-associated fungal primers were manually screened for possible orchid-associated mycorrhizal families based on the information provided in Table 12.1 in Dearnaley et al.71, and only these were retained for further analysis in this dataset.To attempt removing spurious counts due to cross-talk (assignment of reads to a wrong sample) we removed all the zOTUs represented by less than 0.02% of reads in each sample, which is more conservative than previous error estimates72. The datasets were rarefied to the minimum sequencing depth (23,419 for total fungi and 13,074 for orchid-associated fungi), zOTUs present in less than three samples and low abundant zOTUs (with relative abundance  More

1.Falkowski PG, Fenchel T, Delong EF. The microbial engines that drive Earth’s biogeochemical cycles. Science. 2008;320:1034–9.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Arnosti C. Microbial extracellular enzymes and the marine carbon cycle. Ann Rev Mar Sci. 2011;3:401–25.PubMed 
Article 
PubMed Central 

Google Scholar 
3.Arnosti C, Bell C, Moorhead DL, Sinsabaugh RL, Steen AD, Stromberger M, et al. Extracellular enzymes in terrestrial, freshwater, and marine environments: perspectives on system variability and common research needs. Biogeochemistry. 2013;117:5–21.Article 
CAS 

Google Scholar 
4.Moorhead DL, Rinkes ZL, Sinsabaugh RL, Weintraub MN. Dynamic relationships between microbial biomass, respiration, inorganic nutrients and enzyme activities: informing enzyme-based decomposition models. Front Microbiol. 2013;4:223.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Li M, Gao Y, Qian W-J, Shi L, Liu Y, Nelson WC, et al. Targeted quantification of functional enzyme dynamics in environmental samples for microbially mediated biogeochemical processes. Environ Microbiol Rep. 2017;9:512–21.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Song H-S, Thomas DG, Stegen JC, Li M, Liu C, Song X, et al. Regulation-structured dynamic metabolic model provides a potential mechanism for delayed enzyme response in denitrification process. Front Microbiol. 2017;8:1866.PubMed 
PubMed Central 
Article 

Google Scholar 
7.Khersonsky O, Tawfik DS. Enzyme promiscuity: a mechanistic and evolutionary perspective. Annu Rev Biochem. 2010;79:471–505.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Baier F, Copp JN, Tokuriki N. Evolution of enzyme superfamilies: comprehensive exploration of sequence-function relationships. Biochemistry. 2016;55:6375–88.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
9.Sebastián M, Niell FX. Alkaline phosphatase activity in marine oligotrophic environments: implications of single-substrate addition assays for potential activity estimations. Mar Ecol Prog Ser. 2004;277:285–90.Article 

Google Scholar 
10.Catrina I, O’Brien PJ, Purcell J, Nikolic-Hughes I, Zalatan JG, et al. Probing the origin of the compromised catalysis of E. coli alkaline phosphatase in its promiscuous sulfatase reaction. J Am Chem Soc. 2007;129:5760–5.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
11.Sunden F, AlSadhan I, Lyubimov AY, Ressl S, Wiersma-Koch H, Borland J, et al. Mechanistic and evolutionary insights from comparative enzymology of phosphomonoesterases and phosphodiesterases across the alkaline phosphatase superfamily. J Am Chem Soc. 2016;138:14273–87.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
12.Yang K, Metcalf WW. A new activity for an old enzyme: Escherichia coli bacterial alkaline phosphatase is a phosphite-dependent hydrogenase. Proc Natl Acad Sci USA. 2004;101:7919–24.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Copley SD. Shining a light on enzyme promiscuity. Curr Opin Struct Biol. 2017;47:167–75.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Steen AD, Vazin JP, Hagen SM, Mulligan KH, Wilhelm SW. Substrate specificity of aquatic extracellular peptidases assessed by competitive inhibition assays using synthetic substrates. Aquat Micro Ecol. 2015;75:271–81.Article 

Google Scholar 
15.Ivars-Martínez E, D’Auria G, RodrÍGuez-Valera F, SÁNchez-Porro C, Ventosa A, et al. Biogeography of the ubiquitous marine bacterium Alteromonas macleodii determined by multilocus sequence analysis. Mol Ecol. 2008;17:4092–106.PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
16.Tada Y, Taniguchi A, Nagao I, Miki T, Uematsu M, Tsuda A, et al. Differing growth responses of major phylogenetic groups of marine bacteria to natural phytoplankton blooms in the western North Pacific Ocean. Appl Environ Microbiol. 2011;77:4055–65.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Schubert M, Lindgreen S, Orlando L. AdapterRemoval v2: rapid adapter trimming, identification, and read merging. BMC Res Notes. 2016;9:88.PubMed 
PubMed Central 
Article 

Google Scholar 
18.Li D, Liu C, Luo R, Sadakane K, Lam T. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics. 2015;31:1674–6.CAS 
PubMed 
Article 
PubMed Central 

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

Google Scholar 
20.Bushnell B. “BBMap: a fast, accurate, splice-aware aligner,” in Proceedings of the 9th Annual Genomics of Energy & Environment Meeting. Walnut Creek, CA, USA; 2014.21.Scholz J, Besir H, Strasser C, Suppmann S. A new method to customize protein expression vectors for fast, efficient and background free parallel cloning. BMC Biotechnol. 2013;13:12.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
22.McLoughlin SY, Jackson C, Liu JW, Ollis DL. Growth of Escherichia coli coexpressing phosphotriesterase and glycerophosphodiester phosphodiesterase, using paraoxon as the sole phosphorus source. Appl Environ Microbiol. 2004;70:404–12.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
23.Britton J, Dyer RP, Majumdar S, Raston CL, Weiss GA. Ten-minute protein purification and surface tethering for continuous-flow biocatalysis. Angew Chem Int Ed Engl. 2017;56:2296–301.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
24.Ortiz-Tena JG, Rühmann B, Sieber V. Colorimetric determination of sulfate via an enzyme cascade for high-throughput detection of sulfatase activity. Anal Chem. 2018;90:2526–33.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Huitema C, Horsman G. Analyzing enzyme kinetic data using the powerful statistical capabilities of R. 2018. http://biorxiv.org/content/10.1101/316588v1.26.Rainer SF. Soft-bottom benthic communities in Otago Harbour and Blueskin Bay, New Zealand. New Zealand Oceanographic Institute Memoir 80; 1981.27.Grove SL, Probert PK. Sediment macrobenthos of upper Otago Harbour, New Zealand. New Zeal J Mar Fresh. 1999;33:469–80.28.Hoppe HG. Significance of exoenzymatic activities in the ecology of brackish water: measurements by means of methylumbelliferyl-substrates. Mar Ecol Prog Ser. 1983;11:299–308.CAS 
Article 

Google Scholar 
29.Yamaguchi T, Sato M, Hashihama F, Ehama M, Shiozaki T, Takahashi K, et al. Basin‐scale variations in labile dissolved phosphoric monoesters and diesters in the central North Pacific Ocean. J Geophys Res Oceans. 2019;124:3058–72.CAS 
Article 

Google Scholar 
30.Baltar F, Lundin D, Palovaara J, Lekunberri I, Reinthaler T, Herndl GJ, et al. Prokaryotic responses to ammonium and organic carbon reveal alternative CO2 fixation pathways and importance of alkaline phosphatase in the mesopelagic North Atlantic. Front Microbiol. 2016;7:1670.PubMed 
PubMed Central 
Article 

Google Scholar 
31.Yamaguchi H, Arisaka H, Seki M, Adachi M, Kimura K, Tomaru Y. Phosphotriesterase activity in marine bacteria of the genera Phaeobacter, Ruegeria, and Thalassospira. Int Biodeter Biodegr. 2016;115:186–91.CAS 
Article 

Google Scholar 
32.Davidi D, Noor E, Liebermeister W, Bar-Even A, Flamholz A, Tummler K, et al. Global characterization of in vivo enzyme catalytic rates and their correspondence to in vitro kcat measurements. Proc Natl Acad Sci USA. 2016;113:3401–6.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Paytan A, Cade-Menum BJ, McLaughlin K, Faul KL. Selective phosphorus regeneration of sinking marine particles: evidence from 31P-NMR. Mar Chem. 2003;82:55–70.CAS 
Article 

Google Scholar 
34.Wu J, Wang P, Wang Y. Cytotoxic and mutagenic properties of alkyl phosphotriester lesions in Escherichia coli cells. Nucleic Acids Res. 2018;46:4013–21.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.McCarthy JG, Edington BV, Schendel PF. Inducible repair of phosphotriesters in Escherichia coli. Proc Natl Acad Sci USA. 1983;80:7380–4.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
36.Helbert W. Marine polysaccharide sulfatases. Front Mar Sci. 2017;4:6.Article 

Google Scholar 
37.Wegner CE, Richter-Heitmann T, Klindworth A, Klockow C, Richter M, Achstetter T, et al. Expression of sulfatases in Rhodopirellula baltica and the diversity of sulfatases in the genus Rhodopirellula. Mar Genomics. 2013;9:51–61.PubMed 
Article 
PubMed Central 

Google Scholar 
38.Canfield DE, Farquhar J. Animal evolution, bioturbation, and the sulfate concentration of the oceans. Proc Natl Acad Sci USA. 2009;106:8123–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Luo HW, Benner R, Long RA, Hu JJ. Subcellular localization of marine bacterial alkaline phosphatases. Proc Nat Acad Sci USA. 2009;106:21219–23.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
40.Wu J-R, Shien J-H, Shieh HK, Hu C-C, Gong S-R, Chen L-Y, et al. Cloning of the gene and characterization of the enzymatic properties of the monomeric alkaline phosphatase (PhoX) from Pasteurella multocida strain X-73. FEMS Microbiol Lett. 2007;267:113–20.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
41.Kageyama H, Tripathi K, Rai AK, Cha-um S, Waditee-Sirisattha R, Takabe T. An alkaline phosphatase/phosphodiesterase, PhoD, induced by salt stress and secreted out of the cells of Aphanothece halophytica, a halotolerant cyanobacterium. Appl Environ Micro. 2011;77:5178–83.CAS 
Article 

Google Scholar 
42.Rodriguez F, Lillington J, Johnson S, Timmel CR, Lea SM, Berks BC. Crystal structure of the Bacillus subtilis phosphodiesterase PhoD reveals an iron and calcium-containing active site. J Biol Chem. 2014;289:30889–99.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
43.Noskova Y, Likhatskaya G, Terentieva N, Son O, Tekutyeva L, Balabanova L. A novel alkaline phosphatase/phosphodiesterase, CamPhoD, from marine bacterium Cobetia amphilecti KMM 296. Mar Drugs. 2019;17:657.CAS 
PubMed Central 
Article 

Google Scholar 
44.Dyhrman ST, Ammerman JW, Van, Mooy BAS. Microbes and the marine phosphorus cycle. Oceanography. 2007;20:110–6.Article 

Google Scholar 
45.Larson TJ, Ehrmann M, Boos W. Periplasmic glycerophosphodiester phosphodiesterase of Escherichia coli, a new enzyme of the glp regulon. J Biol Chem. 1983;258:5428–32.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.van Veen HW. Phosphate transport in prokaryotes: molecules, mediators and mechanisms. Antonie Van Leeuwenhoek. 1997;72:299–315.PubMed 
Article 
PubMed Central 

Google Scholar 
47.Parthasarathy S, Parapatla H, Nandavaram A, Palmer T, Siddavattam D. Organophosphate hydrolase is a lipoprotein and interacts with Pi-specific transport system to facilitate growth of Brevundimonas diminuta using op insecticide as source of phosphate. J Biol Chem. 2016;291:7774–85.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
48.Hong T, Kong A, Lam J, Young L. Periplasmic alkaline phosphatase activity and abundance in Escherichia coli B23 and C29 during exponential and stationary phase. J Exp Microbiol Immunol. 2007;11:8–13.
Google Scholar 
49.Baltar F, Arístegui J, Gasol J, Yokokawa T, Herndl GJ. Bacterial versus archaeal origin of extracellular enzymatic activity in the Northeast Atlantic deep waters. Micro Ecol. 2013;65:277–88.CAS 
Article 

Google Scholar 
50.Thomson B, Wenley J, Currie K, Hepburn C, Herndl GJ, Baltar F. Resolving the paradox: continuous cell-free alkaline phosphatase activity despite high phosphate concentrations. Mar Chem. 2019;214:103671.CAS 
Article 

Google Scholar 
51.Lei L, Cherukuri KP, Alcolombri U, Meltzer D, Tawfik DS. The dimethylsulfoniopropionate (DMSP) lyase and lyase-like cupin family consists of bona fide DMSP lyases as well as other enzymes with unknown function. Biochemistry. 2018;57:3364–77.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Ferla MP, Brewster JL, Hall KR, Evans GB, Patrick WM. Primordial‐like enzymes from bacteria with reduced genomes. Primordial-like enzymes from bacteria with reduced genomes. Mol Microbiol. 2017;105:508–24.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More

Study siteThis study was conducted at La Selva Biological Station, located in the lowland Atlantic forest of Costa Rica (10°26′ N, 83°59′ W). The mean annual temperature is 26 °C; mean annual precipitation is 4 m and all months have mean precipitation  > 100 mm39. La Selva has undulating topography, with elevation varying between 10 and 140 m above sea level. La Selva Biological Station includes multiple land uses; our analysis includes 103.5 hectares of forest, comprising 33.0 ha of old-growth forest and 70.5 ha of forests with past human disturbance (secondary forests, abandoned agroforestry, abandoned plantation, selectively-logged forests); here, we refer to all areas with past human disturbance as “secondary forests”. This study area does not include the full extent of old-growth or secondary forests at La Selva—we focused our drone data collection on this area because it contained the most severe apparent disturbance from the blowdown. Forests with past human disturbance have been naturally regenerating for a range of time (since 1955–1988); we excluded secondary forests with regeneration starting after 1988.Lidar dataWe use two airborne lidar datasets to quantify dynamics in canopy structure and ACD. Data were collected in 2009 and 2019 (Supplementary Table 2). Data from 2009 were collected by a fixed-wing aircraft over the entire reserve; data from 2019 were collected using the Brown Platform for Autonomous Remote Sensing40. We focused on an area 1.4 km2 in size that includes the region of most severe damage from the blowdown (Supplementary Fig. 1). Both lidar sensors were discrete-return systems. To minimize variation in lidar height estimates from variable laser beam divergence and detector characteristics, we only used data from first returns for all analyses. For the 2019 drone-based lidar with higher native point density and a wider scan angle range40, we limited our analysis to lidar returns with scan angle ± 15 degrees and randomly subsampled data to a homogenous resolution of 10 pts m−2. Previous research demonstrates that lidar data collected above densities of 1 pts m−2 have similar predictive power for determining many forest properties (including tree height, tree density, and basal area)41; both lidar datasets in this study are above this density threshold. All lidar data were projected using EPSG 32,616.For all lidar data, we calculated height above ground using a digital terrain model (DTM) created from lidar data collected in 2006 and validated using 4184 independent measurements within the old-growth forest (intercept =  − 0.406, slope = 0.999, r2 = 0.994, RMSE = 1.85 m; Supplementary Table 2)42. We verified that the horizontal geolocation accuracy with  More

Examining the serum metabolome profiles of bighorn sheep captured by the three primary techniques used to capture wild ungulates revealed significant changes in polar metabolite levels between the different animal groups, and trends that persisted throughout the analyses when directly comparing, in a pairwise fashion, specific capture techniques. Results from PLS-DA modeling and analysis of the top 15 metabolites that contribute most (VIP  > 1.2) to the separation of the three capture groups revealed that amino acid levels of tryptophan, valine, isoleucine, phenylalanine, and proline were highest in animals captured by dart, with intermediate levels in animals capture using dropnets, and lowest in animals captured using the helicopter method (Fig. 3A). One-way ANOVA analyses identified additional amino acids that displayed similar decreasing level trends from dart to dropnet to helicopter capture (dart  > drop net  > helicopter) methods, and included arginine, asparagine, aspartate, cysteine, glutamate, and glutamine, glycine, histidine, leucine, lysine, serine, and tyrosine (Fig. 4). These metabolite level changes suggest a shift in amino acid metabolism, and a potentially higher catabolism of these compounds as a function of increasingly more energetically intense and possibly more stressful capture methods such as helicopter capture.Of these amino acids, aspartate, glycine, and glutamate function as precursors for neurotransmitter synthesis, and may therefore be valuable indicators of the capture techniques’ impacts on animal health and changes to their physiological state. Glutamate is a fundamental component of nitrogen excretion in the urea cycle, and its lower serum levels in animals captured by helicopter support the idea of altered metabolite flow through the urea cycle. In addition to these patterns, decreasing levels of aspartate were observed in samples of dropnet and helicopter captured animals compared to the levels found in the dart-captured animals. The change regarding urea cycle alterations also manifested itself in differential serum urea levels, with fold changes (FC) between the groups decreasing significantly with capture techniques, with a mean FC difference of 1.4 for the dart-captured group, 0.26 for the dropnet-captured group, and − 0.3 for the helicopter-captured animals (Supplementary Table S2). As urea recycling is a prominent feature of ruminant metabolism and urea flux can rapidly change, the urea concentration changes observed between the three capture techniques support an impact on urea cycle intermediates29. While the trend of an overall decrease in urea cycle intermediates parallels a similar trend in amino acid concentrations, the extent to which amino acid metabolism is linked to changes in urea cycle activity is difficult to evaluate due to the nature of nitrogen recycling in the rumen of these ruminants.Other metabolites found in significantly higher concentrations in the serum samples of dart-captured animals compared to the two other techniques included: formate, glucose, 3-hydroxybutyrate, dimethylamine, carnitine (Fig. 3A). Propionate, which was observed to be higher in the dart and dropnet captured animals than that of helicopter captured animals (Fig. 4) is of interest, as it is the main precursor for glucose synthesis in the liver of ruminants30, and potentially reflect a higher dependence of ruminants on gluconeogenesis due to the almost complete conversion of available dietary carbohydrates to volatile fatty acids in the rumen31. As animal capture via nets increases physical activity as the animals struggle to free themselves from entanglement, generally resulting in longer times animals are under physical restraint, as well as the increased physical exertion and stress as they attempt to flee the pursuing helicopter, the observed decrease in serum propionate levels may reflect increased needs to generate glucose de novo via gluconeogenesis.This interpretation of the metabolite data is reinforced by the observation of significantly elevated levels of O-acetylcarnitine in the drop net and helicopter net gun animal capture groups compared to the darted animals (Fig. 4). As an important element of the carnitine/acyl-carnitine shuttle and import of fatty acids into the mitochondria for β-oxidation, acyl-carnitine is a major contributor to the flow of acyl groups into the TCA cycle, and a robust indicator of cardiac output and, by extension, TCA cycle activity levels in mammals32. Additional metabolites that displayed distinctly increasing trends based on capture method (dart  More