Microbiome and Metabolome Analyses in Different Closed-Circulation Aquarium Systems and Their Network Visualization.

Understanding the causes of microbiome formation and its relationship to environmental conditions is important to properly maintain recirculating aquaculture systems (RASs). Although RAS has been applied to numerous fish types and environmental conditions (e.g., loading intensity), the effects of these environmental conditions (i.e., fish type and loading intensity) on microbiome composition are limitedly known. Therefore, we established three experimental aquarium tanks to explore the effects of fish type, loading intensity, filter pore size, and rearing day on microbiome compositions: (1) a tank for Acanthogobius flavimanus, (2) for Girella punctata, and (3) for G. punctata with higher loading intensity. Multivariate analysis showed that the microbial community composition differed among the tanks, indicating that the fish type and loading intensity significantly affected microbiome formation in rearing water. Some microbes, such as Sediminicola and Glaciecola, were detected at a higher loading intensity, indicating that these microbes might be an indicator of eutrophic conditions in the aquacultural systems. In addition, a partial correlation network revealed a connection between microbes and metabolites in the aquarium tanks. Such a microbe-metabolite network might be a clue to control the microbiome by adjusting the molecule abundance in the aquacultural environment.

Introduction

Microbiome plays a crucial role in maintaining the functions and properties of biomediated systems. For example, microbes drive material cycling in their environment through metabolisms by degrading detritus, transforming nitrogen into other chemical forms, and many other processes.[1,2] Simultaneously, microbes sometimes act as pathogens, causing other organisms to fall into the disease.[3,4]

In recent years, aquaculture, which is a biomediated system, has attracted increased attention owing to population growth. In particular, recirculating aquaculture systems (RASs), in which breeding water recirculates within the system through filtering technology, are promising owing to their ability to reduce environmental impact by limiting emissions into outer ecosystems and high resilience to natural disasters.[5] Exploring the microbiome is also important to characterize the system state of RAS because environmental microbiomes can easily affect the health of rearing fish.[6]

Generally, the microbiome community composition can be determined primarily by the source community composition and, secondarily, by environmental conditions—a well-known environmental filtering theory.[7] In RAS with artificial seawater, the potential source is the microbiome associated with the host fish because it supplies microbes on the skin and feces into the surrounding water. Previous studies have demonstrated that such microbiome composition differs according to the fish species. Six fish species in the Gulf of Mexico were sampled throughout a year, where the fish species rather than geographic location and sampling date, was the most critical factor in determining the skin microbial community composition.[8] Different fish species have a unique microbiome in their gut and the associated feces, in which feeding behavior determines the microbiome composition.[9] Therefore, fish type could be the primary factor for the microbiome in RAS. In addition, environmental conditions can alter the microbiome composition. While numerous factors can affect the environmental conditions in RAS, the most critical factor is the loading strength as feeding.[10] Intensive loading supplies excessive organic matter and nitrogen compounds in the tanks, thereby deteriorating water quality and, thus, fish performance, which could select microbes that survive the eutrophic environment. Only a few studies have demonstrated the microbiome in RAS despite the need to control this microbiome.

Generally, each taxonomic group of microbes is likely to have its own niche according to the organic matter type to feed on. A metabolome is a set of small molecules that can be a source for microbial growth or an outcome of microbial metabolism. Therefore, the metabolome in the rearing water may explain the microbial abundance to some extent. While previous studies have revealed such microbe and metabolome connections in environments,[11,12] no studies have adopted such a multi-omics approach for RAS.

In this study, we set up laboratory-scale closed-circulation tanks by considering RAS to analyze the effect of fish types and feed loading on the microbiome and metabolome in rearing water. Furthermore, we used a multi-omics approach (microbiome and metabolome) to visualize the microbe–metabolite connection in these aquarium systems.

Materials and Methods

We set up three 20-L aquariums to rear fish: the first was for Acanthogobius flavimanus (AFL), the second was for Girella punctata (GPL), and the third was for G. punctata with higher loadings (GPH). Three individual fish were in GPL and GPH, and six individual fish were in AFL. The feeding rate was 12 g/day/tank for AFL and GPL, and 36 g/day/tank for GPH. The rearing water circulated within the filtration system, in which a sponge filtered out coarse organic matter, and porous filter media (Eheim, Deizisau, Germany) provide habitats for microbes to perform nitrification. The water sampling was conducted on the 4th, 8th, 13th, 16th, and 38th day after the fish were placed in the artificial seawater tanks. We collected the microbes and metabolites by filtering 1000 mL of the tank water sequentially through 8.0, 1.0, and 0.2 μm filter papers. The filters were stored in a freezer until the analyses, in which half of the filter was used for the microbiome, and the rest was for the metabolome analysis.

