# Disruption of fish gut microbiota composition and holobiont's metabolome during a simulated Microcystis aeruginosa (Cyanobacteria) bloom

Alison Gallet, Sébastien Halary, Charlotte Duval, Hélène Huet, Sébastien Duperron, Benjamin Marie

► To cite this version:

Alison Gallet, Sébastien Halary, Charlotte Duval, Hélène Huet, Sébastien Duperron, et al. Disruption of fish gut microbiota composition and holobiont's metabolome during a simulated Microcystis aeruginosa (Cyanobacteria) bloom. Microbiome, 2023, 11 (1), pp.108. 10.1186/s40168-023-01558-2. mnhn-04102614

HAL Id: mnhn-04102614

https://mnhn.hal.science/mnhn-04102614

Submitted on 22 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
RESEARCH

Open Access

<figure><img src="image_2.png" /></figure>

Check for
updates

# Disruption of fish gut microbiota
composition and holobiont's metabolome
during a simulated *Microcystis aeruginosa*
(Cyanobacteria) bloom

Alison Gallet¹, Sébastien Halary¹, Charlotte Duval¹, Hélène Huet², Sébastien Duperron¹,³*† and Benjamin Marie¹*†

**Abstract**

**Background** Cyanobacterial blooms are one of the most common stressors encountered by metazoans living in freshwater lentic systems such as lakes and ponds. Blooms reportedly impair fish health, notably through oxygen depletion and production of bioactive compounds including cyanotoxins. However, in the times of the "microbiome revolution", it is surprising that so little is still known regarding the influence of blooms on fish microbiota. In this study, an experimental approach is used to demonstrate that blooms affect fish microbiome composition and functions, as well as the metabolome of holobionts. To this end, the model teleost *Oryzias latipes* is exposed to simulated *Microcystis aeruginosa* blooms of various intensities in a microcosm setting, and the response of bacterial gut communities is evaluated in terms of composition and metabolome profiling. Metagenome-encoded functions are compared after 28 days between control individuals and those exposed to highest bloom level.

**Results** The gut bacterial community of *O. latipes* exhibits marked responses to the presence of *M. aeruginosa* blooms in a dose-dependent manner. Notably, abundant gut-associated Firmicutes almost disappear, while potential opportunists increase. The holobiont's gut metabolome displays major changes, while functions encoded in the metagenome of bacterial partners are more marginally affected. Bacterial communities tend to return to original composition after the end of the bloom and remain sensitive in case of a second bloom, reflecting a highly reactive gut community.

**Conclusion** Gut-associated bacterial communities and holobiont functioning are affected by both short and long exposure to *M. aeruginosa*, and show evidence of post-bloom resilience. These findings point to the significance of bloom events to fish health and fitness, including survival and reproduction, through microbiome-related effects. In the context of increasingly frequent and intense blooms worldwide, potential outcomes relevant to conservation biology as well as aquaculture warrant further investigation.

**Keywords** Cyanobacteria, Gut microbiota, Metagenomics, Metabolomics, Medaka fish

<footnote>†Sébastien Duperron and Benjamin Marie equally contributed to this work.</footnote>

<footnote>*Correspondence:</footnote>

<footnote>Sébastien Duperron
sebastien.duperron@mnhn.fr
Benjamin Marie
benjamin.marie@mnhn.fr</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
## Background

Organisms living in lakes and ponds are exposed to cyanobacterial blooms throughout their life [1–4]. Blooms are natural events, yet increased eutrophication and global change associated with human activities make them increasingly frequent, abundant and persistent worldwide [5, 6]. Cyanobacterial blooms affect the whole ecosystem, including the health of teleost fish which occupy higher trophic levels [1, 2, 7]. They can trigger mass mortalities of animals due to oxygen depletion. Cyanobacterial metabolites display a broad range of bioactivities and include various toxins, digestive enzyme inhibitors, antimicrobials, and cytotoxic compounds [8]. As a consequence, cyanobacteria cells, extracts and purified cyanotoxins all induce deleterious effects on teleost fishes [2, 9–12]. One of the most toxic and frequent cyanotoxins is microcystin-LR (MC-LR), a hepatotoxin that accumulates in fish liver mostly following ingestion of contaminated prey or water. Intraperitoneal injection doses of 400 to 1000 µg.kg⁻¹ are lethal to fish [13, 14]. MC-LR adversely affects fish development, reproduction and behavior [15, 16]. It specifically affects the liver in which MC-LR triggers inflammation, oxidative stress, and tissue necrosis [17–19]. Aquatic animals are exposed to cyanobacteria and their toxins through oral ingestion and transfer absorption through the intestine [15, 20, 21]. In this respect, the role of gut-associated microbiota in holobiont's response is currently underestimated [22]. In teleosts, gut microbiome tends to be less diverse than observed in larger vertebrates, and is often dominated by Proteobacteria, Bacteroidetes, and Firmicutes [23, 24]. As for other vertebrates, gut-associated symbionts play important roles in nutrition (polymer degradation, vitamin production), immunity, protection against pathogens, and homeostasis [25–27]. Various host- and environment-related factors are reported to influence gut microbiome structure [28], and the gut microbiota has recently emerged as a primary target for microbiome-aware ecotoxicological concerns [29–33]. However, only a limited number of studies have investigated the effect of cyanobacterial blooms on fish gut microbiota [10, 34–36]. Whole chemical extracts of few *Microcystis* strains were for example shown to influence the composition of gut bacterial communities of medaka fish in microcosm-based experiments, while MC-LR alone did not [35]. This and few other works emphasize that metabolitecocktails and whole cells, rather than toxins alone (microcystins), should be considered for realistic assessment of the microbiome impairs [2, 15, 35, 37], yet the effect of exposure to a range of environmentally relevant levels of cyanobacteria has not been evaluated so far.

In the present study, the impact of a cyanobacterial bloom on the composition and functions of fish

gut-associated bacterial communities, and on the metabolite composition in various host tissues is evaluated using an experimental approach. Due to its resistance to stress and disease, the teleost fish *Oryzias latipes* (medaka) has emerged as a model of choice for aquatic ecotoxicology in general, and a body of literature now specifically documents the effects of cyanotoxins on animal tissue [38–41]. Herein, *O. latipes*, has been exposed to *Microcystis aeruginosa*, the most common bloom-forming cyanobacterium in temperate lentic freshwaters [42]. Three exposure levels were defined within a range of values commonly reported in natural lakes and ponds. The highest exposure level, 100 µg.L⁻¹ *Chla*, corresponds to a high bloom intensity, yet commonly reported in eutrophic lakes [43], while the medium level (10 µg.L⁻¹ *Chla*) is commonly considered a threshold value indicative of a bloom [44]. A first 28-day exposure simulated a long bloom event. Because blooms are highly dynamic events in natural systems, post-bloom resilience was investigated. Then, the hypothesis of a priming effect, translating into a lower impact of a second bloom, was tested. To this end, a post-bloom depuration phase was conducted for 4 days, followed by a second exposure to the highest *M. aeruginosa* concentration for 5 days. Bacterial community compositions were characterized using 16S rRNA gene sequencing, and metabolite contents were profiled by LC-MS/MS. Metagenomes of unexposed control fish gut communities were compared to those of specimens exposed to the highest bloom level to compare their respective annotated functions. Bacterial community and metabolites compositions were then compared using a multi-omics approach to identify correlation networks associated with teleost holobiont response. By testing the dose-dependent effect of a cyanobacterial bloom on teleost gut bacterial microbiota, documenting holobiont post-bloom response, and investigating the effect of a second bloom, this study addresses for the first time the dynamics of holobiont response to cyanobacterial blooms.

## Methods

### Experimental design and sampling

Experimental procedures were carried out in accordance with European legislation on animal experimentation (European Union Directive 2010/63/EU) and were approved for ethical contentment by an independent ethical council (CEEA Cuvier n°68) and authorized by the French government under reference number APAFiS#19316-2019032913284201 v1.

Experiments were performed in 10-liter aquaria (microcosms) with 7-month-old adult male Japanese medaka fish *Oryzias latipes* provided by the AMAGEN platform (Gif-sur-Yvette, France). Before the whole
experiment, five fish were sampled as controls for the histological analyses then fish were pre-acclimatized in clear water (2 weeks) in 15 aquaria, each one containing 8 fish. Prior to the first exposure at day 0, five fish were randomly sampled among the aquaria, and 10 mL of water of each aquarium were pooled, as references of initial fish and water conditions (d0, Fig. 1). Fish were then exposed for 28 days to five treatments: water (control, 0); water containing Z8 medium [45], i.e., the medium used to cultivate cyanobacteria (control for water enrichment with nutrients from Z8); and water containing three environmentally relevant concentrations of live *Microcystis aeruginosa*, i.e., 1, 10, and 100 µg.L⁻¹ Chla (1, 10, 100), respectively (d28, Fig. 1). Each treatment was carried out in three aquaria (labelled a, b, and c). At d28, four fish, and 150 mL of water were sampled in each aquarium. Bottom-growing

biofilms and faeces were also sampled by scraping using sterile scalpel and tweezers, respectively. After sampling, remaining fish from each of the three aquaria exposed to one condition were pooled and transferred to a single aquarium, filled with clear water (treatment 0) for a 4-day depuration period of exactly 110 h (d33, Fig. 1). At d33, three fish and 150 mL of water were sampled in each aquarium. All fish were then exposed in the same aquaria for 5 days (exactly 134 h) to the highest concentration of *M. aeruginosa*, 100 µg.L⁻¹ Chla (d39, Fig. 1). At d39, four fish, 150 mL of water and faeces were sampled in each aquarium. Along the two successive exposures, 10 mL of *M. aeruginosa* culture were sampled every 2 days but only one sample per week was analyzed further. Dataset S1 provides details of sampled individuals and performed analyses.

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Experimental design and *O. latipes* sampling strategy. After a 14-day acclimatization period, adult male *O. latipes* were exposed to two successive exposures. During the first 28 days (d0 to d28), fish were exposed to five different treatments, each carried out in three 10-L aquaria. Treatments consisted of two controls (0, Z8) and three concentrations of *Microcystis aeruginosa* (1, 10, and 100 µg.L⁻¹ Chla). After that, a 4-day depuration phase followed (d28 to d33) performed in clear water (0), then fish were exposed again for 5 days (d33 to d39) to *M. aeruginosa* (100 µg.L⁻¹ Chla). Colors represents the different treatments: blue = water, yellow = Z8 growth medium, light green = 1 µg.L⁻¹ Chla, green = 10 µg.L⁻¹ Chla, dark green = 100 µg.L⁻¹ Chla. Samples were collected at d0 (day 0), d28 (day 28), d33 (day 33), and d39 (day 39). Further analyses were performed on *M. aeruginosa* culture, water, biofilm, fish food, feces samples, and different tissues (Supplementary Note 1 and Figure S1a–d). Dataset S2 provides details on sampling counts.</figcaption></figure>
Microcystis aeruginosa production

