Operational taxon pairs in microbial co-occurrence networks imply complicated bacterial associations

doi:10.21203/rs.3.rs-1535010/v1

Download PDF

Research Article

Operational taxon pairs in microbial co-occurrence networks imply complicated bacterial associations

https://doi.org/10.21203/rs.3.rs-1535010/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

Background

Microbial co-occurrence networks inferred from the abundance data of environmental and human gut samples are widely applied to describe microbial associations and to predict potential microbial interactions. However, the microbial associations and predicted interactions were hardly demonstrated with cultivated microbial isolates. Experimental demonstration of the microbial associations and the predicted interactions needs extensive isolation and cultivation of microbes and matches the isolates to operational taxon units (OTUs) from co-occurrence networks. The high workloads of bacterial isolation and cultivation, the shortage of methods to match cultivated bacterial isolates to co-occurrence network OTUs, and the complexity of the networks themselves were factors that constrained experimental demonstrations.

Results

Here, we integrated droplet microfluidics and bar-coding logistics for high-throughput bacterial isolation and cultivation from samples. The co-occurrence network complexity of samples was reduced by dilution of samples. Phylogenetic tree topology was applied to match bacterial isolates to zero-level OTUs (Zotus). We collected wetland samples from Beijing Olympic Park (BOP) including Potamogeton perfoliatus plants, their roots and root peripheral sediments. Droplets of series diluted homogenates of wetland samples were inoculated into 126 of 96-well plates containing R2A and TSB media. After 10 days of cultivation, 65 plates (each with > 30% wells of the 96 wells showed microbial growth) were elected to reconstruct co-occurrence networks. The prevalent co-occurred zero-level OTUs (Zotus) pairs were selected from the networks and were matched to bacterial isolates. We cultivated 129 bacterial isolates belonged to 15 taxa at species level. Thirty-six bacterial isolate combinations corresponding to 19 Zotus pairs were evaluated for interactions on agar plates. Results suggested that positively associated Zotu pairs in the co-occurrence network implied complicated relations including neutralism, competition, and mutualism. We further tested in broth the bacterial isolate combinations between BOP-1, BOP-5, BOP-11 or BOP-16 and BOP-108 that represented the Zotu1-Zotu10 pair. Results revealed that this Zotu1-Zotu10 pair had also competitive, neutral, or mutual relation, depending on bacterial isolate combination and on cultivation time.

Conclusion

This study established a workflow for experimental demonstration of relations between Zotu pairs from microbial co-occurrence networks. Using a wetland sample, we reconstructed microbial co-occurrence networks and cultivated extensively bacterial isolates, and further demonstrated that positively associated Zotu pairs had neutral, competitive, and mutual relations during bacterial cultivations on solid agar plates and in liquid broth.

Co-occurrence network

matching bacterial isolates and Zotu pairs

microbial interaction

droplet microfluidics

microbial cultivation

Microorganisms in natural or human-engineered ecosystems interact with each other by means of mutualism, commensalism, parasitism/predation, ammenalism, and competition[1]. The sum of all microbial interactions in an ecosystem constitutes a microbial interactive network[2]. So far, the true microbial interactive networks are hardly understood, due to the microbiome complexity arisen from the diverse microbial species and their interactions in the ecosystems. On the other hand, more and more microbiome data are accumulated and those data could be analyzed in silico on multiple platforms, such as parallel-meta suite[3]. Microbial co-occurrence networks were usually reconstructed from metagenomic data or from high-throughput sequencing of 16S rRNA gene amplicons from environmental and host microbiome DNA molecules[4–5]. With the information of abundance data, the microbial co-occurrence networks can be inferred using similarity-based (e.g., Pearson or Spearman correlations)[6–10] or model-based methods (e.g., regression- and rule-based)[11–13]. For examples, microbial co-occurrence networks of activated sludge from wastewater treatment plants[14], human gut[15] and marine environment[10, 16] were inferred from data of high-resolution and extended longitudinal and cross-sectional scales of samples. The inferred microbial co-occurrence networks were widely applied for predicting microbial interactions. But such prediction needs to be carefully examined and cautiously interpreted[5], since most of the predictions were not experimentally demonstrated. A positive relation within a co-occurrence network could be due to cross-feeding, co-colonization, niche overlap, or other random reasons[17–19] and a negative relation could be due to competition, predation, and antagonism[20–21].

Microbial co-occurrence networks compose of centers [operational taxon units (OTUs)] and edges (predicted relations among microbes)[8, 22–23]. The more centers and edges, the more complexity are the microbial co-occurrence networks[24–25]. The microbial co-occurrence network of a real ecosystem usually harbors hundreds of centers and thousands of edges, which would be too complexed to be experimentally demonstrated[1]. We conceived if the complexed co-occurrence networks could be simplified into less complicated networks that were composed of less centers and less edges, and if the microbes representing the centers could be cultivated, the microbial interaction predicted from co-occurrence networks could be tested. Apparently, the simplest microbial interaction would be microbial pairs in the networks. Inspired by this conception, we designed the current study with high-throughput microfluidic microbial cultivation[26], DNA bar-coding[27], zero-level operational taxon units (Zotus) analysis and similarity-based data analysis approaches[28]. We sampled the wetland ecosystem of Beijing Olympic Park (BOP). Microbial co-occurrence networks were reconstructed, prevalent Zotu pairs were identified, and extensive isolation and cultivation of microbes were performed. Zotus pairs were matched to bacterial isolates, and the microbial interactions were experimentally demonstrated.

