Host species and habitat shape fish-associated bacterial communities: phylosymbiosis between fish and their microbiome
Microbiome volume 11, Article number: 258 (2023)
While many studies have reported that the structure of the gut and skin microbiota is driven by both species-specific and habitat-specific factors, the relative importance of host-specific versus environmental factors in wild vertebrates remains poorly understood. The aim of this study was to determine the diversity and composition of fish skin, gut, and surrounding water bacterial communities (hereafter referred to as microbiota) and assess the extent to which host habitat and phylogeny predict microbiota similarity. Skin swabs and gut samples from 334 fish belonging to 17 species were sampled in three Laurentian Great Lakes (LGLs) habitats (Detroit River, Lake Erie, Lake Ontario). We also collected and filtered water samples at the time of fish collection. We analyzed bacterial community composition using 16S metabarcoding and tested for community variation.
We found that the water microbiota was distinct from the fish microbiota, although the skin microbiota more closely resembled the water microbiota. We also found that environmental (sample location), habitat, fish diet, and host species factors shape and promote divergence or convergence of the fish microbiota. Since host species significantly affected both gut and skin microbiota (separately from host species effects), we tested for phylosymbiosis using pairwise host species phylogenetic distance versus bacterial community dissimilarity. We found significant phylogenetic effects on bacterial community dissimilarity, consistent with phylosymbiosis for both the fish skin and gut microbiota, perhaps reflecting the longstanding co-evolutionary relationship between the host species and their microbiomes.
Analyzing the gut and skin mucus microbiota across diverse fish species in complex natural ecosystems such as the LGLs provides insights into the potential for habitat and species-specific effects on the microbiome, and ultimately the health, of the host.
Host-associated microbiomes, specifically the bacterial community (hereafter referred to as microbiota) present inside the host as well as on body surfaces, influence a broad range of host immunological, evolutionary, and ecological processes [1, 2]. The term “microbiome” refers to the entire microbial ecosystem, including the microorganisms (prokaryotes and eukaryotes), their genomes, and their surrounding habitats , and is hypothesized to have co-evolved with its host . Significant research effort has focused on the importance of exogenous abiotic and biotic factors (e.g., habitat, geography, microbial biodiversity, diet) and endogenous host-related factors (e.g., genetics, physiology, immunity) in driving the composition of the microbiome [5,6,7]. At the microbial species level, deterministic (endogenous and exogenous) factors are thought to be the dominant forces shaping species composition of microbiota [8,9,10]. On the other hand, stochastic process-driven microbiome assembly (population growth rates, colonization events, extinction, and speciation) assumes that bacterial species are essentially functionally equivalent (neutral), and community assembly is the result of stochastic dispersal and drift through which organisms are randomly lost and replaced [11, 12]. Recent studies showed that both stochastic and deterministic processes are shaping the microbiome, but still a central question is the extent to which these two processes influence the host, as well as all of their associated microbes [10, 13]. Despite the known effects of habitat and host-specific factors on the microbiome across host taxa, very little is known about the degree of variation in fish-associated microbiota diversity and composition that occurs within and among species in the wild . Moreover, to know how host-associated microbiota may influence host phenotype, and hence contribute to evolutionary processes, quantifying the degree and nature of among-host species microbiome variation and the systematic drivers behind it seems crucial.
We expect that hosts and their microbiomes are linked eco-evolutionarily and the microbiome composition will recapitulate the phylogeny of their host, the basis of phylosymbiosis theory . Essentially, this predicts that hosts that are phylogenetically similar will have microbiomes that are more similar and vice versa . Phylosymbiosis may occur through stochastic and/or deterministic processes , mirroring genetic evolution by drift and/or selection. However, patterns of phylosymbiosis are also expected to be affected by exposure to habitat microbiome(s) . Although most studies in which phylosymbiosis has been identified have focused on microbes inhabiting internal organs, such as the gastrointestinal tract , recent work suggests that external host microbiomes (e.g., skin) can also exhibit a phylosymbiosis signal . In expanding the range of studied vertebrate microbiomes, questions about the range of environmental, ecological, and evolutionary factors that shape gut and skin microbiomes (or more specifically, bacterial communities), and the functions of those communities, still remained unanswered. For example, similar gut microbiome are found among phylogenetically related mammals but also among unrelated mammal species with similar diets [20,21,22].
Most host-microbial interaction studies are mainly focused on mammalian species, with fewer studies focused on other species . Teleosts encompass over half of vertebrate diversity  and are one of the most successful groups of vertebrates on Earth . Teleosts are represented by more than 32,000 species, originated over 600 million years ago, and exhibit a variety of physiologies, natural/life histories, and ecologies . However, their success and species radiation would not have been possible without the help of their microbiome . Additionally, their long history of co-evolution and symbiotic relationships with microbes (compared to mammals which evolved 160 million years ago ) make them good candidates to study host-microbial interactions. However, while most published studies of fish microbiomes include the gut microbiome, few studies have included other key microbial habitats such as skin [28, 29]. Fish skin mucosal immune responses play substantial roles in the course of infection prevention, and healthy skin (mucus) microbial communities are a critical component of that response . In fact, some beneficial bacteria on fish skin can produce antimicrobial compounds which help fish fight pathogenic microorganisms [31, 32]. Moreover, while fish skin is continuously exposed to numerous microorganisms (mainly from water sources), the skin can discriminate between beneficial microbiota and pathogenic microorganisms, although the mechanism is not fully understood . The role of the skin microbiome in host health is under-studied and, importantly for this study, the skin microbiome may thus be under separate selective pressures from the gut microbiome .
The Laurentian Great Lakes (LGLs) in North America form the largest freshwater ecosystem on the planet and have provided valuable ecosystem services for humans for centuries . The LGLs and their associated drainages include diverse ecosystems, complex trophic interactions, and mixed watersheds including forested, wetland, agricultural, and urban areas. The LGLs thus represent a powerful natural ecological laboratory to study important questions concerning microbial ecology and host-microbe interactions and co-evolution. While much work has been done on characterizing microbial communities in LGLs themselves, there is a critical need to examine the roles of host species and habitat on the microbiome of fish at the fish community level. The fish and microbial communities of LGLs provide a novel opportunity to study the intra- and interspecific divergence of host-associated intestinal and skin microbiota among the diverse species of the fish communities.
The aim of this study was to characterize fish bacterial communities among a variety of host species across three locations to (i) identify patterns of fish-associated bacterial communities, (ii) quantify the similarity among water bacterial communities and those of the fishes, and (iii) determine the extent to which microbiota among fish follow a pattern of phylosymbiosis. Our first hypothesis is that host-related deterministic processes (host-based selection pressures on bacteria composition) are the main drivers of bacterial community composition variation, and hence most bacterial community composition variation will be found among species. Our second hypothesis is that the water microbiota will be distinct from the host-associated microbiota, and the gut microbiota will also be different from the skin microbiota, but more similar to the water microbiota than the gut. Finally, we predict that the gut microbiota will be more strongly correlated with host phylogeny than the skin microbiota, as the skin microbiota will be affected by the local environmental microbiome more strongly. Characterizing intraspecific variation in fish microbiomes will help managers, particularly in aquaculture settings, effectively manipulate gut and skin microbiome to promote animal health and well-being. On the other hand, such information may provide fish managers tools to assess fish stock status and health. Specifically, skin swab samples could be used as a non-invasive sampling method to assess the health status of the fish, a valuable option for rare and at-risk fish species. More broadly, demonstrating the possible existence of phylosymbiosis in fish and their microbiomes will provide insight into factors governing the microbial community associated with fish and their co-evolutionary dynamics with teleosts in general.
Materials and methods
The Detroit River of the LGLs is a 51-km channel that comprises the lower portion of the Huron-Erie Corridor, connecting Lake St. Clair to Lake Erie . Samples were collected around Fighting Island (from 42° 10′ 56.5″ N 83° 06′ 39.9″ W up to 42° 14′ 04.2″ N 83° 06′ 32.1″ W) from July 17 to August 29, 2018.
Samples came from the western basin of Lake Erie, a shallow (mean 6 m), warm, and productive basin (area = 3473 km2) fed by the Detroit River and several smaller tributaries draining agricultural watersheds. Samples were collected on October 18, 2018.
Fish were collected from the Bay of Quinte and the eastern basin of Lake Ontario from July 31 to August 1, 2018. The Bay of Quinte is a large (254 km2), Z-shaped embayment with a history of nutrient and anthropogenic stress that feeds into the eastern basin and the upper St Lawrence River. The eastern basin is a mildly eutrophic, bathymetrically complex outlet basin of Lake Ontario with depths averaging ~ 20 m.
Fish were captured using a single anode boat electrofisher (Smith-Root 5.0 GPP) set to use pulsed DC current at 60 Hz using between 30 and 60% of the range to maintain a current of 6–8 A. All fish were euthanized with an overdose of tricaine methanesulfonate (MS222) following the protocols of the Ontario Ministry of Natural Resources and Forestry (MNRF), and Canadian Council on Animal Care. The skin swabs for all fishes were immediately taken after capture by gently rubbing a sterile cotton swab over ~ 50% of the total surface on the right side of each fish. Dorsal, ventral, and pectoral fin areas were swabbed for larger fish. Swab samples were placed into 2-mL tubes and were stored on ice and transported to the Great Lakes Institute for Environmental Research (GLIER) at the University of Windsor and immediately frozen at − 20 °C. The whole fish were transported on ice with the swabs to GLIER, and within 2–4 h after capture, the fish were dissected, and gut content samples were collected and stored frozen at − 20 °C. Water samples (500 mL) were collected at the sampling sites and transported to GLIER on ice for filtration. Water samples were filtered using 0.22-μm pore size, 47-mm-diameter polycarbonate membranes (Isopore™, Millipore, MA), and stored at − 20 °C until DNA extraction.
Fish were captured using bottom-set, graded mesh (32- to 152-mm) stretch monofilament gillnets. After capture, all fish were euthanized in compliance with the protocols of the MNRF and Canadian Council on Animal Care. The skin, fish, and water sampling, as well as the transportation and storage of the samples, were the same as for Detroit River.
Fish were captured using bottom-set, graded mesh (38–152 mm) stretch monofilament gillnets. Upon capture, all fish were euthanized following the protocol described for Detroit River and Lake Erie. Skin swab samples and gut tissue (foregut, midgut, and hindgut) with content were taken at the Glenora Fisheries Station, Ontario, Canada. Samples were stored in a 50 mL falcon tubes filled with 45 mL of a high salt solution (700 g/L ammonium sulfate, 25 mM sodium citrate, 20 mM ethylenediaminetetraacetic acid, pH 5.2) for 48 h to let the salts penetrate the samples (swabs, and gut content) and then stored at − 20 °C until DNA extraction. Water samples (500 mL) were collected and filtered at the Glenora Fisheries Station by using 0.22 μm filters pore size, 47 mm diameter, and the filters kept in high salt solution until delivered to GLIER. All samples were stored at − 20 °C until DNA extraction.
A total of 334 fish from Detroit River (n = 98, 29%), Lake Erie (n = 90, 27%), and Lake Ontario (n = 146, 44%) representing 17 different species were collected for microbiome characterization. The fish species included herbivores, invertivores, piscivore, planktivores, planktivore/invertivore, invertivores/carnivores, invertivores/detritivores, and invertivores/herbivores (Table 1; Supplementary Table S1) [36, 37]. For each of the 334 fish, fork length and total weight were recorded. Gut samples were taken after carefully dissecting the fish with a new razor blade or sterilized scissors to isolate a section comprising midgut and hindgut with both tissue and gut content. For Detroit River, Lake Erie, and Lake Ontario (based on sampling locations), a total of 4, 3, and 5 water samples were collected at the site of capture, respectively.
DNA extraction, library construction, and sequencing
DNA from all samples were extracted using a sucrose lysis buffer protocol as previously described . The V5 (787 F-acctgcctgccg-ATTAGATACCCNGGTAG) and V6 (1046 R-acgccaccgagc-CGACAGCCATGCANCACCT) variable regions of the 16S rRNA were selected for PCR amplification. The first and second PCR conditions were same as previously published methods . Briefly, the PCR protocol for the first round PCR consisted of 95 °C for 3 min followed by 28 cycles of 94 °C for 30 s, 55 °C for 30 s, and 72 °C for 1 m, and a final step at 72 °C for 10 m. For each 96 well PCR plate, one negative control consisting of PCR mix with ultra-pure water instead of DNA template was included. After the first-round PCR amplification, the results were verified by visualizing amplicons on an agarose gel. After checking and verifying first-round PCR products, the PCR was purified using Sera-Mag Magnetic Beads (GE, Healthcare Life Science, UK). A short-cycle second round PCR was conducted on the purified PCR products to ligate the adaptor and barcode (10–12 bp) sequences necessary for sample identification and sequencing. The second PCR was set at 95 °C for 3 min, then 8 cycles of 94 °C for 30 s, 60 °C for 30 s, and 72 °C for 1 min, and final extension at 72 °C for 7 min. The second round PCR amplifications were visualized on an agarose gel and amplicons were combined in proportions based on their estimated concentration (between 1 and 5 μL for samples with strong clear (1 μL) to faint bands (5 μL)). Subsequently, the combined samples were gel extracted from an agarose gel and cleaned and purified using QIAquick Gel Extraction Kit (QIAGEN, Toronto, ON, Canada). In total, 299 (89%), 330 (99%), and 10 (83%) samples for gut, skin swab, and water were amplified successfully and included in the sequencing library. We also included eight PCR blanks (one for each 96 PCR plate) in our library. The concentration of purified PCR product mix (library) was measured on an Agilent 2100 Bioanalyzer with a High Sensitivity DNA chip (Agilent Technologies, Mississauga, ON, Canada). The library concentration was then diluted to 60 pmol·μL−1 and sequenced on an Ion S5™ sequencing system using the Ion S5™ sequencing reagents and an Ion 530™ Chip (Thermo Fisher Scientific, ON, Canada).
Processing of 16S sequences
Two FASTQ files (one for gut and water samples and one for swab samples) were analyzed using the Quantitative Insights Into Microbial Ecology (QIIME2-2020.11) platform . The FASTQ sequence files were demultiplexed using the cutadapt demux-single command to remove sample barcode and primer sequences. Additionally, cutadapt trim-single was used to identify and remove the sequencing adapters for the demultiplexed data . The DADA2 pipeline (dada2 denoise-pyro) was used to denoise single-end sequences, dereplicate, and filter chimeras, followed by Amplicon Sequence Variant (ASV) picking . Chimeric sequences were removed using the removeBimeraDenovo function with the “consensus” method using the default values, except the read truncation length was set to 270 (p-trunc-len 270). The two ASV tables and representative sequences were merged using feature-table merge and feature-table merge-seqs, respectively. Taxonomic classification was performed using the feature-classifier plugin  and the SILVA 138–99 reference database . This plugin supports taxonomic classification of features using the Naive Bayes method. All ASVs were aligned with mafft  (via phylogeny align-to-tree-mafft-fasttree command) and used to construct a phylogeny with fasttree . After quality control, chimera removal and combining the two feature tables, the table was summarized with the feature-table summarize command. A total of 18,147,574 sequences were obtained for 647 samples (299 gut samples, 10 water samples, 330 skin samples plus eight negative controls). The eight negative controls had zero to 10 sequence reads with no consistent taxa present. The decontam  (version 1.8.0) package in the R was used to identify the blank sample ASVs as possible contaminants in our sample microbiota; however, none of the blank sample ASVs were identified as a contaminant. Thus, the negative controls were excluded from the rest of the study and we assumed contamination was not an issue for our data set.
We used the taxa filter-table, to remove ASVs related to mitochondria, chloroplast, and eukaryote sequences (3% of total sequence). We also removed “Unassigned” ASVs (13%), as well as bacteria and archaea ASVs without phylum assignment (1%), resulting in a total of 14,990,847 sequences remaining. The ASV table was rarefied to 3000 reads per sample for the alpha and beta diversity estimation because most of the rarefaction curves plateaued at ~ 3000 reads. Samples with fewer than 3000 reads were deleted. These deleted samples included three gut and six swab samples (Lake Trout (Salvelinus namaycush) (1 gut, 2 swab samples), White Perch (Morone americana) (1 gut), White Bass (Morone chrysops) (1 gut), Blacknose Shiner (Notropis heterolepis) (1 swab), Pumpkinseed (Lepomisgibbosus) (1 swab), White Sucker (Catostomus commersonii) (1 swab), Freshwater Drum (Aplodinotus grunniens) (1 swab)). This decreased the total number of samples to 630 samples (296 gut, 324 swab and 10 water samples) (Table 1). ASVs were retained for further analysis only if they had at least 10 sequence reads in at least two samples. So, the final analysis included 13,884,500 reads (93% of reads were retained) and 6,597,ASVs (15% of total ASVs (43,763)).
Bacteria alpha and beta diversity
Bacteria alpha diversity indices were calculated for each sample using QIIME alpha diversity alpha command. The calculated alpha diversity indices were Chao1 (a metric for species richness), and Faith’s phylogenetic diversity (PD) (a metric that incorporates both species richness and species evenness), Shannon entropy. We estimated beta diversity as the Bray–Curtis dissimilarity matrix among all samples. The rarefied ASV table was used for all subsequent analyses, unless stated. Raw data are available at the sequence read archive of NCBI with PRJNA701818 BioProject accession number.
Statistical analysis of sequence variants
Fish versus environmental (water) microbiota
To test for the effect of the environmental (water) microbiota on gut and skin microbiota, taxonomical compositions of the microbiota from the different sample types (gut, skin and water) were visualized using stacked barplots of the relative abundance of the bacteria at the phylum and family level with the online tool MicrobiomeAnalyst . Subsequently, a non-metric multidimensional scaling (NMDS) plot using the Bray–Curtis distance matrix was generated using vegan (v.2.6–4)  and ggplot2 (v.3.3.5; ) packages in R (v.4.1.0; ) to visualize clustering among the sample types (skin, gut, and water). We then performed permutational analyses of variance (PERMANOVA) using adonis2 in the vegan (v.2.6–4)  package in R  to test for significant differences in beta diversity among the sample types (gut vs water, as well as skin vs water samples). Finally, differences in alpha diversity (species richness and evenness (Shannon entropy, Chao1, PD)) among the different sample types (gut vs water, as well as skin vs water samples) at all locations combined as well as within each location separately were tested using the Kruskal–Wallis (KW) rank test, followed by a post hoc Dunn test with Bonferroni corrected P values in SPSS (IBM SPSS 25).
Fish gut versus skin microbiota
Taxonomic composition of the microbiota of gut and skin samples from seventeen fish species sampled at three locations (Detroit River, Lake Erie, and Lake Ontario) were visualized using stacked barplots of the relative abundance of the bacteria at the family level using phyloseq  and ggplot2  packages built on R . Bacteria families with relative abundance less than 10% (ranging between 0 and 9.99%) in all samples were combined and presented as “Family < 10 percent” in barplots. To test for the differential abundance of bacterial taxa between the gut and skin samples, the ASV table data were aggregated to the family level in R (v.4.1.0) using phyloseq package . Subsequently, the non-normalized ASV table was used for the negative binomial Wald test in DESeq2 . We used the default DESeq2 settings with negative binomial generalized linear model (GLM) fitting with Wald significance tests. P-values were adjusted for multiple testing using the Benjamini–Hochberg false discovery rate correction (FDR) . Bacterial taxa showing differential abundance were defined with the thresholds of FDR < 0.05 and |log2 Fold Change (FC)|> 2. Differences in alpha diversity (Shannon entropy, Chao1, PD) between gut and skin samples were statistically tested using the KW rank test (as described in Fish versus Environmental (Water) microbiota section above). To test for significant beta diversity differences between gut and skin samples, we used PERMANOVAs analysis with adonis2 in vegan (v.2.6–4)  package in R .
Endogenous and exogenous factor analyses
Gut and skin samples clustering based on location and fish species were visualized using an NMDS plot. Observed differences were assessed for significance with PERMANOVA analyses using adonis2 in vegan (v.2.6–4) ; however, we only included samples where fish species occurred in at least two of the three locations. We explored the role of location, fish species, body weight, and the interaction of fish species with sampling location in our PERMANOVA model for both gut and skin samples. Subsequently, we used least absolute shrinkage and selection operator (known as LASSO regression) in R in the glmnet package (v.4.1.7)  to identify the variables and corresponding regression coefficients that lead to a model that minimizes the prediction error and results in variable selection when there is a high likelihood of multicollinearity. We first quantified beta diversity in the gut and skin microbiota using a principal coordinates analysis (PCoA) in R (v.4.1.0) and retained the first five PCoA axes values (gut sample PCoA factors retained = eigenvalue > 2, percent of variation > 6%; skin sample PCoA factors retained = : eigenvalue > 5, percent of variation > 3%). We used Shannon entropy, Chao1, and PD (as described above) to quantify alpha diversity.
Given that we found evidence for significant species, location, and interaction effects on both gut and skin microbiome alpha and beta diversity, we included further analyses of possible mechanisms driving those effects. Specifically, the location effect may reflect habitat preferences among the 17 fish species sampled and the species effect may reflect diet preferences among the sampled species. To explore fish species, diet, and habitat preference effects on gut and skin bacterial community composition (alpha (Shannon entropy, Chao1 and PD) and beta diversity indices (PCoA1-5)), we built several linear mixed-model (LMM) models in R (v.4.1.0) using the lme4  package with fish species, habitat, diet, and body weight as fixed factors and location as a random factor. For selecting the best model for each diversity index, we used AICcmodavg package (v. 2.3.2)  in R (v.4.3.1) using Akaike Information Criterion (AIC) after building the models. We log10 transformed weight, Shannon entropy, Chao1, and PD to meet the normality and heteroscedasticity assumptions of LMMs.
To test if skin and gut microbiota composition were linked with host phylogeny, cytochrome c oxidase I (COX1) and cytochrome b (cytb) sequences for the 17 host fish species were downloaded from the NCBI website (https://www.ncbi.nlm.nih.gov/gene) and combined to create an artificial concatenated sequence. The concatenated sequences were aligned using MUSCLE  on the CIPRES Science Gateway v.3.1 . Subsequently, pair-wise phylogenetic distances between the species were calculated using the Kimura two-parameter model of substitution in MEGA-X (version 10.2.6) . Finally, to evaluate the effect of host phylogeny on microbiome dissimilarity (phylosymbiosis), we performed Mantel tests (method = "pearson", permutations = 999, na.rm = TRUE) in R (version 4.3.1) using the vegan (V2.6–4) package  to compare the bacterial community Bray–Curtis dissimilarity matrixes (skin and gut) and fish species phylogenetic distance matrix. To do so, Bray–Curtis dissimilarity was calculated across all samples within each fish species.
Fish versus environmental (water) microbiota
Taxonomic composition of the skin, gut, and water microbiota were distinct; however, the water microbiota was most divergent. Across all fish species, the fish microbiota were dominated at the phylum level by Pseudomonadota (previously known as Proteobacteria) (78% (skin), 56% (gut)), Fusobacteriota (6% (skin), 19% (gut)), and Bacillota (previously known as Firmicutes) (7% (skin), 18% (gut)). However, the water sample taxa showed a different set of common taxa, with only Pseudomonadota in common among the fish bacteria taxa (Pseudomonadota (38%), Actinobacteriota (27%), Bacteroidota (18%)) (Fig. S1). At the family level, the fish microbiota were dominated by members of Aeromonadaceae (gut (21%), skin (19%)), Fusobacteriaceae (gut (19%), skin (6%)), Enterobacteriaceae (gut (17%), skin (17%)), and Moraxellaceae (gut (1%), skin (20%)) (Fig. S2). However, at the family level, the water microbiota was dominated by Comamonadaceae (23%), Sporichthyaceae (19%), and Chitinophagaceae (5%) (Fig. S2).
Measures of alpha diversity (Shannon entropy, PD, Chao1) showed higher diversity in the water microbiota than in the fish gut and skin microbiota (Shannon entropy: KW 81, P < 0.0001; PD: KW 74, P < 0.0001; Chao1: KW 63, P < 0.0001;) (Fig. 1, Table 2). The post hoc pairwise comparisons revealed that when water samples were compared against gut and skin samples, gut sample microbiota were more divergent (Shannon entropy: test statistic − 320, adj P < 1.38E − 7; PD: test statistic − 316, adj P < 1.97E − 7; Chao1: test statistic − 292, adj P < 2.00E − 06) than the skin microbiota (Shannon entropy: test statistic − 205, adj P = 0.001; PD: test statistic − 208, adj P = 0.001; Chao1: test statistic − 192, adj P = 0.003) (Fig. 1, Table 2). We also tested for differences between gut vs water, and skin vs water separately within each location. The differences between gut vs water compared to skin vs water were more pronounced when we performed the analysis within each location (Table 2).
The NMDS plot also showed clear separation between the fish microbiota and the water microbiota (Fig. 2). The gut and skin microbiota showed considerable overlap in the NMDS; however, the water microbiota clustered separated from the fish microbiota, indicative of different community composition (Fig. 2). PERMANOVA analyses confirmed the statistical significance of the NMDS clusters (PERMANOVA pseudo-F: 12.8, P value < 0.001). Subsequently, individual PERMANOVAs for gut-water and skin-water comparisons showed significant effects for gut microbiota (t-value: 3.02; p-value < 0.001) and skin microbiota (F-value: 3.01; p-value < 0.001) compared to water microbiota.
Fish gut versus skin microbiota
Taxonomic analysis of the microbiota of skin and gut showed that Pseudomonadota (Fig. S1), and specifically members of the Aeromonadaceae and Enterobacteriaceae families generally dominated the fish microbiota (Figs. 3 and 4 and Fig. S2). However, gut and skin samples exhibited different bacterial taxa (Figs. 3 and 4). Differential abundance analysis at the family level was used to acquire more specific insight into differences in microbiota composition between the gut and skin microbiota. Comparing bacterial abundance data at the family level between skin and gut identified 37 bacteria families that had statistically significant differences (Table S2). For example, while members of the Deinococcaceae, Exiguobacteraceae, and Moraxellaceae families were in high abundance in skin samples, they were rare in the gut sample microbiota. On the other hand, Microbacteriaceae and Lachnospiraceae were at higher abundance in the gut microbiota (Table S2).
Moreover, the alpha diversity comparisons between gut and skin samples showed that skin samples had higher diversity relative to gut samples (Shannon entropy: test statistic − 114, adj P < 1.39E − 14; PD: test statistic − 108, adj P < 4.69E − 13; Chao1: test statistic − 99, adj P < 2.94E − 11) (Table 2, Fig. 1). Although gut and skin samples showed considerable overlap in the NMDS plot (Fig. 2), our PERMANOVA analysis showed that gut and skin samples had significantly different microbiota based on Bray–Curtis dissimilarity (t-value: 4.08; p-value < 0.001).
Endogenous and exogenous factor analyses
The microbiota composition in the gut and skin samples of different fish species sampled at different locations often showed high levels of variation among and within host species and among locations, highlighting the potential effect of location and fish species on fish microbiota (Figs. 3 and 4). Furthermore, NMDS analysis showed that microbiota composition clustered based on the host fish location (Lake Erie, Lake Ontario, Detroit River) (Fig. 5A, B) and host fish species within the location for both gut and skin microbiota (Fig. 5C, D). PERMANOVA analyses supported those clustering patterns and showed highly significant effects of location, host fish species, and the interaction between host fish species and location on both the gut and skin microbiota (Table 3). Moreover, body weight also showed a significant effect on both skin and gut samples.
Our LASSO regression analysis showed that the best predictor variables for gut samples were different based on the diversity indices (alpha (Shannon entropy, PD, Chao1) and beta diversity (PCoA 1- 5)), with fish species identified as the least predictor for alpha diversity indices and diet and location as the best predictor (Fig. S3). On the other hand, fish species and locations were the best predictors for beta diversity indices. For skin samples, both alpha and beta diversity indices showed fish species followed by weight as the best predictors (Fig. S4). Our LMM analysis (with fish species, habitat, and diet included as fixed factors and location as a random factor) of the effects on bacterial diversity indices (alpha (Chao1, PD) and beta (PCoA 1–5)) showed that fish species, diet, and habitat had significant effects on alpha diversity indices for gut samples but not for skin samples (Table 4). Moreover, fish species, diet, and habitat had significant effects for both gut and skin samples for the beta diversity indices. Despite including diet and habitat in our model, fish species remained as a significant effect, indicating the species effect is likely more complex than simple species-specific diet or habitat preferences.
Mantel tests of pairwise correlations between host phylogenetic distances and pairwise Bray–Curtis bacterial community dissimilarity (across all locations) values revealed a significant, but weak, positive relationship between bacterial community dissimilarity and host evolutionary distance for both gut (r = 0.18, P < 0.05) skin samples (r = 0.26, P < 0.01) supporting the phylosymbiosis for fish gut and skin samples (Fig. 6). However, as seven of the data points were separated from the rest (due to specific pylogenetic relationships—see Fig. 6), those seven data points were removed and the Mantel tests of pairwise correlations were re-calculated. By removing the seven comparisons from the data, the Mantel test revealed a significant positive relationship only between skin bacterial community dissimilarity and host evolutionary distance (Mantel statistic r: 0.25, P value < 0.05) (Fig. S5).
Diverse exogenous and endogenous factors have been reported to contribute to determining the composition of fish microbiomes [5, 14, 29]. However, the relative contributions of those factors in determining the composition of the teleost fish microbiome remain poorly characterized [7, 29]. Our study evaluated how the bacterial component of the gut and skin microbiota in 17 wild freshwater fish species sampled at three locations is shaped according to the surrounding water, trophic guild, diet, body weight, host species, and host phylogeny. While the fish microbiota (skin and gut) were different from the surrounding water microbiota, sample location (Detroit River, Lake Erie, Lake Ontario) had a strong effect on the composition of both the gut and skin microbiota across all species of fish. Moreover, host species as well as host species-by-location interactions were significant but did not explain as much variation in the fish bacterial community composition as the location’s main effect. Importantly, when we included fish diet and habitat preferences in our analysis, diet and habitat were important, but host species was still significant—indicating that the fish species effect was phylogenetically consistent and independent of diet and habitat, consistent with long-term co-evolutionary effects. When we tested the relationship between host phylogenetic distance and bacterial community dissimilarity, we found a significant correlation, consistent with phylosymbiosis. Overall, our results indicate that the host’s environment has a greater role than host-specific selection in the assembly and composition of both host-associated microbiota (gut and skin) and must be accounted for in any assessment of host-specific effects on the microbiome; however, evolutionary relationships between the host and its associated microbiomes are also important contributors.
The aquatic microbial community is thought to be the main source of the bacterial component of fish microbiomes, for both the gut and skin . Nevertheless, even with the ongoing and constant exposure to the surrounding aquatic microbiome, studies indicate that fish harbor microbiome that are distinct from the water microbiome [62, 63]. Our work supports those previous findings; we showed that the microbiota in the fish gut and skin were highly distinct from the surrounding water microbiota. Similar to other studies, we found Pseudomonadota, Actinobacteriota, and Bacteroidota were the most abundant phyla in the aquatic environment [28, 29, 64]. Other studies report that Pseudomonadota, Fusobacteriota, Bacillota, and Bacteroidota often comprise up to 95% of the fish bacterial communities [6, 29, 64, 65], also consistent with our work. These divergent patterns of bacterial abundance are expected, as Pseudomonadota play an important role in the growth of fishes through nutrient cycling and the mineralization of organic compounds , while Bacillota and Fusobacteriota have roles in fatty acid absorption, lipid metabolism, fermentative process, and degradation of oligosaccharides in fish . While, in our study, the dominant bacterial phyla in the fish had some overlap with the water microbiota, major differences were apparent, likely due to differences in the functional roles of the fish-associated versus aquatic microbiome. While differences between the fish microbiomes and the surrounding water microbiome are expected, we were surprised by the high level of divergence of the fish skin microbiota—clearly, the Teleost skin mucous microbiota plays important functional roles in the host’s performance and does not simply reflect the surrounding water’s microbiota.
While broad comparisons at the phylum level are valuable, a more detailed differential abundance analysis at the family level provides more specific insights into microbiota variation between gut and skin microbiota. Our differential abundance analysis showed fundamental differences between the microbiota of skin and gut across a diverse array of host fish species. For example, Deinococcaceae and Exiguobacteraceae were generally more abundant in the skin mucus microbiota, relative to the gut-content microbiota. On the other hand, members of Microbacteriaceae and Lachnospiraceae were more common in the gut microbiota. Bacteria in the family Deinococcaceae are obligate aerobes and have a high resistance to ionizing radiation (gamma- and/or ultra-violet (UV) radiation) . This makes it unsurprising that Deinococcaceae was more abundant in the skin microbiome, as skin experiences more exposure to solar radiation than the gut. Their resistance to ionizing radiation may also contribute to Deinococcaceae being reported in habitats associated with high levels of solar radiation, including fish skin , hot springs , and rivers . In our study, we found a general pattern of elevated levels of Deinococcaceae in pelagic species (such as Yellow Perch, Walleye) while lower levels were observed in deep water or benthic species (such as Brown Trout, Lake Trout, and Round Goby). The prevalence of Exiguobacterium in the skin microbiota in this study is consistent with its reported occurrence under a wide range of environmental conditions , including freshwater  and skin in humans . Perhaps one reason for the elevated abundance of Exiguobacterium in fish skin samples relative to the gut samples in our study is the ability of Exiguobacterium to thrive under highly variable conditions , an important consideration for fish skin microbiome that are exposed to considerable environmental variation relative to the gut habitat. Our detection of Microbacteriaceae at higher levels in the gut microbiota is consistent with previous work that reported them in various terrestrial and aquatic ecosystems , as well as associated with fish at various life stages [76,77,78]. The association of Lachnospiraceae with the fish gut microbiota likely reflects their important functional role as intestinal symbionts of vertebrates , acting as butyrate producers in the intestinal microbiome. Members of Lachnospiraceae have been previously reported in the fish gut [80, 81] and indeed have been shown to have a symbiotic link to their host in surgeonfish . Finally, we found substantially elevated levels of Enterobacteriaceae in the fish skin and gut samples, relative to the water samples (Fig. S1). Dominance of Enterobacteriaceae taxa in fish microbiomes has been widely reported in freshwater fish microbiomes [82,83,84] likely reflecting the importance of these bacteria for fish health and homeostasis [28, 85].
A large volume of published work shows that both endogenous and exogenous factors contribute to the Teleost microbiome composition [4, 5, 7, 8]. Endogenous factors can act at the individual level to drive variation in the microbiome (e.g., via life history, host health status, or host gene expression), as well as at the population level (e.g., via adaptation to local selection pressures) or even the species-level (e.g., via genomic variation and phylogenetic ancestry) . On the other hand, exogenous abiotic (climate, water chemistry, geography, etc.) and biotic (such as water microbiome) factors can contribute to variation in microbiome composition as well . Our analyses of alpha and beta diversity indices for both the skin and gut microbiota showed significant effects of location, habitat (exogenous), diet, and host fish species (endogenous). A variety of studies have reported both environmental and host species effects on the gut and skin microbiome composition in fishes [5, 29, 88, 89]. Generally, the skin microbiota is reported to be more affected by environmental factors than the gut microbiome [8, 29], which is not surprising, given the close contact between the aquatic environment and fish skin habitat . Surrounding environments such as water and sediment are thought to be major sources of skin and gut microbiome [91, 92]. Moreover, we observed significant differences in bacterial diversity between skin, gut, and the surrounding aquatic environment, with fish skin and water having higher diversity than gut, similar to other studies [29, 62, 93]. This is in contrast with other studies that included humans  and fish  that showed lower diversity for the skin microbiome (although this depended on the specific measure of diversity). One reason for this could be that our fish are from wild populations and other studies often use captive fish—captivity generally reduces alpha diversity and changes the composition/structure of vertebrate skin and gut microbiomes . Similarly, in humans, westernization significantly affects human microbiome diversity by lowering the diversity of bacteria . In this study, sample location dominated host species effects, with, as expected, a stronger effect for the skin microbiota (R2 = 0.17) than for the gut microbiota (R2 = 0.13), although host species effects were still substantial (skin microbiota (R2 = 0.05); gut microbiota (R2 = 0.05)). This pattern of effects on skin versus gut microbiota was despite the strong microbiota divergence between both fish microbiota and the water microbiota. However, our Kruskal–Wallis results also showed that skin microbiota were more similar to that found in the water than the gut microbiota. As the skin is in constant direct contact with the surrounding water microbiome, the skin microbiome would be expected to reflect at least part of the bacteriological composition of the surrounding water . By contrast, the gut microbiome is known to be strongly influenced by host-related factors and diet . In this study, we included 17 different fish species sampled at three different locations, and while we found a strong host species and location effects for both gut and skin, we recognized that those relationships may actually reflect local dietary/habitat preference variation among the study species [89, 99]. Thus, our reported exogenous factor (location) may actually include a component of endogenous effects. This is further supported by the significant species-by-location interaction effect; such an effect reflects species-specific sample location effects, which would logically be due to differences in host diet, at least for the gut microbiome. Moreover, our study showed approximately equal effects of host species effects on the gut and skin microbiota, perhaps due to the inclusion of multiple host species that utilized the aquatic habitat quite differently. Our observed host fish species effects on both gut and skin microbiome microbiota agree with previous research showing interindividual, population, and species variation for microbial community composition , all of which were interpreted as due to host endogenous factors. Given that the gut microbiome habitat is highly controlled by the host’s physiology, only bacterial species adapted to that environment would be expected to thrive in the gut; hence, the host should have a considerable effect on the gut bacterial community composition . The specific mechanisms driving variation in the fish bacterial community are still unclear, despite substantial research (including this work), likely due to complex interactions among possible mechanisms.
Microbiome similarity among species in a community is predicted to decrease with increasing evolutionary divergence of the host organisms , this is the basis for phylosymbiosis. However, phylosymbiotic patterns can result from diverse factors, including phenotypic divergence among fishes that are phylogenetically distant , co-evolution between the individual bacteria in the microbiome and the host , and even patterns of host behavior or life history that may be correlated with phylogeny but also indirectly affect the microbiome (e.g., feeding preferences) . Additionally, evolutionary processes such as selection and drift can also shape the species relatedness and their associated bacterial community, thus resulting in phylosymbiosis [16, 102]. Numerous studies have documented an effect of the host fish species on microbiome [4, 85, 88]. We showed correlations between bacterial community composition divergence and host fish taxonomic divergence for both the gut and skin microbiota. While specifies-specific diet differences do contribute to bacterial community similarity, the species effects were still present after correcting for diet, indicating that the mechanism behind our observed phylosymbiosis cannot be due to diet alone. Previous studies showed that host-specific microbiomes are a widespread pattern in nature, occurring in many host organisms , including mammals , birds , insects , and fish [4, 94, 105]. Although some bacterial lineages may still co-diversify with hosts, it is important to note that phylosymbiosis by itself is not an indicator of host–microbiome adaptive co-evolution. Evidence for phylosymbiosis in non-mammalian vertebrate animals is incomplete and inconsistent . For example, some studies in fish showed evidence for phylosymbiosis [4, 94], whereas others report mixed or weak evidence [7, 29]. Curiously, despite the fact that multiple studies (including ours) have shown that the fish skin microbiome is generally more affected by the environment than the gut microbiome [29, 62], we found a stronger signal of phylosymbiosis for the fish skin microbiota than for the gut microbiota. This was not expected, as the fish skin microhabitat is much more affected by environmental physicochemical parameters , not to mention the water microbiome. Clearly, the host species still has a substantial effect on even their peripheral microbiomes, highlighting the functional importance of all host-associated microbiomes. Given that our sampled fishes included 17 species representing 7 orders, our phylosymbiosis analysis is powerful and provides a solid foundation for future evaluations of the role of co-evolution between host organisms and their microbiome communities.
Overall, our findings contribute to the characterization of the modulators of microbiome composition and diversity across fish taxa. While many studies have characterized fish microbiome, yet few of those studies included multiple species sampled in the wild across multiple locations. Our study design provided a robust test of the relative effects of habitat, host diet, and host species on the bacterial communities of two key microbiomes associated with fish health and fitness. Not surprisingly, we found that the fish microbiota are distinct from the aquatic environmental microbiota, but that sampling location had a strong effect on microbiota composition, nevertheless. Curiously, we found strong host species-by-location interaction effects for both skin and gut microbiota, indicating that the species effects varied among the three sampled locations, possibly due to local fish diet and/or habitat-use differences. As expected, we also found a significant, but less strong effect of host fish species on both the gut and skin microbiota. Based on the host fish species effect that persisted after correcting for diet and habitat preferences, we tested for, and identified, significant phylosymbiotic signals between host phylogeny and both the gut and skin microbiome. This suggests that both the gut and skin microbiota co-evolved with their host species, although ecological covariation also contributes substantially since the variance in microbiota similarity explained was modest. Investigations of the nature of fish-microbe associations and, whether they are sustained, functional relationships or transient effects of fish and habitat associations are critical to further our understanding of the potential beneficial interactions between hosts and their microbiomes.
Availability of data and materials
The raw 16S rRNA gene sequencing data are available at the Sequence Read Archive of NCBI with PRJNA701818 BioProject accession number. For a complete list of packages and code for microbiome analyses, see https://github.com/javad30/Host-Species-and-Habitat-Shape-Fish--associated-Bacterial-Communities.
Bordenstein SR, Theis KR. Host biology in light of the microbiome: ten principles of holobionts and hologenomes. PLoS Biol. 2015;13(8):e1002226.
Woodhams DC, Bletz MC, Becker CG, Bender HA, Buitrago-Rosas D, Diebboll H, et al. Host-associated microbiomes are predicted by immune system complexity and climate. Genome Biol. 2020;21(1):23.
Marchesi JR, Ravel J. The vocabulary of microbiome research: a proposal. Microbiome. 2015;3:31.
Doane MP, Morris MM, Papudeshi B, Allen L, Pande D, Haggerty JM, et al. The skin microbiome of elasmobranchs follows phylosymbiosis, but in teleost fishes, the microbiomes converge. Microbiome. 2020;8(1):93.
Minich JJ, Petrus S, Michael JD, Michael TP, Knight R, Allen EE. Temporal, environmental, and biological drivers of the mucosal microbiome in a wild marine fish, Scomber japonicus. mSphere. 2020;5(3):e00401–20.
Minich JJ, Poore GD, Jantawongsri K, Johnston C, Bowie K, Bowman J, et al. Microbial ecology of Atlantic salmon (Salmo salar) hatcheries: impacts of the built environment on fish mucosal microbiota. Appl Environ Microbiol. 2020;86(12):e00411–20.
Riiser ES, Haverkamp THA, Varadharajan S, Borgan O, Jakobsen KS, Jentoft S, et al. Metagenomic shotgun analyses reveal complex patterns of intra- and interspecific variation in the intestinal microbiomes of codfishes. Appl Environ Microbiol. 2020;86(6):e02788–19.
Chiarello M, Paz-Vinas I, Veyssiere C, Santoul F, Loot G, Ferriol J, et al. Environmental conditions and neutral processes shape the skin microbiome of European catfish (Silurus glanis) populations of Southwestern France. Environ Microbiol Rep. 2019;11(4):605–14.
Rothschild D, Weissbrod O, Barkan E, Kurilshikov A, Korem T, Zeevi D, et al. Environment dominates over host genetics in shaping human gut microbiota. Nature. 2018;555(7695):210–5.
Zhou J, Ning D. Stochastic community assembly: does it matter in microbial ecology? Microbiol Mol Biol Rev. 2017;81(4):10–1128.
Heys C, Cheaib B, Busetti A, Kazlauskaite R, Maier L, Sloan WT, et al. Neutral processes dominate microbial community assembly in Atlantic salmon, Salmo salar. Appl Environ Microbiol. 2020;86(8):e02283–19.
Sadeghi J, Chaganti SR, Shahraki AH, Heath DD. Microbial community and abiotic effects on aquatic bacterial communities in north temperate lakes. Sci Total Environ. 2021;781:146771.
Kohl KD, Dearing MD, Bordenstein SR. Microbial communities exhibit host species distinguishability and phylosymbiosis along the length of the gastrointestinal tract. Mol Ecol. 2018;27(8):1874–83.
Uren Webster TM, Consuegra S, Hitchings M, Garcia de Leaniz C. Interpopulation variation in the Atlantic salmon microbiome reflects environmental and genetic diversity. Appl Environ Microbiol. 2018;84(16):e00691–18.
Brooks AW, Kohl KD, Brucker RM, van Opstal EJ, Bordenstein SR. Phylosymbiosis: relationships and functional effects of microbial communities across host evolutionary history. PLoS Biol. 2016;14(11):e2000225.
Mallott EK, Amato KR. Host specificity of the gut microbiome. Nat Rev Microbiol. 2021;19(10):639–53.
Lim SJ, Bordenstein SR. An introduction to phylosymbiosis. Proc Biol Sci. 1922;2020(287):20192900.
Lutz HL, Jackson EW, Webala PW, Babyesiza WS, Kerbis Peterhans JC, Demos TC, et al. Ecology and host identity outweigh evolutionary history in shaping the bat microbiome. mSystems. 2019;4(6):10–1128.
Ross AA, Muller KM, Weese JS, Neufeld JD. Comprehensive skin microbiome analysis reveals the uniqueness of human skin and evidence for phylosymbiosis within the class Mammalia. Proc Natl Acad Sci U S A. 2018;115(25):E5786–95.
Ley RE, Lozupone CA, Hamady M, Knight R, Gordon JI. Worlds within worlds: evolution of the vertebrate gut microbiota. Nat Rev Microbiol. 2008;6(10):776–88.
Ley RE, Hamady M, Lozupone C, Turnbaugh PJ, Ramey RR, Bircher JS, et al. Evolution of mammals and their gut microbes. Science. 2008;320(5883):1647–51.
Muegge BD, Kuczynski J, Knights D, Clemente JC, Gonzalez A, Fontana L, et al. Diet drives convergence in gut microbiome functions across mammalian phylogeny and within humans. Science. 2011;332(6032):970–4.
Pascoe EL, Hauffe HC, Marchesi JR, Perkins SE. Network analysis of gut microbiota literature: an overview of the research landscape in non-human animal studies. ISME J. 2017;11(12):2644–51.
Nelson JS, Grande TC, Wilson MV. Fishes of the World. Wiley; 2016.
Colston TJ, Jackson CR. Microbiome evolution along divergent branches of the vertebrate tree of life: what is known and unknown. Mol Ecol. 2016;25(16):3776–800.
Wang F, Men X, Zhang G, Liang K, Xin Y, Wang J, et al. Assessment of 16S rRNA gene primers for studying bacterial community structure and function of aging flue-cured tobaccos. AMB Express. 2018;8(1):182.
Ghanbari M, Kneifel W, Domig KJ. A new view of the fish gut microbiome: advances from next-generation sequencing. Aquaculture. 2015;448:464–75.
Krotman Y, Yergaliyev TM, Alexander Shani R, Avrahami Y, Szitenberg A. Dissecting the factors shaping fish skin microbiomes in a heterogeneous inland water system. Microbiome. 2020;8(1):9.
Sylvain FE, Holland A, Bouslama S, Audet-Gilbert E, Lavoie C, Val AL, et al. Fish skin and gut microbiomes show contrasting signatures of host species and habitat. Appl Environ Microbiol. 2020;86(16).
Salinas I, Fernandez-Montero A, Ding Y, Sunyer JO. Mucosal immunoglobulins of teleost fish: a decade of advances. Dev Comp Immunol. 2021;121:104079.
Ross AA, Rodrigues Hoffmann A, Neufeld JD. The skin microbiome of vertebrates. Microbiome. 2019;7(1):79.
Hussain A, Sachan SG. Fish epidermal mucus as a source of diverse therapeutical compounds. Int J Pept Res Ther. 2023;29(3):36.
Grice EA, Segre JA. The skin microbiome. Nat Rev Microbiol. 2011;9(4):244–53.
Ozersky T, Bramburger AJ, Elgin AK, Vanderploeg HA, Wang J, Austin JA, et al. The changing face of winter: lessons and questions from the Laurentian Great Lakes. Wiley Online Library; 2021.
Lapointe NW. Effects of shoreline type, riparian zone and instream microhabitat on fish species richness and abundance in the Detroit River. J Great Lakes Res. 2014;40:62–8.
Eakins RJ. Ontario Freshwater Fishes Life History Database. Version 5.04. 2020. https://www.ontariofishes.ca.
Heuvel CE, Haffner GD, Zhao Y, Colborne SF, Despenic A, Fisk AT. The influence of body size and season on the feeding ecology of three freshwater fishes with different diets in Lake Erie. J Great Lakes Res. 2019;45(4):795–804.
Shahraki AH, Chaganti SR, Heath D. Assessing high-throughput environmental DNA extraction methods for meta-barcode characterization of aquatic microbial communities. J Water Health. 2019;17(1):37–49.
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(8):852–7.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet journal. 2011;17(1):10–2.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJ, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13(7):581–3.
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(1):90.
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(Database issue):D590-6.
Katoh K, Misawa K, Kuma K, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30(14):3059–66.
Price MN, Dehal PS, Arkin AP. FastTree 2–approximately maximum-likelihood trees for large alignments. PLoS ONE. 2010;5(3):e9490.
Davis NM, Proctor DM, Holmes SP, Relman DA, Callahan BJ. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome. 2018;6(1):226.
Chong J, Liu P, Zhou G, Xia J. Using MicrobiomeAnalyst for comprehensive statistical, functional, and meta-analysis of microbiome data. Nat Protoc. 2020;15(3):799–821.
Oksanen J, Blanchet FG, Kindt R, Legendre P, Minchin P, O’hara R, et al. Community ecology package. R package version. 2013;2(0):103–32.
Wickham H. ggplot2: elegant graphics for data analysis: springer; 2016:203–20.
Team RC. R: A language and environment for statistical computing. 2013:275–86.
McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE. 2013;8(4):e61217.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc: Ser B (Methodol). 1995;57(1):289–300.
Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1.
Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. arXiv preprint arXiv:14065823. 2014;67(1):201–210.
Mazerolle MJ, Mazerolle MMJ. Package ‘AICcmodavg’. R package. 2017;281:1–220.
Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32(5):1792–7.
Miller MA, Pfeiffer W, Schwartz T, editors. The CIPRES science gateway: a community resource for phylogenetic analyses. Proceedings of the 2011 TeraGrid Conference: extreme digital discovery; 2011:1–8.
Tamura K, Dudley J, Nei M, Kumar S. MEGA4: molecular evolutionary genetics analysis (MEGA) software version 4.0. Mol Biol Evol. 2007;24(8):1596–9.
Hammer Ø, Harper DA, Ryan PD. PAST: Paleontological statistics software package for education and data analysis. Palaeontol Electron. 2001;4(1):9.
Galbraith H, Iwanowicz D, Spooner D, Iwanowicz L, Keller D, Zelanko P, et al. Exposure to synthetic hydraulic fracturing waste influences the mucosal bacterial community structure of the brook trout (Salvelinus fontinalis) epidermis. AIMS Microbiol. 2018;4(3):413–27.
Reinhart EM, Korry BJ, Rowan-Nash AD, Belenky P. Defining the distinct skin and gut microbiomes of the northern pike (Esox lucius). Front Microbiol. 2019;10:2118.
Chiarello M, Auguet JC, Bettarel Y, Bouvier C, Claverie T, Graham NAJ, et al. Skin microbiome of coral reef fish is highly variable and driven by host phylogeny and diet. Microbiome. 2018;6(1):147.
Uren Webster TM, Rodriguez-Barreto D, Castaldo G, Gough P, Consuegra S, de Garcia Leaniz C. Environmental plasticity and colonisation history in the Atlantic salmon microbiome: a translocation experiment. Mol Ecol. 2020;29(5):886–98.
Li X, Zhou L, Yu Y, Ni J, Xu W, Yan Q. Composition of gut microbiota in the gibel carp (Carassius auratus gibelio) varies with host development. Microb Ecol. 2017;74(1):239–49.
Kirchman DL. The ecology of Cytophaga-Flavobacteria in aquatic environments. FEMS Microbiol Ecol. 2002;39(2):91–100.
Slade D, Radman M. Oxidative stress resistance in Deinococcus radiodurans. Microbiol Mol Biol Rev. 2011;75(1):133–91.
Hu C, Huang Z, Liu M, Sun B, Tang L, Chen L. Shift in skin microbiota and immune functions of zebrafish after combined exposure to perfluorobutanesulfonate and probiotic Lactobacillus rhamnosus. Ecotoxicol Environ Saf. 2021;218:112310.
Ferreira AC, Nobre MF, Rainey FA, Silva MT, Wait R, Burghardt J, et al. Deinococcus geothermalis sp. nov. and Deinococcus murrayi sp. nov., two extremely radiation-resistant and slightly thermophilic species from hot springs. Int J Syst Bacteriol. 1997;47(4):939–47.
Lee JJ, Lee YH, Park SJ, Lee SY, Park S, Lee DS. Deinococcus knuensis sp. nov., a bacterium isolated from river water. Antonie Van Leeuwenhoek. 2017;110(3):407–14.
Dastager SG, Mawlankar R, Sonalkar VV, Thorat MN, Mual P, Verma A, et al. Exiguobacterium enclense sp. nov., isolated from sediment. Int J Syst Evol Microbiol. 2015;65(Pt 5):1611–6.
White RA 3rd, Soles SA, Gavelis G, Gosselin E, Slater GF, Lim DSS, et al. The complete genome and physiological analysis of the eurythermal firmicute Exiguobacterium chiriqhucha strain RW2 isolated from a freshwater microbialite, widely adaptable to broad thermal, pH, and salinity ranges. Front Microbiol. 2018;9:3189.
Tena D, Martinez NM, Casanova J, Garcia JL, Roman E, Medina MJ, et al. Possible Exiguobacterium sibiricum skin infection in human. Emerg Infect Dis. 2014;20(12):2178–9.
Vishnivetskaya TA, Kathariou S, Tiedje JM. The Exiguobacterium genus: biodiversity and biogeography. Extremophiles. 2009;13(3):541–55.
Evtushenko LI, Takeuchi M. The family microbacteriaceae The prokaryotes. 2006;3:1020–98.
Liu Y, de Bruijn I, Jack AL, Drynan K, van den Berg AH, Thoen E, et al. Deciphering microbial landscapes of fish eggs to mitigate emerging diseases. ISME J. 2014;8(10):2002–14.
Jami M, Ghanbari M, Kneifel W, Domig KJ. Phylogenetic diversity and biological activity of culturable Actinobacteria isolated from freshwater fish gut microbiota. Microbiol Res. 2015;175:6–15.
Nguyen CDH, Amoroso G, Ventura T, Minich JJ, Elizur A. Atlantic salmon (Salmo salar L., 1758) gut microbiota profile correlates with flesh pigmentation: cause or effect? Mar Biotechnol (NY). 2020;22(6):786–804.
Arroyo FA, Pawlowska TE, Choat JH, Clements KD, Angert ER. Recombination contributes to population diversification in the polyploid intestinal symbiont Epulopiscium sp. type B. ISME J. 2019;13(4):1084–97.
Ricaud K, Rey M, Plagnes-Juan E, Larroquet L, Even M, Quillet E, et al. Composition of intestinal microbiota in two lines of rainbow trout (Oncorhynchus mykiss) divergently selected for muscle fat content. Open Microbiol J. 2018;12:308–20.
Escalas A, Auguet J-C, Avouac A, Seguin R, Gradel A, Borrossi L, et al. Ecological specialization within a carnivorous fish family is supported by a herbivorous microbiome shaped by a combination of gut traits and specific diet. Front Mar Sci. 2021;8:91.
Egerton S, Culloty S, Whooley J, Stanton C, Ross RP. The gut microbiota of marine fish. Front Microbiol. 2018;9:873.
Khurana H, Singh DN, Singh A, Singh Y, Lal R, Negi RK. Gut microbiome of endangered Tor putitora (Ham.) as a reservoir of antibiotic resistance genes and pathogens associated with fish health. BMC Microbiol. 2020;20(1):249.
Gajardo K, Rodiles A, Kortner TM, Krogdahl A, Bakke AM, Merrifield DL, et al. A high-resolution map of the gut microbiota in Atlantic salmon (Salmo salar): a basis for comparative gut microbial research. Sci Rep. 2016;6:30893.
Huang Q, Sham RC, Deng Y, Mao Y, Wang C, Zhang T, et al. Diversity of gut microbiomes in marine fishes is shaped by host-related factors. Mol Ecol. 2020;29(24):5019–34.
Wong S, Rawls JF. Intestinal microbiota composition in fishes is influenced by host ecology and environment. Mol Ecol. 2012;21(13):3100–2.
Llewellyn MS, McGinnity P, Dionne M, Letourneau J, Thonier F, Carvalho GR, et al. The biogeography of the atlantic salmon (Salmo salar) gut microbiome. ISME J. 2016;10(5):1280–4.
Fu H, Zhang L, Fan C, Liu C, Li W, Cheng Q, et al. Environment and host species identity shape gut microbiota diversity in sympatric herbivorous mammals. Microb Biotechnol. 2020;14(4):1300–15.
Kim PS, Shin NR, Lee JB, Kim MS, Whon TW, Hyun DW, et al. Host habitat is the major determinant of the gut microbiome of fish. Microbiome. 2021;9(1):166.
Guivier E, Pech N, Chappaz R, Gilles A. Microbiota associated with the skin, gills, and gut of the fish Parachondrostoma toxostoma from the Rhône basin. Freshw Biol. 2020;65(3):446–59.
Xing M, Hou Z, Yuan J, Liu Y, Qu Y, Liu B. Taxonomic and functional metagenomic profiling of gastrointestinal tract microbiome of the farmed adult turbot (Scophthalmus maximus). FEMS Microbiol Ecol. 2013;86(3):432–43.
Wu S, Wang G, Angert ER, Wang W, Li W, Zou H. Composition, diversity, and origin of the bacterial community in grass carp intestine. PLoS ONE. 2012;7(2):e30440.
Ruiz-Rodriguez M, Scheifler M, Sanchez-Brosseau S, Magnanou E, West N, Suzuki M, et al. Host species and body site explain the variation in the microbiota associated to wild sympatric Mediterranean teleost fishes. Microb Ecol. 2020;80(1):212–22.
Minich JJ, Harer A, Vechinski J, Frable BW, Skelton ZR, Kunselman E, et al. Host biology, ecology and the environment influence microbial biomass and diversity in 101 marine fish species. Nat Commun. 2022;13(1):6978.
Dallas JW, Warne RW. Captivity and animal microbiomes: potential roles of microbiota for influencing animal conservation. Microb Ecol. 2023;85(3):820–38.
Clemente JC, Pehrsson EC, Blaser MJ, Sandhu K, Gao Z, Wang B, et al. The microbiome of uncontacted Amerindians. Sci Adv. 2015;1(3):e1500183.
Steiner K, Heasman K, Laroche O, Pochon X, Preece M, Bowman JP, et al. The microbiome of Chinook salmon (Oncorhynchus tshawytscha) in a recirculation aquaculture system. Aquaculture. 2021;534:736227.
Xiao F, Zhu W, Yu Y, He Z, Wu B, Wang C, et al. Host development overwhelms environmental dispersal in governing the ecological succession of zebrafish gut microbiota. NPJ Biofilms Microbiomes. 2021;7(1):5.
Sevellec M, Laporte M, Bernatchez A, Derome N, Bernatchez L. Evidence for host effect on the intestinal microbiota of whitefish (Coregonus sp.) species pairs and their hybrids. Ecol Evol. 2019;9(20):11762–74.
Boutin S, Sauvage C, Bernatchez L, Audet C, Derome N. Inter individual variations of the fish skin microbiota: host genetics basis of mutualism? PLoS ONE. 2014;9(7):e102649.
Miyake S, Ngugi DK, Stingl U. Phylogenetic diversity, distribution, and cophylogeny of giant bacteria (Epulopiscium) with their surgeonfish hosts in the Red Sea. Front Microbiol. 2016;7:285.
O’Brien PA, Tan S, Yang C, Frade PR, Andreakis N, Smith HA, et al. Diverse coral reef invertebrates exhibit patterns of phylosymbiosis. ISME J. 2020;14(9):2211–22.
Trevelline BK, Sosa J, Hartup BK, Kohl KD. A bird’s-eye view of phylosymbiosis: weak signatures of phylosymbiosis among all 15 species of cranes. Proc Biol Sci. 1923;2020(287):20192988.
Leigh BA, Bordenstein SR, Brooks AW, Mikaelyan A, Bordenstein SR. Finer-scale phylosymbiosis: insights from insect viromes. mSystems. 2018;3(6):10–1128.
Pollock FJ, McMinds R, Smith S, Bourne DG, Willis BL, Medina M, et al. Coral-associated bacteria demonstrate phylosymbiosis and cophylogeny. Nat Commun. 2018;9(1):4921.
We thank the Ontario Ministry of Natural Resources and Forestry (MNRF), especially staff at the Glenora Fisheries Station for helping during the sampling. We also thank Brent Nawrocki and Amy Weinz for helping with Detroit River sampling.
This research received financial support from Canada’s Natural Sciences and Engineering Research Council (NSERC; 819047) and Genome Canada (GEN-FISH) to DDH and JS received an Ontario Trillium Scholarship (OTS).
Ethics approval and consent to participate
All fish handling was performed by other agencies; however, we also had Canadian Council for Animal Care permitting for handling and humanely euthanizing fish.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Relative abundance of bacterial community composition presented at the phylum level for gut, skin, and water microbiota (samples are combined across sample types). Phyla with less than 0.01% of relative abundance are combined and presented as “others”). Figure S2. Bacterial community composition (relative abundance at the family level) for gut, skin, and water microbiomes across all fish species collected at three sites in the Great Lakes (Lake Erie, Lake Ontario and Detroit River). Bacterial families with less than 0.1% relative abundance are combined and presented as “others”. Fig S3. LASSO regression analysis was used to indentify the best predictor variables for gut samples based on their alpha (Shannon entropy, PD, Chao1) and beta diversity (PCoA 1- 5) indices. Inside each bar is showing the coeffitient value and X axis is showing the importance of the predictor variable for alpha and beta diversity indices. Diet, location and fish species was identified as the best predictor for most of the diversity indecies. Fig S4. LASSO regression analysis was used to indentify the best predictor variables for fish skin samples based on their alpha (Shannon entropy, PD, Chao1) and beta diversity (PCoA 1- 5) indices. Inside each bar is showing the coeffitient value and X axis is showing the importance of the predictor variable for alpha and beta diversity indices. Location and fish species was identified as the best predictor for most of the diversity indecies. Figure S5. Scatterplot of pairwise host phylogenetic distance vs pairwise Bray Curtis dissimilarity for both gut (a) and skin (b) samples. Samples were combined within host species. Host phylogenetic distance was estimated using of CO1 and CytB mitochondrial gene sequences. Table S1. Summary of Great Lakes fish species sampled for gut and skin microbiome. We provide a description of the sample locations and total sample size. Table S2. Comparison of differentially abundant bacterial taxa at the family level for gut and skin microbiomes across all fish species and sample locations using DESeq2 method (Benjamini-Hochberg false-discovery rate [BH FDR] 0.05, |log2fold change| > 2). Positive log2 FC indicate higher abundance in skin samples and negative log2 FC specify higher abundance in gut samples.
About this article
Cite this article
Sadeghi, J., Chaganti, S.R., Johnson, T.B. et al. Host species and habitat shape fish-associated bacterial communities: phylosymbiosis between fish and their microbiome. Microbiome 11, 258 (2023). https://doi.org/10.1186/s40168-023-01697-6