Front. Toxicol. Frontiers in Toxicology Front. Toxicol. 2673-3080 Frontiers Media S.A. 1634241 10.3389/ftox.2025.1634241 Toxicology Original Research A short-term sublethal oral exposure to microcystin-LR disrupts cecal microbiome homeostasis in mallard Drovetski et al. 10.3389/ftox.2025.1634241 Drovetski Serguei V. 1 * Shearn-Bochsler Valerie I. 2 Hofmeister Erik K. 2 Karouna-Renier Natalie K. 1 Dusek Robert J. 2 1 U.S. Geological Survey, Eastern Ecological Science Center at the Patuxent Research Refuge, Laurel, MD, United States 2 U. S. Geological Survey, National Wildlife Health Center, Madison, WI, United States

Edited by: Alena Sevcu, Technical University of Liberec, Czechia

Reviewed by: Yu Zhang, Jinan University, China

Rene Shahmohamadloo, Washington State University, United States

 *Correspondence: Serguei V. Drovetski, sdrovetski@usgs.gov 15 10 2025 2025 7 1634241 23 05 2025 15 09 2025 Copyright © 2025 Drovetski, Shearn-Bochsler, Hofmeister, Karouna-Renier and Dusek. 2025 Drovetski, Shearn-Bochsler, Hofmeister, Karouna-Renier and Dusek

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

 Introduction

The frequency of cyanobacterial blooms seems to have increased globally in recent decades due to human induced eutrophication and climate change. Cyanobacterial blooms can produce several groups of toxins, among which microcystin-LR (MC-LR) is one of the most abundant. Effects of MC-LR on avian microbiome have not been studied and studies in laboratory murines have been limited to metabarcoding of prokaryotes.

 Methods

Using RNA shotgun sequencing, we compared the richness and composition of metabolically active prokaryotes, expressed virulence factors, antimicrobial resistance genes, metabolic pathways, Gene Ontology terms, enzymes, and proteins in mallards (Anas platyrhynchos) that were orally exposed to a sublethal dose of MC-LR for one week and unexposed birds.

 Results

Richness and composition of all compared features did not differ between exposed and control birds and none were differentially expressed between exposure groups. However, richness and/or composition of all features except virulence factors and Carbohydrate Active enzymes had multiple-fold greater dispersion in exposed birds than in controls. This effect was especially pronounced in expressed metabolic (MetaCyc) pathways.

 Discussion

Our results suggest that MC-LR exposure had a stochastic (rather than deterministic) effect on cecal microbiota, especially its function. Observed disturbance of the microbiota homeostasis is consistent with the Anna Karenina Principle. This principle has been documented in a wide range of eukaryotes using primarily microbial community metabarcoding. Although stochastic disturbance of microbiota function has been hypothesized, our study seems to be the first to demonstrate this in an experimental study.

 microcystin-LR mallard cecum microbiota function Anna Karenina principle metatranscriptomics section-at-acceptance Environmental Toxicology 1 Introduction

During recent decades, the frequency of harmful cyanobacterial blooms seems to have increased in fresh and brackish waterbodies globally (Paerl and Paul, 2012; Buratti et al., 2017; Anderson et al., 2021). Increasing anthropogenic eutrophication and global climate change synergistically facilitate expansion of cyanobacterial blooms (Paerl et al., 2011; Paerl and Paul, 2012; Michalak et al., 2013). Microcystin (MC) is the most frequently detected group of toxins produced by cyanobacteria (Campos and Vasconcelos, 2010; Li et al., 2017). They are potent hepatotoxins and among the >250 currently known congeners (Rattner, 2022), MC-LR is the most frequently reported (Gupta et al., 2003).

MC-LR toxicity has been well studied in laboratory mice (House mouse, Mus musculus). The intraperitoneal median lethal dose (LD50) for MC-LR in adult VAF/plus CD-1 and BULB/c mice was estimated at 0.065 mg/kg (Robinson et al., 1989; Yoshida et al., 1997), whereas the 24-h oral LD50 was estimated to be over two orders of magnitude higher - 10.9 mg/kg (Yoshida et al., 1997). Avian toxicity of MC-LR is less well known and LD50 concentrations have been estimated only for intraperitoneal injections but not for oral exposure. Intraperitoneal 24-h LD50 estimates for adult mallards (Anas platyrhynchos; 0.085 mg/kg body weight; mix of 96.4% MC-LR and 3.6% MC-YR) and week-old ducklings (0.065 mg/kg body weight) were similar to those of mice (Li et al., 2012). By more environmentally realistic oral administration, an acute outcome has not been achieved regardless of the dosage used. The highest concentration of MC-LR in these failed attempts to induce acute toxicity in mallards is reported at 17.5 mg/kg (Bong, 2021). Notwithstanding the failure to induce an acute outcome in oral exposure experiments, studies universally reported sub-lethal effects of MCs on different avian organs and systems, including liver, kidney, spleen, gut, and reproductive and immune systems (Pašková et al., 2008; Damkova et al., 2009; Damkova et al., 2 011; Peckova et al., 2009; Kral et al., 2012; Ondracek et al., 2012).

The great difference in outcomes between intraperitoneal injection of MCs and their oral administration suggests that the avian digestive system and microbiota may play an important role in MCs metabolism, especially in species that have well-developed ceca, such as galliforms and herbivorous waterfowl. The paired intestinal ceca play a crucial role in fowl adaptation to extracting nutrients from otherwise indigestible plant foods (Remington, 1989). They serve as microbial fermentation chambers in which non-starch polysaccharides (dietary fiber) are broken into volatile fatty acids (Clench and Mathias, 1995; Svihus et al., 2013) that can be utilized by the host as an energy source (Bergman, 1990). Furthermore, cecal microbiota are involved in detoxification of plant secondary compounds (Kohl et al., 2016) and environmental toxins (Jeong et al., 2013). In a study of the greater sage-grouse (Centrocercus urophasianus), the greatest seasonal differences in microbiota richness and composition related to dietary shifts from primarily insects and forbs in the summer to chemically defended sagebrush (Artemisia sp.) leaves in the winter were observed in the cecum (Drovetski et al., 2019). Likewise, during a chronic MC-LR exposure experiment with male BALB/c mice, the greatest changes in microbiota richness and abundance of sentinel bacteria were observed in the cecum (Chen et al., 2015). Given that gut microbiota plays a crucial role in its host’s health, immunity, and protection against pathogens (Hooper et al., 2012; Liang et al., 2018), a full understanding of MCs’ effects on the host is unlikely to be achieved without elucidating their effects on the gut microbiota. The authors are not aware of studies reporting effects of MC-LR exposure on avian microbiome.