16S rRNA amplicon sequencing was performed for the microbiome analysis. DNA was extracted from the filter by ethanol precipitation and then amplified by PCR using forward and reverse primers (Eurofins Genomics K.K., Tokyo, Japan), dNTP, and ExTaq (Takara Bio Inc., Shiga, Japan). The PCR products were purified using magnetic beads of AMpure (Beckman Coulter, Inc., Brea, USA) and then re-dissolved in pure water. The purified DNAs were diluted to the same concentration after measuring the concentration of the purified DNA with a Qubit Fluorometer (Thermo Fisher Scientific, Waltham, USA) and then mixed into one library. The pooled library and PhiX Control v3 (Illumina, Inc., San Diego, USA) were denatured under NaOH. The library was diluted to a final concentration of 3.5 pM with the HT1 hybridization buffer and then spiked with the denatured PhiX solution. The final solution was loaded into MiSeq Reagent Kits V2 (Illumina, Inc., San Diego, USA) and analyzed using MiSeq (Illumina, Inc., San Diego, USA). The sequence data in FASTQ format were processed using QIIME. The forward and backward reads were joined, and the chimeras were then filtered out. The reads were clustered into OTUs based on 97% similarity and annotated with an RDP database.

Nuclear magnetic resonance (NMR) analysis was performed for the metabolome analysis of the collected filter. The filter was placed in 15 mL of pure water and then sonicated using a Bioruptor water bath sonicator (BR2024A, Sonic Bio Co., Ltd., Kanagawa, Japan) to remove metabolites. The suspension was frozen and lyophilized, which was extracted in potassium phosphate (KPi) buffer solution in D2O containing 1 mmol/L sodium 2,2-dimethyl-2-silapentane-5-sulfonate (DSS) at 65 °C for 15 min under shaking (1400 rpm). After centrifugation at 14,000 rpm for 5 min, the supernatants were transferred into NMR tubes. NMR spectra of 2D J-resolved (2D J-res) were recorded using an AVANCE II 700 MHz NMR spectrometer (Bruker BioSpin, Rheinstetten, Germany). Acquisition parameters of NMR were as follows: time domain data size was 16,384 for f2 and 32 for f1, spectral width was 12,500 Hz for f2 and 50 Hz for f1, and number of scans was 32 with 2 s recycle delay. Using software Topspin (Bruker Biospin, Rheinstetten, Germany), we processed the spectra with sine-bell window function, zero-fillings to 128 points, tilt correction, and symmetrization. For the following analysis, the 2D J-res spectra were projected to 1D along a 1H chemical shift. Using rNMR software,[13] we defined 78 regions of interest (ROI) of the projected 2D J-res spectra and calculated the maximum intensity within each ROI. The peak intensities were preprocessed to convert the integral intensity to 1, which generated relative proportion data. Peak annotation was performed by superimposing the standard material spectra of 2D J-res and 1H–13C heteronuclear single quantum coherence (HSQC) after narrowing down the signal candidates using SpinAssign[14] and SpinCouple.[15]

