Skip to main content
Download PDF

Core species and interactions prominent in fish-associated microbiome dynamics

Microbiome volume 11, Article number: 53 (2023) Cite this article

Abstract

Background

In aquatic ecosystems, the health and performance of fish depend greatly on the dynamics of microbial community structure in the background environment. Nonetheless, finding microbes with profound impacts on fish’s performance out of thousands of candidate species remains a major challenge.

Methods

We examined whether time-series analyses of microbial population dynamics could illuminate core components and structure of fish-associated microbiomes in the background (environmental) water. By targeting eel-aquaculture-tank microbiomes as model systems, we reconstructed the population dynamics of the 9605 bacterial and 303 archaeal species/strains across 128 days.

Results

Due to the remarkable increase/decrease of constituent microbial population densities, the taxonomic compositions of the microbiome changed drastically through time. We then found that some specific microbial taxa showed a positive relationship with eels’ activity levels even after excluding confounding effects of environmental parameters (pH and dissolved oxygen level) on population dynamics. In particular, a vitamin-B12-producing bacteria, Cetobacterium somerae, consistently showed strong positive associations with eels’ activity levels across the replicate time series of the five aquaculture tanks analyzed. Network theoretical and metabolic modeling analyses further suggested that the highlighted bacterium and some other closely-associated bacteria formed “core microbiomes” with potentially positive impacts on eels.

Conclusions

Overall, these results suggest that the integration of microbiology, ecological theory, and network science allows us to explore core species and interactions embedded within complex dynamics of fish-associated microbiomes. 

Video Abstract

Background

Microbial communities are essential factors of the life of vertebrates [1,2,3,4], playing key roles in the development and homeostasis of their hosts [5,6,7]. Gut microbiomes, for example, play key roles in the nutrition and disease prevention of human and other mammal species [8, 9]. Such physiological and ecological effects of gut microbes on hosts have been reported as well for fish [6, 10, 11]. Meanwhile, because fish are continuously exposed to numerous pathogenic and non-pathogenic microbial species in the water, their performance (or fitness) depends not only on gut-associated microbes [6, 10] but also on the microbiomes of the background environment [12,13,14]. Therefore, finding key microbiome components whose dynamics determine fish’s health or performance is of interdisciplinary interest spanning from microbiology to zoology and environmental science. However, due to the tremendous diversity of bacteria and archaea in aquatic ecosystems [15, 16], exploring such core microbial species associated with fish health remains a challenge.

A starting point for finding fish-health-associated microbes in aquatic ecosystems is to track the dynamics of microbial community compositions [17]. Nonetheless, we still have limited knowledge of the extent to which the structure of fish-associated microbiomes changes through time (but see [18]). Although time-series data of microbiomes have become available in large-scale projects of human-associated microbes [19, 20], few attempts have been made to monitor microbiomes associated with other animals over tens of time points. Moreover, continuous sampling of fecal or body-part samples of targeted vertebrate individuals is generally much harder in aquatic environments than in terrestrial environments. Thus, developing model systems for time-series analyses of microbe–fish ecological interactions is a demanding but essential step for exploring core bacteria and/or archaea out of thousands of candidate species in microbial communities.

Despite the hardship in gaining time-series microbiome samples at the individual level, fish-associated microbiome dynamics can be monitored at the population or community level by sampling environmental water samples [21,22,23]. Because excrements of fish are released to water, the samples of background water are expected to reflect the gut microbiomes of fish populations or communities. Furthermore, as individual fish are continuously exposed to the background microbiomes, analyses of water samples provide essential insights into surrounding environmental conditions, decomposition processes of non-ingested feed, and potential sources of gut microbiomes [12,13,14]. In this respect, time-series analyses of aquaculture or aquarium systems offer an ideal opportunity for investigating the relationship between microbial community structure, core microbial species, and vertebrate health.

By targeting a recirculating aquaculture system of the Japanese eel (Anguilla japonica), we herein integrate microbiology, community ecology, and network science for detecting key species and structures within fish-rearing-water microbiomes. Based on the DNA metabarcoding of prokaryote (bacterial and archaeal) communities for the 128-day time series, we revealed to what extent the compositions of aquaculture microbiomes fluctuate through time. We then reconstructed the population dynamics (i.e., increase/decrease) of the 9908 microbial amplicon sequencing variants (ASVs) constituting the aquaculture microbiomes, screening bacteria or archaea whose abundance was tightly linked with the health condition of eels. We then found that several microbial ASVs showed positive associations with eel health consistently across the five replicate aquaculture tanks examined, even after controlling the effects of these microbes’ environmental preference (e.g., preference to pH and dissolved oxygen level). With the approaches of network science and metabolic modeling, we further examined potential interactions between the core microbes, uncovering the structure of “core microbiomes” potentially determining fish health/performance. Overall, this study illustrates how core species and interactions are detected based on time series datasets of microbiome dynamics.

Results

Microbiome dynamics

Monitoring of microbiome dynamics was conducted by targeting the five water tanks of an aquaculture farm of the Japanese eel. In each water tank (diameter = 5 m; height = 1 m; volume = 20 m3), 1400–4300 eel individuals (average weight = 80–130 g) had been kept. The pH and dissolved oxygen (DO) concentrations were recorded for each tank every day. In addition, as a measure of ecosystem-level functions of microbiomes, the health condition of eels was evaluated based on eight criteria, yielding eel activity scores on a scale of 0 to 40 (see the “Methods” section). For the analyses of microbiome dynamics, water was sampled from each aquaculture tank every 24 h during 128 days. By applying a quantitative amplicon sequencing approach for estimating 16S ribosomal RNA gene (16S rRNA gene) copy concentrations of respective microbes [24], we obtained time-series datasets representing the increase/decrease of 9605 bacterial and 303 archaeal ASVs representing 618 genera and 325 families (Fig. 1a). Thus, our data offered a novel opportunity to test synchronizations among microbial population dynamics, environmental factors (pH and DO), and vertebrate performance (eel activity level).

Fig. 1
figure 1

Microbiome dynamics in the eel aquaculture system. a Dynamics of absolute abundance. For each water sample of each aquaculture tank, the absolute abundance of prokaryotes was inferred as 16S rRNA gene copy concentration based on the quantitative amplicon sequencing approach with standard DNA gradients. b Dynamics of relative abundance. The time series of the family-level taxonomic compositions are shown for each aquaculture tank. See Extended Data Figures 1–3 for phylum-, order-, and genus-level taxonomic compositions

Full size image

At the community level, drastic taxonomic turnover was observed in the time series of each aquaculture tank (Fig. 1; Additional file 1: Fig. S1, Additional file 2: Fig. S2, Additional file 3: Fig. S3). In Tanks 1 and 2, for example, the community structure characterized by the predominance of Fusobacteriaceae and Microbacteriaceae was suddenly altered by a Flavobacteriaceae-dominated state around day 45 (Fig. 1). Meanwhile, microbiomes of Tanks 3–5 displayed more complex dynamics represented by frequent shifts between Flavobacteriaceae-dominated and Chitinophagaceae-dominated states, although a clear classification of community states was difficult (Fig. 1).

A multivariate analysis of the prokaryote community structure suggested that the community state characterized by the dominance of Fusobacteriaceae and Microbacteriaceae was potentially associated with high eels’ activity (Fig. 2). In contrast, the Flavobacteriaceae-dominated and Chitinophagaceae-dominated states, which were observed in high-pH conditions, were associated with low eels’ activity (Fig. 2). At the genus level, the high-eel-activity-related state of dominance by Fusobacteriaceae and Microbacteriaceae was characterized by the high relative abundance of Cetobacterium, which includes species potentially contribute to fish physiological homeostasis [25]. On the other hand, the Flavobacteriaceae-dominated and Chitinophagaceae-dominated states associated with low eels’ activity were represented by Flavobacterium and Edaphobaculum, respectively (Additional file 4: Fig. S4). Among the genera, Flavobacterium includes fish pathogens [26], while Edaphobaculum [27] has been poorly investigated in terms of their effects on fish physiology. These results suggest the potential impacts of environmental microbiome dynamics on fish health/behavior in aquaculture systems.