Several experimental studies reported effects of MC-LR exposure on prokaryotic microbiota of laboratory mice (Chen et al., 2015; Sarkar et al., 2019; Guilin et al., 2020; Lee et al., 2020; Mills et al., 2021; Zhuang et al., 2021; Yang et al., 2022; Song et al., 2024). All of these studies used 16S rRNA gene sequences to assess changes in microbiota richness and composition. In one study, MC-LR was administered through intraperitoneal injection of 0.01 mg/kg body mass five times a week for 2 weeks (Sarkar et al., 2019), whereas in the other six studies MC-LR was administered orally through gavage or drinking water. Despite the substantial variation in MC-LR concentration (0.001–0.120 mg/L in drinking water or 0.04–4 mg/kg of body mass every 24–48 h through gavage) and duration of the administration period (1 week to 1 year), all but one study indicated significant changes in abundance of sentinel bacteria (e.g., decline in Firmicutes/Bacteroidetes ratio) in fecal or colonic microbiota. The only study that failed to detect compositional changes in fecal microbiota (Mills et al., 2021), exposed mice to the highest concentration of MC-LR (daily gavage of 4 mg/kg) for the shortest time period (1 week), suggesting that time of exposure may affect detection of microbiota changes. However, an equally short study with lower daily dose (1 mg/kg) reported significant compositional changes in Sprague–Dawley rat (Rattus norvegicus domestica) fecal microbiota (Lin et al., 2015). Only two studies sampled cecal microbiota of mice during experimental MC-LR exposure (Chen et al., 2015; Zhuang et al., 2021). In addition to changes in abundance of some bacterial taxa, both studies reported increase in cecal microbiota richness. Increased prokaryotic richness, likely due to colonization by opportunistic pathogens, has been reported for cecal microbiota of prairie grouse (Tympanuchus spp.) exposed to crop production (Drovetski et al., 2022) and may represent a common outcome of exposure to environmental stressors.

In the present study, we evaluate effects of an experimental week-long sublethal exposure to purified MC-LR administered through oral gavage on metabolically active prokaryotic microbiota, expression of virulence factors, antimicrobial resistance genes, and metabolic pathways in the ceca of the mallard. To achieve this goal, we employ deep shotgun metatranscriptomics (median of 121,129,328 sequences/sample) of the cecal content of fully-grown juvenile female mallards divided into two exposure categories: birds exposed to a daily dose of 0.75 mg MC-LR (n = 8) and control birds with no MC-LR exposure (n = 8).

 2 Materials and methods 2.1 Experimental design and sampling

We used 16 fully grown juvenile (16–20 weeks old) female mallards raised together on an outdoor farm in Wisconsin, United States (Table 1). During the experiment, all birds were kept in individual cages in the same environmental conditions and fed the same diet. MC-LR doses were prepared by diluting 1 mg/mL MC in >99.95% ethanol, dispensing 0.75 mL of the mixture into a gelatin caplet, and allowing the ethanol to evaporate. Mallards were randomly assigned to two groups: those that received a caplet daily with 0.75 mg MC-LR for 7 days orally (hereafter exposed birds; n = 8) and birds that received a similarly prepared caplet with no MC-LR (hereafter control birds; n = 8). The dosage was confirmed and liver MC-LR concentration was measured in exposed birds (Jenkins et al., 2025). All birds were euthanized approximately 24 h after their last dose.

Mallard (Anas platyrhynchos) cecal content sample IDs, sequencing depth, National Center for Biotechnology Information (NCBI) accession numbers, exposure categories (Supplementary Table 1 in Drovetski et al., 2025).

ID NWHC case number-ID Number of pair-ended reads NCBI accession Exposure category
PI202111 46882-01 69269328 SRX22659392 Control
PI202112 46882-02 59914980 SRX22659393 Control
PI202113 46882-03 74399926 SRX22659400 MC-LR
PI202114 46882-04 63006094 SRX22659401 MC-LR
PI202115 46882-05 61214348 SRX22659402 Control
PI202116 46882-06 65352132 SRX22659403 Control
PI202117 46882-07 73182852 SRX22659404 MC-LR
PI202118 46882-08 57385111 SRX22659405 MC-LR
PI202119 46882-09 64965961 SRX22659406 Control
PI202120 46882-10 49218452 SRX22659407 Control
PI202121 46882-11 54965880 SRX22659394 MC-LR
PI202122 46882-12 49817301 SRX22659395 MC-LR
PI202123 46882-13 53046416 SRX22659396 Control
PI202124 46882-14 59180630 SRX22659397 Control
PI202125 46882-15 58290957 SRX22659398 MC-LR
PI202126 46882-16 66641048 SRX22659399 MC-LR

Immediately after each bird was euthanized, its abdominal cavity was opened with a sterile scalpel blade to expose ceca. One cecum was sliced across with a sterile scalpel blade and approximately 100 mg of cecal content was expelled directly into a 2 mL screw cap tube pre-filled with 1 mL of DNA/RNA Shield and a mix of 0.5 mm and 0.1 mm ultra-high density BashingBeads (Zymo Research, Irvine, CA, USA). Samples were kept at 4 °C until DNA extraction within a maximum of 30 days after collection. Age of birds was confirmed by bursa Fabricii presence (well developed, fleshy in juvenile birds, atrophies in adults), and sex was confirmed by gonad examination.

 2.2 Molecular procedures

Total RNA from the cecal content was extracted using the ZymoBIOMICS DNA/RNA Miniprep Kit (Zymo Research, Irvine, California, USA) following the manufacturer’s instructions. RNA quantity was measured on a Qubit 4.0 using the High-Sensitivity RNA Assay Kit (ThermoFisher Scientific, Waltham, Massachusetts, USA). To test the collection tubes and the extraction kit for possible contamination, we included a single blank sample (a collection tube with 1 mL of DNA/RNA Shield and a mix of 0.5 mm and 0.1 mm ultra-high density BashingBeads but without sample) into our RNA extraction process. The Qubit assay did not detect any RNA using 10 μL of the blank sample (lower detection limit is 0.2 ng/μL), whereas the lowest RNA concentration estimated using 2 μL of our samples was 37.8 ng/μL (median 53.5 ng/μL). We did not sequence this negative control sample as our previous re-sequencing of the negative controls produced inconsistent sets of microbial strains reflecting other samples in the sequencing lane and did not bias microbial composition of real samples containing multiple orders of magnitude greater quantities of nucleic acids.

 2.3 Shotgun metatranscriptomic sequencing and bioinformatics

RNA extracts were shipped on dry ice to GeneWiz/Azenta Life Sciences (South Plainfield, New Jersey, USA) for quality control and sequencing. RNA libraries were prepared using the NEB Ultra II RNA kit (New England Biolabs, Inc., Ipswich, Massachusetts, USA) and QIAseq FastSelect -rRNA HMR Kit (QIAGEN, Germantown Maryland) for rRNA depletion. To control for possible batch effects, all libraries were sequenced simultaneously on an Illumina HiSeq 4000 platform (2 × 150bp; Illumina, Inc., San Diego, California, USA). The total number of resulting sequences was 1,959,702,832 or 979,851,416 pair-ended reads from our 16 samples. The number of reads obtained from individual samples ranged from 49,218,452 to 74,399,926 with a median of 60,564,664 (Table 1).