2.1 Sample collection and treatments

In this study, submerged Potamogeton perfoliatus plant and root and sediment were collected from the dragon-shape wetland ecosystem of BOP (Longitude E116°23′2.98″; Latitude N40°01′3.00″) on November 5, 2019. Three sites at distance of 40 m were selected, and five P. perfoliatus plants, their roots and root peripheralsediment were collected. Samples were transferred quickly to laboratory and were processed according to Bai et al. ^[25]. Briefly, the aboveground (leaves and stems) and underground parts (roots) of each P. perfoliatus were separated with sterile scissors. The aboveground parts were washed with 150 mL of sterilized PBS solution (containing NaCl 136 mM, Na₂HPO₄ 8 mM, KH₂PO₄ 2 mM, KCl 2.6 mM, pH 7.4) in a 50 mL sterile tube. The roots were washed into a 250 mL conical flask with 100 mL sterilized PBS (shaken at 160 rpm for 30 min, twice). To cultivate the microbial communities associated with samples, 0.5 g of aboveground parts and 0.5 g of roots were individually homogenized with 1 mL MgCl₂ solution (10 mM) and then 5 mL sterile PBS buffer were added. Sediments (0.5 g) were suspended in 5 mL sterile PBS. The plant homogenates and sediment suspension were filtered twice using 40-µm cell strainer (BD Falcon, USA) and diluted into 10^− 1 ~ 10^− 7 series.

2.2 Microbial community inoculation and cultivation with microfluidics

The 10^− 1 ~ 10^− 7 dilution series were used for inoculation and subsequently microbial community cultivation. A automatic liquid dispenser (MultiDrop Combin nL, Thermo Fisher Scientific Inc., USA) was applied. The microfluidic infusion tubes of liquid dispenser were washed with 50 mL 75% ethanol followed by sterile water. Then, 10 mL of diluted samples were pipetted into separation bottles, and 1 µL homogenate of each diluted sample was inoculated into each well containing 150 µL R2A[29] or Tryptone Soy Broth (TSB) (Hope Bio-Technology, China) of 96-well plates in triplicates, with the MultiDrop Combin nL liquid dispenser. Totally, 126 of 96-well plates were inoculated. The inoculated plates were cultivated at 30 ℃ for 10 days. Depending on the dilutions and cell densities in samples, the numbers of wells that showed microbial growth in the 96-well plates were varied. We maintained those plates that more than 30% of wells showed microbial growth for further analysis. With this cutoff value, 96 of the 126 plates were selected and were stored in 20% glycerol stocks (by adding 50 µL 60% glycerol in each well) at -80 ℃ until use.

2.3 DNA extraction, PCR amplification, bar-coding and sequencing

The biomasses in the 96 plates were harvested and used for DNA extraction. The V4 regions of 16S rRNA genes were amplified from DNA molecules with 515F and 806R primers (Supplementary Table S1) by using 2x Dream Taq Green PCR Master Mix (Thermo Fisher Scientific Inc., USA), and the PCR products were purified with AMPure XP (Beckman Coulter GmbH, Krefeld, Germany). The PCR products were then used for DNA library construction[30–31] and sequenced on a Novaseq 6000 platform (Illumina Inc., USA) (250 bp paired-end reads) by Guangdong Magigene Biotechnology Co., Ltd, China. We adopted a two-sided barcode polymerase chain reaction (PCR) system, as described in the references [27, 32]. Barcode primers used in this study are listed in Supplementary Table S1.

2.4 Data analysis, reconstruction of co-occurrence networks, and identification of Zotu pairs

The raw DNA sequence data analysis was processed using UPARSE pipeline (http://drive5.com/usearch/manual/uparse_pipeline.html) with assemble paired reads (-fastq_mergepairs). Sequences were orientated by Silva 132 database(-orient)[28, 33], length- trimmed for 250bp (-fastx_truncate), and the singleton, chimeric and plant DNA sequences were removed. After quality control, the sequences were extracted into feature table (-otutab). Assignment for Zotu was performed according to the Silva 132 database (https://www.arb-silva.de/) and National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/). The pairwise Spearman’s rank correlations between detected Zotus in 96-well plates were calculated by using R (version 3.6.2; https://www.r-project.org/) package psych (version 2.1.6). The Zotus showing occurrence frequencies < 30% of wells in a given 96-well plate were filtered. Subsequently, the edges and nodes information of co-occurrence networks were calculated using the R package igraph (version 1.2.6). In order to simplify the downstream experiments, we summarized the co-occurrence networks with degree centrality < 5 by using Pajek (version 5.11)[34] and Gephi (version 0.9.2)[35] and the frequency of occurrence of a given Zotu pair was calculated as \(\text{F}\text{O}=\text{N}\text{p}/\text{N}\text{g}\), where FO was the frequency of occurrence of the given Zotu pair, Np was the number of the given Zotu pair and Ng was the number of 96-well plate in each group. Finally, 65 microbial co-occurrence networks were reconstructed. Zotu pairs were identified from these co-occurrence networks, and the robust and prevalent Zotu pairs were defined as those Zotu pairs showing Spearman’s |ρ| >0.6 and P < 0.01 and occurrence frequencies > 30% from the in-silico analysis.

2.5 Isolation of bacteria and matching the robust and prevalent Zotu pairs

To obtain bacterial isolates representing the robust and prevalent Zotu pairs in the networks, we examined 16S rRNA gene abundances in the 96-well plates, and selected the wells with high abundances of targeted Zotus for bacterial isolation. The isolation was performed as follows: 1) Transferred 20 µL cell suspension from the selected well to a new sterile 1.5 mL EP tube; 2) The cell suspension was diluted into 10^− 1 ~ 10^− 7 series; 3) 100 µL of each dilution were spread onto the R2A or TSB agar plates; 4) Single colonies on the agar plates after incubation for 5 and 10 days were picked up and checked for purity by 16S rRNA gene sequencing. The taxonomy of bacterial isolates was determined following a protocol as previously described [36]. The bacterial isolates that matched the targeted Zotus were determined based on phylogenetic tree topology and 16S rRNA genes sequence similarity between the isolates and the targeted Zotus, and were used for demonstration of Zotu pairs relation on agar plates or liquid broth.