Fig. 2
figure 2

Multivariate analysis of community structure. a Community state space. Community compositions of the samples are plotted on the two-dimensional surface defined with non-metric multidimensional scaling (NMDS). The NMDS was performed based on the Bray-Curtis β-diversity of family-level taxonomic compositions. The projections of the data points onto the vectors have a maximum correlation with the variables examined (pH, DO, and eels’ activity level). See Extended Data Figure 4 for an additional analysis based on genus-level taxonomic compositions. b Examples of community structure in the NMDS surface. For several points within the NMDS surface (panel a), family-level taxonomic compositions are shown. The example points are ordered along the vector representing high eels’ activity level

Full size image

Exploring microbes with key roles

To reveal the relationship between microbiome structure and eels’ activity level, we next evaluated how the population dynamics of each microbial ASV were associated with environmental variables and eels’ activity level. Specifically, we examined how the population size (absolute abundance) of each ASV varied with pH, DO, and eels’ activity level (Fig. 3a) based on correlation analyses with twin-surrogate permutations [28] (Fig. 3b-c). The ASVs varied in their environmental preference for pH and DO conditions as well as in their associations with eels’ activity levels (Additional file 5: Fig. S5). We also found that ASVs’ relationship with eels’ activity level displayed tank-dependent complex associations with pH or DO preference (Fig. 3d–e). Thus, for each microbial ASV in each aquaculture tank, we calculated a partial correlation between absolute abundance and eels’ activity scores through the time series by controlling the effects of pH or DO. Because partial correlation coefficients were consistent between the pH-controlled and DO-controlled calculations (Fig. 3f), the pH-controlled partial correlation coefficients were used in the following analyses.

Fig. 3
figure 3

Microbes associated with eels’ activity. a Time series of pH, dissolved oxygen (DO) level, and eels’ activity score are shown for each aquaculture tank. b Example of the correlation analysis. For each variable shown in panel b, Spearman’s correlation with the absolute abundance of each ASV in each aquaculture tank was examined. c Randomization analysis of correlation. The significance of correlation coefficients was examined based on a twin-surrogate randomization analysis of time-series data (100,000 permutations). Coefficients less than −0.3 and those larger than 0.3 roughly represent significant negative and positive correlations, respectively. d Each ASV’s correlation with pH and eels’ activity level. e Each ASV’s correlation with DO and eels’ activity level. f Partial correlation with eels’ activity level. To control the effects of pH or DO, a partial correlation between absolute abundance and eels’ activity scores was calculated for each ASV in each tank. g Taxonomic comparison of relationship with eels’ activity level. A partial correlation with eels’ activity level is shown for the genera that appeared in all the aquaculture tanks (shown in the decreasing order of mean values). h Time-lag analysis of correlations. In calculating the partial correlation between eels’ activity level and the absolute abundance of the Cetobacterium ASV (X_0002), a defined time-lag was introduced to the eels’ activity variable

Full size image

The partial correlation coefficients with eels’ activity levels varied greatly depending on prokaryote taxa (Fig. 2e). Nonetheless, ASVs belonging to some bacterial genera showed a consistently positive correlation with eels’ activity scores across the five tanks (Fig. 2g; Additional file 5: Fig. S5). The list of those ASVs included bacteria belonging to the genera Cetobacterium (Fusobacteriaceae; Fusobacteriia; ASV ID = X_0002), Plesiomonas (Enterobacteriaceae; Gammaproteobacteria; X_0020), Turicibacter (Erysipelotrichaceae; Bacilli; X_0041), Paraclostridium (Clostridiaceae; Clostridia; X_0014), Romboutsia (Peptostreptococcaceae; Clostridia; X_0028), Edwardsiella (Hafniaceae; Gammaproteobacteria; X_0027), Clostridium (Clostridiaceae; Clostridia; X_0029), and an ASV belonging to Barnesiellaceae (Bacteroidia; X_0064) (Fig. 3g; Additional file 5: Fig. S5d).

An additional database search of the 16S rRNA sequences suggested that some of the ASVs with positive associations with eels’ activity levels belonged to bacterial species with potential physiological impacts on fish. For example, the Cetobacterium ASV, which showed the strongest positive partial correlation with eels’ activity level, was represented by the 16S rRNA sequences completely matching that of Cetobacterium somerae (formerly recognized as “Bacteroides type A”) in the NCBI nucleotide database. This Cetobacterium species has been known to produce high concentrations of vitamin B12, and hence, their potential contributions to fish’s physiology have been anticipated [25]. Meanwhile, the Edwardsiella ASV listed above was allied to the notorious fish pathogen E. tarda [29], illuminating paradoxical relationships with eels’ health. However, our supplementary phylogenetic analysis based on the sodB gene marker [30] indicated that 95.1 % of Edwardsiella bacteria detected in the focal eel aquaculture system belonged to non-pathogenic clades [30, 31] within the genus Edwardsiella (Additional file 6: Fig. S6).

In terms of negative impacts on eels’ activity level, bacteria in the genera Aeromonas (Aeromonadaceae; Gammaproteobacteria), Methylobacterium (alternatively, Methylorubrum; Beijerinckiaceae; Alphaproteobacteria), and Acinetobactor (Moraxellaceae; Gammaproteobacteria) were highlighted (Fig. 3g). Among them, Aeromonas and Acinetobactor have been known to include fish pathogens [32, 33]. At the ASV level, an ASV allied to the cvE6 clade within the order Chlamydiales (Chlamydiae; Verrucomicrobiota) showed the strongest negative correlation with eels’ activity scores (Additional file 5: Fig. S5d).

Although the above analysis controlling environmental preferences of respective bacteria allows high-throughput screening for species with potential positive/negative impacts on target biological functions, the simple statistical approach with partial correlation analyses precludes insights into the direction of causality. Specifically, it is important to consider the possibility that the high/low abundance of an ASV is a consequence but not a cause of eels’ high/low activity. Therefore, we performed an additional analysis introducing time lags into eels’ activity scores throughout the time series. We then found that the abundance of the Cetobacterium ASV was positively correlated with eels’ activity scores of the next day, while correlations between Cetobacterium abundance and past eels’ activity scores were much lower than those with no time lags (Fig. 3h). Meanwhile, the high correlation between 5-day-ago eels’ activity level and present-day Cetobacterium abundance was observed in some tanks (Tanks 1 and 4; Fig. 3h), illuminating the importance of carefully interpreting the results of the time series analysis.

Networks of interactions

We then reconstructed webs of potential microbe-to-microbe interactions to illuminate microbial groups or interactions positively associated with eels’ health. We first applied the Meinshausen-Bühlmann (MB) method, which was designed to evaluate patterns of coexistence realized by the effects of microbe–microbe interactions as well as those of niche sharing between microbes. For each aquaculture tank, the reconstructed network of microbe–microbe coexistence (Additional file 7: Fig. S7, Additional file 8: Fig. S8; Additional file 9: Table S1) was compartmentalized into several modules, which differed in mean partial correlations with eels’ activity scores (Fig. 4). We then found that each of the five networks included a module constituted by the abovementioned Cetobacterium ASV and several other ASVs with consistently positive associations with eels’ activity level (Fig. 4; Additional file 10: Fig. S9). The bacteria consistently formed network modules of coexistence with the Cetobacterium ASV were Plesiomonas (X_0020), Turicibacter (X_0041), Paraclostridium (X_0014), Romboutsia (X_0028), Edwardsiella (X_0027), and Clostridium (X_0029) (Fig. 4).

Fig. 4
figure 4