Blooms were simulated in lab using the non-axenic and
easy-to-cultivate *M. aeruginosa* mono-clonal strain PMC
728.11 maintained in the Paris Museum Collection which
can produce diverse variants of bioactive metabolites
(Supplementary Note 2). The strain was cultivated in Z8
medium at 25 ± 1 °C with a 16-h:8-h light/dark cycle (at
14 µmol.m⁻².s⁻¹) in 2-L bottles all along the experiment
(LED strips TR-LD-24VDC-SMD5050-NW and modula-
tor). The concentrations of *M. aeruginosa* were estimated
using Chlorophyll *a* extraction [46] and absorbance
measurements as a proxy using a spectrophotometer
(Cary 60 UV-Vis, Agilent). Every 2 days, after water
renewal, *M. aeruginosa* was measured in each aquarium
using a fluorometer (FluoroProbe III, bbe Moldaenke),
and appropriate volume of the culture was added as per
the desired final concentrations (1, 10, or 100 µg.L⁻¹
Chla).

Monitoring of experimental parameters

Water was aerated using a pump (Air-Flow 4, SuperFish,
Netherlands). Every second day throughout the experi-
ment, water parameters were monitored (pH, tempera-
ture, conductivity, nitrates, and nitrites), aquaria were
cleaned (faeces removed by aspiration), half of the water
was replaced with fresh water composed of 2/3 osmosis
(RiOs 5, Merck Millipore) and 1/3 filtered. On the same
timeline, *M. aeruginosa* concentrations were measured
and adjusted to maintain bloom levels in exposed treat-
ments, and 1 mL of sterile Z8 medium was added in
the Z8 treatment. Fish were exposed to constant tem-
perature (23 ± 1 °C), pH (7.5 ± 0.1), and conductivity
(234 ± 22 µS.cm⁻¹), to low levels of nitrates (≤ 1 mg.L⁻¹)
and nitrites (≤ 4 mg.L⁻¹), to a 12-h:12-h light/dark cycle
(Philips MASTER TL-D 36W/840). They were fed twice
daily (~ 3–5% of the fish biomass per day) with Nutra
HP 0.3 (Crude protein 57, Crude fat 17, N.F.E 7.5, Ash
10, Crude fiber 0.5, Phosphorus 1.7, Vitamins A, D3, E;
Skretting, Norway). Microcystin (MC) concentration
was monitored on a regular basis on water samples in *M.
aeruginosa*-containing treatments and quantified using
enzyme-linked immunosorbent assay (ELISA) analyses
(Microcystins-ADDA SAES ELISA, Eurofins Abraxis).
Each sample was analyzed in duplicates and MC con-
centrations were determined according to the MC-LR
response standard curve (Dataset S2).

Fish, *M. aeruginosa* culture, water, biofilm, food, and feces
processing

Fish were anesthetized in 0.1% tricaine methanesulfonate
(MS-222; Sigma, St. Louis, MO) buffered with 0.1%
NaHCO₃ and sacrificed. Whole guts (including content,
due to small size), muscles and livers were dissected,

flash-frozen in liquid nitrogen and stored at – 80 °C. For
histopathological examinations, livers from fish sam-
pled before the whole experiment and at d28 and d39
(one fish per aquarium) were dissected, fixed in David-
son fixative as previously described [15], maintained for
24 h at 4 °C then dehydrated in 70% ethanol and con-
served at 4 °C. Livers samples were then embedded in
paraffin and blocks were cut into 4-µm thick sections,
stained with hematoxylin-eosin-saffron (HES), periodic
acid schiff (PAS), and Perls Prussian blue, and observed
under photonic microscope (Zeiss, Germany). Aquarium
water samples were filtered on a 0.22-µm filter (Nucleo-
pore Track-Etch Membrane) and frozen. *M. aeruginosa*
culture (4 mL) and biofilm samples were centrifugated
(10 min, 10 °C, 3220×g) and pellets were frozen. Feces
pellets were directly frozen. A food sample was kept for
DNA extraction.

Metabolites extraction

Metabolite contents were extracted from fish livers,
whole guts (including content) and muscles, *M. aerugi-
nosa* cultures, and biofilms. Prior to sonication, muscles
were freeze-dried then ground using a bead beater (Tis-
sueLyser II, Qiagen) while cultures and biofilms were
only freeze-dried. Samples were weighted, then sonicated
in 75% methanol (1 mL per 100 mg of tissue, 3x, on ice)
and centrifuged (10 min, 4 °C, 15,300×g). Supernatants
containing metabolite extracts were kept at – 20 °C for
mass spectrometry analyses. All pellets were discarded,
except gut pellets dried and kept at – 80 °C to perform a
subsequent DNA extraction on the same gut tissue.

Mass spectrometry data processing and analysis

Each metabolite extract from fish livers, guts and mus-
cles was analyzed by ultra high-performance liquid chro-
matography (UHPLC; ELUTE, Bruker) coupled with a
high-resolution mass spectrometer (ESI-Qq-TOF Com-
pact, Bruker) at 2 Hz speed, on simple MS mode then
on broad-band Collision Ion Dissociation (bbCID) or
autoMS/MS mode on the 50–1500 *m/z* range. Three fea-
ture peak lists were generated from MS spectra within
a retention time window of 1–15 min and a filtering of
5000 counts using MetaboScape 4.0 software (Bruker).
The three peak lists consisted of the area-under-the-
peaks of extracted analytes from the three tissues (gut,
liver, muscle) sampled at d28, d33, and d39, resulting in
1672, 909, and 3127 analytes, respectively. The genuine
metabolite content of the culture was investigated on
metabolite extracts using LC-MS/MS approach, com-
bined with molecular network analysis and metabolite
annotation using a cyanobacterial metabolite reference
database, as previously described [47]. Prior to analyses,
Pareto scaling was applied on the datasets [48]. Principal
component analyses (PCA) were performed to compare the metabolite composition among groups using the *mixOmics* package [49] in R 4.1.0 (R Core Team, 2021). The variance among groups was compared conducting PERMANOVA (999 permutations) based on euclidean distance with *vegan* [50] followed by Bonferroni-adjusted pairwise comparisons with *RVAideMemoire* [51].

## DNA extraction

DNA was extracted from gut pellets resulting from the metabolites extraction, culture, biofilm, faeces, water filters and food using the ZymoBIOMICS DNA Miniprep kit (Zymo Research, California). Prior to DNA extraction, all pellets were re-suspended in Eppendorf tubes with 750 µL of the ZymoBIOMICS™ lysis solution, then the contents were transferred to the ZR BashingBead™ lysis tubes. Water filters were cut into pieces then transferred to the ZR BashingBead™ lysis tubes. All steps were conducted following the manufacturer's instructions except for mechanical lysis, achieved on a bead beater (Tissue-Lyser II, Qiagen) during 6 × 1 min. An extraction blank was sequenced as a control and the 8 ASVs that were > 1% of reads were removed from the whole dataset as potential contaminants. Together, they represented 0 to 2.8% of reads in the different samples (mean 0.46%). The quality and quantity of the extracted DNA was tested on Q-bit (Thermo).

## Bacterial 16S rRNA gene sequencing and analyses

The V4–V5 variable region of the 16S rRNA gene was amplified using 479F (5'-CAGCMGCYGCNG-TAANAC-3') and 888R primers (5'-CCGYCAATTC-MTTTRAGT-3') [52], and sequenced (Illumina MiSeq paired-end, 2 × 250 bp, GenoScreen, France). Paired-end reads were demultiplexed, quality controlled, trimmed and assembled with FLASH [53]. Sequence analysis was performed using the QIIME 2 2020.11 pipeline [54]. Chimeras were removed and sequences were trimmed to 367 pb then denoised using the *DADA2* plugin, resulting in Amplicon Sequence Variants (ASVs) [55]. ASVs were affiliated to taxa using the SILVA database release 138 [56] using the *feature-classifier* plugin and *classify-sklearn* module [57, 58]. Sequences assigned as Eukaryota, Archaea, Mitochondria, Chloroplast and unassigned were removed from the dataset then the sample dataset was rarefied to a list of 6978 sequences.

Alpha- and beta-diversity analyses were performed using the *phyloseq* [59], *vegan* and *RVAideMemoire* packages in R. Linear mixed models (LMMs) were used to compare species richness among the five treatments and the three replicate aquaria within each treatment at d28, using the *MuMln* [60] and *lmerTest* [61] R packages. We adapted the LMMs to the non-independency of