Raw sequences were uploaded to the CosmosID Metagenomics Cloud app.cosmosid.com (CosmosID Inc., Germantown, Maryland, USA, www.cosmosid.com) for taxonomic identification, virulence and antimicrobial resistance profiling, and functional analyses. The CosmosID Metagenomics Cloud uses the KEPLER pipeline (CosmosID Hab, 2025c) for host-agnostic microbial taxonomic profiling. KEPLER utilizes k-mer exact-matching and probabilistic alignment to identify and estimate normalized abundance of microbial taxa. Antimicrobial resistance (AMR) genes and virulence factors (VF) are classified using the Resfinder (Florensa et al., 2022) and VFDB (Liu et al., 2022) databases respectively, and the KEPLER-AMR/VF Profiler (CosmosID Hub, 2025a). For metatranscriptomic functional profiling, the quality-controlled and adaptor-trimmed in BBduk (Joint Genome Institute, 2023) reads are translated and searched against a comprehensive and non-redundant protein sequence database, UniRef 90 provided by UniProt (The UniProt Consortium, 2017). The mapping of the reads to genes is weighted by mapping quality, coverage and gene length to estimate community-wide weighted gene family abundances following Franzosa et al. (Franzosa et al., 2018). Gene families are then annotated to MetaCyc (Caspi et al., 2008; Caspi et al., 2018) metabolic pathways (Franzosa et al., 2018) and grouped into gene ontology (GO) terms, carbohydrate active (CAZy) and Enzyme Commission enzymes, and protein domains (Pfam).

We used the normalized read frequency of prokaryotic strains provided by the CosmosID Metagenomics Cloud for statistical analyses. The normalized read frequency is a probabilistic estimate of the number of reads aligning to a reference genome normalized by its size and is suitable for comparative or differential abundance analyses. The normalized read frequencies were cumulative sum scaled (CSS; McMurdie and Holmes, 2014) and log2 transformed to account for unequal sequencing depth among samples and non-normal distribution of strains in the samples. Because VFs and AMRs databases are gene-based rather than organism-based, we used the percent total matches (the number of total k-mers identified out of all possible k-mers), which approximates gene abundance in the sample and is directly comparable among samples without normalization. We multiplied these percentages by 106 to produce “copies per million” which we CSS + log2 transformed for analyses. For the functional classification, we used copies per million values that represent the Total-Sum Scaled (TSS) abundance values normalized for comparisons across samples with unequal sequencing depth (CosmosID Hub, 2025b). Detailed description of the CosmosID Metagenomics Cloud methodology is provided on the CosmosID Methods page (CosmosID Hab, 2025c).

 2.4 Data analyses

We compared richness and Shannon Index (Shannon, 1948) values of metabolically active microbial strains, expressed VFs, AMR genes, MetaCyc pathways, Gene Ontology (GO) terms, CAZy and Enzyme Commission enzymes, and Protein families (Pfam) between MC-LR exposure categories (exposed birds and controls; Table 1) using linear model regression (lm function) implemented in Lme4 v. 1.1–36 (Bates et al., 2015) R package and asymptotic test for the equality of coefficients of variation (Feltz and Miller, 1996) implemented in cvequality v.0.2.0 (Marwick and Krishnamoorthy, 2019) R package. To compare composition of microbial strains, VFs, AMR genes, MetaCyc pathways, GO Terms, CAZy and Enzyme Commission enzymes, and Pfam between exposure categories, we calculated pairwise Bray-Curtis dissimilarities (Bray and Curtis, 1957) among samples and used them for permutational multivariate analysis of variance PERMANOVA (Anderson, 2001), tests of multivariate dispersion, and principal coordinate analyses (PCoA) implemented in vegan v. 2.4–4 (Oksanen et al., 2019) R package. In all analyses listed in this paragraph, we used p ≤ 0.05 as indicative of the null hypothesis rejection.

To identify the microbial strains, VFs, AMR genes, MetaCyc pathways, GO Terms, CAZy and Enzyme Commission enzymes, and Pfam that best explain differences between exposed and control birds, we used the Microbiome Multivariate Associations with Linear Models MaAsLin 3 (Mallick et al., 2021) implemented in R (R Core Team, 2025). Although MaAsLin 3 includes a variety of statistical models for identifying differentially abundant features, we chose the log-transformed linear model on Total Sum Scaled (TSS) strain abundance data but did not normalize other data because they were normalized during processing in the CosmosID Metagenomics Cloud. We used minimum prevalence of 0.25 to exclude features that were found in less than half of the individuals sampled in a single exposure category. We used the maximum false discovery rate corrected p-value (q-value) of 0.05.

 3 Results 3.1 Prokaryotes

The Kepler pipeline identified a total of 263 metabolically active prokaryotic strains in our dataset - two archaeons and 261 bacteria (Supplementary Table 2 in Drovetski et al., 2025). Archaea were represented by a mesophilic ammonia-oxidizing chemolithoautotroph Nitrososphaera viennensis (Thermoproteota, 3.2% of total abundance, detected in all 16 birds) (Stieglmeier et al., 2014) and a methanogen Methanobrevibacter_A woesei (Methanobacteriota, 0.7%, detected in two exposed and two control birds) (Miller and Lin, 2002). Bacteria were represented by 11 phyla: Bacteroidota (37.9% of total abundance; detected in all 16 birds), Bacillota_C (29.1%, n = 16), Spirochaetota (10.80%, n = 16), Actinomycetota (5.92%, n = 16), Bacillota_A (5.90%, n = 15), Bacillota (2.37%, n = 14), Desulfobacterota (1.87%, n = 15), Pseudomonadota (1.64%, n = 16), Fusobacteriota (0.49% control n = 4, exposed n = 3), Campylobacterota (0.06%, n = 1 in each exposure category), and Deferribacterota (0.05%, control n = 2, exposed n = 3).

Prokaryotic richness in individual samples varied from 21 to 126 strains with a median of 84. Control birds had higher median richness (89.5) than exposed birds (75.5) but this difference was not statistically supported (Adj. R 2 = 0.041, F(1, 14) = 1.638, p = 0.222; Figure 1a), apparently due to 2.5 times greater coefficient of richness variation in exposed birds (0.430) than that in controls (0.171; D’AD = 4.405, p = 0.036). Shannon Index values varied from 1.611 to 3.565, with a median of 2.829. Neither Shannon Index values nor their coefficients of variation differed between exposed and control birds (Adj. R 2 = −0.052, F(1, 14) = 0.256, p = 0.256; D’AD = 0.565, p = 0.452).

Boxplot depicting richness (a) and principal coordinates analysis (PCoA) plot based on Bray-Curtis dissimilarities (b) of metabolically active prokaryotic strains in ceca of mallards (Anas platyrhynchos) experimentally exposed to microcystin-LR (MC-LR) versus controls. Boxes in the panel (a) represent interquartile ranges (from lower quartile to the upper quartile), whiskers represent minimum and maximum values, and the thick horizontal lines represent medians. Ellipses in panel (b) represent one standard deviation around the centroid. The equality of richness coefficients of variation (a) was rejected, whereas equality of multivariate dispersion (b) was not. The value of D′AD test statistic measures how far each sample coefficient of variation is from estimated population coefficient of variation.

Two-part graphic comparing prokaryote richness in control and MC-LR conditions. Panel (a) shows a box plot with higher richness variance for MC-LR, significant at p=0.036. Panel (b) displays a PCoA plot indicating sample clustering; the ellipses for control and MC-LR overlap and showing no significant difference in multivariate dispersions (p=0.315). Control is in white, MC-LR in gray.