Microbe-to-microbe coexistence networks. For each aquaculture tank, patterns of coexistence were analyzed based on the sparse inverse covariance estimation for ecological associations with the Meinshausen-Bühlmann (MB) model. Only the ASVs that appeared in 30 or more samples were targeted in the analysis of each tank. Within the networks, pairs of microbial ASVs that may interact with each other in facilitative ways and/or those potentially sharing environmental preferences are linked with each other. Network modules, which represent groups of densely linked ASVs, are shown for each network. The color of nodes indicates a partial correlation between ASV abundance and eels’ activity level (controlled variable = pH). The inferred network modules are shown by colors for each tank in a box. The ASVs that consistently displayed positive or negative correlation with eels’ activity level (Additional file 5: Fig. S5) are highlighted with the defined symbols. See Additional files 7–9 for additional information of the nodes (ASVs) and modules within the network. ASVs included in minor sub-networks (number of nodes < 5) are not shown

Full size image

To infer the presence/absence of direct interactions between these bacteria with a positive relationship with eels’ activity level, we conducted an additional network analysis based on the sparse and low-rank (SLR) decomposition method, which allowed us to remove latent effects of environmental conditions. In the networks reconstructed with the SLR method (Fig. 5), the potential effects of niche sharing were controlled and hence the links between bacterial ASVs were expected to represent potential positive interactions. The estimated interaction coefficients were highly correlated between the MB and SLR methods (Additional file 11: Fig. S10). Meanwhile, in the SLR-based network, removing the effects of potential niche sharing (sharing of environmental preference) resulted in the simplification of network structure, in which estimated direct interactions between microbes were focused (Fig. 5). Despite the considerable difference between MB- and SLR-based network topology, the Cetobacterium ASV with the strongest associations with eels’ activity level was, again, linked with the Plesiomonas, Turicibacter, Paraclostridium, Romboutsia, Edwardsiella, and Clostridium ASVs within the SLR network (Fig. 5), suggesting positive interactions with these bacteria.

Fig. 5
figure 5

Inferred direct interactions between microbes. Based on the “sparse and low-rank” (SLR) model, direct interactions between microbial ASVs were inferred by controlling the effects of shared environmental preference. Only the ASVs that appeared in 30 or more samples were targeted in the analysis of each tank. The links between nodes represent potentially positive interactions between ASVs. The color of nodes indicates a partial correlation between ASV abundance and eels’ activity level (controlled variable = pH). ASVs included in minor sub-networks (number of nodes < 5) are not shown

Full size image

Potential metabolic interactions

To estimate functional interactions between microbes, we focused on the genomic compositions of respective microbes within the aquaculture microbiomes. After retrieving the information of genomic compositions from reference databases, we analyzed the inferred gene repertoires (KEGG metabolic pathway/process profiles) of the microbial ASVs based on multivariate analysis. Along the principal component axes, Cetobacterium, which showed consistent correlations with eels’ activity (Fig. 3g; Additional file 5: Fig. S5d), was located distantly from Edwardsiella, Plesiomonas, and Turicibacter (Fig. 6a). In contrast, Romboutsia, Paraclostridium, and Clostridium displayed similar metabolic gene repertories with Cetobacterium (Fig. 6a).

Fig. 6
figure 6

Metabolic interactions between microbes. a Metagenomic niche space. Microbial ASVs are plotted on a two-dimensional surface of PCoA based on their KEGG metabolic pathway/process profiles inferred with a phylogenetic prediction of genomes. Microbial ASVs plotted closely within the surface are expected to have similar gene repertoires. The ASVs highlighted in Figs. 4 and 5 are shown with large symbols. b Potential competitive and facilitative interactions. Based on the NCBI RefSeq genome information, potential metabolic interactions between each pair of ASVs were inferred in terms of metabolic resource overlap (MRO) and metabolic interaction potential (MIP). Histograms of MRO and MIP are shown on the horizontal and vertical axes, respectively. ASV pairs including the Cetobacterium ASV, whose abundance was positively associated with eels’ activity level in all the five water tanks (Fig. 3g; Additional file 5: Fig. S5), are shown in pink. Relationships between the Cetobacterium ASV and the ASVs highlighted in Figs. 4 and 5 are indicated as well

Full size image

We next evaluated potential competitive and facilitative interactions between microbes based on a genome-scale metabolic modeling approach. In the analysis, reference genomic information was used to infer competition for available resources and exchanges of metabolites, yielding metabolic resource overlap and metabolic interaction potential scores for each pair of microbial ASVs. We then found that the Romboutsia, Edwardsiella, and Plesiomonas ASVs had relatively low metabolic resource overlap and relatively high metabolic interaction potential with the Cetobacterium ASV among the prokaryotes examined (Fig. 6b).

Discussion

Through the 128-day monitoring of thousands of microbial species/strains, we here found that aquatic microbiomes associated with fish could show drastic shifts of community structure through time. Such the dynamical nature of community processes has been intensively investigated in human-associated microbiomes in light of the potential influence on host status [19, 20]. In particular, shifts (collapse) of microbial community structure to disease-related states (i.e., dysbiosis) have been considered as essential mechanisms determining human health [34, 35]. Given the growing literature on microbiome dynamics in medical science, knowledge of shifts between alternative states of fish-related microbiomes [14] is expected to shed new light on the physiological and ecological processes of vertebrates.

The aquaculture microbiome dynamics were described as shifts among Fusobacteriaceae-abundant states, Flavobacteriaceae-dominated states, and Chitinophagaceae-dominated states, although intermediate states existed through the time series (Fig. 1). Among them, Fusobacteriaceae-abundant states, which were characterized by high abundance of Cetobacterium, were designated as microbiome compositions positively associated with eels’ activity level (Fig. 2; Additional file 4: Fig. S4). In fact, among the 9,908 microbial ASVs examined, the ASV representing Cetobacterium somerae showed the strongest associations with eel’s activity level through the time series even after controlling the effects of environmental preference (Fig. 3; Additional file 5: Fig. S5). This Cetobacterium species has been reported from a broad taxonomic range of freshwater fish [36,37,38,39], especially from the intestines of species that do not require dietary vitamin B12 [25]. Although vitamin B12 (cobalamin) plays essential roles in animal physiology (e.g., normal functioning of nervous systems and the maturation of red blood cells), they can be synthesized only by specific clades of bacteria and archaea [40, 41]. Genomic studies have shown that C. somerae has a series of genes for anaerobic vitamin B12 biosynthesis [42]. Indeed, the bacterium produces the highest concentrations of vitamin B12 compared to other culturable bacteria within freshwater fish-associated microbiomes [25, 43]. Given the prevalence of Cetobacterium in freshwater fish species [36,37,38], our results suggest that maintaining microbiomes at Cetobacterium-abundant states is the key to build general platforms for stably keeping freshwater aquaculture/aquarium systems.

Further analyses based on network theory and metabolic modeling indicated the possibility that the Cetobacterium species form facilitative interactions with some other microbial species/ASVs (Figs. 4, 5 and 6). Among the bacteria for which interactions with Cetobacterium were inferred from multiple analyses, Edwardsiella tarda has been known to include notorious pathogens of broad taxonomic groups of fish including eels [29, 30, 44]. However, we found that the E. tarda population of the investigated aquaculture system was dominated by non-pathogenic strains [30, 31] of the species (Additional file 6: Fig. S6). Thus, the presence of microbial species/strains belonging to well-known taxa of pathogens do not necessarily result in negative impacts on fish. Rather, our analyses suggested that “seemingly pathogenic” microbes could be involved in core microbiome components (network modules) constituted by microbes contributing to the maintenance of fish health. Further studies are awaited to explore potential mechanisms such as the competitive exclusion of pathogenic strains by non-pathogenic strains [45, 46] or indirect negative impacts on pathogenic strains through the activation of fish immune systems [6, 47] by non-pathogenic strains. In contrast to E. tarda, Romboutsia, and Plesiomonas, which were inferred as microbes with facilitative interactions with C. somerae, too (Figs. 4, 5 and 6), have been poorly investigated in terms of their functions. Their potential roles in the competitive exclusion of pathogens or activation of host immune systems deserve further investigations.