The differences in the whole microbiome and metabolome between samples were visualized by non-metric multidimensional scaling (nMDS) based on the Bray–Curtis dissimilarity index. Two-way PERMANOVA and pairwise two-way PERMANOVA were performed to test whether the microbiome and metabolome profiles differed according to the aquarium tank and filter pore size. Because our data set had a layered structure, we constructed a generalized linear model (GLM) for factor analysis.[16] The effects of the experimental conditions (i.e., fish type, loading intensity, filter pore size, and rearing day) on each microbe and metabolite were assessed using the GLM.where β0 is the intercept, β1, β2, β3, and β4 are the slopes of the dummy variables (x1; A. flavimanus = 0 and G. punctata = 1, x2; low loading = 0 and high loading = 1, x3; 0.2 μm = 0 and 1 μm = 1, x4; 0.2 μm = 0 and 8 μm = 1), and β5 is the slope of the incubation day (x5). The multi-omics analysis for each microbe–microbe, metabolite–metabolite, and microbe–metabolite pair was performed using a partial correlation network based on the following GLMwhere β0–β5 and x1–x5 are same as the above model, and β6 is the slope for the proportion of another microbe/metabolite (x6). For the simplification of the network, the microbial taxa with a maximum abundance exceeding 1% and the annotated metabolites were selected for constructing the GLM. The slope of β6 and its significance for microbe–metabolite pairs were used for the partial network visualization: the significant partial correlations[17] (p > 0.05) were connected by coloring with a positive or negative slope. The network layout was calculated using the Kamada–Kawai algorithm,[18] and the community structure of the network was inferred by the leading eigenvector of the modularity matrix.[19] For comparison, a normal correlation network was also constructed for the significant (p < 0.05) microbe–metabolite pairs. The network visualization was conducted using Tidygraph and Ggraph packages on R.

Results and Discussion

Microbiome Analysis

Microbiome in the rearing water comprised three main phyla: Proteobacteria, Bacteroidetes, and Verrucomicrobia (maximum of 95.1, 75.8, and 58.4%, respectively, Figure S2). Such a bacterial composition has been commonly observed, in which a recent review reported that Proteobacteria and Bacteroidetes are the top two phyla in their abundance for RAS regardless of the rearing fish species.[20] In addition, some studies have reported the dominance of Verrucomicrobia followed by Proteobacteria and Bacteroidetes in aquaculture systems.[21,22] The GLM suggested that the proportion of Proteobacteria and Bacteroidetes did not differ among fish types, whereas that of Verrucomicrobia in the G. punctata tank was larger (Table S1), suggesting that Verrucomicrobia might be a more host-specific phylum in the rearing water. At the genus level, Sediminicola, Glaciecola, unassigned Saprospiraceae, and Rubritalea exceeded 50% of relative abundance, whereas such dominance varied with sample (Figure A).

Figure 1

Water microbiome and metabolome in the three experimental tanks; A. flavimanus with low loading intensity, G. punctata with low loading intensity, and G. punctata with high loading intensity. (A) Microbial composition at the genus level for three tanks collected via the three filter types. (B) nMDS score plot for the genus-level microbiome. (C) Metabolite composition for three tanks collected via the three filter types. (D) nMDS score plot for metabolome.

nMDS revealed the β-diversity of the microbiome at the genus level. Microbiome compositions were clearly separated according to the tank and filter pore size on the nMDS score plot, where the microbiomes in the A. flavimanus tank were plotted on the higher nMDS axis 1 and those on the coarse filter were plotted on the higher nMDS axis 2 (Figure B). Statistically, two-way PERMANOVA demonstrated that the microbial community compositions varied with tank and filter factors (p < 0.05). In addition, pairwise two-way PERMANOVA showed that each tank pair (AFL vs GPL, AFL vs GPH, and GPL vs GPH) significantly differed in the microbiome. Because our rearing experiment started with artificial seawater, the microbiome difference according to tank would be ascribed to fish type and eutrophic conditions due to the excess feeding. Importantly, although GPL and GPH were different tanks, the microbiomes were similar to each other compared with the AFL tank (Figure B). Such discrepancy between A. flavimanus and G. punctata tanks suggests that the microbiome composition in the rearing water was primarily influenced by the original supply source of microbes, probably originating from the host surface and intestinal bacteria.[23,24]

To detect some taxa that changed their abundance with each of the experimental conditions, a GLM was applied to each genus, where the proportion was explained by four factors: fish type, loading intensity, filter pore size, and rearing day. The GLM extracted some taxa whose proportions varied with these experimental conditions. For example, Rubritalea and Vibrio were greater in G. punctata tanks, whereas Sediminicola and Pseudoalteromonas were greater in the A. flavimanus tank (Table S2). These trends were consistent with the fish mucocutaneous microbes for some of the genera: Vibrio was abundant on the G. punctata skin and Sediminicola was abundant on the A. flavimanus skin (Figure S3). This suggests that the fish mucocutaneous microbes are associated with the water microbiome as a direct contact site between them and that the microbiome on the fish skin might influence the water microbiome as a microbial source. In addition, Sediminicola and Glaciecola, for example, increased their proportions in the high loading tank. Because these genera increased their abundance under intensive feeding conditions, the abundance of these genera might be an indicator of high environmental loading in aquarium tanks. Likewise, Glaciecola, for example, was found to have greater abundance in conventional intensive mariculture systems compared with RAS.[25]