Exposure to MC-LR did not affect the composition of the metabolically active prokaryotic community (PERMANOVA R 2 = 0.042, F(1, 14) = 0.617, p = 0.954), and the two exposure groups overlapped substantially on the PCoA plot (Figure 1b). Although the mean distance to centroid was higher in the MC-LR exposed birds (0.392) than in controls (0.336), the null hypothesis of the homogeneity of multivariate dispersions between exposure groups was not rejected (F(1, 14) = 1.081, p = 0.315). MaAsLin3 did not identify any prokaryotic strains differentially active between control and exposed birds.

 3.2 Virulence factors

KEPLER-AMR/VF Profiler identified a total of 102 expressed VFs in our dataset (Supplementary Table 3 in Drovetski et al., 2025). Richness of expressed VFs varied from 1 to 48 with a median of 7.5. No relationship between exposure to MC-LR and VF richness (Adj. R 2 = −0.049, F(1, 14) = 0.293, p = 0.597) or difference in coefficient of richness variation (D’AD = 0.003, p = 0.953) between exposed and control birds was observed. Shannon Index values varied from 0 to 3.779, with a median of 1.908. Shannon Index values did not differ between exposed and control birds (Adj. R 2 = −0.071, F(1, 14) = 0.006, p = 0.942) and the null hypothesis of equality of their coefficients of variation was not rejected (D’AD = 1.044, p = 0.307). VFs load (the sum of all VFs total matches in the sample) varied from 0.105 to 23.907, with a median of 1.585. Neither the VF load (Adj. R 2 = −0.055, F(1, 14) = 0.216, p = 0.649) nor its coefficients of variation (D’AD = 0.023, p = 0.881) differed between exposed and control birds. Exposure to MC-LR did not affect the composition of expressed VFs (PERMANOVA R 2 = 0.053, F(1, 14) = 0.776, p = 0.808), and the two exposure groups overlapped substantially on the PCoA plot. Although the mean distance to centroid was higher in the MC-LR exposed birds (0.622) than in controls (0.592), the null hypothesis of the homogeneity of multivariate dispersions between exposure groups was not rejected (F(1, 14) = 0.443, p = 0.477). MaAsLin3 identified no VFs whose expression differed between exposed and control birds.

 3.3 Antimicrobial resistance genes

KEPLER-AMR/VF Profiler identified a total of 67 expressed AMRs from nine classes in our dataset (Supplementary Table 4 in Drovetski et al., 2025). Richness of expressed AMRs varied from 6 to 25 with a median of 15.0. There was no relationship between exposure to MC-LR and richness of AMRs (Adj. R 2 = −0.022, F(1, 14) = 1.334, p = 0.267). There was no difference in coefficient of richness variation between exposed birds and controls (D’AD = 1.956, p = 0.162). Shannon Index values varied from 1.592 to 3.077, with a median of 2.491. Shannon Index values did not differ between exposed and control birds (Adj. R 2 = 0.073, F(1, 14) = 2.186, p = 0.161), whereas the null hypothesis of equality of their coefficients of variation was rejected (control cv = 0.093, exposed cv = 0.209; D’AD = 3.908, p = 0.048; Figure 2a). The load of expressed AMRs varied from 2.634 to 14.009, with a median of 6.452. Neither the load (Adj. R 2 = −0.025, F(1, 14) = 0. 0.641, p = 0.437) nor its coefficients of variation (D’AD = 0.813, p = 0.367) differed between exposed and control birds. Exposure to MC-LR did not affect the composition of expressed AMRs (PERMANOVA R 2 = 0.082, F(1, 14) = 1.246, p = 0.216), and the two exposure groups overlapped completely on the PCoA plot (Figure 2b). The mean distance to centroid was higher in the MC-LR exposed birds (0.362) than in controls (0.266; F(1, 14) = 5.053, p = 0.035) rejecting the null hypothesis of the homogeneity of multivariate dispersions between exposure groups. MaAsLin3 identified no AMRs differentially expressed between exposed and control birds.

Boxplot depicting the Shannon Index (a) and principal coordinates analysis (PCoA) plot based on Bray-Curtis dissimilarities (b) of expressed antimicrobial resistance genes (AMRs) in ceca of mallards (Anas platyrhynchos) experimentally exposed to microcystin-LR (MC-LR) versus controls. Boxes in the panel (a) represent interquartile ranges (from lower quartile to the upper quartile), whiskers represent minimum and maximum values, and the thick horizontal lines represent medians. Ellipses in panel (b) represent one standard deviation around the centroid. Both the equality of richness coefficients of variation (a) and equality of multivariate dispersion (b) were rejected. The value of D′AD test statistic measures how far each sample coefficient of variation is from estimated population coefficient of variation.

Chart a is a box plot comparing AMR Shannon Index between Control and MC-LR groups, showing higher varience in MC-LR. Chart b is a PCoA plot illustrating microbial community differences between the groups, with a significant difference in multivariate dispersions evident. 3.4 Metabolic pathways

CosmosID Metagenomics Cloud identified a total of 242 expressed MetaCyc pathways in our dataset (Supplementary Table 5 in Drovetski et al., 2025). Their richness varied from 32 to 167 with a median of 124.5. The richness of expressed MetaCyc pathways was greater in controls (median = 130.5) than exposed (118.5) birds but this differences was not statistically supported (Adj. R 2 = 0.072, F(1, 14) = 2.157, p = 0.164), likely due to three fold greater coefficient of richness variation in exposed (cv = 0.409) than in control birds (cv = 0.138; D’AD = 5.957, p = 0.015; Figure 3a). Shannon index values varied from 2.474 to 4.764. Similar to differences in richness, Shannon Index values were greater in controls (median = 4.213) than in exposed birds (median = 4.127) but this difference was marginally supported (Adj. R 2 = 0.120, F(1, 14) = 3.049, p = 0.103) also likely due to over three fold greater coefficients of variation in exposed birds (cv = 0.172) than in controls (cv = 0.051; D’AD = 8.073, p = 0.004; Figure 3b). Exposure to MC-LR did not alter composition of expressed MetaCyc pathways (PERMANOVA R 2 = 0.087, F(1, 14) = 1.339, p = 0.203), and the exposure groups overlapped completely on the PCoA plot (Figure 3c). The mean distance to centroid was greater in the MC-LR exposed birds (0.265) than in controls (0.171), but the null hypothesis of the homogeneity of multivariate dispersions between exposure groups was marginally rejected (F(1, 14) = 2.707, p = 0.059). MaAsLin3 identified no MetaCyc pathways that were differentially expressed between exposed and control birds.

Boxplots depicting richness (a) Shannon Index (b) and principal coordinates analysis (PCoA) plot based on Bray-Curtis dissimilarities of expressed metabolic (MetaCyc) pathways (c) in ceca of mallards (Anas platyrhynchos) experimentally exposed to microcystin-LR (MC-LR) versus controls. Boxes in the panel (a, b) represent interquartile ranges (from lower quartile to the upper quartile), whiskers represent minimum and maximum values, and the thick horizontal lines represent medians. Ellipses in panel (c) represent one standard deviation around the centroid. The equality of coefficients of variation of richness and Shannon Index (a,b) and equality of multivariate dispersion (c) were rejected, although marginally for the latter. The value of D′AD test statistic measures how far each sample coefficient of variation is from estimated population coefficient of variation.