While the time series dataset allowed us to highlight core species and interactions within microbial communities, more sophisticated statistical platforms beyond simple correlational approaches are necessary for confirming causative relationships between microbiome dynamics and vertebrate health/performance. In this respect, methods based on information theory or nonlinear mechanics, such as transfer entropy [48] and empirical dynamic modeling [49, 50], are expected to help us infer causative interactions among microbial population dynamics, environmental factors, and vertebrate performance. Albeit promising, these methods require substantial computational resources when we try to analyze microbiomes consisting of thousands of ASVs. Further methodological advances will deepen our understanding of the mechanisms by which microbiome dynamics and vertebrate performance are linked with each other.

Conclusions

We showed that the structure of fish-associated microbiomes could drastically change through time as has been reported in studies on dysbiosis of human gut microbiomes. As analyses of microbiome dynamics are extended from medical science to researches targeting other vertebrates, we will be more and more aware of overlooked roles of microbes in both terrestrial and aquatic ecosystems. Feedback between intestine and environmental microbiomes, for example, deserves future intensive research in terms of potential great impacts on animal population/community dynamics. In particular, given that aquatic vertebrates are continuously exposed to the excrements of other individuals or species, their gut microbiome dynamics (and related health conditions) may be more likely to be synchronized at the population or community levels than those of terrestrial vertebrates. Therefore, simultaneous monitoring of the intestine and background environmental microbiomes will provide platforms for uncovering such feedback and synchronization processes. Further insights into fish-associated microbiome dynamics will reorganize our basic understanding of aquatic ecosystem dynamics, advancing technologies for sustainable food production through stable aquaculture systems [17, 51,52,53].

Methods

Sampling

Monitoring of microbiome dynamics was conducted by targeting the five water tanks of the indoor eel-aquaculture system of A-Zero Inc. (Nishiawakura, Okayama Prefecture, Japan). In each water tank (diameter = 5 m; height = 1 m; volume = 20 m3), 1400–4300 eel individuals (average weight = 80–130 g) had been kept. In the eel aquaculture system, the introduction of purchased young fish and shipment of a fraction of grown fish occur irregularly. The benefit of targeting such a commercial aquaculture system rather than fully controlled experimental systems is that we can analyze complex temporal dynamics of real aquaculture systems and thereby explore diverse microbial species potentially having great impacts on fish health. The eels were fed mainly with fishmeal, cereals, pollack oil, and krill. The amount of the feed had been manually optimized for each tank everyday by an expert keeper. The rearing water had been kept stirred with the rigorous swimming of eels and water flow from the filtration system, and the water temperature in the tanks was kept at around 30 °C. About 10 % of tank water was replaced with warmed fresh well water every day. The drainage from the five tanks was mixed and processed in a series of filtration equipment. The filtered drainage was returned to each tank after being processed in another filtration equipment adjacent to each tank. The eels were fed with a mixture of commercial artificial diets. The pH, dissolved oxygen (DO), and eels’ activity levels were recorded for each tank every day. The eels’ activity level was evaluated based on the sum of the scores of the following eight criteria: initial responses to feeders, the proportion of eels responding to feeders, sharpness of movement, the proportion of eels eating the artificial diet, the level of splashes, the amount of scattered diet, the time to consume the diet, and the proportion of foraging eels at the end of feeding. For each of the criteria, scoring was done on a five-point scale (maximum point = 5) by an expert of eel aquaculture maintenance: thus, 40 (5 × 8 criteria) is the maximum point of the eels’ activity score. Albeit subjective, the criteria evaluated continuously by a professional provide inferences of eel’s health conditions throughout the time series. The water in the tanks was continuously mixed by the movement of eels.

From each aquaculture tank (Tanks 1–5), ca. 1.5 mL of water was sampled in the morning every day during the 128 days from March 25 to July 30, 2020, except for 9 days (days 102, 103, 120, 121, 122, 123, 124, 125, and 126), i.e., the samples of 119 days were available. In Tank 4, the samples were unavailable on additional 3 days (days 67–69) due to the cleaning and the entire replacement of water. Note that such maintenance events do not necessarily result in the complete replacement of rearing water microbiomes because the five tanks in the aquaculture system share filtration equipment as mentioned above. Consequently, the number of collected samples was 592 (119 days × 5 tanks–3 days in Tank 4). Water sample was collected in a 2.0-mL microtube, and they were immediately stored at −20 °C in a freezer until DNA extraction.

Quantitative 16S rRNA sequencing

To extract DNA from each sample, 250 μL of the collected water was mixed with 400 μL lysis buffer (0.0025 % SDS, 20 mM Tris (pH 8.0), 2.5 mM EDTA, and 0.4 M NaCl) and 250 μL 0.5 mm zirconium beads in a 2.0 mL microtube. The microtubes were then shaken at 25 Hz for 5 min using TissueLyser II (Qiagen, Venlo). After centrifugation, the aliquot was mixed with proteinase K solution (×1/100 of the total volume) within a sterilized laminar-flow cabinet, being incubated at 40 °C for 60 min followed by 95 °C for 5 min.

We then performed PCR by applying a quantitative amplicon sequencing method [24, 54]. Although most existing microbiome studies were designed to infer a “relative” abundance of microbial amplicon sequence variants (ASVs) or operational taxonomic units (OTUs), information of “absolute” abundance provides additional insights into microbiome dynamics, i.e., insights into increase/decrease of the population size of each prokaryote ASV/OTU within a microbiome throughout a time series [24]. The quantitative amplicon sequencing approach is based on the addition of artificial (standard) DNA sequences with defined concentrations into PCR master solutions. Therefore, even if compositions or concentrations of PCR inhibitor molecules in DNA extracts vary among time-series samples, potential bias caused by such inhibitors can be corrected based on the use of the internal standards (i.e., standard DNAs within PCR master solutions).

The V4 hypervariable region of the prokaryote 16S rRNA gene was PCR-amplified with the forward primer 515f [55] fused with 3–6-mer Ns for improved Illumina sequencing quality and the forward Illumina sequencing primer (5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG- [3–6-mer Ns] – [515f]-3′) and the reverse primer 806rB [56] fused with 3–6-mer Ns for improved Illumina sequencing quality [57] and the reverse sequencing primer (5′- GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA G [3–6-mer Ns] - [806rB] -3′) (0.4 μM each). To apply the quantitative amplicon sequencing, five standard DNA sequence variants with different concentrations of artificial 16S rRNA sequences (0.1, 0.05, 0.02, 0.01, and 0.005 nM) were added to the PCR master mix solution [24], which was based on the buffer and polymerase system of KOD One (Toyobo). To prevent contamination of external DNA, the PCR master mix solution was automatically dispensed with Mantis liquid dispenser (Formulatrix) in a sterilized laminar-flow cabinet. The temperature profile of the PCR was set as to 35 cycles at 98 °C for 10 s, 55 °C for 30 s, and 68 °C for 30 s: the ramp rate through the thermal cycles was set to 1 °C/s to prevent the generation of chimeric sequences [58]. Illumina sequencing adaptors were then added to respective samples in the supplemental PCR using the forward fusion primers consisting of the P5 Illumina adaptor, 8-mer indexes for sample identification [59], and a partial sequence of the sequencing primer (5′-AAT GAT ACG GCG ACC ACC GAG ATC TAC AC - [8-mer index] - TCG TCG GCA GCG TC-3′) and the reverse fusion primers consisting of the P7 adaptor, 8-mer indexes, and a partial sequence of the sequencing primer (5′-CAA GCA GAA GAC GGC ATA CGA GAT - [8-mer index] - GTC TCG TGG GCT CGG-3′). KOD One was used with a temperature profile: followed by 8 cycles at 98 °C for 10 s, 55 °C for 5 s, and 68 °C for 30 s (ramp rate = 1°C/s). PCR-negative control samples were introduced into the above reactions to check the effectiveness of the sterilized experimental environment.

