Droplet microfluidics-based high-throughput bacterial cultivation for validation of taxon pairs in microbial co-occurrence networks

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

Microbial co-occurrence networks are composed of nodes and edges, which usually represent microbes and statistically significant associations between microbes, respectively. We hypothesized that the microbial associations could be validated if the topological properties of networks are simplified, and if the microbes representing the nodes can be cultivated. To test this hypothesis, we designed a workflow as shown in Fig. 1. A total of 12,096 wells from 126 96-well plates were inoculated with droplets of series diluted environmental samples, wells from each 96-well plate represented the same combination of given culture condition, sample type (plants, roots, and sediments) and dilution rate (from 10^–1 to 10^–7). After being cultivated at 30 °C for 10 days, 69 effective (Supplementary Table S3) plates with > 30% wells showing microbial growth were retained for downstream microbial community analysis. Microbial DNA in each well was extracted, bar-coded, and sequenced for the inference of co-occurrence networks. The wells of plates showing high abundances of target Zotus were targeted for microbial isolations. Lastly, the cultivated microbial isolates were matched to Zotus in the network and used for demonstration of microbial interactions.

Prevalent Zotu pairs in the co-occurrence networks

Depending on the microbial density in samples, the 96-well plates harbored different numbers of wells with microbial growth. We obtained 65 96-well plates (6,091 wells) that were effective with microbial growth and data analysis for co-occurrence network reconstruction. After quality control and denoise, we obtained 130 Gbp sequence data. A total of 14,377 Zotus were annotated (Supplementary Table S4). 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. A total of 29,805 unique Zotu pairs were identified from the 65 co-occurrence networks (Fig. 2). The Jaccard similarity matrices estimated based on edge presence-absence patterns of these sub-networks revealed that networks constructed for the same medium and dilution level of a given sample were more similar to each other than those of the other media and dilution levels (Supplementary Figure S1). As expected, the edge numbers of sub-networks decreased with the increase in dilution levels (Supplementary Figure S2). 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 S5), suggesting that Zotu pairs were most positively associated. 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 the genus level are listed in Table 1. Comparison between the co-occurrence networks constructed with Spearman rank correlation (Supplementary Table S5) and FlashWeave (Supplementary Table S6) revealed that an average of 51.2% Zotu pairs across the networks were shared by these two different computational methods.

Cultivation of bacterial strains and matching to Zotu pairs

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 focused on 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 S7) 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. There are 32 isolated strains (i.e., BOP-1 to BOP-32) that shared identical 16S rRNA gene V4 region and clustered with Zotu 1 and Zotu12259 into the Aeromonas lineage in the phylogenetic tree (Fig. 3), we, therefore, randomly selected four bacterial isolates, BOP-1, BOP-5, BOP-11, and BOP-16, members of the genus Aeromonas, as representative isolates of either Zotu1 or Zotu12259. Similarly, BOP-61, BOP-73, BOP-74, and BOP-80 were selected from 34 bacterial strains (i.e., BOP-60 to BOP-94) sharing identical 16S rRNA gene V4 region as representative isolates of either Zotu7 or Zotu49. In addition, 20 (i.e., BOP-108 to BOP-128) and 8 bacterial strains (i.e., BOP-98 to BOP-106) that respectively shared identical 16S rRNA gene v4 region were taxonomically annotated as members of the genus Pseudomonas and clustered into the Pseudomonas lineage in Fig. 3. The isolates of BOP-102 and BOP-108 were selected as representative bacteria of Zotu10. 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 were paired in the co-occurrence networks but not on the list of top 3 prevalent Zotu pairs.

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 four 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. The selected Zotu pairs for experimental validation were confirmed by FlashWeave, except Zotu12259–Zotu49. All 36 cultivated bacterial isolate combinations were tested for their interactions on agar plates, and the results (Fig. 4) showed their mutualistic, competitive relations, or neutral (no interaction). 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 Zotu pair of bacteria from the genera Aeromonas (Zotu1) and Pseudomonas (Zotu10) were co-cultured in broth. Given that bacterial isolates BOP-1, BOP-5, BOP-11, and BOP-16 shared identical 16S rRNA gene V4 region with Zotu1, we only selected the BOP-108 that showed a close phylogenetic relationship with Zotu10 as the representative isolate (Supplementary Figure S3). These bacterial isolates were cultivated in axenic or co-culture, and their growths (cell densities) were monitored with qPCR. Results showed that the axenic growths of bacterial isolates (Fig. 5A) were very different from co-culture. As shown in Fig. 5B, C and E, 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 us 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.

Source: Ecology - nature.com