Three-panel data visualization. Panel a: Box plots comparing MetaCyc richness between Control and MC-LR groups. The Control group shows higher variance, D’AD equals 6.0, p-value equals 0.015. Panel b: Box plots of MetaCyc Shannon index, with the Control group showing higher variance, D’ AD equals 8.1, p-value equals 0.004. Panel c: Principal coordinates analysis (PCoA) plot with ellipses representing Control and MC-LR groups. PCo one accounts for 50 percent and PCo two for 12.5 percent of variance. MC-LR birds had marginally greater multivariate dispersion. Fstatistic equals 2.7, p-value equals 0.069. 3.5 Gene ontology terms (GO terms), protein families (Pfam), Enzyme Commission enzymes, and carbohydrate active enzymes (CAZy)

There were no differences in α-diversity measurements, their coefficients of variation, or in composition of GO terms, CAZy and Enzyme Commission enzymes, and Pfam (Supplementary Tables 6–9, respectively in Drovetski et al., 2025) between control and MC-LR exposed birds. MaAsLin 3 did not identify any features that were differentially expressed between the experimental groups. However, β-diversity multivariate dispersion was greater in the MC-LR exposed birds than in controls for Pfam (F(1,14) = 4.4, N permutations = 999, p = 0.033; Figure 4a) and Enzyme Commission enzymes (F(1, 14) = 2.8, N permutations = 999, 0.049; Figure 4b), and the null hypothesis of the equality of multivariate dispersion was marginally rejected for GO terms (F(1, 14) = 3.3, N permutations = 999, p = 0.064; Figure 4c). Only for carbohydrate active enzymes, the equality of multivariate dispersion was not rejected (F(1, 14) = 2.3, n permutations = 999, p = 0.132).

Principal coordinates analysis (PCoA) plots based on Bray-Curtis dissimilarities of expressed protein families Pfam; (a) Enzyme Commission enzymes (b) and gene ontology (GO) terms (c). Ellipses represent one standard deviation around centroid. In all cases the equality of multivariate dispersion was rejected, although marginally for the latter.

Three panels (a, b, c) showing Principal Coordinates Analysis (PCoA) plots comparing control and MC-LR treatments. Ellipses indicate the confidence regions for each group. Each panel has specific PCo1, PCo2 percentages, F-value, permutations, and p-values. Panel a: PCo1 at 67.7%, PCo2 at 10.9%, F=4.4, p=0.033. Panel b: PCo1 at 49.5%, PCo2 at 12.8%, F=2.8, p=0.049. Panel c: PCo1 at 75.4%, PCo2 at 8.1%, F=3.3, p=0.064. In all three cases, MC-LR birds have greater multivariate dispersion. Control is marked with white and MC-LR with gray. 4 Discussion

The goal of our study was to investigate changes in mallard cecal microbiota activity and function following a week-long, orally administered exposure to a relatively low, sublethal dose of purified MC-LR. We chose mallard for this experimental study because this species is widely used in toxicity experiments and is one of the most frequently affected during wildlife mortality events coinciding with cyanobacterial blooms (Rattner, 2022). Mallard’s abundance, ecology, and foraging habits routinely put them in contact with algal blooms and facilitate ingestion of cyanotoxins. To the authors’ knowledge, this is the first study of MC-LR exposure effects on avian microbiome and the first to use shotgun RNA sequencing in terrestrial wildlife.

Neither richness, Shannon index values, nor composition of metabolically active prokaryotes, expressed VFs, AMRs, MetaCyc pathways, GO terms, Pfam, Enzyme Commission nor CAZy enzymes differed between exposed and control birds suggesting a lack of deterministic shifts in composition and function of mallard cecal microbiota in response to short-term, low-dose MC-LR exposure. Likewise, there were no metabolically active strains, VFs, AMRs, metabolic pathways, GO terms, enzymes or proteins that were differentially active between exposed and control birds. This appears to contradict findings of both earlier murine studies of MC-LR effects on cecal prokaryotic microbiota based on denaturing gradient gel electrophoresis (Chen et al., 2015) or short amplicon sequencing (Zhuang et al., 2021) of 16S rRNA gene. Both studies reported increases in prokaryotic microbiota richness, Shannon index values, and changes in abundance of sentential cecal bacteria. However, during library preparation in our study, the 16S rRNA gene fragments were depleted. It is plausible that an increase in operational taxonomic units (OTU) richness identified using substitutions and indels in hypervariable segments of 16S rRNA gene does not necessarily translate into proportional increases in strain richness calculated using variation in expressed non-ribosomal mRNA.

The lack of deterministic response, however, does not mean there was no effect of MC-LR exposure on the mallard cecal microbiota and its function. We observed increase in variance of α-diversity indicators and/or β-diversity multivariate dispersion of metabolically active prokaryotes and their expressed AMRs, MetaCyc pathways, GO terms, enzymes and proteins in exposed birds. Only VFs and carbohydrate active enzymes did not display elevated stochasticity in the birds exposed to MC-LR. The elevated stochasticity in the cecal microbiota and its function of exposed birds relative to controls, observed in our data, is consistent with the Anna Karenina Principle (Zaneveld et al., 2017) that has been documented in a variety of animal and plant study systems (Zaneveld et al., 2017; Ahmed et al., 2019; Ma, 2020; Arnault et al., 2023; Li and Yang, 2025). This principle, when applied to host-associated microbiota, suggests that in response to stressors, symbiotic microbial community composition may not shift to a new stable state but rather become unstable leading to its greater variation among stressor-exposed than unexposed hosts. Stressor-related loss of the ability to regulate microbial community composition by the host and/or its microbiota and resulting dysbiosis, the disruption of microbiota homeostasis, have been linked to many host diseases, adverse outcomes, and reduced fitness across multicellular eukaryotes from plants and corals to humans (Zaneveld et al., 2017; Ahmed et al., 2019; Ma, 2020; Arnault et al., 2023; Li and Yang, 2025). Unfortunately, we cannot determine whether the observed disruption of microbiome homeostasis is a direct effect of microbiota contact with the toxin or might be combined with indirect effect(s) through modulation of some host organ(s)’ function because our understanding of MC-LR toxicity in birds is still in its infancy (Rattner, 2022). Even in much better studied mammals, we have little information on MC-LR toxicokinetics and toxicodynamics besides liver and male reproductive system (Ma et al., 2021).

Disruption of cecal microbiota homeostasis is likely to be particularly consequential to the host. In predominantly herbivorous waterfowl, including mallard and gamebirds, paired ceca function as microbial fermentation chambers that break up indigestible complex carbohydrates (i.e., dietary fiber) into short chain fatty acids that serve as the primary source of energy to the host and convert uric acid to ammonia (Clench and Mathias, 1995; Svihus et al., 2013; Svihus, 2014). Strong functional selection and long digesta retention time in the cecum greatly reduce interindividual α- and β-diversity variance in cecal microbiota relative to microbiotas of all other gut regions (Drovetski et al., 2018; Drovetski et al., 2019). Anna Karenina effects, i.e., stochastic responses to sublethal MC-LR exposure, which we observed in the mallard cecal microbiota, have the lowest tolerance limits but greatest effects on host energy balance and, ultimately, fitness than similar responses in other digestive tract regions.

