Conception of the workflow to demonstrate the microbial associations from co-occurrence networks with microbial cultivation
Microbial co-occurrence networks compose of nodes (OTUs) and edges (microbe-microbe associations). We hypothesized that the microbial associations could be demonstrated if the complex co-occurrence networks are simplified and composed of fewer centers and fewer edges, and if the microbes representing the centers be cultivated. To work on this hypothesis, we designed a workflow chart as shown in Figure 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 and microbial cultivation. Microbial DNAs in each plate were extracted, bar-coded, and sequenced for the inference 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.
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 considered those plates that 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 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 (Figure 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 the genus level are listed in Table 1.
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 shown in Figure 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 Figure 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.
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 (Figure 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 (Figure 5A) were very different from co-culture. As shown in Figures 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 (Figure 5D). The growth of BOP-11 and BOP-108 apparently synchronized, suggesting that they were neutral or mutual to each other.