The PCR amplicons of the samples were then pooled after a purification/equalization process with the AMPureXP Kit (Beckman Coulter). Primer dimers, which were shorter than 200 bp, were removed from the pooled library by supplemental purification with AMpureXP: the ratio of AMPureXP reagent to the pooled library was set to 0.6 (v/v) in this process. Because the quality of forward sequences is generally higher than that of reverse sequences in Illumina sequencing, we optimized the MiSeq run setting in order to use only forward sequences. Specifically, the run length was set at 271 forward (R1) and 31 reverse (R4) cycles to enhance forward sequencing data: the reverse sequences were used only for screening 16S rRNA sequences in the following bioinformatic pipeline.

Bioinformatics

In total, 16,298,203 sequencing reads were obtained in the Illumina sequencing. The raw sequencing data were converted into FASTQ files using the program bcl2fastq 1.8.4 distributed by Illumina. The output FASTQ files were demultiplexed using Claident v0.2. 2018.05.29 [60]. The removal of low-quality sequences and ASV inferences was done using DADA2 [61] v.1.22.0 of R 4.1.2 [62]. Singletons and possibly chimeric sequences were excluded from the dataset within the DADA2 pipeline. The taxonomy of the output ASVs was inferred based on the naive Bayesian classifier method [63] using the SILVA v.138 database [64]. The number of ASVs recached plateaus with an increasing number of sequencing reads in each sample (Additional file 11: Fig. S10). On average, 15,308 reads (SD = 8990) were obtained per sample. Based on the calibration with the concentration gradients of the five standard DNAs [24, 50], concentrations of respective ASVs were obtained for each sample (16S rRNA gene copy numbers per unit volume of tank water samples; copies/μL). As the number of 16S rRNA gene copies per genome generally varies among prokaryotic taxa [65], 16S rRNA gene copy concentration is not directly the optimal proxy of cell or biomass concentration. Meanwhile, in this study, estimates of 16S rRNA gene copy concentrations were used to observe increase/decrease of abundance (i.e., population dynamics) within the time series of respective microbial ASVs. Thus, variation in the number 16S rRNA gene copy numbers among microbial taxa had no qualitative effects on the subsequent population- and community-ecological analyses. The samples in which Pearson’s coefficients of correlations between sequencing read numbers and standard DNA copy numbers (i.e., correlation coefficients representing calibration curves) were less than 0.8 were removed as those with unreliable estimates. Samples with less than 1000 reads were discarded as well. In total, microbiome data were successfully obtained from 577 out of 592 samples. For each aquaculture tank, we then obtained a sample × ASV matrix, in which a cell entry depicted the concentration of 16S rRNA gene copies of an ASV in a sample.

Community structure

For each aquaculture tank, Bray-Curtis β-diversity was calculated for all pairs of time points based on the matrix describing the relative abundance of prokaryote families using the vegan 2.6.2 package [66] of R. Based on the β-diversity estimates, the community structure of all the samples across the five water tanks were visualized on the surface of non-metric multidimensional scaling (NMDS). The vectors representing the environmental variables (pH and DO) and eels’ activity level were calculated with the “envfit” function of R, and they were shown on the NMDS surface. The analysis was conducted as well based on the matrix describing the relative abundance of genera.

Within the dataset, we explored the presence of bacterial ASVs belonging to potentially fish-pathogenic taxa. The ASVs belonging to the genera including notorious fish pathogens were screened in light of review papers of fish pathogenic prokaryotes [67,68,69,70]. The total absolute abundance (16S rRNA gene copy concentrations) of bacteria belonging to the genera Aeromonas, Edwardsiella, Flavobacterium, Mycobacterium, Pseudomonas, and Renibacterium was then displayed through the time series to roughly inspect the microbiome dynamics (Additional file 13: Fig. S12).

Environmental preference of ASVs

To evaluate the environmental preference of each microbial ASV, Spearman’s correlation between absolute abundance (in the metric of DNA copy numbers of the 16S rRNA gene) and pH was calculated. For each tank, the ASVs that appeared in 30 or more samples were subjected to the analysis. For each ASV in each water tank, the statistical significance of the obtained correlation coefficient was examined with a randomization analysis obtained based on a twin-surrogate method for time-series data [28] (100,000 permutations). Correlation coefficients less than −0.3 and those larger than 0.3 tended to show statistically significant negative and positive correlations with pH, respectively, after Benjamini-Hochberg adjustment of P values in multiple testing [i.e., false discovery rate (FDR)]. Likewise, Pearson’s correlation coefficients between respective ASVs’ absolute abundance and DO concentrations were calculated.

ASV abundance and eel’s activity

We explored microbial ASVs that potentially have profound impacts on eels’ health. For each water tank, Spearman’s correlation between absolute abundance and eels’ activity score was calculated for the ASVs that appeared in 30 or more samples. However, because ASV abundance could be affected by pH or DO concentration, the use of such simple correlation coefficients might be misleading. Therefore, we controlled potential effects by environmental factors/conditions based on a partial correlation approach as follows:

$${r}_{xy\bullet z}=\frac{r_{xy}-{r}_{xz}{r}_{yz}}{\sqrt{1-{r}_{xz}^2}\ \sqrt{1-{r}_{yz}^2}},$$

where rxy, rxz, and ryz were a correlation between ASV abundance and eels’ activity level, that between ASV abundance and an environmental factor (pH or dissolved oxygen concentration), and that between eels’ activity level and an environmental factor, respectively. For each ASV, a randomization analysis was performed with the twin-surrogate method (100,000 permutations).

Time-lag analysis

We extended the analysis of the partial correlation between microbial abundance and eels’ activity level by introducing time lags between the two variables. Specifically, a partial correlation between an ASV’s abundance on Day x and eels’ activity score on Day x + l was calculated. The time lag l ranged from −5 to 5 in the analysis (l = 0 means no delay introduced to eels’ activity level).

Pathogenic and non-pathogenic Edwardsiella

We performed an additional analysis to infer the proportion of opportunistically pathogenic and potentially non-pathogenic clades [30, 31] of Edwardsiella bacteria in the aquaculture system. In a previous phylogenetic study based on an internal fragment of iron-cofactored superoxide dismutase gene (sodB), Edwardsiella species and strains have been classified into two major clades, which differ in the presence of pathogenicity to fish [30] (hereafter, “pathogenic” and “non-pathogenic” clades). Therefore, we characterized Edwardsiella bacteria in the aquaculture tanks based on the illumina sequencing of the Edwardsiella sodB gene sequences. The fragment of the sodB region was PCR-amplified with the forward primer E1F [30] fused with 3–6-mer Ns for improved Illumina sequencing quality and the forward Illumina sequencing primer (5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG- [3–6-mer Ns] – [E1F]-3′) and the reverse primer 497R [30] fused with 3–6-mer Ns for improved Illumina sequencing quality [57] and the reverse sequencing primer (5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA G [3–6-mer Ns] - [497R]-3′) (0.2 μM each). The buffer and polymerase system of KOD One (Toyobo) was used with the temperature profile of 35 cycles at 98 °C for 10 s, 55 °C for 5 s, and 68 °C for 30 s (ramp rate = 1°C/s). The sequencing adaptors and sample identifier indexed were added to the amplicons, and the purification of the library and sequencing was performed as detailed above.

The output sequencing data were demultiplexed and processed with DADA2. The ASVs that were not aligned to the sodB sequences of Edwardsiella [30] were discarded. The remaining ASVs and previously reported Edwardsiella sequences [30] were aligned using mafft v.7.475 [71]. A neighbor-joining tree was reconstructed based on the maximum composite likelihood method with a bootstrap test using MEGA v.11 [72] (1000 permutations). The ASVs belonging to the pathogenic clade and those belonging to the non-pathogenic clade of Edwardsiella were distinguished within the phylogeny.