In conclusion, our study revealed that short-term sublethal oral exposure to purified MC-LR affects cecal microbiome in the semi-domestic mallard. These effects are consistent with the Anna Karenina Principle postulating that microbiota of the host exposed to a stressor, may have a stochastic rather than deterministic response. Indeed, the coefficient of variation of the metabolically active microbiota richness was greater in the birds exposed to MC-LR than in controls. At the same time there were neither shifts in cecal microbiota richness, Shannon Index, composition nor differentially active strains between exposure groups. However, manifestation of the Anna Karenina Principle was even more pronounced in microbiota function. Increased stochasticity among exposed birds was detected in richness, Shannon Index, and composition of expressed metabolic pathways, resulting in increased multidimensional dispersion of expressed AMRs, GO terms, proteins and enzymes. To the authors’ knowledge, this is the first documented extension of the Anna Karenina Principle beyond community composition to its function, which has been implied previously (Zaneveld et al.,2017; Arnault et al., 2023) but not tested.

 Data availability statement

Publicly available datasets were analyzed in this study. These data can be found here: https://doi.org/10.5066/P1OCUFM6. Raw sequence archives were deposited to the National Center for Biotechnology Information as the bioproject PRJNA1045981: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1045981/.

 Ethics statement

The animal study was approved by the U.S. Geological Survey, National Wildlife Health Center Institutional Animal Care and Use Committee (tracking number EP210811). The study was conducted in accordance with the local legislation and institutional requirements.

 Author contributions

SD: Writing – original draft, Methodology, Visualization, Data curation, Resources, Formal Analysis, Funding acquisition, Conceptualization, Project administration, Validation, Supervision, Investigation. VS-B: Validation, Methodology, Writing – review and editing, Investigation. EH: Investigation, Methodology, Validation, Writing – review and editing. NK-R: Supervision, Writing – review and editing, Resources. RD: Conceptualization, Investigation, Resources, Data curation, Funding acquisition, Project administration, Supervision, Methodology, Writing – review and editing.

 Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This research was supported by the U.S. Geological Survey Ecosystems Mission Area Environmental Health Program (Toxic Substances Hydrology and Contaminant Biology) through the Toxins and Harmful Algal Blooms Science Team. The funding body played no role in the design of the study, collection, analysis, and interpretation of data, or in writing the manuscript.