2.6 Demonstration of microbial interactions with agar plate and liquid broth

The interactions between the bacterial isolates that matched Zotu pairs were validated by cultivation of the representative isolates on agar plates and in liquid broth. For agar plate experiments, 3 µL of the representative bacterial isolate cultures were dripped onto a TSB agar plate surface, and then 3 µL of the other matched bacterial isolate were dripped at external tangency to the firstly dripped representative isolate. The plates were cultivated at 30 ℃ for 24 h and photograph was recorded.

In addition to the experiment performed on agar plates, we also evaluated the interactions of bacterial isolates matching Zotu pairs in liquid broth. The bacterial isolates were cultivated in 5 mL TSB medium at 30 ℃ and 160 rpm overnight. The growing cells were harvested at 8,000 rpm for 5 min and washed twice, and resuspended in sterile PBS to OD₆₀₀ = 0.5 (Ultramicro spectrophotometer B-500Nabi; Yuanxi Instrument Co., Ltd, China). The mono- and co-culture experiments were inoculated with the same volumes (25 µL) of cell suspensions and cultured in 5 mL sterile TSB buffer. The inoculated TSB medium was cultured at 30 ℃ and 160 rpm. Experiments were done in triplicates. The growths were monitored with an Automated Microbiology Growth Curve Analysis System Screen C (OY Growth Curves AB Ltd., Finland) and quantified using SYBR Green quantitative real-time PCR (qPCR) (qTOWER3/G, Analytik Jena, Germany), with primer sets specifically targeting genus Pseudomonas (Forward: CGTAGGTGGTTTGTTAAGTTGGATGT; Reverse: GCACCTCAGTGTCAGTATCAGT) and genus Aeromonas (Forward: GATTTGGAGGCTGTGTCCTTGAGAC; Reverse: AGGATTCCAGACATGTCAAGGCCA). qPCR procedures are detailed in the next paragraphs.

2.7 qPCR for quantification of Pseudomonas and Aeromonas cells in liquid cultures