Microbe–microbe interactions

Potential positive/negative interactions between microbial ASVs were inferred based on the framework of sparse inverse covariance estimation for ecological associations (SPIEC-EASI [73]). For each water tank, patterns in the coexistence (co-occurrence) were examined with the Meinshausen-Bühlmann (MB) method as implemented in the SpiecEasi package [73] of R. The network inference based on coexistence patterns allowed us to detect pairs of microbial ASVs that potentially interact with each other in facilitative ways and/or those potentially sharing environmental preferences. Because the estimation of coexistence patterns was not feasible for rare nodes, the microbial ASVs that appeared in less than 30 samples were excluded from the input matrices of the network analysis. Network modules, within which closely associated ASVs were interlinked with each other, were identified with the algorithm based on edge betweenness using the igraph package [74] of R. For each network module in each water tank, mean partial correlation with eels’ activity level was calculated across ASVs constituting the module.

In addition to the networks representing whole coexistence patterns, we reconstructed networks depicting direct interactions between microbial ASVs. To separate the effects of direct microbe–microbe interactions from those of shared environmental preferences between microbes (i.e., shared niches), 10 latent components (latent variables) were included in the analysis based on the “sparse and low-rank” (SLR) model [75].

KEGG pathway/process profiles

To infer metabolic interactions between microbial ASVs, we performed a series of analyses based on reference genome information. We performed phylogenetic prediction of gene repertoires using PICRUSt2 v2.3.0-b [76] in order to gain the overview of the niche space defined with metagenomic information [77, 78]. ASVs that appeared in 30 or more samples across the five tanks were subjected to the analysis. Based on the inferred KEGG metabolic pathway/process profiles [79], microbial ASVs were plotted on a two-dimensional surface of a principal coordinate analysis (PCoA) based on Bray-Curtis β-diversity of KEGG metabolic pathway/process profiles.

Metabolic modeling