We are grateful to Jonas Lai for RNA extraction from cecal content samples. We thank E. Falendysz, H. Lamb, and the rest of the animal care staff for providing veterinary and daily care of the animals and S. Steinfeldt and B. Wolfe for assistance with necropsies and sample processing. The reviewers were Yu Zhang, Rene Shahmohamadloo and Viviana Mazzei provided much appreciated comments on the original draft. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Ahmed H. I. Herrera M. Liew Y. J. Aranda M. (2019). Long-term temperature stress in the coral model aiptasia supports the ‘Anna Karenina Principle’ for bacterial microbiomes. Front. Microbiol. 10, 975. 10.3389/fmicb.2019.00975 31139158 Anderson M. J. (2001). A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26 (1), 32–46. 10.1111/j.1442-9993.2001.01070.pp.x Anderson D. M. Fensin E. Gobler C. J. Hoeglund A. E. Hubbard K. A. Kulis D. M. (2021). Marine harmful algal blooms (HABs) in the United States: history, current status and future trends. Harmful Algae 102, 101975. 10.1016/j.hal.2021.101975 33875183 Arnault G. Mony C. Vandenkoornhuyse P. (2023). Plant microbiota dysbiosis and the Anna Karenina principle. Trends Plant Sci. 28 (1), 18–30. 10.1016/j.tplants.2022.08.012 36127241 Bates D. (2015). Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48. 10.18637/jss.v067.i01 Bergman E. N. (1990). Energy contributions of volatile fatty acids from the gastrointestinal tract in various species. Physiol. Rev. 70 (2), 567–590. 10.1152/physrev.1990.70.2.567 2181501 Bong Y. T. (2021). Microcystin-LR toxicity in mallards and delta smelt - ProQuest. University of California Davis. Available online at: https://www.proquest.com/openview/fe00e06261368154adf06b6faa291ee3/1?cbl=18750&diss=y&pq-origsite=gscholar (Accessed December 14, 2022). Bray J. R. Curtis J. T. (1957). An ordination of the upland forest communities of southern Wisconsin. Ecol. Monogr. 27 (4), 325–349. 10.2307/1942268 Buratti F. M. Manganelli M. Vichi S. Stefanelli M. Scardala S. Testai E. (2017). Cyanotoxins: producing organisms, occurrence, toxicity, mechanism of action and human health toxicological risk evaluation. Archives Toxicol. 91 (3), 1049–1130. 10.1007/s00204-016-1913-6 28110405 Campos A. Vasconcelos V. (2010). Molecular mechanisms of microcystin toxicity in animal cells. Int. J. Mol. Sci. 11 (1), 268–287. 10.3390/ijms11010268 20162015 Caspi R. Foerster H. Fulcher C. A. Kaipa P. Krummenacker M. Latendresse M. (2008). The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases. Nucleic Acids Res. 36 (Suppl. l_1), D623–D631. 10.1093/nar/gkm900 17965431 Caspi R. Billington R. Fulcher C. A. Keseler I. M. Kothari A. Krummenacker M. (2018). The MetaCyc database of metabolic pathways and enzymes. Nucleic Acids Res. 46 (D1), D633–D639. 10.1093/nar/gkx935 29059334 Chen J. Xie P. Lin J. Zeng C. (2015). Effects of microcystin-LR on gut microflora in different gut regions of mice. J. Toxicol. Sci. 40 (4), 485–494. 10.2131/jts.40.485 26165645 Clench M. H. Mathias J. R. (1995). “The avian cecum: a review,” 107(1), p. 29. CosmosID Hab (2025c). KeplerTM host-agnostic taxonomic profiling, cosmos-hub. Available online at: https://docs.cosmosid.com/docs/kepler-microbiome-profiler (Accessed August 13, 2025). CosmosID Hub (2025a). AMR/VF profiling, cosmos-hub. Available online at: https://docs.cosmosid.com/docs/amr-vf-pipeline (Accessed August 13, 2025). CosmosID Hub (2025b). Functional host-agnostic profiling, cosmos-hub. Available online at: https://docs.cosmosid.com/docs/functional-analysis (Accessed August 13, 2025). Damkova V. Sedlackova J. Bandouchova H. Peckova L. Vitula F. Hilscherova K. (2009). Effects of cyanobacterial biomass on avian reproduction: a Japanese quail model. Neuro Endocrinol. Lett. 30 (1), 205–210.20027172 Damkova V. Paskova V. Bandouchova H. Hilscherova K. Sedlackova J. Novotny L. (2011). Testicular toxicity of cyanobacterial biomass in Japanese quails. Harmful Algae 10 (6), 612–618. 10.1016/j.hal.2011.04.013 Drovetski S. V. O'Mahoney M. Ransome E. J. Matterson K. O. Lim H. C. Chesser R. T. (2018). Spatial organization of the gastrointestinal microbiota in urban Canada geese. Sci. Rep. 8 (1), 3713. 10.1038/s41598-018-21892-y 29487373 Drovetski S. V. O'Mahoney M. J. V. Matterson K. O. Schmidt B. K. Graves G. R. (2019). Distinct microbiotas of anatomical gut regions display idiosyncratic seasonal variation in an avian folivore. Anim. Microbiome 1 (1), 2. 10.1186/s42523-019-0002-6 33499946 Drovetski S. V. Schmidt B. K. Lai J. E. Gross M. S. Hladik M. L. Matterson K. O. (2022). Exposure to crop production alters cecal prokaryotic microbiota, inflates virulome and resistome in wild prairie grouse. Environ. Pollut. 306, 119418. 10.1016/j.envpol.2022.119418 35526643 Drovetski S. V. Dusek R. J. Karouna-Renier N. K. Hofmeister E. K. Shearn-Bochsler V. I. (2025). “Functional changes in cecal microbiome of mallards (Anas platyrhynchos) experimentally exposed to sublethal doses of Microcystin-LR,”. U.S. Geological Survey: ScienceBase. 10.5066/P1OCUFM6 Feltz C. J. Miller G. E. (1996). An asymptotic test for the equality of coefficients of variation from k populations. Statistics Med. 15 (6), 647–658. 10.1002/(sici)1097-0258(19960330)15:6<647::aid-sim184>3.0.co;2-p 8731006 Florensa A. F. Kaas R. S. Clausen P. T. L. C. Aytan-Aktug D. Aarestrup F. M. (2022). ResFinder – an open online resource for identification of antimicrobial resistance genes in next-generation sequencing data and prediction of phenotypes from genotypes. Microb. Genomics 8 (1), 000748. 10.1099/mgen.0.000748 35072601 Franzosa E. A. McIver L. J. Rahnavard G. Thompson L. R. Schirmer M. Weingart G. (2018). Species-level functional profiling of metagenomes and metatranscriptomes. Nat. Methods 15 (11), 962–968. 10.1038/s41592-018-0176-y 30377376 Guilin Z. Pengyu Z. Wei L. Fengqi H. Chen F. Yu Y. (2020). Reduction of gut microbial diversity and short chain fatty acids in BALB/c mice exposure to microcystin-LR. Ecotoxicology 29 (9), 1347–1357. 10.1007/s10646-020-02254-9 32772242 Gupta N. Pant S. C. Vijayaraghavan R. Rao P. V. L. (2003). Comparative toxicity evaluation of cyanobacterial cyclic peptide toxin microcystin variants (LR, RR, YR) in mice. Toxicology 188 (2), 285–296. 10.1016/S0300-483X(03)00112-4 12767698 Hooper L. V. Littman D. R. Macpherson A. J. (2012). Interactions between the microbiota and the immune system. Science 336 (6086), 1268–1273. 10.1126/science.1223490 22674334 Jenkins J. A. Wood K. Baudoin B. A. Foss A. J. Shearn-Bochsler V. I. Laughrey Z. A. (2025). “Laboratory analysis of liver tissue and cells from mallard ducks (Anas platyrhynchos) after oral exposure to microcystin-LR, 2021: U.S. Geological Survey data release.” U.S. Geological Survey: Science Base. 10.5066/P9EF18W0 Jeong H. G. Kang M. J. Kim H. G. Oh D. G. Kim J. S. Lee S. K. (2013). Role of intestinal microflora in xenobiotic-induced toxicity. Mol. Nutr. and Food Res. 57 (1), 84–99. 10.1002/mnfr.201200461 23166009 Joint Genome Institute (2023). “BBDuk guide,” in DOE joint genome Institute. Available online at: https://jgi.doe.gov/data-and-tools/software-tools/bbtools/bb-tools-user-guide/bbduk-guide/(Accessed April 20, 2023). Kohl K. D. Connelly J. W. Dearing M. D. Forbey J. S. (2016). Microbial detoxification in the gut of a specialist avian herbivore, the Greater Sage-Grouse. FEMS Microbiol. Lett. 363 (14), fnw144. 10.1093/femsle/fnw144 27242374 Kral J. Pikula J. Bandouchova H. Damkova V. Hilscherova K. Misik J. (2012). Avian high-dose toxicity of cyanobacterial biomass. Neuro Endocrinol. Lett. 33 (Suppl. 3), 161–165.23353861 Lee J. Lee S. Mayta A. Mrdjen I. Weghorst C. Knobloch T. (2020). Microcystis toxin-mediated tumor promotion and toxicity lead to shifts in mouse gut microbiome. Ecotoxicol. Environ. Saf. 206, 111204. 10.1016/j.ecoenv.2020.111204 32871519 Li W. Yang J. (2025). Investigating the Anna Karenina principle of the breast microbiome. BMC Microbiol. 25 (1), 81. 10.1186/s12866-024-03738-y 39979818 Li X.-Y. Wang S.-H. Wang C.-Y. Bai X. Ma J.-G. (2012). Exposure to crude microcystins via intraperitoneal injection, but not oral gavage, causes hepatotoxicity in ducks. Afr. J. Biotechnol. 11 (48), 10894–10898. 10.5897/ajb12.121 Li J. Li R. Li Ji (2017). Current research scenario for microcystins biodegradation – a review on fundamental knowledge, application prospects and challenges. Sci. Total Environ. 595, 615–632. 10.1016/j.scitotenv.2017.03.285 28407581 Liang D. Leung R. K. K. Guan W. Au W. W. (2018). Involvement of gut microbiome in human health and disease: brief overview, knowledge gaps and research opportunities. Gut Pathog. 10, 3. 10.1186/s13099-018-0230-4 29416567 Lin J. Chen J. He J. Chen J. Yan Q. Zhou J. (2015). Effects of microcystin-LR on bacterial and fungal functional genes profile in rat gut. Toxicon 96, 50–56. 10.1016/j.toxicon.2015.01.011 25617596 Liu B. Zheng D. Zhou S. Chen L. Yang J. (2022). VFDB 2022: a general classification scheme for bacterial virulence factors. Nucleic Acids Res. 50 (1), D912–D917. 10.1093/nar/gkab1107 34850947 Ma Z. (2020). Testing the Anna Karenina principle in human microbiome-associated diseases. iScience 23 (4), 101007. 10.1016/j.isci.2020.101007 32305861 Ma Y. Liu H. Du X. Shi Z. Liu X. Wang R. (2021). Advances in the toxicology research of microcystins based on Omics approaches. Environ. Int. 154, 106661. 10.1016/j.envint.2021.106661 34077854 Mallick H. Rahnavard A. McIver L. J. Ma S. Zhang Y. Nguyen L. H. (2021). Multivariable association discovery in population-scale meta-omics studies. PLOS Comput. Biol. 17 (11), e1009442. 10.1371/journal.pcbi.1009442 34784344 Marwick B. Krishnamoorthy K. (2019). Cvequality: tests for the equality of coefficients of variation from multiple groups. Available online at: https://cran.r-project.org/web/packages/cvequality/index.html (Accessed April 28, 2023). McMurdie P. J. Holmes S. (2014). Waste not, want not: why rarefying microbiome data is inadmissible. PLOS Comput. Biol. 10 (4), e1003531. 10.1371/journal.pcbi.1003531 24699258 Michalak A. M. Anderson E. J. Beletsky D. Boland S. Bosch N. S. Bridgeman T. B. (2013). Record-setting algal bloom in Lake Erie caused by agricultural and meteorological trends consistent with expected future conditions. Proc. Natl. Acad. Sci. 110 (16), 6448–6452. 10.1073/pnas.1216006110 23576718 Miller T. L. Lin C. (2002). Description of Methanobrevibacter gottschalkii sp. nov., Methanobrevibacter thaueri sp. nov., Methanobrevibacter woesei sp. nov. and Methanobrevibacter wolinii sp. nov. Int. J. Syst. Evol. Microbiol. 52 (3), 819–822. 10.1099/00207713-52-3-819 12054244 Mills M. C. Evans M. V. Lee S. Knobloch T. Weghorst C. Lee J. (2021). Acute cyanotoxin poisoning reveals a marginal effect on mouse gut microbiome composition but indicates metabolic shifts related to liver and gut inflammation. Ecotoxicol. Environ. Saf. 215, 112126. 10.1016/j.ecoenv.2021.112126 33721663 Oksanen J. Simpson G. L. Blanchet F. G. Kindt R. Legendre L. Minchin P. R. (2025). “vegan: R package for community ecologists”. vegandevs. Available online at: https://github.com/vegandevs/vegan (Accessed: March 10, 2025). Ondracek K. Bandouchova H. Damkova V. Hilscherova K. Kral J. Osickova J. (2012). Risk of combined exposure of birds to cyanobacterial biomass containing microcystins, acetylcholinesterase inhibitor and anticoagulant. Neuro Endocrinol. Lett. 33 (Suppl. 3), 155–160.23353860 Paerl H. W. Paul V. J. (2012). Climate change: links to global expansion of harmful cyanobacteria. Water Res. 46 (5), 1349–1363. 10.1016/j.watres.2011.08.002 21893330 Paerl H. W. Hall N. S. Calandrino E. S. (2011). Controlling harmful cyanobacterial blooms in a world experiencing anthropogenic and climatic-induced change. Sci. Total Environ. 409 (10), 1739–1745. 10.1016/j.scitotenv.2011.02.001 21345482 Pašková V. Adamovský O. Pikula J. Skocovská B. Band'ouchová H. Horáková J. (2008). Detoxification and oxidative stress responses along with microcystins accumulation in Japanese quail exposed to cyanobacterial biomass. Sci. Total Environ. 398 (1), 34–47. 10.1016/j.scitotenv.2008.03.001 18423819 Peckova L. Bandouchova H. Hilscherova K. Damkova V. Sedlackova J. Vitula F. (2009). Biochemical responses of juvenile and adult Japanese quails to cyanobacterial biomass. Neuro Endocrinol. Lett. 30 (Suppl. 1), 199–204.20027171 R Core Team (2025). R: the R Project for statistical computing v 4.5.1. Vienna, Austria: R Foundation for Statistical Computing. Available online at: https://www.r-project.org/(Accessed August 13, 2025). Rattner B. A. Wazniak C. E. Lankton J. S. McGowan P. C. Drovetski S. V. Egerton T. A. (2022). Review of harmful algal bloom effects on birds with implications for avian wildlife in the Chesapeake Bay region. Harmful Algae 120, 102319. 10.1016/j.hal.2022.102319 36470599 Remington T. E. (1989). Why do grouse have ceca? A test of the fiber digestion theory. J. Exp. Zoology 252 (S3), 87–94. 10.1002/jez.1402520515 2575134 Robinson N. A. Miura G. A. Matson C. F. Dinterman R. E. Pace J. G. (1989). Characterization of chemically tritiated microcystin-LR and its distribution in mice. Toxicon 27 (9), 1035–1042. 10.1016/0041-0101(89)90154-2 2508271 Sarkar S. Kimono D. Albadrani M. Seth R. K. Busbee P. Alghetaa H. (2019). Environmental microcystin targets the microbiome and increases the risk of intestinal inflammatory pathology via NOX2 in underlying murine model of Nonalcoholic Fatty Liver Disease. Sci. Rep. 9 (1), 8742. 10.1038/s41598-019-45009-1 31217465 Shannon C. E. (1948). A mathematical theory of communication. Bell Syst. Tech. J. 27 (3), 379–423. 10.1002/j.1538-7305.1948.tb01338.x Song Y. Wang X. Lu X. Wang T. (2024). Exposure to microcystin-LR promotes colorectal cancer progression by altering gut microbiota and associated metabolites in APCmin/+ mice. Toxins 16 (5), 212. 10.3390/toxins16050212 38787064 Stieglmeier M. Klingl A. Alves R. J. E. Rittmann S. K. M. R. Melcher M. Leisch N. (2014). Nitrososphaera viennensis gen. nov., sp. nov., an aerobic and mesophilic, ammonia-oxidizing archaeon from soil and a member of the archaeal phylum Thaumarchaeota. Int. J. Syst. Evol. Microbiol. 64 (Pt_8), 2738–2752. 10.1099/ijs.0.063172-0 24907263 Svihus B. (2014). Function of the digestive system. J. Appl. Poult. Res. 23 (2), 306–314. 10.3382/japr.2014-00937 Svihus B. Choct M. Classen H. L. (2013). Function and nutritional roles of the avian caeca: a review. World’s Poult. Sci. J. 69 (2), 249–264. 10.1017/S0043933913000287 The UniProt Consortium (2017). UniProt: the universal protein knowledgebase. Nucleic Acids Res. 45 (D1), D158–D169. 10.1093/nar/gkw1099 27899622 Yang Y. Wang H. Wang X. Chen L. Liu W. Cai D. (2022). Long-term environmental levels of microcystin-LR exposure induces colorectal chronic inflammation, fibrosis and barrier disruption via CSF1R/Rap1b signaling pathway. J. Hazard. Mater. 440, 129793. 10.1016/j.jhazmat.2022.129793 36029734 Yoshida T. Makita Y. Nagata S. Tsutsumi T. Yoshida F. Sekijima M. (1997). Acute oral toxicity of microcystin-LR, a cyanobacterial hepatotoxin, in mice. Nat. Toxins 5 (3), 91–95. 10.1002/1522-7189(1997)5:3<91::AID-NT1>3.0.CO;2-H 9285911 Zaneveld J. R. McMinds R. Vega Thurber R. (2017). Stress and stability: applying the Anna Karenina principle to animal microbiomes. Nat. Microbiol. 2 (9), 17121–17128. 10.1038/nmicrobiol.2017.121 28836573 Zhuang L. Jin Z. Li H. Wu S. Tong X. Wang H. (2021). Effects of chronic exposure to microcystin-LR on the gut microbiota of male mice. Int. J. Toxicol. 40 (2), 171–177. 10.1177/1091581820972311 33307919