Samples and sequencing results
Five different gilthead sea bream families selected by growth derived from the PROGENSA (Spanish selection program of gilthead sea bream) broodstock were used in this study (F3 generation) [27,28,31]. The five families were then grouped in three sets, named suprafamilies, by their different growth trajectories as previously described [31]: fast (families e5e2 and e6e2, constituting suprafamily e5e6), intermediate (family c2c7) and slow (families c4c2 and e4e1, constituting suprafamily c4e4). It is important to bear in mind that these families are not of clonal origin, they have substantial genomic heterogeneity, and were grouped considering phenotypic characters. These animals fed a control (D1) or a well-balanced plant-based diet (D2) for 9 months, were kept in common garden in order to eliminate the disturbing effects that could appear by rearing the analyzed families in different tanks. Six fish of each family and diet were used to sample the adherent microbiota of the anterior intestine. In total, 60 samples were sequenced: 24 samples for the fast- and slow-growth suprafamilies and 12 samples for the intermediate-growth group. More details on the rearing of the fish and sampling can be found in the methods section.
After Illumina sequencing of the 60 samples, three samples were eliminated from further analysis due to low quality of the reads (one from suprafamily c4e4 fed with D1 and two from suprafamily e5e6 fed with D1). The remaining 57 samples yielded 8,803,202 high quality reads, with a mean of 154,442 reads per sample, ranging from 99,040 to 303,044 (Additional File 1: Table S1). The reads were assigned to 2327 OTUs at 97% identity threshold. Almost half of these OTUs (49.8%) were classified up to the level of species, 89.3% to the level of genus and more than 95% to the levels of family (95.1%), order (97.5%), class (98.8%) and phylum (99.9%). Rarefaction analysis showed curves that approximated saturation (horizontal asymptote), thus a good coverage of the bacterial community was achieved and the number of sequences for analysis was considered appropriate (Additional File 2: Figure S1).
Microbiota diversity and composition
When comparing the bacterial diversity and composition of the three groups of families, regardless of the diet, no significant differences were found in richness (Chao1 and ACE) and diversity (Shannon and Simpson) indexes (Table 1). The overall bacterial composition at the phylum level was also not different among groups (Kruskal-Wallis test, Dunn’s post-test), indicating that all these animals harbor the typical microbiota expected in gilthead sea bream intestines (Figure 1). In all groups of families, Proteobacteria was the most abundant phylum, constituting ³ 50% of the total resident bacteria in the anterior intestine. Vibrionaceae and Oceanospirillaceae families were the most abundant Proteobacteria in all groups, representing 18.9% and 11.6%, 19.5% and 14.5%, and 22.3% and 8.2% of the bacterial communities from e5e6, c2c7 and c4e4, respectively (Additional File 3: Figure S2). The second most abundant phylum was Firmicutes (³ 27% in all groups), with the families Clostridiaceae, Listeriaceae and Staphylococcaceae representing ~20% of the total microbiota in all groups. The phylum Actinobacteria (³ 7%, particularly families Corynebacteriaceae, Micrococcaceae and Propionibacteriaceae), followed by Cyanobacteria (³ 1.5%) and Bacteroidetes (³ 1%), were also abundant in the three groups (Additional File 3: Figure S2).
Diet and family effects on microbiota
A PERMANOVA test was used to study differences in bacterial composition by diet, family and the interaction diet × family. Not taking into account the diet, the family variable showed no statistically significant effect on the microbiota(P = 0.169, F = 1.308, R2 = 0.044). However, statistical differences were detected when comparing animals fed different diets, regardless of the genetic background (P = 0.041, F = 1.913, R2 = 0.032). Although R2 values detected were quite low, they were in line with what was reported in other microbiota studies [61]. This is due to the complexity and variability of microbiota samples. To corroborate and study in more detail these differences, a PLS-DA model was constructed and statistically validated. In Figure 2, the PLS-DA model shows a clear separation of fish fed D1 or D2 along component 1 (79.69%). The PLS-DA model was successfully validated with a permutation test discarding the possibility of over-fitting of the supervised model (Additional File 4: Figure S3A and B). These results highlight that the diet has a significant impact on the composition of the adherent bacterial communities of the anterior intestine, masking differences among families.
The PERMANOVA results of the interaction showed that significant differences appeared when taking into account both, genetic background (suprafamilies) and diet (P = 0.017, F = 1.843, R2 = 0.062). The PLS-DA models of the different suprafamilies showed a separation of fish fed D1 from those fed D2 (Figure 3 A, B and C). The permutation tests (Additional File 4: Figure S3 C-H) showed that only the models for c2c7 and c4e4 could be validated, while the model for e5e6 was over-fitted. This means that significant differences in bacterial communities due to diet are mainly occurring in the groups of families with intermediate (c2c7) and slow (c4e4) growth.
To determine which groups of bacteria were driving these separations with the diet changes, a more detailed analysis of the variable importance in projection (VIP) was performed throughout a heatmap representation. Hierarchical clustering of samples was applied and the minimum VIP values significantly driving the separation of the groups in the model were calculated being the OTUs within these values selected for further analysis. Differences in diet were mainly changing 128 OTUs (VIP > 1.4) in suprafamily e5e6, 158 OTUs (VIP > 1) in family c2c7 and 232 OTUs (VIP > 1) in suprafamily c4e4 (Figure 3 D, E and F). A detailed list of the VIPs can be found in Additional File 5: Table S2.
Microbiota changes due to diet depend on the genetic background
A detailed study of the different OTUs driving the separation due to diet in the different suprafamilies revealed that most were family exclusive, and only 10 of these OTUs were changing in all groups (Figure 4). The number of OTUs exclusively changing in groups e5e6, c2c7 and c4e4 were 89, 91 and 170, respectively. The greatest difference was found when studying the abundance of these exclusive OTUs in the different groups. Exclusive OTUs changing in groups c2c7 and c4e4 accounted for 25.8% and 30.4%, respectively, of the total bacterial composition. Interestingly, the 89 OTUs changing in group e5e6 only accounted for 9.2% of the total bacterial communities in these animals.
Figure 5 shows the most abundant bacteria of those that exclusively drive the separation by diet in the different groups of families. In suprafamily e5e6 Actinobacteria (particularly Kokuria and Micrococcus) and Cyanobacteria decreased and some Firmicutes (Bacillus sp.) disappeared when animals were fed D2. Regarding Proteobacteria, some genera increased (Afipia, Bradyrhizobium) and others decreased (Pseudoalteromonas, Photobacterium, Vibrio) with the plant-based diet (Figure 5 A). In the intermediate growth family (c2c7), some Actinobacteria (Arthrobacter sp.) appeared, and different species of the genus Staphylococcus (Firmicutes) increased with the D2 diet. Proteobacteria again appeared more variable, with Novosphingobium, Ralstonia, and Marinomonas decreasing, and Sphingomonas, Haemophilus, Photobacterium, and Stenotrophomonas increasing with D2 (Figure 5 B). In suprafamily c4e4, the Actinobacteria genera Corynebacterium and Zhihengliuella, and the Firmicutes Bacillus increased, while the Firmicutes genera Brochothrix decreased with D2. Within Proteobacteria, the genera Caulobacter, Vibrio, and some Photobacterium (P. damselae and P. leiognathi) decreased and Mesorhizobium and other Photobacterium sp. increased with D2 (Figure 5 C).
Predicted metabolic differences are larger in groups with fewer bacteria changes
In an attempt to evaluate the biological significance of the differences induced by diets in the microbiota of the different groups of families, pathway analysis was performed using the inferred metagenomes of the OTUs driving the separation by diet. The results showed that in e5e6 (fast growth) 59 pathways could be significantly changing with the different diets. Whereas, in c2c7 (intermediate growth), 84 pathways showed differences, and only 15 pathways were predicted to be changing in the slow growth group (c4e4) (Additional File 6: Table S3). Very few overlaps were found when comparing the pathways predicted to be changing in each group of families when fed D2 in comparison to D1. In group c2c7, the number of over- and under-represented pathways was balanced, whereas most of the changing pathways in group e5e6 were under-represented, and in group c4e4 were over-represented in D2 (Figure 6 A and B). Among the most under-represented pathways in e5e6 there were several pathways related to infection, inflammation or activation of the immune system (Bacterial invasion of epithelial cells, Staphylococcus aureus infection, Arabinogalactan biosynthesis or RIG-I-like receptor signaling pathway). The only highly over-represented pathways (logFC > 1) were Flavone and flavonol biosynthesis and Glycosphingolipid biosynthesis (Figure 6 C). In group c2c7, most of the highly differentially abundant pathways (logFC > |1|) were over-represented, highlighting Arabinogalactan biosynthesis, Sesquiterpenoid and triterpenoid biosynthesis and Complement and coagulation cascades. In c4e4 very few pathways were differentially abundant, with only 4 pathways with logFC > |1|, which were highly over-represented: Indole alkaloid biosynthesis, Biosynthesis of enediyne antibiotics, Furfural degradation and RIG-I like receptor signaling pathway (Figure 6 C).
Genetic background effects on disease resistance
Taking into account the differences detected in intestinal morphology, plasticity and microbial composition between fish selected for fast and slow growth, the different response upon an intestinal pathogen was also evaluated. Two families were selected for this study, e5e2 (fast growth) and c4c3 (slow growth). These animals were kept together in the same tanks, fed D2 diet and challenged with the intestinal myxozoan parasite Enteromyxum leei. This parasite lives and divides in the intestinal epithelium inducing intestinal damage, impaired nutrient absorption and anorexia. Therefore, among the most common disease signs induced by this parasite are diminished growth and weight loss [62]. After 70 days post-infection (dpi) biometric parameters showed that, as expected, the weight gain, relative to the weight at 0 dpi, of the control uninfected fish was significantly higher in the fast-growth family (12.6%) than in the low-growth group (4.7%). Upon infection, the fast-growth family showed a significantly lower weight loss (5%) than the slow growth group (11%) (Figure 7 A). The prevalence of infection was 78% and 87% for groups e5e2 and c4c3, respectively. Moreover, the fast-growth group showed significantly lower intensity of infection and parasite abundance values when compared to c4c3, the low-growth family (Figure 7 B and C).