To infer potential metabolic interactions between microbes, we performed the species metabolic interaction analysis [80]. For the ASVs that appeared in 30 or more samples (day) in at least one aquaculture tank, we explored NCBI RefSeq genome sequences (https://www.ncbi.nlm.nih.gov/refseq/) whose 16S rRNA sequences matched those of query ASVs with ≥ 99 % identity. In the database exploration, reference genome information was available for 181 out of 417 ASVs examined. The reference genome information was subjected to genome-scale metabolic modeling as implemented in CarveMe [81] 1.5.0. Metabolic resource overlap (MRO) and metabolic interaction potential (MIP) were then estimated for each pair of microbial ASVs as implemented in SMETANA [80] 1.0.0.

Availability of data and materials

The 16S rRNA sequencing data are available from the DNA Data Bank of Japan (DDBJ) with the accession number PRJDB14313. The microbial community data are deposited at the GitHub repository (https://github.com/hiro-toju/EelMicrobiome128). All the R scripts used to analyze the data are available at the GitHub repository (https://github.com/hiro-toju/EelMicrobiome128).

Abbreviations

ASV:

Amplicon sequence variants

DDBJ:

DNA Data Bank of Japan

DO:

Dissolved oxygen

FDR:

False discovery rate

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MB:

Meinshausen-Bühlmann

MIP:

Metabolic interaction potential

MRO:

Metabolic resource overlap

NMDS:

Nonmetric multidimensional scaling

PCoA:

Principal coordinate analysis

rRNA:

Ribosomal ribonucleic acid

SLR:

Sparse and low-rank

References

  1. Ley RE, Hamady M, Lozupone C, Turnbaugh PJ, Ramey RR, Bircher JS, et al. Evolution of mammals and their gut microbes. Science. 1979;2008(320):1647–51.

    Google Scholar 

  2. Youngblut ND, Reischer GH, Walters W, Schuster N, Walzer C, Stalder G, et al. Host diet and evolutionary history explain different aspects of gut microbiome diversity among vertebrate clades. Nat Commun. 2019;10:1–15.

    Article  CAS  Google Scholar 

  3. Arumugam M, Raes J, Pelletier E, le Paslier D, Yamada T, Mende DR, et al. Enterotypes of the human gut microbiome. Nature. 2011;473:174–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Ley RE, Lozupone CA, Hamady M, Knight R, Gordon JI. Worlds within worlds: evolution of the vertebrate gut microbiota. Nat Rev Microbiol. 2008;6:776–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. McFall-Ngai M, Hadfield MG, Bosch TCG, Carey HV, Domazet-Lošo T, Douglas AE, et al. Animals in a bacterial world, a new imperative for the life sciences. Proc Natl Acad Sci U S A. 2013;110:3229–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Yukgehnaish K, Kumar P, Sivachandran P, Marimuthu K, Arshad A, Paray BA, et al. Gut microbiota metagenomics in aquaculture: factors influencing gut microbiome and its physiological role in fish. Rev Aquac. 2020;12:1903–27.

    Google Scholar 

  7. Huttenhower C, Gevers D, Knight R, Abubucker S, Badger JH, Chinwalla AT, et al. Structure, function and diversity of the healthy human microbiome. Nature. 2012;486:207–14.

    Article  CAS  Google Scholar 

  8. Barko PC, McMichael MA, Swanson KS, Williams DA. The gastrointestinal microbiome: a review. J Vet Intern Med. 2018;32:9–25.

    Article  CAS  PubMed  Google Scholar 

  9. Cho I, Blaser MJ. The human microbiome: at the interface of health and disease. Nat Rev Genet. 2012;13:260–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Legrand TPRA, Wynne JW, Weyrich LS, Oxley APA. A microbial sea of possibilities: current knowledge and prospects for an improved understanding of the fish microbiome. Rev Aquac. 2020;12:1101–34.

    Article  Google Scholar 

  11. Tarnecki AM, Burgos FA, Ray CL, Arias CR. Fish intestinal microbiome: diversity and symbiosis unravelled by metagenomics. J Appl Microbiol. 2017;123:2–17.

    Article  CAS  PubMed  Google Scholar 

  12. Bartelme RP, Smith MC, Sepulveda-Villet OJ, Newton RJ. Component microenvironments and system biogeography structure microorganism distributions in recirculating aquaculture and aquaponic systems. mSphere. 2019;4:e00143–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Rud I, Kolarevic J, Holan AB, Berget I, Calabrese S, Terjesen BF. Deep-sequencing of the bacterial microbiota in commercial-scale recirculating and semi-closed aquaculture systems for Atlantic salmon post-smolt production. Aquac Eng. 2017;78:50–62.

    Article  Google Scholar 

  14. Infante-Villamil S, Huerlimann R, Jerry DR. Microbiome diversity and dysbiosis in aquaculture. Rev Aquac. 2021;13:1077–96.

    Article  Google Scholar 

  15. Venter JC, Remington K, Heidelberg JF, Halpern AL, Rusch D, Eisen JA, et al. Environmental genome shotgun sequencing of the Sargasso Sea. Science. 1979;2004(304):66–74.

    Google Scholar 

  16. Wang Y, Sheng HF, He Y, Wu JY, Jiang YX, Tam NFY, et al. Comparison of the levels of bacterial diversity in freshwater, intertidal wetland, and marine sediments by using millions of illumina tags. Appl Environ Microbiol. 2012;78:8264–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Borges N, Keller-Costa T, Sanches-Fernandes GMM, Louvado A, Gomes NCM, Costa R. Bacteriome structure, function, and probiotics in fish larviculture: the good, the bad, and the gaps. Annu Rev Anim Biosci. 2021;9:423–52.

    Article  CAS  PubMed  Google Scholar 

  18. Nikouli E, Meziti A, Antonopoulou E, Mente E, Kormas KA. Host-associated bacterial succession during the early embryonic stages and first feeding in farmed gilthead sea bream (Sparus aurata). Genes (Basel). 2019;10:483.

    Article  CAS  PubMed  Google Scholar 

  19. Ravel J, Brotman RM, Gajer P, Ma B, Nandy M, Fadrosh DW, et al. Daily temporal dynamics of vaginal microbiota before, during and after episodes of bacterial vaginosis. Microbiome. 2013;1:29.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Halfvarson J, Brislawn CJ, Lamendella R, Vázquez-Baeza Y, Walters WA, Bramer LM, et al. Dynamics of the human gut microbiome in inflammatory bowel disease. Nat Microbiol. 2017;2:1–7.

    Article  Google Scholar 

  21. Stoeck T, Frühe L, Forster D, Cordier T, Martins CIM, Pawlowski J. Environmental DNA metabarcoding of benthic bacterial communities indicates the benthic footprint of salmon aquaculture. Mar Pollut Bull. 2018;127:139–49.

    Article  CAS  PubMed  Google Scholar 

  22. Klase G, Lee S, Liang S, Kim J, Zo YG, Lee J. The microbiome and antibiotic resistance in integrated fishfarm water: implications of environmental public health. Sci Total Environ. 2019;649:1491–501.

    Article  CAS  PubMed  Google Scholar 

  23. Djurhuus A, Closek CJ, Kelly RP, Pitz KJ, Michisaki RP, Starks HA, et al. Environmental DNA reveals seasonal shifts and potential interactions in a marine community. Nat Commun. 2020;11:254.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Ushio M. Interaction capacity as a potential driver of community diversity. Proc R Soc B Biol Sci. 2022;289:20212690.

    Article  Google Scholar 

  25. Sugita H, Miyajima C, Deguchi Y. The vitamin B12-producing ability of the intestinal microflora of freshwater fish. Aquaculture. 1991;92 C:267–76.

    Article  Google Scholar 

  26. Cai W, de La Fuente L, Arias CR. Biofilm formation by the fish pathogen flavobacterium columnare: development and parameters affecting surface attachment. Appl Environ Microbiol. 2013;79:5633–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Cao M, Huang J, Li J, Qiao Z, Wang G. Edaphobaculum flavum gen. Nov., sp. nov., a member of family Chitinophagaceae, isolated from grassland soil. Int J Syst Evol Microbiol. 2017;67:4475–81.

    Article  CAS  PubMed  Google Scholar 

  28. Thiel M, Romano MC, Kurths J, Rolfs M, Kliegl R. Twin surrogates to test for complex synchronisation. Europhys Lett. 2006;75:535.

    Article  CAS  Google Scholar 

  29. Park S, bin, Aoki T, Jung TS. Pathogenesis of and strategies for preventing Edwardsiella tarda infection in fish. Vet Res. 2012;43:67.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Yamada Y, Wakabayashi H. Identification of fish-pathogenic strains belonging to the genus Edwardsiella by sequence analysis of sodB. Fish Pathol. 1999;34:145–50.

    Article  Google Scholar 

  31. Li GY, Li J, Xiao P, Guo YH, Mo ZL. Detection of type III secretion gene as an indicator for pathogenic Edwardsiella tarda. Lett Appl Microbiol. 2011;52:213–9.

    Article  CAS  PubMed  Google Scholar 

  32. Esteve C, Biosca E, Amaro C. Virulence of Aeromonas hydrophita and some other bacteria isolated from European eels Anguilla anguilla reared in fresh water. Dis Aquat Org. 1993;16:15–20.

    Article  Google Scholar 

  33. Kozińska A, Paździor E, Pȩkala A, Niemczuk W. Acinetobacter johnsonii and Acinetobacter lwoffii - the emerging fish pathogens. Bull Vet Inst Pulawy. 2014;58:193–9.

    Article  Google Scholar 

  34. Carding S, Verbeke K, Vipond DT, Corfe BM, Owen LJ. Dysbiosis of the gut microbiota in disease. Microb Ecol Health Dis. 2015;26:26191.

    PubMed  Google Scholar 

  35. Levy M, Kolodziejczyk AA, Thaiss CA, Elinav E. Dysbiosis and the immune system. Nat Rev Immunol. 2017;17:219–32.

    Article  CAS  PubMed  Google Scholar 

  36. Kim PS, Shin NR, Lee JB, Kim MS, Whon TW, Hyun DW, et al. Host habitat is the major determinant of the gut microbiome of fish. Microbiome. 2021;9:166.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Peng M, Luo H, Kumar V, Kajbaf K, Hu Y, Yang G. Dysbiosis of intestinal microbiota induced by dietary oxidized fish oil and recovery of diet-induced dysbiosis via taurine supplementation in rice field eel (Monopterus albus). Aquaculture. 2019;512:734288.

    Article  CAS  Google Scholar 

  38. Huang W, Cheng Z, Lei S, Liu L, Lv X, Chen L, et al. Community composition, diversity, and metabolism of intestinal microbiota in cultivated European eel (Anguilla anguilla). Appl Microbiol Biotechnol. 2018;102:4143–57.

    Article  CAS  PubMed  Google Scholar 

  39. Hsu HY, Chang FC, Bin WY, Chen SH, Lin YP, Lin CY, et al. Revealing the compositions of the intestinal microbiota of three Anguillid eel species using 16S rDNA sequencing. Aquac Res. 2018;49:2404–15.

    Article  CAS  Google Scholar 

  40. Degnan PH, Taga ME, Goodman AL. Vitamin B12 as a modulator of gut microbial ecology. Cell Metab. 2014;20:769–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Hazra AB, Han AW, Mehta AP, Mok KC, Osadchiy V, Begley TP, et al. Anaerobic biosynthesis of the lower ligand of vitamin B12. Proc Natl Acad Sci U S A. 2015;112:10792–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. LaFrentz BR, LaFrentz SA, Beck BH, Arias CR. Draft genome sequences of Cetobacterium somerae 2G large and two novel Cetobacterium Isolates from intestines of channel catfish (Ictalurus punctatus). Microbiol Resour Announc. 2020;9:e01006–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Tsuchiya C, Sakata T, Sugita H. Novel ecological niche of Cetobacterium somerae, an anaerobic bacterium in the intestinal tracts of freshwater fish. Lett Appl Microbiol. 2008;46:43–8.

    CAS  PubMed  Google Scholar 

  44. Wakabayashi H, Egusa S, Yabe K. Edwardsiella tarda (Paracolobactrum anguillimortiferum) associated with pond-cultured eel disease. Nippon Suisan Gakkaishi. 1973;39:931–6.

    Article  Google Scholar 

  45. De BC, Meena DK, Behera BK, Das P, das Mohapatra PK, Sharma AP. Probiotics in fish and shellfish culture: Immunomodulatory and ecophysiological responses. Fish Physiol Biochem. 2014;40:921–71.

    CAS  Google Scholar 

  46. Moriarty DJW. Control of luminous Vibrio species in penaeid aquaculture ponds. Aquaculture. 1998;164:351–8.

    Article  Google Scholar 

  47. Gomez D, Sunyer JO, Salinas I. The mucosal immune system of fish: the evolution of tolerating commensals while fighting pathogens. Fish Shellfish Immunol. 2013;35:1729–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Schreiber T. Measuring information transfer. Phys Rev Lett. 2000;85:461–4.

    Article  CAS  PubMed  Google Scholar 

  49. Sugihara G, May R, Ye H, Hsieh C, Deyle E, Fogarty M, et al. Detecting causality in complex ecosystems. Science. 2012;338:496–500.

    Article  CAS  PubMed  Google Scholar 

  50. Fujita H, Ushio M, Suzuki K, Abe M, Yamamichi M, Iwayama K, et al. Alternative stable states, nonlinear behavior, and predictability of microbiome dynamics. Microbiome. 2023. (in press). https://doi.org/10.1101/2022.08.23.505041.

  51. Klinger D, Naylor R. Searching for solutions in aquaculture: charting a sustainable course. Annu Rev Environ Resour. 2012;37:247–76.

    Article  Google Scholar 

  52. Bostock J, McAndrew B, Richards R, Jauncey K, Telfer T, Lorenzen K, et al. Aquaculture: global status and trends. Philos Transact R Soc B Biol Sci. 2010;365:2897–912.

    Article  Google Scholar 

  53. Dawood MAO, Koshio S, Abdel-Daim MM, van Doan H. Probiotic application for sustainable aquaculture. Rev Aquac. 2019;11:907–24.

    Article  Google Scholar 

  54. Ushio M, Murakami H, Masuda R, Sado T, Miya M, Sakurai S, et al. Quantitative monitoring of multispecies fish environmental DNA using high-throughput sequencing. Metabarcoding Metagenom. 2018;2:1–15.

    Google Scholar 

  55. Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Lozupone CA, Turnbaugh PJ, et al. Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc Natl Acad Sci U S A. 2011;108(SUPPL. 1):4516–22.

    Article  CAS  PubMed  Google Scholar 

  56. Apprill A, Mcnally S, Parsons R, Weber L. Minor revision to V4 region SSU rRNA 806R gene primer greatly increases detection of SAR11 bacterioplankton. Aquat Microb Ecol. 2015;75:129–37.

    Article  Google Scholar 

  57. Lundberg DS, Yourstone S, Mieczkowski P, Jones CD, Dangl JL. Practical innovations for high-throughput amplicon sequencing. Nat Methods. 2013;10:999–1002.

    Article  CAS  PubMed  Google Scholar 

  58. Stevens JL, Jackson RL, Olson JB. Slowing PCR ramp speed reduces chimera formation from environmental samples. J Microbiol Methods. 2013;93:203–5.

    Article  CAS  PubMed  Google Scholar 

  59. Hamady M, Walker JJ, Harris JK, Gold NJ, Knight R. Error-correcting barcoded primers for pyrosequencing hundreds of samples in multiplex. Nat Methods. 2008;5:235–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Tanabe A. Claident v0.2.2018.05.29, a software distributed by author at http://www.fifthdimension.jp/. 2018.

  61. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. R Core Team. R: A language and environment for statistical computing. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2020. https://www.R-project.org/

    Google Scholar 

  63. Wang Q, Garrity GM, Tiedje JM, Cole JR. Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol. 2007;73:5261–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. 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 Res. 2013;41:D590–6.

    Article  CAS  PubMed  Google Scholar 

  65. Klappenbach JA, Saxman PR, Cole JR, Schmidt TM. Rrndb: the ribosomal RNA operon copy number database. Nucleic Acids Res. 2001;29:181–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Oksanen J. The vegan package available at https://cran.r-project.org/web/packages/vegan/index.html. 2007.

  67. Sakai M. Current research status of fish immunostimulants. Aquaculture. 1999;172:63–92.

    Article  CAS  Google Scholar 

  68. Magnadottir B. Immunological control of fish diseases. Mar Biotechnol. 2010;12:361–79.

    Article  CAS  Google Scholar 

  69. Toranzo AE, Magariños B, Romalde JL. A review of the main bacterial fish diseases in mariculture systems. Aquaculture. 2005;246:37–61.

    Article  Google Scholar 

  70. Loch TP, Faisal M. Emerging flavobacterial infections in fish: A review. J Adv Res. 2015;6:283–300.

    Article  CAS  PubMed  Google Scholar 

  71. Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30:772–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Tamura K, Stecher G, Kumar S. MEGA11: molecular evolutionary genetics analysis version 11. Mol Biol Evol. 2021;38:3022–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Kurtz ZD, Müller CL, Miraldi ER, Littman DR, Blaser MJ, Bonneau RA. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput Biol. 2015;11:e1004226.

    Article  PubMed  PubMed Central  Google Scholar 

  74. Csardi G, Nepusz T. The igraph software package for complex network research. InterJournal Complex Syst. 2006;1695:1–9.

    Google Scholar 

  75. Kurtz ZD, Bonneau R, Müller CL. Disentangling microbial associations from hidden environmental and technical factors via latent graphical models. bioRxiv. 2019. https://doi.org/10.1101/2019.12.21.885889.

  76. Douglas GM, Maffei VJ, Zaneveld JR, Yurgel SN, Brown JR, Taylor CM, et al. PICRUSt2 for prediction of metagenome functions. Nat Biotechnol. 2020;38:685–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Fahimipour AK, Gross T. Mapping the bacterial metabolic niche space. Nat Commun. 2020;11:4887.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Alneberg J, Bennke C, Beier S, Bunse C, Quince C, Ininbergs K, et al. Ecosystem-wide metagenomic binning enables prediction of ecological niches from genomes. Commun Biol. 2020;3:119.

    Article  PubMed  PubMed Central  Google Scholar 

  79. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017;45:D353–61.

    Article  CAS  PubMed  Google Scholar 

  80. Zelezniak A, Andrejev S, Ponomarova O, Mende DR, Bork P, Patil KR. Metabolic dependencies drive species co-occurrence in diverse microbial communities. Proc Natl Acad Sci U S A. 2015;112:6449–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Machado D, Andrejev S, Tramontano M, Patil KR. Fast automated reconstruction of genome-scale metabolic models for microbial species and communities. Nucleic Acids Res. 2018;46:7542–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Yuta Nogi and A Zero Inc. for their support in sampling. We are grateful to the reviewers for their comments that improved the manuscript.

Funding

This work was financially supported by JSPS Grant-in-Aid for Scientific Research (20K20586), JST FOREST (JPMJFR2048) to H.T., JSPS Grant-in-Aid for Scientific Research (20K06820 and 20H03010) to K.S., and JSPS Fellowship to H.F.

Author information

Authors and Affiliations

  1. Center for Ecological Research, Kyoto University, Otsu, Shiga, 520-2133, Japan

    Daii Yajima, Hiroaki Fujita, Ibuki Hayashi, Genta Shima & Hirokazu Toju

  2. Integrated Bioresource Information Division, BioResource Research Center, RIKEN, Tsukuba, Ibaraki, 305-0074, Japan

    Kenta Suzuki

Authors
  1. Daii Yajima
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Hiroaki Fujita
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ibuki Hayashi
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Genta Shima
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Kenta Suzuki
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Hirokazu Toju
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

H.T. designed the work with D.Y.. D.Y., G.S., and H.F. performed the experiments. D.Y., H.F., I.H., K.S., and H.T. analyzed the data. H.T. and D.Y. wrote the paper with all the authors. The authors read and approved the final manuscript.

Corresponding author

Correspondence to Hirokazu Toju.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Phylum-level community structure.

Additional file 2: Figure S2.

Order-level community structure.

Additional file 3: Figure S3.

Genus-level community structure.

Additional file 4: Figure S4.

Multivariate analysis of community structure (genus level).

Additional file 5: Figure S5.

ASV-level comparison of correlation with environmental variables and eels’ activity level.

Additional file 6: Figure S6.

Phylogenetic analysis of Edwardsiella.

Additional file 7: Figure S7.

Taxonomy of the nodes within the coexistence networks.

Additional file 8: Figure S8.

Correlations with environmental variables.

Additional file 9: Table S1.

Properties of the ASVs within the networks.

Additional file 10: Figure S9.

Properties of network modules.

Additional file 11: Figure S10.

Comparison of network reconstruction methods.

Additional file 12: Figure S11.

Rarefaction curves of the samples.

Additional file 13: Figure S12.

Temporal dynamics of potential fish pathogens.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Yajima, D., Fujita, H., Hayashi, I. et al. Core species and interactions prominent in fish-associated microbiome dynamics. Microbiome 11, 53 (2023). https://doi.org/10.1186/s40168-023-01498-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40168-023-01498-x

Keywords

Microbiome

ISSN: 2049-2618

Contact us