We amplified the DNA fragments specific to Pseudomonas and Aeromonas with the above primers by using 2x Dream Taq Green PCR Master Mix (Thermo Fisher Scientific Inc., USA) and purified by Universal DNA Purification Kit (Tiangen Biotech Co., Ltd, China). Thus, we obtained PCR products featuring the genus Pseudomonas (185 bp) or the genus Aeromonas (186 bp). The PCR products were linked with pGM-T plasmid using pGM-T Cloning Kit (Tiangen Biotech Co., Ltd, China). The ligated DNA product was transformed into E. coli DH5α competent cells and plated onto LB agar plate (containing 50 µg/mL Ampicillin) at 37 ℃, overnight. The constructed plasmids were extracted using TIANprep Mini Plasmid Kit (Tiangen Biotech Co., Ltd, China), sequenced, and named p-Pseu and p-Aero, respectively. The reconstructed plasmid concentration was measured using Qubit 4.0 (Thermo Fisher Scientific, USA). The templates were diluted with sterile ddH₂O in 8 folds (10^− 1 ~ 10^− 8 copies/ µL)[37]. The standard curve and quantitation of co-cultured bacterial isolates were amplified by using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) and normalized to the copy numbers of the 16S rRNA gene of each isolate (https://rrndb.umms.med.umich.edu/).

3.1 Conception of the workflow to demonstrate the microbial associations from co-occurrence networks with microbial cultivation

Microbial co-occurrence networks compose of centers (OTUs) and edges (microbe-microbe associations). We hypothesized that the microbial associations could be demonstrated if the complexed co-occurrence networks are simplified and composed of less centers and less edges, and if the microbes representing the centers be cultivated. To work on this hypothesis, we designed a workflow chart as showed in Fig. 1. In this workflow chart, samples were first processed, serially diluted in buffer and dilution parts were used for co-occurrence network reconstruction and subsequently targeted bacterial isolation. A high-through droplet microfluidic system was applied for inoculation of 96-well plates. After cultivation 10 days, microbial DNAs in each plate were extracted, bar-coded, and sequenced for reconstruction of co-occurrence networks. The Zotu pairs in networks from plates were extracted according to statistics, and prevalent Zotus were elected. The wells of plates showing high abundances of prevalent Zotus were targeted for subsequent microbial isolation and cultivation. Lastly, the cultivated microbial isolates were matched to Zotus in the network and used for demonstration of microbial interactions.

3.2 Prevalent Zotu pairs in the co-occurrence networks

Following the workflow chart we developed in this study, 126 of 96-well plates were inoculated with samples from the BOP wetland system with dilutions of 10^− 1 to 10^− 7. Depending on the microbial density in samples, the 96-well plates harbored different numbers of wells with microbial growth. We considered that those plates had less than 30% wells with microbial growth would significantly decrease the quality of co-occurrence network reconstruction, and those plates were ruled out for subsequent data analysis and co-occurrence network reconstruction. With this cutoff value, we obtained 65 of 96-well plates (totally 6,091 wells) that were effective with microbial growth and were used for data analysis and co-occurrence network reconstruction. The biomasses were harvested from these plates and used for DNA extraction, library construction and high throughput 16S rRNA amplicon sequencing. After quality control and denoise, we obtained 136 Gbp sequence data. A total of 14,377 Zotus were annotated (Supplementary Table S2). There were 217 ± 94 (average ± standard deviation) prevalent Zotus, i.e., these Zotus appeared at frequencies ≥ 30% of wells in a given 96-well plate.

Next, we analyzed Zotus compositions and abundances in each well of the 65 plates. Accordingly, we reconstructed 65 independent microbial co-occurrence networks and further retrieved the robust (Spearman’s |ρ| >0.6 and P < 0.01) and prevalent Zotu pairs from these microbial networks. The Spearman rank correlation coefficients of Zotu pairs in the co-occurrence networks were all positive except Zotu11-Zotu12 pair from sediment sample (Supplementary Table S3), suggesting that Zotu pairs were most positively associated. Using Pajek (v5.11), we reconstructed co-occurrence networks with Zotu nodes of degree centrality < 5 (Fig. 2) and retrieved Zotu pairs. Table 1 shows the top 3 prevalent Zotu pairs from co-occurrence sub-networks of samples (plants, roots and sediments) and conditions (R2A or TSB medium). Using the 16S rRNA gene sequences, we tried to identify the phylogenetically related taxa of these Zotus, and the most closely related taxa at genus level are listed in Table 1.

Table 1

Top 3 prevalent Zotu pairs, their frequencies of occurrence networks and phylogenetically related taxa. The origins of those pairs (plants, roots, sediments) and the media used (R2A or TSB) for microfluidic 96-well plate cultivation are indicated.
Zotu Pairs
Zotu1-Zotu15673	0.62	Zotu1	Aeromonas	Roots/R2A
		Zotu15673	Aeromonas
Zotu18445-Zotu10	0.38	Zotu18445	Aeromonas
		Zotu10	Pseudomonas
Zotu8404-Zotu6	0.31	Zotu8404	Acinetobacter
		Zotu6	Acinetobacter
Zotu1-Zotu15673	0.81	Zotu1	Aeromonas	Roots/TSB
		Zotu15673	Aeromonas
Zotu1-Zotu15942	0.72	Zotu1	Aeromonas
		Zotu15942	Aeromonas
Zotu12259-Zotu1	0.36	Zotu12259	Aeromonas
		Zotu1	Aeromonas
Zotu1-Zotu15673	0.57	Zotu1	Aeromonas	Sediments/R2A
		Zotu15673	Aeromonas
Zotu1-Zotu15942	0.5	Zotu1	Aeromonas
		Zotu15942	Aeromonas
Zotu12707-Zotu7	0.42	Zotu12707	Aeromonas
		Zotu7	Enterobacter
Zotu5008-Zotu11151	0.78	Zotu5008	Aeromonas	Sediments/TSB
		Zotu11151	Aeromonas
Zotu1-Zotu15673	0.78	Zotu1	Aeromonas
		Zotu15673	Aeromonas
Zotu16-Zotu46	0.78	Zotu16	Citrobacter
		Zotu46	Methylobacter
Zotu8404-Zotu6	0.78	Zotu8404	Acinetobacter	Plants/R2A
		Zotu6	Acinetobacter
Zotu140-Zotu91	0.56	Zotu140	Azorhizobium
		Zotu91	Azorhizobium
Zotu8404-Zotu12231	0.44	Zotu8404	Acinetobacter
		Zotu12231	Aeromonas
Zotu1-Zotu15673	0.89	Zotu1	Aeromonas	Plants/TSB
		Zotu15673	Aeromonas
Zotu7-Zotu11295	0.78	Zotu7	Enterobacter
		Zotu11295	Aeromonas
Zotu1-Zotu15942	0.67	Zotu1	Aeromonas
		Zotu15942	Aeromonas

3.3 Cultivation of bacterial strains and matching to prevalent Zotu pairs

In order to experimentally demonstrate microbial interactions, it was essential to cultivate bacteria that corresponded to the Zotus from the co-occurrence networks. Based on the partial 16S rRNA genes (V4 regions) representing the Zotu pairs in Table 1, they were related to 8 bacterial genera (Aeromonas, Acinetobacter, Citrobacter, Methylobacter, Azorhizobium, Enterobacter, Pseudomonas, and an unidentified bacterial group). Referring to the Zotu abundances in plates, we located the wells and plates that showed high abundances of the Zotus in Table 1, and these wells were selected for bacterial isolation. We successfully obtained 129 bacterial strains (supplementary Table S4) and they were phylogenetically close to 15 bacterial species based on 16S rRNA gene identities, including Aeromonas caviae (3 isolates), Aeromonas hydrophila (7 isolates), Aeromonas media (5 isolates), Aeromonas rivipollensis (17 isolates), Elizabethkingia anopheles (27 isolates), Enterobacter ludwigii (8 isolates), Enterobacter soli (6 isolates), Klebsiella aerogenes (1 isolate), Microbacterium oxydans (2 isolates), Pantoea agglomerans (18 isolates), Pectobacterium aroidearum (2 isolates), Pleomorphomonas oryzae (1 isolate), Pleomorphomonas plecoglossicida (10 isolates), Pseudomonas protegens (21 isolates), and Raoultella ornithinolytica (1 isolate).

Next, we were trying to match the cultivated bacterial strains to the Zotus of the co-occurrence networks, and paid special attention to the Zotu pairs in Table 1. Based on the topology of the phylogenetic tree with the 16S rRNA genes of the bacterial isolates and the V4 regions of Zotus, we observed that 96 of the 129 bacterial isolates, representing 10 of the 15 bacterial species, matched 5 Zotus from co-occurrence networks. Specifically, four bacterial isolates, BOP-1, BOP-5, BOP-11 and BOP-16, members of the genus Aeromonas, matched either to Zotu1 or Zotu12259. As showed in Fig. 3, they all clustered into the Aeromonas lineage in the phylogenetic tree. Similarly, the isolates BOP-102 and BOP-108 matched Zotu10, and they were taxonomically annotated as members of the genus Pseudomonas and clustered into the Pseudomonas lineage in Fig. 3. Thirty-three of the 129 isolates, representing 5 of the 15 bacterial species, were not able to match any Zotus in Table 1. We also observed that Zotu15673, Zotu18445, Zotu8404, Zotu6, Zotu15942, Zotu12707, Zotu5008, Zotu11151, Zotu16, Zotu46, Zotu140, Zotu91, Zotu12231, Zotu11295 in Table 1 and the co-occurrence networks, respectively, did not match any cultivated bacterial isolates. However, we found that isolates BOP-61, BOP-73, BOP-74 and BOP-80 matched with Zotu7 and Zotu49 that was paired in the co-occurrence networks but not on the top 3 pairs list.

3.4 Bacterial associations on agar plates and in broth

Using the cultivated bacterial isolates and Zotu pairs, we tried to experimentally demonstrate the microbial interactions predicted from in silico co-occurrence networking analysis. Taking all cultivated bacterial isolates and the Zotu pairs from in silico analysis into consideration, we had 36 cultivated bacterial isolate combinations that represented 4 Zotu pairs in the co-occurrence networks. They were 6 cultivated bacterial isolate combinations for Zotu1-Zotu12259, 6 cultivated bacterial isolate combinations for Zotu7-Zotu49, 8 cultivated bacterial isolate combinations for Zotu1-Zotu10, and 16 cultivated bacterial isolate combinations for Zotu12259-Zotu49. All 36 cultivated bacterial isolate combinations were tested for their interactions on agar plates, and the results (Fig. 4) showed that their relations were neutral, mutualistic, or competitive. These results suggested that positively associated Zotu pairs from co-occurrence networks implied complicated bacterial associations.

In addition to the observations on TSB agar plates, the bacterial isolates BOP-1, BOP-5, BOP-11, and BOP-16 for Zotu1 and bacterial isolate BOP-108 for Zotu10 were cultivated in axenic or co-culture, and their individual growth (cell densities) were monitored with qPCR. Results showed that the axenic growths of bacterial isolates (Fig. 5A) were very different from co-culture. As showed in Fig. 5B, 5C and 5E, BOP-1, BOP-5, and BOP-16 for Zotu1 showed exponential growth in the first 12 h after inoculation, and then their cell densities decreased sharply. It was noteworthy that the BOP-108 for Zotu10 started growth right on 12 h after inoculation. This observation reminded that BOP-108 for Zotu10 was possibly competitive or even inhibitive to BOP-1, BOP-5, and BOP-16 for Zotu1. The isolate BOP-11 for Zotu1 showed very differently from axenic culture when co-cultivated with BOP-108 for Zotu10 (Fig. 5D). The growth of BOP-11 and BOP-108 apparently synchronized, suggesting that they were neutral or mutual to each other.

Co-occurrence network analyses based on metagenomic and amplicon sequences have been frequently used to infer microbial relations in environmental or host microbiomes, and further to predict microbial interactions (positive or negative) in these environments[8, 22–23, 38–39], which are informatic for discovering the true microbial interactions attributing to either microbial properties (such as metabolic coupling, nutrient competition, etc.) or from environmental filtration (such as responding similarly to environmental factors or due to physical separation)[40–41]. In the present study, we established an experimental workflow for demonstrating the interactions of Zotu pairs in co-occurrence networks with bacterial isolates. The workflow consisted of high-throughput microbial cultivation with droplet microfluidics, DNA bar-coding and sequencing, and matching microbial isolates to Zotu pairs and testing interactions of bacterial isolates on solid and broth media. The application of dilution cultivation during high-throughput microbial cultivation with droplet microfluidics significantly simplified the microbial community compositions, which is important for the identification of robust and reliable correlations in ecosystems with extremely high microbial diversity. This workflow also took the advances of droplet microfluidics, which enables high-throughput cultivation and enumeration in an array of nanoliter droplets[42]. This droplet microfluidics also holds the biggest promise to obtain less complex microbial consortia in microliter droplets via a high-throughput manner. Although biases of selected microbial growth might be accumulated due to differences in growth rates and to culture conditions (medium, temperature, aeration, etc.), thousands of droplets containing microbial sub-communities from a single sample allowed us to infer the co-occurrence networks using limited numbers of biological samples and to identify the robust associated Zotu pairs. In this study, we applied the workflow and demonstrated the Zotu pairs from co-occurrence network implied complicated relations in BOP wetland ecosystem. We observed neutral, competitive and mutual relations for positively associated Zotu pairs in co-occurrence network, suggesting that microbial interactions might be more complicated in real-world ecosystems. We observed in broth culture that the bacterial isolates interacted differently from that in agar plates and the interactions changed over time (e.g., Zotu1-Zotu10), which further complicated the interpretation of Zotu associations from co-occurrence networks.

It was also important that bar-coding of DNAs in this study from different plates for tracing prevalent Zotus pairs and subsequently guiding microbial isolation and cultivation. Co-occurrence networks can be inferred from the abundances of Zotus (usually based on 16S rRNA genes) [43–46] or microbial taxa (here refers to those based on metagenomes) or functional genes[30–31, 47]. During matching bacterial isolates to Zotus, we found that Zotus assignment to lower taxonomic levels using 16S rRNA genes from high-throughput amplicon sequencing (i.e., V4 or V3V4 regions of 16S rRNA genes) [43–46] was not ideal for the established workflow. One Zotu might be a sum of several bacterial strains sharing high identity of V4 regions[48], and this happened to the Zotu1 and Zotu49 that each matched four bacterial isolates (Fig. 3). Taking Zotu49 as an example, the different combinations for cultivated bacterial isolates showed inconsistent interactions on agar plates. The fact that one Zotus matched more than one bacterial isolates complicated experimental design and significantly increased workload of experimental demonstration. Due to the limited resolution of 16S rRNA genes for taxon assignments, it was hard to distinguish such observation of one Zotu matching more than one bacterial isolates was the true-fact in the BOP wetland ecosystem or resulted from the logistics of matching bacterial isolates to Zotus. Nevertheless, we considered whole genomic information from metagenomic data (such as MAGs) would increase the resolution for taxon assignments and would improve our workflow, and we are currently working on this idea with different microbial ecosystems.

A workflow for demonstration with cultivated bacterial strains of microbial associations from co-occurrence networks has been established and was applied for investigation of BOP wetland ecosystem. Results demonstrated that positively associated Zotus pairs displayed neutral, competitive and mutual relations, and suggested that the microbial interactions behind the Zotu pairs in the co-occurrence networks were complicated. This study supports the conception that experimental demonstration is indispensable for interpreting microbial association and provided an example of how to carry out such a demonstration. Limitations of the current study are discussed, and further improvement of the workflow is proposed.

Ethics approval and consent to participate: This study does not contain materials associated with human or animals.

Consent for publication: All authors have read and proved this submission.

Availability of data and material: All 16S rRNA sequence data are available as supplementary tables with this submission.

Competing interests: Authors declare not conflict of interest.

Funding: This study is supported by grants from National Key R&D Program of China (2019YFA0905500) and NSFC-EU Environmental Biotechnology joint program (NSFC grant No. 31861133002).

Authors' contributions: MZJ, HZZ, NZ and SJL performed sampling and the experiments. CL, CYJ and SJL conceived the study. MZJ, YLW and SJL analyzed the data wrote the manuscript.

Acknowledgements: Authors acknowledge the help from Dr. Xuejun Wu at Beijing Iloveearth Ecology Compony for his assistance of sampling.

Faust K, Raes J. Microbial interactions: from networks to models. Nature Reviews Microbiology. 2012; doi:10.1038/nrmicro2832.
Ings TC, Hawes J. The history of ecological networks. In: Dáttilo W, Rico-Gray V, editors. Ecological Networks in the Tropics. Springer; 2018. p. 15–28.
Chen YZ, Jian Li J, Zhang YF, Zhang MQ, Sun Z, Jing GC, et al. Parallel-Meta Suite: Interactive and rapid microbiome data analysis on multiple platforms. iMeta. 2022; doi: 10.1002/imt2.1.
Liu X, Shi Y, Yang T, Gao GF, Zhang L, Xu R, et al. Distinct Co-occurrence Relationships and Assembly Processes of Active Methane-Oxidizing Bacterial Communities Between Paddy and Natural Wetlands of Northeast China. Frontiers In Microbiology. 2022; doi:10.3389/fmicb.2022.809074.
Goberna M, Verdú M. Cautionary notes on the use of co-occurrence networks in soil ecology. Soil Biology and Biochemistry. 2022; doi: 10.1016/j.soilbio.2021.108534.
Barberan A, Bates ST, Casamayor EO, Fierer N. Using network analysis to explore co-occurrence patterns in soil microbial communities. ISME Journal. 2012; doi:10.1038/ismej.2011.119.
Zhou J, Deng Y, Luo F, He Z, Tu Q, Zhi X. Functional molecular ecological networks. mBio. 2010; doi:10.1128/mBio.00169-10.
Xun W, Liu Y, Li W, Ren Y, Xiong W, Xu Z, et al. Specialized metabolic functions of keystone taxa sustain soil microbiome stability. Microbiome. 2021; doi:10.1186/s40168-020-00985-9.
Liu F, Wang Z, Wu B, Bjerg JT, Hu W, Guo X, et al. Cable bacteria extend the impacts of elevated dissolved oxygen into anoxic sediments. ISME Journal. 2021; doi:10.1038/s41396-020-00869-8.
Gao P, Du G, Zhao D, Wei Q, Zhang X, Qu L, et al. Influences of Seasonal Monsoons on the Taxonomic Composition and Diversity of Bacterial Community in the Eastern Tropical Indian Ocean. Frontiers in Microbiology. 2020; doi:10.3389/fmicb.2020.615221.
Tsai KN, Lin SH, Liu WC, Wang DY. Inferring microbial interaction network from microbiome data using RMN algorithm. Bmc Systems Biology. 2015; doi:10.1186/s12918-015-0199-2.
Kurtz ZD, Muller CL, Miraldi ER, Littman DR, Blaser MJ, Bonneau RA. Sparse and Compositionally Robust Inference of Microbial Ecological Networks. Plos Computational Biology. 2015; doi:10.1371/journal.pcbi.1004226.
Steinway SN, Biggs MB, Loughran TP, Papin JA, Albert R. Inference of Network Dynamics and Metabolic Interactions in the Gut Microbiome. Plos Computational Biology. 2015; doi:10.1371/journal.pcbi.1004338.
Werner JJ, Knights D, Garcia ML, Scalfone NB, Smith S, Yarasheski K, et al. Bacterial community structures are unique and resilient in full-scale bioenergy systems. Proceedings of the National Academy of Sciences of the United States of America. 2011; doi:10.1073/pnas.1015676108.
Halfvarson J, Brislawn CJ, Lamendella R, Vazquez-Baeza Y, Walters WA, Bramer LM, et al. Dynamics of the human gut microbiome in inflammatory bowel disease. Nature Microbiology. 2017; doi:10.1038/nmicrobiol.2017.4.
Milici M, Deng ZL, Tomasch J, Decelle J, Wos-Oxley ML, Wang H, et al. Co-occurrence Analysis of Microbial Taxa in the Atlantic Ocean Reveals High Connectivity in the Free-Living Bacterioplankton. Frontiers in Microbiology. 2016; doi:10.3389/fmicb.2016.00649.
Woyke T, Teeling H, Ivanova NN, Huntemann M, Richter M, Gloeckner FO, et al. Symbiosis insights through metagenomic analysis of a microbial consortium. Nature. 2006; doi:10.1038/nature05192.
Ribeck N, Lenski RE. Modeling and quantifying frequency-dependent fitness in microbial populations with cross-feeding interactions. Evolution. 2015; doi:10.1111/evo.12645.
Wang X, Lu L, Zhou X, Tang X, Kuang L, Chen J, et al. Niche Differentiation of Comammox Nitrospira in the Mudflat and Reclaimed Agricultural Soils Along the North Branch of Yangtze River Estuary. Frontiers in Microbiology. 2020; doi:10.3389/fmicb.2020.618287.
Williams RJ, Howe A, Hofmockel KS. Demonstrating microbial co-occurrence pattern analyses within and between ecosystems. Frontiers in Microbiology. 2014; doi:10.3389/fmicb.2014.00358.
Frey-Klett P, Burlinson P, Deveau A, Barret M, Tarkka M, Sarniguet A. Bacterial-Fungal Interactions: Hyphens between Agricultural, Clinical, Environmental, and Food Microbiologists. Microbiology and Molecular Biology Reviews. 2011; doi:10.1128/Mmbr.00020-11.
Forster D, Qu Z, Pitsch G, Bruni EP, Kammerlander B, Proschold T, et al. Lake Ecosystem Robustness and Resilience Inferred from a Climate-Stressed Protistan Plankton Network. Microorganisms. 2021; doi:10.3390/microorganisms9030549.
Brandon-Mong GJ, Shaw GTW, Chen WH, Chen CC, Wang DR. A network approach to investigating the key microbes and stability of gut microbial communities in a mouse neuropathic pain model. Bmc Microbiology. 2020; doi:10.1186/s12866-020-01981-7.
Guimera R, Nunes Amaral LA. Functional cartography of complex metabolic networks. Nature. 2005; doi:10.1038/nature03288.
Guimera R, Amaral LA. Cartography of complex networks: modules and universal roles. Journal of statistical mechanics. 2005; doi:10.1088/1742-5468/2005/02/P02001.
Jiang CY, Dong LB, Zhao JK, Hu XF, Shen CH, Qiao YX, et al. High-Throughput Single-Cell Cultivation on Microfluidic Streak Plates. Applied and Environmental Microbiology. 2016; doi:10.1128/Aem.03588-15.
Zhang J, Liu YX, Guo X, Qin Y, Garrido-Oter R, Schulze-Lefert P, et al. High-throughput cultivation and identification of bacteria from the plant root microbiota. Nature Protocols. 2021; doi:10.1038/s41596-020-00444-7.
Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010; doi:10.1093/bioinformatics/btq461.
Reasoner DJ, Geldreich EE. A new medium for the enumeration and subculture of bacteria from potable water. Applied and Environmental Microbiology. 1985; doi:10.1128/aem.49.1.1-7.1985.
Zhang JY, Liu YX, Zhang N, Hu B, Jin T, Xu HR, et al. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nature Biotechnology. 2019; doi:10.1038/s41587-019-0104-4.
Bai Y, Muller DB, Srinivas G, Garrido-Oter R, Potthoff E, Rott M, et al. Functional overlap of the Arabidopsis leaf and root microbiota. Nature. 2015; doi:10.1038/nature16192.
Walters WA, Caporaso JG, Lauber CL, Berg-Lyons D, Fierer N, Knight R. PrimerProspector: de novo design and taxonomic analysis of barcoded polymerase chain reaction primers. Bioinformatics. 2011; doi:10.1093/bioinformatics/btr087.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Research. 2013; doi:10.1093/nar/gks1219.
Batagelj V, Mrvar A. Pajek - Analysis and Visualization of Large Networks. In: Jünger M, Mutzel P, editors. Graph Drawing Software. Springer; 2003. p. 77–103
Bastian M, Heymann S, Jacomy M. Gephi: an open source software for exploring and manipulating networks. In: International AAAI Conference on Web and Social Media. 2009. p. 361–362
Liu C, Du MX, Abuduaini R, Yu HY, Li DH, Wang YJ, et al. Enlightening the taxonomy darkness of human gut microbiomes with a cultured biobank. Microbiome. 2021; doi:10.1186/s40168-021-01064-3.
Green MR, Sambrook J. Constructing a Standard Curve for Real-Time Polymerase Chain Reaction (PCR) Experiments. Cold Spring Harbor Protocols. 2018; doi:10.1101/pdb.prot095026.
Wang L, Xu YN, Chu CC, Jing Z, Chen Y, Zhang J, et al. Facial Skin Microbiota-Mediated Host Response to Pollution Stress Revealed by Microbiome Networks of Individual. mSystems. 2021; doi:10.1128/mSystems.00319-21.
Szoboszlay M, Tebbe CC. Hidden heterogeneity and co-occurrence networks of soil prokaryotic communities revealed at the scale of individual soil aggregates. Microbiologyopen. 2021; doi:10.1002/mbo3.1144.
Levy R, Borenstein E. Metabolic modeling of species interaction in the human microbiome elucidates community-level assembly rules. Proceedings of the National Academy of Sciences of the United States of America. 2013; doi:10.1073/pnas.1300926110.
Liu R, Wang Z, Wang L, Li Z, Fang J, Wei X, et al. Bulk and Active Sediment Prokaryotic Communities in the Mariana and Mussau Trenches. Frontiers in Microbiology. 2020; doi:10.3389/fmicb.2020.01521.
Villa MM, Bloom RJ, Silverman JD, Durand HK, Jiang S, Wu A, et al. Interindividual Variation in Dietary Carbohydrate Metabolism by Gut Bacteria Revealed with Droplet Microfluidic Culture. mSystems. 2020; doi:10.1128/mSystems.00864-19.
Lane DJ, Pace B, Olsen GJ, Stahl DA, Sogin ML, Pace NR. Rapid-Determination of 16s Ribosomal-Rna Sequences for Phylogenetic Analyses. Proceedings of the National Academy of Sciences of the United States of America. 1985; doi:DOI 10.1073/pnas.82.20.6955.
Hamady M, Knight R. Microbial community profiling for human microbiome projects: Tools, techniques, and challenges. Genome Research. 2009; doi:10.1101/gr.085464.108.
Tringe SG, von Mering C, Kobayashi A, Salamov AA, Chen K, Chang HW, et al. Comparative metagenomics of microbial communities. Science. 2005; doi:10.1126/science.1107851.
Qin JJ, Li RQ, Raes J, Arumugam M, Burgdorf KS, Manichanh C, et al. A human gut microbial gene catalogue established by metagenomic sequencing. Nature. 2010; doi:10.1038/nature08821.
Zhalnina K, Louie KB, Hao Z, Mansoori N, da Rocha UN, Shi S, et al. Dynamic root exudate chemistry and microbial substrate preferences drive patterns in rhizosphere microbial community assembly. Nature Microbiology. 2018; doi:10.1038/s41564-018-0129-3.
Baxter NT, Schmidt AW, Venkataraman A, Kim KS, Waldron C, Schmidt TM. Dynamics of Human Gut Microbiota and Short-Chain Fatty Acids in Response to Dietary Interventions with Three Fermentable Fibers. mBio. 2019; doi:10.1128/mBio.02566-18.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

Operational taxon pairs in microbial co-occurrence networks imply complicated bacterial associations

Archived Versions:

Version 1

Abstract

Background

Results

Conclusion

Figures

1. Background

2. Methods And Materials

2.1 Sample collection and treatments

2.2 Microbial community inoculation and cultivation with microfluidics

2.3 DNA extraction, PCR amplification, bar-coding and sequencing

2.4 Data analysis, reconstruction of co-occurrence networks, and identification of Zotu pairs

2.5 Isolation of bacteria and matching the robust and prevalent Zotu pairs

2.6 Demonstration of microbial interactions with agar plate and liquid broth

2.7 qPCR for quantification of Pseudomonas and Aeromonas cells in liquid cultures

3 Result

3.1 Conception of the workflow to demonstrate the microbial associations from co-occurrence networks with microbial cultivation

3.2 Prevalent Zotu pairs in the co-occurrence networks

3.3 Cultivation of bacterial strains and matching to prevalent Zotu pairs

3.4 Bacterial associations on agar plates and in broth

4. Discussion

5. Conclusion

Declarations

References

Additional Declarations

Supplementary Files

Archived Versions:

Version 1