individuals within each replicate, and defined the replicates as random effects and the five treatments as fixed effects, according to the formula $Y \sim \text{treatment} + (1 | \text{treatment} : \text{aquarium})$. Principal coordinates analyses (PCoA) based on weighted and unweighted UniFrac distances were performed to examine the dissimilarity of bacterial composition between groups. Among- and within-group variance levels were compared using PERMANOVA (999 permutations) and PERMDISP (999 permutations), respectively. Differentially abundant taxa across groups were identified using the linear discriminant analysis (LDA) effect size (LEfSe) tool [62] in the Galaxy workspace [63] (http://huttenhower.sph.harvard.edu/lefse/). Default parameters were applied using a LDA score threshold of 3.5 and the multi-class strategy (one-against-all).

## Microbiome-metabolome integrative analysis

The integration of datasets, i.e., the area-under-the-peaks in metabolite profiles and the ASV counts describing the bacterial communities in the same sample, was performed using the *mixOmics* package in R. Pareto scaling was applied on the metabolome data, and a centered log-ratio transformation then a pre-filtering keeping only abundant ASVs, (i.e., representing at least 1% of the reads in at least one sample), were applied on the microbiome data. Following unsupervised analyses on each dataset, completed to explore and visualize any similar changes according to treatments, the integration was carried out. A supervised Projection to Latent Structures Discriminant Analysis (PLS-DA) was performed using DIABLO (Data Integration Analysis for Biomarker discovery using Latent cOmponent) [64], enabling to identify highly-correlated variables (metabolites and ASVs) also discriminating the different treatments. The integration of both datasets was realized using the full-weighted design matrix and the *block.plsda* function implemented in *mixOmics*. The *plotDiablo* function enabled to check the well maximized covariation between datasets by displaying a Pearson correlation score. Then, relevance networks displaying the most discriminant covariates (metabolites, ASVs) were produced using the *network* function with Pearson correlation cut-offs [65].

## Metagenomic sequencing and analysis

Shotgun metagenome sequencing was performed on DNA from 10 gut samples collected at d28, 5 in treatment d28_0 and 5 from d28_100 (Illumina HiSeq, 2 × 150 bp, GenoScreen, France). Reads corresponding to animal sequences were identified by aligning each dataset against *Oryzias latipes* available at the NCBI, using *BBMap/bbsplit* [66], and discarded. Remaining reads from each sample were assembled using *metaSPAdes* with default
parameters [67]. Scaffolds were first taxonomically annotated using Contig Annotation Tool (CAT) [68] and Kaiju [69] allowing to detect sequences from *O. latipes* retroviruses and *Microcystis* genome which were then discarded. All scaffolds were clustered using MyCC [70] (k-mer size = 4, minimal sequence size = 1000) and bins were taxonomically annotated using Bin Annotation Tool [68]. Completeness of bins was assessed using CheckM [71]. Relative abundance of bins in each sample were also determined using BBMap. Significant bins between the two treatments were determined using Wilcoxon rank-sum test. Finally, coding sequences predicted by Prodigal [72] were functionally annotated using eggNOG-emapper [73]. Resulting KEGG annotations were used as input to MinPath [74] in order to obtain the MetaCyc pathway information, and to Omixer-RPM [75] to specifically characterize gut metabolic modules (GMMs) as firstly described in [76], with pathway coverage threshold ≥ 80%.

## Results

### Monitoring of experiments

*Oryzias latipes* fish were maintained in suitable and stable conditions, and no fish died during the whole experimentation (Fig. 1 and Dataset S2). During the first 28 days long exposure (d0–d28), concentrations of *Microcystis aeruginosa* in the three treatments (d28_1, d28_10, and d28_100) corresponded to expected levels, 1.0 ± 0.2, 10.0 ± 0.8, and 100.2 ± 7.8 µg.L⁻¹ *Chla*, respectively. Microcystin levels in water were 0.4 ± 0.0 and 10.4 ± 2.1 µg MC-LR eq.L⁻¹ in the d28_10 and d28_100 treatments, respectively, while microcystin was below detection level (< 0.15 µg.L⁻¹) in d28_1. During the second *M. aeruginosa* exposure (d33–d39), fish were exposed to 102.8 ± 2.9 µg.L⁻¹ *Chla* and 11.2 ± 2.8 µg MC-LR eq.L⁻¹. The metabolic content of *Microcystis aeruginosa* cultures was examined (see Supplementary Note 2, Dataset S3, and Table S1). Histopathological analyses did not reveal noticeable visual differences in fish liver tissue, with little to no carbohydrate reserves, and no lipofuscin and macrophagic hemosiderin.

### Diversity and composition of the gut bacterial microbiota after 28 days of exposure

At d28, gut communities display between 42 and 219 ASVs with higher average bacterial richness (136 ± 56

ASVs) and evenness (0.556) reported in fish guts exposed to d28_Z8, and lower average bacterial richness (82 ± 25 ASVs) and evenness (0.391) in the treatment d28_0 (Dataset S4). Gut-associated species richness is highest in the d28_Z8 treatment (LMM, *p* < 0.05), while no differences are observed among replicates within each treatment (LMM, *p* > 0.05). The Shannon index increases slightly with *Microcystis* concentration (d28_1: 2.09, d28_10: 2.14, d28_100: 2.28). Visual comparisons on individual plots of principal coordinates analyses (PCoA) based on the weighted or unweighted UniFrac distances suggest changes occur in terms of both abundances as well as community membership when fish are exposed to *M. aeruginosa* (d28_1, d28_10, d28_100) or d28_Z8 compared to d28_0 (Fig. 2a, b). Bacterial community composition appears different among treatments (PERMANOVA, weighted UniFrac, *p* < 0.001), notably between d28_100 and the other four treatments (*p* < 0.01), as well as between d28_0 and d28_Z8, d28_1, and d28_100 (*p* < 0.02). Differences between d28_Z8, d28_1, and d28_10 are not significant (*p* > 0.23) indicating comparable effect of a Z8 nutrients enrichment and lower *M. aeruginosa* doses. In addition, levels of variance in the different treatments are not significantly different (PERMDISP, *p* > 0.19).

In treatment d28_0 (water-exposed control), Fusobacteriota (50 ± 16%), Firmicutes (28 ± 15%), Proteobacteria (10 ± 6%), and Bacteroidota (5 ± 3%) dominate the microbiota, altogether representing 92.8 ± 6.9% of the reads. Relative abundances of these phyla are different among other treatments. Notably, Firmicutes are far less abundant in d28_100 (3 ± 2%) compared to other treatments (d28_0: 28 ± 15%, d28_Z8: 28 ± 15%, d28_1: 41 ± 14%, d28_10: 36 ± 12%) (Fig. 2c).

Many ASVs display significant differences across treatments in their relative abundances (Fig. 2e and Figure S2a). We focus on the six dominant ASVs, i.e., representing at least 10% of reads in at least one sample, accounting from 11 to 72% of the reads. ASV1349, ASV1564, and ASV2363 exhibit lower relative abundances in d28_0 (average below 0.7%). ASV1349 (*Flavobacterium*) and ASV1564 (*Aeromonas*) are more abundant in d28_Z8 (7.6 ± 10% and 3.8 ± 3.5%, respectively) and d28_100 (4.2 ± 4.1% and 6.2 ± 2.9%, respectively). ASV2363 (*Reyranella*) is also more abundant in

(See figure on next page.)

**Fig. 2** Changes in the composition of gut bacterial communities after 28 days of exposure. **a, b** PCoA using the weighted (a) or unweighted (b) UniFrac distances on fish gut bacterial community in the five treatments. **c** Relative abundance of bacterial phyla across treatments. **d** Relative abundance of Firmicutes members (family level). Firmicutes are mainly represented by the Erysipelotrichaceae family, and a single ASV (ASV1662). Colored horizontal bars represent the different treatments (see Fig. 1). **e** Significant ASVs from the linear discriminant analysis (LDA) effect size (LEfSe) with a LDA score above 3.5, and their relative abundances across d28 treatments. The colored boxes on the y-axis represent the treatment where each ASV is most abundant; dots represent the average relative abundance; line spreads over the range of observed values. Only dominant ASVs, i.e., representing at least 10% of reads in at least one sample were further considered
a

PCoA (wU) on the gut microbiota at d28

<figure><img src="image_4.png" /></figure>

b

PCoA (unwU) on the gut microbiota at d28

<figure><img src="image_5.png" /></figure>

c

<figure><img src="image_6.png" /></figure>

Phylum

<figure><img src="image_7.png" /></figure>

d

<figure><img src="image_8.png" /></figure>

Family

<figure><img src="image_9.png" /></figure>

e

<figure><img src="image_10.png" /></figure>

Treatment

<figure><img src="image_11.png" /><figcaption>**Fig. 2** (See legend on previous page.)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 3** Heatmap of GMM presence/absence in bins obtained from the medaka microbiome. GMM label colors represent the metabolic categories. Bin label colors correspond to their abundance variation between experimental conditions. Dendrograms are based on hierarchical clustering</figcaption></figure>

d28_100 (5.8 ± 11.1%). ASV1662 affiliated to the genus ZOR0006, is the main Firmicutes (96–100% of all Firmicutes reads), and as previously mentioned, is least abundant in d28_100 compared to other treatments (Fig. 2d). ASV2042 (Romboutsia), dominant in a single sample (17.8%) while below 3.5% in all others, was not further considered. Finally, ASV1644 (unassigned Bacteria) displays similar abundances in treatments d28_0 (5.6 ± 6.0%) and d28_100 (7.2 ± 9.7%) but is much less abundant in treatment d28_1 (0.7 ± 1.2%). Aside from the dominant ASVs that vary, ASV1620 (Cetobacterium) displays non-significant differences in abundance among treatments, representing 48 ± 16% of reads in d28_0 and average 27% to 46% of reads in other treatments.

## Metagenome-based comparison of gut communities in d28_0 and d28_100

Metagenome datasets obtained from 5 fish guts from d28_0 and 5 others from d28_100 yielded negligible amounts of Archaea and non-fish Eukaryota sequences. Among the 32 bacterial bins obtained, 8 were found to be more abundant after treatment (*p* < 0.05, Wilcoxon rank-sum test) (Figure S2b), including seven which vary

concordantly with ASVs sharing the same taxonomic annotation (Flavobacterium, Aeromonas, Gemmobacter, Rhizobiales, Figure S2a, b). Three bins showed a decreasing abundance after treatment, including bin24, affiliated to the Firmicutes, with a majority of coding sequences associated to ZOR0006 (58.06%), the genus corresponding to aforementioned ASV1662.

Gut Metabolic Modules were also investigated, with a total of 93 found within the 32 bins (Fig. 3). Each bin possesses a specific set of GMM involved in central metabolism and degradation of organic molecules or cobalamin (vitamin B₁₂) biosynthesis. A hierarchical clustering based on the GMM presence/absence failed to gather bins following their abundance variation behavior between treatments. Interestingly, cobalamin biosynthesis is a largely shared metabolic feature in the medaka gut microbiome, since 17 bins harboured this GMM, including bin29 associated with *Cetobacterium*. ZOR0006 related bin24 carries a limited set of GMM mostly involved in Carbohydrates degradation. It also harbours the lactate pyruvate interconversion function that is only shared with the *Cetobacterium*-affiliated bin29. When considering
<figure><img src="image_13.png" /><figcaption>**Fig. 4** Modification of the gut metabolome composition associated with bacteria changes. a Principal component analysis (PCA) representing gut metabolite profiles at d28. b Relevance network analysis illustrating the most correlated metabolites (in white) and ASVs (in orange) discriminating among the five treatments. Only variables (ASVs, metabolites) with correlation values above ± 0.695 are displayed. Pearson correlations between covariate metabolites and ASVs are represented by coloured segments (blue: negative association, red: positive association)</figcaption></figure>

genomic investment in the three major organic
substrate types degradation, namely amino-acids,
carbohydrates and lipids, bins were split into two sub-
populations (Figure S3): a first group gathering bac-
teria with a generalist strategy (in the center of the
ternary plot), and a second one that is highly special-
ized in carbohydrates degradation, the latter includ-
ing the three bins showing an abundance decrease in
exposure condition (bin11, bin24, and bin22). Beyond
GMM, only 0.3% (6) of all enzymes are specific to
these three bins (Figure S2c), including one involved
in lactose degradation (KO2788) but not essential for
this pathway. On the other hand, 4.3% (95) of enzymes
are specifically found in the 8 bins that are enriched
in d28_100 (namely bin9, bin10, bin12, bin16, bin19,
bin23, bin25, bin26). These represent 95 specific
enzymes, including 14 involved in the biosynthesis of
secondary metabolites, 10 in porphyrin and chloro-
phyll metabolism, and 6 in the biosynthesis of cofac-
tors. The analysis of clusters of orthologous groups
(COGs), regardless of their annotation status, reveals
that 963 of them (4.4% of the total COGs) are specific
to the three bins that are the most abundant in d28_0,
and that 2316 others (10.7% of the total gene families)
are specific to the 8 bins that are the most abundant

in d28_100. However, one should notice that most of
these COGs carry no known functional annotation,
compromising our understanding of the influence of
gene content variations on microbiota functioning
based on metagenomic assumption.

Gut metabolome variation and integration with gut
microbiota composition after 28 days of exposure

A total of 1674 metabolites were detected across gut
samples at d28. The PCA analysis separates the metab-
olite profiles of fish exposed to the different *M. aerugi-
nosa* levels (d28_1, d28_10, d28_100) along the first axis,
while the second axis mostly separates d28_0 from other
treatments (Fig. 4a). The gut metabolite composition is
different among treatments (PERMANOVA, euclidean
distance, *p* < 0.001), between d28_0 and the three *M.
aeruginosa* treatments (*p* = 0.01), and between d28_100
and all treatments (*p* < 0.03) but d28_10 (*p* = 0.19). Treat-
ments d28_Z8, d28_1, and d28_10 do not display dif-
ferences (*p* > 0.3). No significant variation is observed
among treatments in livers (*p* = 0.13) (Figure S4a). In
muscles, differences occur (*p* < 0.001) especially between
d28_0 and both d28_Z8 and d28_1 (*p* = 0.01), and
between d28_100 and all other treatments (*p* = 0.01)
except d28_0 (*p* = 0.33) (Figure S4b).
A joint analysis of gut metabolites and bacterial communities was performed, as described by Singh and colleagues [64]. The resulting combined dataset consists of two matrices of the 70 abundant ASVs, i.e., representing at least 1% of reads in at least one sample, and the 1674 metabolites, sufficiently well correlated together (76%) to explore associations between ASVs and metabolites. The relevance network illustrates correlations between the most highly associated ASVs and metabolites discriminating among the five treatments. Four ASVs (ASV1349, ASV1564, ASV475, ASV1995) appear negatively correlated with several metabolites and are the least abundant in d28_0 (Fig. 2e and Fig. 4b). Differently, the two other ASVs, namely ASV1543 (*Epulopiscium*) and ASV2350 (Barnesiellaceae), are positively correlated with few metabolites and are the most abundant in d28_0 (Fig. 4b). Unfortunately, due to the general lack of proper metabolite databases [77], only few of the metabolites presenting high correlation with abundant ASVs could be annotated, and these mainly corresponded to amino acids or small peptides (Dataset S5).

## Diversity and composition of the gut bacterial microbiome and metabolome after depuration (d33) and a second exposure (d39)

After d28, fish from all 5 treatments were transferred to clean water (Fig. 1). At d33, some specimens were sampled, while others were transferred to a second exposure to *M. aeruginosa* (100 µg.L⁻¹ Chla), then sampled at d39. The species richness of gut-associated communities decreases between d28 and d33 (over 82 versus 50 ± 24 ASVs), then increases at d39 (76 ± 22 ASV). The PCoA discriminates between gut communities from fish exposed to 100 µg.L⁻¹ Chla (d28_100 and d39) and other treatments (Figure S5a). Community compositions differ among treatments (PERMANOVA, weighted UniFrac, *p* < 0.001). Treatment d33 differs from all other treatments (*p* < 0.05) except for d28_0 (*p* = 0.063); in other words, bacterial communities at d33 are mostly similar to those observed in the d28_0 treatment. Interestingly, *Microcystis*-affiliated ASVs are almost absent (below 0.2% of reads) from the water and absent from gut samples at d33. Communities at d39 differ from all other treatments (*p* < 0.05), except for d28_100 (*p* = 0.063); indicating that bacterial communities at d39 are similar to those from the d28 treatment exposed to the highest bloom intensity. Interestingly, *Microcystis*-affiliated ASVs are abundant at d39 in the water (9.8 to 18.1% of reads) and occur in the gut samples (up to 4.6%). Variance differs among treatments (PERMDISP, *p* < 0.01), particularly between d33 and all other treatments (*p* < 0.05), and between d39 and treatments d33 and d28_Z8 (*p* < 0.05), with d33 and d39 displaying lower variance. Some phyla are present in

similar abundances at d33 and d39, including Fusobacteriota (62 ± 8% and 66 ± 14%, respectively), Proteobacteria (12 ± 6% and 19 ± 9%) and Bacteroidota (9 ± 3% and 7 ± 3%, Fig. 5a, b). Firmicutes, again mostly consisting of ASV1662, are present at d33 and almost absent at d39 (15.4 ± 6.4% versus 0.8 ± 1.8%). Other significant changes have been observed in specific ASVs between d33 and d39 (Figure S6a) among which seven were already observed to vary between the different d28 treatments (Fig. 2e and Figure S2a). Interestingly, ASV1620 (*Cetobacterium*), that was found abundant and stable among the different d28 treatments, is still remarkably stable at d33 and d39 (58.8 ± 7.4% and 61.7 ± 13.2%).

The gut metabolite profiles are significantly different among treatments (PERMANOVA, euclidean distance, *p* < 0.001, Figure S5b), especially between d33 and all treatments (*p* < 0.05) except d28_Z8 and d28_1, as well as between d39 and all others (*p* < 0.05) except d28_Z8 and d28_10. Contrary to bacterial microbiota composition, the d33 metabolomes are overall different from d28_0, and those of d39 are different from d28_100. A good correlation (81%) was observed between the gut ASVs and the metabolites discriminating between d33 and d39 when investigating the two matrices, containing 31 abundant ASVs and 1674 metabolites, respectively. The relevance network notably revealed two ASVs, belonging to genus Vibrio (ASV1406, ASV2161), that are negatively correlated with numerous metabolites, while two ASVs, belonging to genus *Reyranella* (ASV1030, ASV2363) appear positively correlated with numerous metabolites (Figure S6b). No differences occur among treatments in liver samples (*p* > 0.05, Figure S5c). In muscles, metabolomes at d33, d39, and d28_100 appear similar (*p* > 0.10), but different from all other treatment groups (Figure S5d).

## Occurrence of dominant ASVs in compartments other than gut

Dominant ASVs were searched for in fish food, *M. aeruginosa* culture, water, biofilm, and faeces (Figure S7). ASV1620 (*Cetobacterium*) is most abundant in guts (44.3 ± 20.4%) then in faeces (28.6 ± 15.3%), congruent with its best matches in the database which are fish gastrointestinal bacteria (zebrafish, Nile Tilapia). ASV1349 (*Flavobacterium*) and ASV1564 (*Aeromonas*) are most abundant in faeces (8.8 ± 9.6% and 7.8 ± 2.8%, respectively), then guts (2.1 ± 4.4% and 3.4 ± 2.6%, respectively) and biofilms (2.0 ± 3.3% and 3.4 ± 2.5%, respectively). This is congruent with their respective BLASTN hits with bacteria from fish intestinal tracts or aquatic environments. ASV1644 (unassigned Bacteria) is also most abundant in faeces (10.7 ± 7.4%), then in guts (3.7 ± 6.6%), but could not be assigned taxonomically.
<figure><img src="image_14.png" /><figcaption>**Fig. 5** Depuration and a second *Microcystis* bloom impact the gut microbiota community. Relative abundance of gut bacterial phyla at d33, after 4 days in water following the d28 treatment (a), and at d39 after 5 additional days exposed to *M. aeruginosa* (b). Colored horizontal bars represent the different imposed treatments: blue = water, dark green = 100 µg.L⁻¹ Chla. Individuals sampled from a given aquarium (blue then green box below histograms) are organized in columns</figcaption></figure>

The Firmicutes ASV1662 (ZOR0006) mainly occurs in guts (19.3 ± 16.9%) and is slightly abundant in faeces (1.2 ± 1.4%), in accordance with the habitat of its closest matches, namely fish intestinal bacteria. Finally, ASV2363 (Reyranella) is most abundant in biofilms (4.4 ± 8.7%) compared to other compartments (water 2.2 ± 6.3%; feces 1.3 ± 2.7%; guts 1.2 ± 4.1%), and is related to various environmental bacteria.

Discussion

Results from the 28-day exposure indicate that *Microcystis aeruginosa* blooms modify the fish gut bacterial microbiota compositions and have different effects on different taxa. Firmicutes, a phylum commonly found in gut bacterial communities of vertebrates and very likely implied in host metabolism processes [78–80], appear particularly sensitive as they decrease sharply upon exposure to the highest bloom intensity. Firmicutes were largely represented by a single bacterium (Erysipelotrichaceae\_ZOR0006) that was almost absent

outside of gut samples, suggesting being an indigenous and resident symbiont of *O. latipes* gut [81]. Kaakoush and colleagues [82] have previously discussed the central role of Erysipelotrichaceae in host metabolism and health in relation to diet specificity. According to our results based on the functional classification proposed by Vieira-Silva, ZOR0006 displays a limited repertoire of metabolic modules compared to other bins such as *Cetobacterium* (16 vs. 54 GMM, respectively). Furthermore, most of the functions present in ZOR0006 are also encoded by other bacteria in the microbiome, with the notable exception of the lactate-pyruvate interconversion module, a rare yet potentially important function otherwise only found in *Cetobacterium*. Indeed, lactate homeostasis is thought to be essential in the gut [83]. Although high concentrations are associated with inflammatory bowel disease in humans [84, 85], lactate has been shown to inhibit the growth of pathogenic bacteria [86], whereas pyruvate stimulates it [87]. In addition, lactate has recently been shown to be involved in the repair of the gut epithelium
[88]. Symbiont abundance changes could thus affect the gut epithelium function as a barrier. Interestingly, the degradation activities of ZOR0006, like those of the other two bacteria whose abundance decreased after exposure, are totally oriented towards carbohydrates, suggesting a specialized degradation metabolism. The ZOR0006 genus has already been observed to decrease in fish gut microbiota after exposing zebrafish to high concentrations of antibiotics or fungicides [89, 90]. Interestingly, the drop of ZOR0006 abundance observed in this study occurs at a level between 10 and 100 µg.L⁻¹ Chla. Together with the observation of the greater influence of the highest bloom condition on the whole gut community, this suggests the existence of a cyanobacterial bloom threshold above which the gut microbiota composition is particularly altered. On the other hand, some abundant gut bacteria appear stable throughout the bloom, and even during the depuration and second exposure, the most remarkable being the bacterium affiliated to *Cetobacterium* (Fusobacteria). This genus is generally associated with healthy fish microbiota and notably contributes to host health as a cobalamin (vitamin B12) producer [91–93]. In addition to the confirmation of cobalamin biosynthesis ability, our analysis proved *Cetobacterium* to possess the largest set of GMM, a functional feature that could suggest a highly significant role as a symbiont. *Cetobacterium* is most likely a fish gut resident and has previously been reported stable in medaka fish upon exposure to pure microcystin-LR and cell extract of the *Microcystis* strain [35], supporting its maintenance through the exposure to cyanobacterial bloom and respective metabolites. Finally, relative abundances of some other bacteria increase during the *Microcystis* bloom. These could include transient gut bacteria that originate from the environment and proliferate once established in the gut, one example being *Reyranella* which also occurs in biofilms. These may also include rare gut resident taxa that can take advantage of peculiar conditions to proliferate. In our study, some of these bacteria, corresponding to *Flavobacterium*, *Aeromonas*, and *Shewanella*, are common inhabitants of fish guts or the environment, or potential fish gut pathogens according to the literature [94–96]. A similar increase of opportunistic bacteria was recently documented in guts of zebrafish exposed 96h to *M. aeruginosa* [10]. As previously shown for several metabolite mixtures from cyanobacterial cell extracts [35], exposure to whole cyanobacterial cells thus has a major impact on gut community compositions, indicating that fish microbiota might be impacted during the bloom as well as after bloom senescence which causes the release of the cyanobacterial cell contents into the surrounding water [97]. However, only the highest concentration of the bloom and/or its cell extract, that was explored in the present

study, could induce most evident microbiota changes, suggesting that the microbiota responsiveness might be dose dependent. The cyanobacterial strain PMC 728.11 contains various secondary metabolites that may be responsible for variations in fish gut microbiota. Among those already identified and specifically those potentially produced by the strain [35, 98], cyanopeptides, such as aerucyclamides and bacteriocins, are thought to exhibit potent antimicrobial or cytotoxic bioactivities [99, 100], and could directly impact the microbiota during *Microcystis* cell digestion into the intestine lumen.

Metagenome-based investigation confirms the variations of taxa abundances, but very few taxa-specific known (i.e., annotated) functions were identified for the bacteria that displayed the highest abundance variations, including for the Firmicutes ZOR0006 which displays the smallest set of GMM. This is at first consistent with the common claim that a change in microbiota composition does not necessarily imply a change in the functions as estimated by gene content [101, 102]. However, similar genes are not necessarily expressed in the same way in different bacterial taxa, and thus phenotypes might still change dramatically despite the potential occurrence of similar gene contents. Indeed, variations in metabolite profiles reveal obvious functional variations induced by the different treatments, in particular changes in gut metabolite composition associated with increasing bloom concentrations. These changes could thus be related to unannotated COGs, much more numerous among the differentially abundant bacteria compared to annotated KO, or could result from variations in the expression of metabolic pathways in the holobiont. Whatsoever, the correlation observed between some dominant bacteria and many metabolites supports that the response of gut microbiota composition and that of the holobiont's gut metabolome are linked. Gut community disruption is associated with metabolic changes, especially in the presence of higher bloom levels. This could lead to a potential dysbiotic state, however defining eubiotic versus dysbiotic state is not straightforward, as discussed elsewhere [103, 104], and will require exploring various additional holobiont physiological variables. Altogether, the crosstalk between gut bacteria and metabolites remains difficult to characterize. For example, the presence of pathogens enhanced by bloom-induced toxicological impairs of fish physiology could secondarily affect the holobiont gut metabolism. Alternatively, the variation of the quantity of some gut metabolites could allow the development of opportunistic pathogens, as reported in mammals [105, 106]. Interestingly, variations of metabolite composition in muscles and livers remains limited in comparison to those observed in the gut, even after 28 days, implying that the gut metabolome compartment is more responsive to the *Microcystis aeruginosa* exposure. This observation is congruent with the fact that the gut is
more directly exposed to the surrounding environment, through ingestion of water that contains various organisms and their associated compounds, and emphasizes the relevance of investigating the gut microbiota, as the gut epithelium sits at the interface between the host and its environment [32, 33]. Changes in microbiome and overall holobiont functions seem to be deeply linked, despite that the respective contribution of the host and the bacteria community cannot yet be disentangled due to the current impossibility to link a metabolite to its organism of origin within the holobiont and will require further dedicated functional analysis [107].

After the first 28-day bloom, depuration in clear water seems to have restored, to a certain extent, gut microbiota compositions. Indeed, compositions at d33 tended to resemble those of naïve specimens, not previously exposed to *M. aeruginosa* or to Z8 (d28_0). This indicates that the gut community composition is resilient, even over a relatively short time. The following second high-intensity bloom on the other hand yielded bacterial communities resembling those of specimens exposed to the first high-intensity bloom (d28_100), including the disappearance of Firmicutes and the stability of *Cetobacterium*. This indicates that a short duration bloom (5 days) already has a strong effect [37]. Gut communities thus quickly respond to the presence or absence of *M. aeruginosa* and associated bioactive compounds. Effects of the first and the second bloom are however not strictly identical. Firmicutes for example show reproducible behavior, disappearing upon first and second bloom, and seem to quickly re-establish post-bloom, suggesting they are not resistant, but remarkably resilient. However, the relative abundance of other taxa (including *Flavobacterium*, *Reynanella*, *Shewanella*), whose abundance increased upon the 28 days bloom, did not increase during the second 5-day exposure. It remains here difficult to compare variations of most non-dominant bacteria between d33 and d39 because of the substantial inter-individual variation of microbiota composition, combined with the limited sample size in terms of individual number per condition, compared to d28 conditions. Stress associated with specimen handling, the transfer of fish to newly cleaned aquaria after 28 days, has previously been referenced as a moderate, but possible, source of stress [108] and also could be involved to a certain extent in some of the differences. Contrary to gut community compositions, the metabolite composition in fish guts after the depuration and the second exposure did not tend towards those observed in unexposed and highest bloom at d28, respectively. So, metabolite compositions do not linearly follow the trends found in community compositions. This also suggests that the depuration and the second exposure may have been too prompt to induce major shifts in the

gut metabolome, as these compartments may present different kinetics. Thus, the dynamics of gut microbiota and metabolome are very likely not identical, and the microbiota composition may somehow respond to changes faster than the metabolome. This would suggest that longer blooms occurring in nature may have more effect and functional consequences for the holobiont homeostasis compared to shorter blooms [109, 110].

## Conclusions

Overall, this study emphasizes that cyanobacterial blooms have the potential to alter fish gut microbiota composition and holobiont functions. Additionally, the drop of Firmicutes abundance induced by *M. aeruginosa* bloom which threshold level would be comprised between 10 and 100 µg.L⁻¹ *Chla*, together with the observation of a greater influence of higher bloom condition (100) on the overall gut community, support the existence of a notable tipping point for the responsiveness of gut microbiome to cyanobacterial bloom intensity. This finding could have important eco(toxico)logical consequences, as these levels are commonly reached in natural ecosystems during typical *Microcystis* bloom episodes worldwide, suggesting that destabilization of fish gut communities might be a very common event [42]. In nature, cyanobacterial blooms are nowadays becoming increasingly frequent [3], and sometimes persistent among seasons. Freshwater fish, especially those living in shallow water ponds, face numerous bloom episodes of varying durations during their lifetime. Consequences of iterative (up to chronic) exposures on the holobiont should be explored. Indeed, successive blooms without enough recovery time could induce a cumulative drift of the gut microbiota and metabolome, leading to suboptimal states that may lead to host health impairment [111]. This phenomenon could be of major significance for fish health in eutrophic natural ecosystems as well as aquaculture ponds where cyanobacteria often proliferate. Testing the effects of a single live *M. aeruginosa* strain administrated by simple balneation on a model teleost fish is a first step towards understanding bloom effects on fish microbiota and holobiont health. In the future, particular attention should also be paid to the natural diversity of cyanobacterial within blooms, as successive blooms often involve different species that produce different metabolite cocktails [112].

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s40168-023-01558-2.

<footnote>Additional file 1: Supplementary Note 1. Dataset pre-processing and bacterial community comparison within the entire sample set. **Supplementary Note 2.** Metabolic content of *Microcystis aeruginosa* cultures. **Figure S1.** (a-c) Alpha-diversity metrics (species richness, Shannon,</footnote>
evenness) in the different sample types. Means are represented by a red dot. (d) Principal coordinates analysis (PCoA) representing bacterial communities from the different compartments using the weighted Uni-Frac distance. **Figure S2.** (a-b) Differentially abundant fish gut bacteria between treatments d28_0 (blue) and d28_100 (dark green) based on 16S rRNA reads (a) and shotgun metagenome sequencing (b). The absence of 16S rRNA sequences in metagenome bins does not allow to directly connect ASVs with bins, however the coloured taxon names represent bacteria with the same taxonomic affiliation. Lowest assigned taxonomic levels are displayed with the level specified by a single letter. (a) Discriminant ASVs based on the Linear discriminant analysis (LDA) effect size (LEfSe) with a LDA score above 3.5 and their variations in relative abundance in fish guts. The coloured boxes represent the treatment where each ASV is found the most abundant. The dots represent average relative abundance, the lines spread over the range of observed values. Only dominant ASVs (see text), were considered and further investigated. (b) Taxonomic affiliations (using the BAT method) of bins displaying significantly different relative abundances based on the Wilcoxon rank-sum test (p < 0.05) and their relative abundance across fish gut samples. The CAT method was used to affiliate bin10 (Rhodospirillales) at a lower taxonomic level, i.e. *Reyranella massiliensis*. The treatment where each bin is significantly more abundant is represented by the coloured box. (c) Overlap of KO counts among three groups: KO from bins significantly more abundant in d28_0 (blue) or in d28_100 (dark green) (p <0.05, Wilcoxon rank-sum test), or KO from bins non-significantly different between the two treatments (beige) (p > 0.05). **Figure S3.** Triplet representation of bin contributions in the degradation potential of the whole microbiome, defined as the fraction of GMM coded by a bin in each of the three major degradation types, amino acids, lipids and carbohydrates. Bin colors correspond to their abundance variation between experimental conditions. **Figure S4.** Principal component analyses (PCA) illustrating the metabolite composition in fish livers (a) and muscles (b) from the five different treatments during 28 days. **Figure S5.** Comparisons of the composition of bacterial communities (a) and metabolite profiles (b) in fish guts, and the composition of metabolite profiles between the three different sampled organs, guts (b), livers (c) and muscles (d). The three organs were either exposed long-term during 28 (d28) days (full circles or diamonds) or short-term during 4 (d33) or 5 (d39) days (open circles or diamonds). (a) Principal coordinates analysis (PCoA) on weighted UniFrac distance illustrating bacterial composition in fish guts. (b-d) Principal component analyses (PCA) representing metabolite profiles in fish guts, livers and muscles. **Figure S6.** (a) Significant ASVs from the LEfSe analysis, differentially abundant between d33 gut samples (in water) and d39 gut samples (in 100 µg.L⁻¹ Chla). The coloured boxes represent the treatment where ASVs are found more abundant. Dots and lines indicate mean and range of values, respectively. The displayed taxonomic affiliations correspond to the lowest assigned level, displayed using the first taxon letter. The two ASVs (ASV2363, ASV1662) underlined are also found differentially abundant between d28_0 and d28_100. (b) Relevance network analysis representing the most correlated ASVs and metabolites discriminating d33 and d39 gut samples. Only ASVs (in orange) and metabolites (in white) associated with Pearson correlation scores above ±0.714 are displayed. Coloured segments represent Pearson correlation values, either positive (red) or negative (blue). **Figure S7.** Relative abundance of dominant ASVs across compartments other than gut (culture, water, biofilm, gut and faeces). ASVs were searched for in the fish food sample discarded at the rarefaction step, but any ASVs were found (only 1% of reads for ASV1620). **Table S1.** Composition of annotated metabolites of the PMC 728.11 *Microcystis aeruginosa* strain. Percentages represent the proportions of each annotated cluster on the total annotated clusters. Abbreviation: AA = Amino Acids.

**Additional file 2: Dataset S1.** Sampling counts during the whole experiment.

**Additional file 3: Dataset S2.** Monitoring of abiotic and biotic parameters during the experiment.

**Additional file 4: Dataset S3.** Annotations of metabolites from the *Microcystis aeruginosa* strain PMC 728.11.

**Additional file 5: Dataset S4.** 16S rRNA gene sequencing and quality filtered read counts, SRA database accession numbers, and measures of alpha-diversity metrics (species richness, Shannon, evenness).

**Additional file 6: Dataset S5.** Annotations of the most correlated metabolites with bacteria from fish guts at day 28.

## Acknowledgements

Kandiah Santhirakumar managed fish maintenance, Claude Yéprémian advised on cyanobacterial cultivation and Arthur Escalas helped with statistical analyses. We thank the Amagen platform for providing medaka fish, ENVA for histological data processing, and the PtSMB platform of MNHN for metabolomics.

## Authors' contributions

S.D. and B.M. conceived the study. A.G., C.D., S.D., and B.M. conceived the experiment. A.G. and C.D. conducted the experiment. S.D. and B.M. took part in the experiment. H.H. conducted histological data processing and analyses. A.G. and C.D. conducted molecular data processing. A.G., S.H., S.D., and B.M. analyzed data. A.G., S.D., and B.M. wrote the manuscript. S.H. edited the manuscript. All authors contributed and agreed on the contents.

## Funding

Research was funded by ATM 3M and 3M2 from the MNHN and by AcSymb project from the Institut Universitaire de France. AG's PhD fellowship was funded by the "Ecole doctorale ED227" (MNHN-SU).

## Availability of data and materials

16S rRNA gene sequencing and shotgun raw data are deposited into the GenBank SRA database under the BioProject PRJNA746242 (samples SAMN20080275 to SAMN20080438). The full R code and QIIME 2 script are available on Zenodo (https://doi.org/10.5281/zenodo.5922616).

## Declarations

### Ethics approval and consent to participate

Experiments were conducted according to best practices within the framework of an authorized program (APAFiS#19316-2019032913284201 v1) approved by Ethical Committee of the National Museum of Natural History (Paris) CEEA n°68.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

### Author details

¹UMR7245 Molécules de Communication et Adaptation des Micro-organismes. Muséum National d'Histoire Naturelle, CNRS, Paris, France. ²UMR1161 Virologie, École Nationale Vétérinaire d'Alfort, INRA - ANSES - ENVA, Maisons-Alfort, France. ³Institut Universitaire de France, Paris, France.

Received: 18 February 2022 Accepted: 26 April 2023

Published online: 16 May 2023

## References

1. Havens KE. Chapter 33: Cyanobacteria blooms: effects on aquatic ecosystems. In: Hudnell HK, editor. Cyanobacterial Harmful Algal Blooms State Sci Res Needs. New York: Springer; 2008. 733–47. https://doi.org/10.1007/978-0-387-75865-7_33

2. Pavagadhi S, Balasubramanian R. Toxicological evaluation of microcystins in aquatic fish species: Current knowledge and future directions. Aquat Toxicol. 2013;142–143:1–16.

3. Huisman J, Codd GA, Paerl HW, Ibelings BW, Verspagen JMH, Visser PM. Cyanobacterial blooms. Nat Rev Microbiol. 2018;16:471–83.
4. Escalas A, Catherine A, Maloufi S, Cellamare M, Hamlaoui S, Yéprémian C, et al. Drivers and ecological consequences of dominance in periurban phytoplankton communities using networks approaches. Water Res. 2019;163: 114893.

5. O'Neil JM, Davis TW, Burford MA, Gobler CJ. The rise of harmful cyanobacteria blooms: The potential roles of eutrophication and climate change. Harmful Algae. 2012;14:313-34.

6. Paerl HW, Otten TG. Harmful cyanobacterial blooms: causes, consequences, and controls. Microb Ecol. 2013;65:995-1010.

7. Malbrouck C, Kestemont P. Effects of microcystins on fish. Environ Toxicol Chem. 2006;25:72.

8. Singh S, Kate BN, Banerjee UC. Bioactive compounds from cyanobacteria and microalgae: an overview. Crit Rev Biotechnol. 2005;25:73-95 (Taylor & Francis).

9. Marie B, Huet H, Marie A, Djediat C, Puiseux-Dao S, Catherine A, et al. Effects of a toxic cyanobacterial bloom (*Planktothrix agardhii*) on fish: Insights from histopathological and quantitative proteomic assessments following the oral exposure of medaka fish (*Oryzias latipes*). Aquat Toxicol. 2012;114-115:39-48.

10. Qian H, Zhang M, Liu G, Lu T, Sun L, Pan X. Effects of different concentrations of *Microcystis aeruginosa* on the intestinal microbiota and immunity of zebrafish (*Danio rerio*). Chemosphere. 2019;214:579-86.

11. Saraf SR, Frenkel A, Harke MJ, Jankowiak JG, Gobler CJ, McElroy AE. Effects of *Microcystis* on development of early life stage Japanese medaka (*Oryzias latipes*): comparative toxicity of natural blooms, cultured *Microcystis* and microcystin-LR. Aquat Toxicol. 2018;194:18-26.

12. Sotton B, Paris A, Le Manach S, Blond A, Lacroix G, Millot A, et al. Global metabolome changes induced by cyanobacterial blooms in three representative fish species. Sci Total Environ. 2017;590-591:333-42.

13. Chorus I, Bartram J, editors. Toxic cyanobacteria in water: a guide to their public health consequences, monitoring, and management. London ; New York: E & FN Spon; 1999.

14. Ibelings BW, Havens KE. Cyanobacterial toxins: a qualitative meta-analysis of concentrations, dosage and effects in freshwater, estuarine and marine biota. In: Hudnell HK, editor. Cyanobacterial Harmful Algal Blooms State Sci Res Needs. New York: Springer; 2008. 675-732. https://doi.org/10.1007/978-0-387-75865-7_32.

15. Le Manach S, Sotton B, Huet H, Duval C, Paris A, Marie A, et al. Physiological effects caused by microcystin-producing and non-microcystin producing *Microcystis aeruginosa* on medaka fish: a proteomic and metabolomic study on liver. Environ Pollut. 2018;234:523-37.

16. Qiao Q, Le Manach S, Huet H, Duvernois-Berthet E, Chaouch S, Duval C, et al. An integrated omic analysis of hepatic alteration in medaka fish chronically exposed to cyanotoxins with possible mechanisms of reproductive toxicity. Environ Pollut. 2016;219:119-31.

17. Ibelings BW, Bruning K, de Jonge J, Wolfstein K, Pires LMD, Postma J, et al. Distribution of microcystins in a Lake Foodweb: no evidence for biomagnification. Microb Ecol. 2005;49:487-500.

18. Deng D-F, Zheng K, Teh F-C, Lehman PW, Teh SJ. Toxic threshold of dietary microcystin (-LR) for quart medaka. Toxicon. 2010;55:787-94.

19. Amado LL, Monserrat JM. Oxidative stress generation by microcystins in aquatic animals: Why and how. Environ Int. 2010;36:226-35.

20. Ernst B, Hitzfeld B, Dietrich D. Presence of *Planktothrix* sp. and cyanobacterial toxins in Lake Ammersee, Germany and their impact on whitefish (*Coregonus lavaretus* L.). Environ Toxicol. 2001;16:483-8.

21. Ibelings BW, Chorus I. Accumulation of cyanobacterial toxins in freshwater "seafood" and its consequences for public health: A review. Environ Pollut. 2007;150:177-92.

22. Wang AR, Ran C, Ringø E, Zhou ZG. Progress in fish gastrointestinal microbiota research. Rev Aquac. 2018;10:626-40.

23. Llewellyn MS, Boutin S, Hoseinifar SH, Derome N. Teleost microbes: the state of the art in their characterization, manipulation and importance in aquaculture and fisheries. Front Microbiol. 2014;5. Cited 2018 Oct 17. Available from: http://journal.frontiersin.org/article/10.3389/fmicb.2014.00207/abstract.

24. Egerton S, Culloty S, Whooley J, Stanton C, Ross RP. The gut microbiota of marine fish. Front Microbiol. Frontiers; 2018;9. Cited 2021 Jul 2. Available from: https://www.frontiersin.org/articles/https://doi.org/10.3389/fmicb.2018.00873/full.

25. Sullam KE, Essinger SD, Lozupone CA, O'Connor MP, Rosen GL, Knight R, et al. Environmental and ecological factors that shape the gut bacterial communities of fish: a meta-analysis. Mol Ecol. 2012;21:3363-78.

26. de Bruijn I, Liu Y, Wiegertjes GF, Raaijmakers JM. Exploring fish microbial communities to mitigate emerging diseases in aquaculture. FEMS Microbiol Ecol. 2018;94. Cited 2018 Nov 22. Available from: https://academic.oup.com/femsec/article/doi/https://doi.org/10.1093/femsec/fix161/4675208.

27. Nayak SK. Role of gastrointestinal microbiota in fish. Aquac Res. 2010;41:1553-73.

28. Ringø E, Zhou Z, Vecino JLG, Wadsworth S, Romero J, Krogdahl Å, et al. Effect of dietary components on the gut microbiota of aquatic animals. A never-ending story? Aquac Nutr. 2016;22:219-82.

29. Claus SP, Guillou H, Ellero-Simatos S. The gut microbiota: a major player in the toxicity of environmental pollutants? NPJ Biofilms Microbiomes. 2016;2:16003. https://doi.org/10.1038/npjbiofilms.2016.3.

30. Jin Y, Wu S, Zeng Z, Fu Z. Effects of environmental pollutants on gut microbiota. Environ Pollut. 2017;222:1-9.

31. Rosenfeld CS. Gut Dysbiosis in Animals Due to Environmental Chemical Exposures. Front Cell Infect Microbiol. 2017;7. Cited 2018 Oct 17. Available from: Available from:http://journal.frontiersin.org/article/10.3389/fcimb.2017.00396/full.

32. Evariste L, Barret M, Mottier A, Mouchet F, Gauthier L, Pinelli E. Gut microbiota of aquatic organisms: a key endpoint for ecotoxicological studies. Environ Pollut. 2019;248:989-99.

33. Duperron S, Halary S, Gallet A, Marie B. Microbiome-aware ecotoxicology of organisms: relevance, pitfalls, and challenges. Front Public Health. 2020;8. Cited 2021 Jan 28. Available from: https://www.frontiersin.org/article/10.3389/fpubh.2020.00407/full.

34. Duan Y, Xiong D, Wang Y, Dong H, Huang J, Zhang J. Effects of *Microcystis aeruginosa* and microcystin-LR on intestinal histology, immune response, and microbial community in *Litopenaeus vannamei*. Environ Pollut. 2020;265: 114774.

35. Duperron S, Halary S, Habiballah M, Gallet A, Huet H, Duval C, et al. Response of fish gut microbiota to toxin-containing cyanobacterial extracts: a microcosm study on the Medaka (*Oryzias latipes*). Environ Sci Technol. Lett. 2019;6:341-7.

36. Li J, Chen C, Zhang T, Liu W, Wang L, Chen Y, et al. µEvaluation of microcystin-LR absorption using an *in vivo* intestine model and its effect on zebrafish intestine. Aquat Toxicol. 2019;206:186-94.

37. Foucault P, Gallet A, Duval C, Marie B, Duperron S. Gut microbiota and holobiont metabolome composition of the medaka fish (*Oryzias latipes*) are affected by a short exposure to the cyanobacterium *Microcystis aeruginosa*. Aquat Toxicol. 2022;253: 106329.

38. Metcalfe CD, Gray MA, Kiparissis Y. The Japanese Medaka (*Oryzias latipes*): An *In Vivo* Model For Assessing The Impacts Of Aquatic Contaminants On The Reproductive Success Of Fish. In Impact Assessment of Hazardous Aquatic Contaminants Concept and Approaches. ed. S. S. Rao. Boca Raton: FL. Lewis Publishers. 1999;29-52.

39. Patyanna PJ, Davi RA, Parkerton TF, Brown RP, Cooper KR. A proposed multigeneration protocol for Japanese medaka (*Oryzias latipes*) to evaluate effects of endocrine disruptors. Sci Total Environ. 1999;233:211-20.

40. Sotton B, Paris A, Le Manach S, Blond A, Lacroix G, Millot A, et al. Metabolic changes in Medaka fish induced by cyanobacterial exposures in mesocosms: an integrative approach combining proteomic and metabolomic analyses. Sci Rep. 2017;7. Cited 2018 Oct 17. Available from: http://www.nature.com/articles/s41598-017-04423-z.

41. Wittbrodt J, Shima A, Schartl M. Medaka — a model organism from the far east. Nat Rev Genet. 2002;3:53-64.

42. Harke MJ, Steffen MM, Gobler CJ, Otten TG, Wilhelm SW, Wood SA, et al. A review of the global ecology, genomics, and biogeography of the toxic cyanobacterium. Microcystis spp. Harmful Algae. 2016;54:4-20.

43. Reichwaldt ES, Song H, Ghadouani A. Effects of the Distribution of a *Toxic Microcystis* Bloom on the Small Scale Patchiness of Zooplankton. PLoS One. 2013;8: e66674.

44. Binding CE, Pizzolato L, Zeng C. EOLakeWatch; delivering a comprehensive suite of remote sensing algal bloom indices for enhanced monitoring of Canadian eutrophic lakes. Ecol Indic. 2021;121: 106999.

45. Rippka R. [1] Isolation and purification of cyanobacteria. Methods Enzymol. Elsevier; 1988. 3-27. Cited 2021 Apr 14. Available from: https://linkin.nghub.elsevier.com/retrieve/pii/0076687988670042.
46. Yéprémian C, Catherine A, Bernard C, Congestri R, Elersek T, Pilkaityte R. Chlorophyll a Extraction and Determination. Handb Cyanobacterial Monit Cyanotoxin Anal. John Wiley & Sons, Ltd; 2016. 331–4. Cited 2021 Apr 14. Available from: http://onlinelibrary.wiley.com/doi/abs/10.1002/ 9781119068761.ch34.

47. Le Manach S, Duval C, Marie A, Djediat C, Catherine A, Edery M, et al. Global metabolomic characterizations of *Microcystis* spp. highlights clonal diversity in natural bloom-forming populations and expands metabolite structural diversity. Front Microbiol. Frontiers; 2019;0. Cited 2021 Jul 17. Available from: https://www.frontiersin.org/articles/10. 3389/fmicb.2019.00791/full.

48. Tugizimana F, Steenkamp PA, Piater LA, Dubery IA. A conversation on data mining strategies in LC-MS untargeted metabolomics: pre-processing and pre-treatment steps. Metabolites. 2016;6:40.

49. Rohart F, Gautier B, Singh A, Cao K-AL. mixOmics: An R package for 'omics feature selection and multiple data integration. PLOS Comput Biol. 2017;13:e1005752 (Public Library of Science).

50. Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: Community Ecology Package. R package version 2.5-7. 2020; Available from: https://CRAN.R-project.org/package=vegan.

51. Hervé M. RVAideMemoire: Testing and Plotting Procedures for Biostatistics. R package version 0.9-79. 2021; Available from: https://CRAN.R- project.org/package=RAideMemoire.

52. Terrat S, Christen R, Dequiedt S, Lelièvre M, Nowak V, Regnier T, et al. Molecular biomass and MetaTaxogenomic assessment of soil microbial communities as influenced by soil DNA extraction procedure. Microb Biotechnol. 2012;5:135–41.

53. Magoč T, Salzberg SL. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics. 2011;27:2957–63.

54. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7 (Nature Publishing Group).

55. Callahan BJ, McMurdie PJ, Holmes SP. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 2017;11:2639–43.

56. 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.

57. Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al. Scikit-learn: Machine Learning in Python. J Mach Learn Res. 2011;12:2825–30.

58. Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2's q2-feature-classifier plugin. Microbiome. 2018;6:90.

59. McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLOS One. 2013;8:e61217 (Public Library of Science).

60. Barton K. MuMln: Multi-model inference. R package version 1.43.17. 2020; Available from: https://CRAN.R-project.org/package=MuMln.

61. Kuznetsova A, Brockhoff PB, Christensen RHB. ImerTest Package: Tests in Linear Mixed Effects Models. J Stat Softw. 2017;82:1–26.

62. Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, et al. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12:R60.

63. Goecks J, Nekrutenko A, Taylor J. The Galaxy Team. Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome Biol. 2010;11:R86.

64. Singh A, Shannon CP, Gautier B, Rohart F, Vacher M, Tebbutt SJ, et al. DIABLO: an integrative approach for identifying key molecular drivers from multi-omics assays. Bioinformatics. 2019;35:3055–62.

65. González I, Cao K-AL, Davis MJ, Déjean S. Visualising associations between paired 'omics' data sets. BioData Min. 2012;5:19.

66. Bushnell B. BBMap: A Fast, Accurate, Splice-Aware Aligner [Internet]. Lawrence Berkeley National Lab. (LBNL), Berkeley, CA (United States); 2014 Mar. Report No.: LBNL-7065E. Available from: https://www.osti. gov/biblio/1241166-bbmmap-fast-accurate-splice-aware-aligner.

67. Bankevich A, Nurk S, Antipov D, Gurevich A, Dvorkin M, Kulikov A, et al. SPAdes: a new genome assembly algorithm and its applications

to single-cell sequencing. J Comput Biol J Comput Mol Cell Biol. 2012;19:455–77.

68. von Meijenfeldt FAB, Arkhipova K, Cambuy DD, Coutinho FH, Dutilh BE. Robust taxonomic classification of uncharted microbial sequences and bins with CAT and BAT. Genome Biol. 2019;20:217.

69. Menzel P, Ng KL, Krogh A. Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat Commun. 2016;7:11257 (Nature Publishing Group).

70. Lin H-H, Liao Y-C. Accurate binning of metagenomic contigs via automated clustering sequences using information of genomic signatures and marker genes. Sci Rep. 2016;6:24175 (Nature Publishing Group).

71. Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 2015;25:1043–55.

72. Hyatt D, Chen G-L, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:1–11 (BioMed Central).

73. Huerta-Cepas J, Szklarczyk D, Heller D, Hernández-Plaza A, Forslund SK, Cook H, et al. eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 2019;47:D309–14.

74. Ye Y, Doak TG. A Parsimony approach to biological pathway reconstruction/inference for genomes and metagenomes. PLOS Comput Biol. 2009;5:e1000465 (Public Library of Science).

75. Darzi Y, Falony G, Vieira-Silva S, Raes J. Towards biome-specific analysis of meta-omics data. ISME J. 2016;10:1025–8 (Nature Publishing Group).

76. Vieira-Silva S, Falony G, Darzi Y, Lima-Mendez G, Garcia Yunta R, Okuda S, et al. Species-function relationships shape ecological properties of the human gut microbiome. Nat Microbiol. 2016;1:16088.

77. Gertsman I, Barshop BA. Promises and pitfalls of untargeted metabolomics. J Inherit Metab Dis. 2018;41:355–66.

78. Ley RE, Bäckhed F, Turnbaugh P, Lozupone CA, Knight RD, Gordon JI. Obesity alters gut microbial ecology. Proc Natl Acad Sci. 2005;102:11070–5 (National Academy of Sciences).

79. Ley RE, Turnbaugh PJ, Klein S, Gordon JI. Human gut microbes associated with obesity. Nature. 2006;444:1022–3 (Nature Publishing Group).

80. Li X, Yan Q, Xie S, Hu W, Yu Y, Hu Z. Gut microbiota contributes to the growth of fast-growing transgenic common carp (*Cyprinus carpio* L.). PLOS One. 2013;8:e64577.

81. Duval C, Marie B, Foucault P, Duperron S. Establishment of the bacterial microbiota in a lab-reared model teleost fish, the Medaka *Oryzias latipes*. Microorganisms. 2022;10:2280 (Multidisciplinary Digital Publishing Institute).

82. Kaakoush NO. Insights into the role of Erysipelotrichaceae in the human host. Front Cell Infect Microbiol. 2015;5. Cited 2021 Mar 31. Available from: http://journal.frontiersin.org/Article/https://doi.org/10.3389/ fcimb.2015.00084/abstract

83. Wang SP, Rubio LA, Duncan SH, Donachie GE, Holtrop G, Lo G, et al. Pivotal roles for pH, lactate, and lactate-utilizing bacteria in the stability of a human colonic microbial ecosystem. Cotter PD, editor. mSystems. 2020;5. Cited 2022 Feb 17. Available from: https://journals.asm.org/ doi/https://doi.org/10.1128/mSystems.00645-20

84. Vernia P, Caprilli R, Latella G, Barbetti F, Magliocca FM, Cittadini M. Fecal lactate and ulcerative colitis. Gastroenterology. 1988;95:1564–8.

85. Guo F-F, Yu T-C, Hong J, Fang J-Y. Emerging roles of hydrogen sulfide in inflammatory and neoplastic colonic diseases. Front Physiol. 2016;7. Cited 2022 Feb 17. Available from: https://www.frontiersin.org/ article/https://doi.org/10.3389/fphys.2016.00156

86. McWilliam Leitch EC, Stewart CS. Escherichia coli O157 and non-O157 isolates are more susceptible to L-lactate than to D-lactate. Appl Environ Microbiol. 2002;68:4676–8.

87. Anderson CJ, Medina CB, Barron BJ, Karvelyte L, Aaes TL, Lambertz I, et al. Microbes exploit death-induced nutrient release by gut epithelial cells. Nature. 2021;596:262–7 (Nature Publishing Group).

88. Lee Y-S, Kim T-Y, Kim Y, Lee S-H, Kim S, Kang SW, et al. Microbiota-derived lactate accelerates intestinal stem-cell-mediated epithelial development. Cell Host Microbe. 2018;24:833-846.e6.

89. Qian M, Wang J, Ji X, Yang H, Tang B, Zhang H, et al. Sub-chronic exposure to antibiotics doxycycline, oxytetracycline or florfenicol impacts
gut barrier and induces gut microbiota dysbiosis in adult zebrafish (*Daino rerio*). *Ecotoxicol Environ Saf.* 2021;221: 112464.

90. Bao Z, Zhao Y, Wu A, Lou Z, Lu H, Yu Q, et al. Sub-chronic carbendazim exposure induces hepatic glycolipid metabolism disorder accompanied by gut microbiota dysbiosis in adult zebrafish (*Daino rerio*). *Sci Total Environ.* 2020;739: 140081.

91. 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.

92. Li T, Li H, Gatesoupe F-J, She R, Lin Q, Yan X, et al. Bacterial signatures of "red-operculum" disease in the gut of crucian carp (*Carassius auratus*). *Microb Ecol.* 2017;74:510–21.

93. Ma C, Chen C, Jia L, He X, Zhang B. Comparison of the intestinal microbiota composition and function in healthy and diseased Yunlong Grouper. AMB Express. 2019;9:187.

94. Hayatgheib N, Moreau E, Calvez S, Lepelletier D, Pouliquen H. A review of functional feeds and the control of *Aeromonas* infections in freshwater fish. *Aquac Int.* 2020;28:1083–123.

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

96. Austin B, Austin DA. Bacterial Fish Pathogens. Springer International Publishing; 2016. Cited 2021 Jun 9. Available from: http://link.springer.com/https://doi.org/10.1007/978-3-319-32674-0

97. Merel S, Walker D, Chicana R, Snyder S, Baurès E, Thomas O. State of knowledge and concerns on cyanobacterial blooms and cyanotoxins. *Environ Int.* 2013;59:303–27.

98. Halary S, Duval C, Gallet A, Duperron S, Piquet B, Demay J, et al. Draft genome sequence of the toxic freshwater *Microcystis aeruginosa* strain PMC 728.11 (*Cyanobacteria, Chroococcales*). Roux S, editor. *Microbiol Resour Announc.* 2020;9. Cited 2021 Mar 31. Available from: https://mra.asm.org/content/9/48/e01096-20

99. Riley MA, Wertz JE. Bacteriocin diversity: ecological and evolutionary perspectives. *Biochimie.* 2002;84:357–64.

100. Martins J, Vasconcelos V. Cyanobactins from cyanobacteria: current genetic and chemical state of knowledge. *Mar Drugs.* 2015;13:6910–46 (Multidisciplinary Digital Publishing Institute).

101. Moya A, Ferrer M. Functional redundancy-induced stability of gut microbiota subjected to disturbance. *Trends Microbiol.* 2016;24:402–13.

102. Louca S, Polz MF, Mazel F, Albright MBN, Huber JA, O'Connor MI, et al. Function and functional redundancy in microbial systems. *Nat Ecol Evol.* 2018;2:936–43 (Nature Publishing Group).

103. Iebba V, Totino V, Gagliardi A, Santangelo F, Cacciotti F, Trancassini M, et al. Eubiosis and dysbiosis: the two sides of the microbiota. *New Microbiol.* 2016;39:1–12.

104. Zaneveld JR, McMinds R, Vega Thurber R. Stress and stability: applying the Anna Karenina principle to animal microbiomes. *Nat Microbiol.* 2017;2:17121.

105. Vogt SL, Peña-Díaz J, Finlay BB. Chemical communication in the gut: effects of microbiota-generated metabolites on gastrointestinal bacterial pathogens. *Anaerobe.* 2015;34:106–15.

106. Li Z, Quan G, Jiang X, Yang Y, Ding X, Zhang D, et al. Effects of metabolites derived from gut microbiota and hosts on pathogens. *Front Cell Infect Microbiol.* 2018;8. Cited 2021 Jun 7. Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6152485/

107. Fontaine S, Kohl KD. Gut microbiota of invasive bullfrog tadpoles responds more rapidly to temperature than a noninvasive congener. *Mol Ecol.* 2020;29:2449–62.

108. Portz DE, Woodley CM, Cech JJ. Stress-associated impacts of short-term holding on fishes. *Rev Fish Biol Fish.* 2006;16:125–70.

109. Gerphagnon M, Macarthur DJ, Latour D, Gachon CMM, Ogtrop FV, Gleason FH, et al. Microbial players involved in the decline of filamentous and colonial cyanobacterial blooms with a focus on fungal parasitism. *Environ Microbiol.* 2015;17:2573–87.

110. Moreno-Ostos E, Cruz-Pizarro L, Basanta A, George DG. The influence of wind-induced mixing on the vertical distribution of buoyant and sinking phytoplankton species. *Aquat Ecol.* 2009;43:271–84.

111. Blaser MJ. Antibiotic use and its consequences for the normal microbiome. *Science.* 2016;352:544–5 (American Association for the Advancement of Science).

112. Reynolds CS. The ecology of phytoplankton. Cambridge: Cambridge University Press; 2006.

## Publisher's Note

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

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_15.png" /></figure>