Metabolome Analysis

We found 78 signals in the 2D J-res spectra, among which 35 signals were annotated (Figure S1). The peak signal intensities for the same annotated substances were summed together to calculate the relative proportion (Figure C). In rearing water, betaine and acetate frequently accounted for a large proportion. In functional group units, the amino acid groups (Ala, Leu, Val, Phe, Ile, Tyr, Glu, and Met) were found in a large proportion.

nMDS was performed based on the relative intensity of the ROIs for the 2D J-res spectra (Figure D). According to the two-way PERMANOVA, metabolome profiles varied with aquarium tanks and filter types (p < 0.05). In pairwise PERMANOVA, significant differences appeared between GPL and GPH, and between AFL and GPH, whereas no significant difference was observed between GPL and AFL, suggesting that intensive loading rather than the host fish type affected the metabolite profiles in these aquarium tanks. In addition, the metabolome discrepancy between the high-loading and low-loading tanks was evident, especially in the coarser filter. This indicates that the compositional changes in the metabolome were caused by the external input of coarse materials, such as feeds and feces, rather than small molecular microbial metabolites in the aquarium. Focusing on individual metabolites, the GLM suggested that intensive loading increased the proportion of betaine and methionine (Table S3). Because betaine and methionine are typically used as food additives to promote fish growth[26−28] and are also detected in our feed and fecal samples, these accumulations would be ascribed to leakage from feeds and digested feces. In contrast, acetate, formate, and ethanol decreased their proportion in the high-loading tank (Table S3). These molecules are likely to be the biproducts of biological metabolism: the proportion of these substances in the high-loading tank was reduced, probably because the influence of the feed components was too large.

Network Visualization

Because our data set had a nested structure, the normal correlation network analysis between the microbiome and metabolome (Figure A) would be influenced by the spurious correlation of other experimental factors (i.e., fish type, loading intensity, filter type, and rearing day in this case). Therefore, the GLM-based partial correlation network was constructed to infer the microbe–metabolite correspondence (Figure B). For the visualization, we extracted only the edges between the microbes and metabolites, and each node was depicted with degree centrality for the microbe–metabolite links.

Figure 2

Network visualization. (A) Normal correlation network and (B) partial correlation network for genus-level microbiome and metabolome data set. Each node shows the metabolite or microbe, and each link shows a significant correlation. Node size shows degree centrality, node shape shows the node type (microbe/metabolite), node color shows the community inferred by leading eigenvalue, and edge color shows whether positive or negative correlation.

Importantly, the partial correlation network rather than the normal correlation network explicitly showed the relationships with microbes specific to functional metabolite groups. The layout algorithm (i.e., Kamada–Kawai method in this study) evenly distributed the microbe–metabolite nodes for the normal correlation network (Figure A) but localized them for the partial correlation network (Figure B), indicating that the partial correlation network extracted several microbe–metabolite nodes with similar links. Interestingly, community structure based on the leading eigenvector of the modularity matrix indicated that the metabolite nodes with similar chemical functions connected their edges with similar microbes. For example, the amino acids of Ala, Val, Leu, Ile, and Glu belonged to the same community, which were shared with Marinicella, Thalassomonas, Thalassospira, Glaciecola, and so forth (Figure B), indicating that these microbes might be more actively involved in amino acid metabolisms. Organic acids of lactate and formate, as well as alcohols of ethanol and glycerol, belonged to the same community (Figure B). Among the microbes within this community, Roseivirga and unannotated Alteromonadaceae had positive relationships with these organic acids and alcohols and a negative relationship with glucose. These metabolites of organic acids and alcohols are the products of the glycolysis, and the contrasting relationships with glucose might explain why these microbes play a central role in the glycometabolism in the aquarium tanks.

Conclusions

In this study, we revealed that microbiome compositions in aquarium tanks were differentiated according to fish type and feeding intensity. Such differences were attributed to the source origin of microbes (fish type) and the rearing conditions due to intense loading. In addition, the partial correlation network visualized the microbes–metabolite links, providing insights into the control of microbiomes in biomediated systems.