Ecological Specialization Within a Carnivorous Fish Family Is Supported by a Herbivorous Microbiome Shaped by a Combination of Gut Traits and Specific Diet

Animals have been developing key associations with micro-organisms through evolutionary processes and ecological diversification. Hence, in some host clades, phylogenetic distance between hosts is correlated to dissimilarity in microbiomes, a pattern called phylosymbiosis. Teleost fishes, despite being the most diverse and ancient group of vertebrates, have received little attention from the microbiome perspective and our understanding of its determinants is currently limited. In this study, we assessed the gut microbiome of 12 co-occurring species of teleost representing a large breadth of ecological diversity and originating from a single family (i.e., the Sparidae). We tested how host evolutionary history, diet composition and morphological traits are related to fish gut microbiome. Despite fish species having different microbiomes, there is no phylosymbiosis signal in this fish family, but gut length and diet had a strong influence on the microbiome. We revealed that the only species with a specialized herbivorous diet, Sarpa salpa had a 3.3 times longer gut than carnivorous species and such a long gut favor the presence of anaerobic bacteria typical of herbivorous gut microbiomes. Hence, dietary uniqueness is paired with both unique gut anatomy and unique microbiome.


INTRODUCTION
All animals evolved and diversified within a microbial world, and consequently developed a myriad of associations with micro-organisms (McFall-Ngai et al., 2013). The abundant and diverse microbial communities present within and on the surface of animal bodies, i.e., the microbiomes, play a central role in the fitness of their host, providing benefits that include successful development of the immune system and protection against pathogens, but also enhanced food processing and nutrient absorption (McKenney et al., 2018a;Moran et al., 2019). These long histories of animalmicrobes interactions resulted for many clades in phylosymbiosis, i.e., the eco-evolutionary pattern where phylogenetically related hosts tend to have more similar microbiomes than distantly related ones (Brooks et al., 2016;Lim and Bordenstein, 2019). Further, phylosymbiosis signal in the gut microbiome was shown to be stronger when hosts share a recent common ancestor (i.e., intra clade) than when shared ancestry is more ancient (i.e., inter clades) (Groussin et al., 2017). Phylosymbiosis has been evidenced in many animals, notably those for which vertical transmission of the microbiome has been reported, such as in the gut of insects (Brucker and Bordenstein, 2012;Haddow et al., 2013;Sanders et al., 2014;Brooks et al., 2016), and mammals (Phillips et al., 2012;Sanders et al., 2014), or on the skin of mammals  and elasmobranchs (Doane et al., 2020). However, evidences are far less pronounced in animals with horizontal microbiome transmission, with the report of weak to no support for skin phylosymbiosis in amphibians (Bletz et al., 2017) and marine fishes (Chiarello et al., 2018(Chiarello et al., , 2019Doane et al., 2020).
Marine teleost fishes gather >10,000 species and >100 families, and a large range of ecological strategies, with for instance diet ranging from herbivory to top-predators even within some families (e.g., Labridae, Sparidae). However, to date only a few studies have assessed the microbiome of several fish species from a single ecosystem using sequencingbased approaches (Givens et al., 2015;Miyake et al., 2015;Sullam et al., 2015;Sylvain et al., 2019;Ruiz-Rodríguez et al., 2020). Consequently, the presence or absence of gut phylosymbiosis in the gut microbiome of marine fishes has never been tested. However, a case of co-phylogeny, i.e., congruent phylogeny between hosts and individual members of their microbiome, has been reported between Acanthuridae (i.e., surgeonfishes) and the bacteria Epulopiscium spp., a genus from the class Clostridia in the Firmicutes phylum (Miyake et al., 2015(Miyake et al., , 2016. This specific host-microbes association between a bacterial clade and a fish family could still be coupled with low association between the whole microbiome and host phylogeny, especially when considering fishes from different families. Beside hosts identity and long term association with microorganisms that result in phylosymbiosis, the composition of the gut microbiome in marine teleosts is influenced by many factors (Sullam et al., 2012;Nikouli et al., 2020). Environment plays a role by influencing bacterioplankton communities, notably for marine fishes who continuously drink seawater to maintain osmolarity (Baldisserotto et al., 2019). This flow of environmental micro-organisms contributes to the early stages of gut microbial colonization and will influence microbiome composition at later stages (Egerton et al., 2018;Nikouli et al., 2019). Hosts diet also shapes the gut microbiome of fishes as food composition influences the type of resources available to gut micro-organisms and places the hosts in contact with the microbiomes of their prey (Sullam et al., 2012). Additionally, fishes have specific gut morphological traits that depend on their diet with, for instance, the presence of a gizzard in some species that consume sediment or a differentiated hindgut in herbivores which result in anaerobic conditions that favor certain micro-organisms (Mountfort et al., 2002;Clements et al., 2009;Egerton et al., 2018). Despite these recent advances, the overall role of diet and anatomy in shaping the gut microbiome still remains to be unfolded, especially within families where ecological diversification occurred.
In this study, we assessed the gut microbiome of 12 species from a single family (Sparidae) representing a large breadth of diet and sampled from a single location (i.e., the North-Western Mediterranean Sea). We then assessed how host evolutionary history, diet and gut anatomy are related to gut microbiome.

Samples Collection
Fishes were collected in the shallow habitats (2-10 m) of the North-Western coast of Gulf of Lions (Hérault, France, Mediterranean Sea) between June 2018 and February 2019. Fishes were caught at night by artisanal fishermen at the end of the night using gillnets, stored on ice on the boat and were processed within 3 h after being caught. Fish were dissected using tools cleaned with 70 • ethanol, the gut was stretched and a picture was taken for morphometric measurements. The last third of the gut, before the cloaca, was squeezed out on a piece of parafilm paper using a sterile pipette. If this gut content was made of digested food (i.e., not only mucus) it was then homogenized and collected in a 3 mL cryotube before storage at −80 • C until DNA extraction. Gut microbiome from a total of 48 adult individuals from 12 Sparidae species were sampled (2 to 8 individuals per species, Figure 1 and Supplementary Table 1).

Phylogenetic, Dietary and Morphological Dissimilarity Between Hosts
Phylogenetic distances between fish species were estimated as the cophenetic distance extracted from the only multi-locus timecalibrated phylogeny of Sparidae (Santini et al., 2014; Figure 1).
Diet of all species was recovered from food items listed in FishBase (Froese and Pauly, 2003), selecting only reports from the Mediterranean Sea and only included adults records. FishBase uses different levels of food items classification and we used the level that allows discriminating between groups of plants and invertebrates. This included the following food items: benthic algae, finfish, planktonic crustaceans, planktonic invertebrates, small benthic crustaceans, large benthic crustaceans, sessile invertebrates, echinoderms, annelids, cnidarians, mollusks, and shelled mollusks (Figure 1). The records were averaged for each species and resulted in a matrix containing the proportions of each item in the diet of each species. Species were classified into three groups according to their main type of prey (plants, plankton, and macrofauna). The prey proportion matrix was used to compute diet dissimilarity between species using the Bray-Curtis index.
Gut morphology was characterized using gut length (mm) and gut length relative to fish standard length (from snout to basis of the caudal fin). Fishes morphological features were measured at the individual level and species-level traits were estimated as the average across individuals from that species (Figure 1). Inter-species morphological dissimilarity was computed using Euclidean distance. Phylogenetic conservatism of gut traits was tested using Abouheif test (Abouheif, 1999), which is designed to detect phylogenetic autocorrelation in a quantitative trait, as implemented in the adephylo package (Jombart et al., 2010).  (from Santini et al., 2014). N is the number of individuals analyzed per species. Diet classification and proportions of food items in the diet were collected in FishBase.

DNA Extraction
DNA extractions were performed in the molecular biology platforms of the MARBEC laboratory (UMR 9190) 1 and at the GenSeq platform 2 , using the Qiagen MagAttract PowerSoil DNA KF Kit, selected for its compliance with the Earth Microbiome Project (Marotz et al., 2017). Extractions were performed in 96 wells plates in which three wells were left empty to serve as negative controls and three wells were loaded using standard mock communities (ZymoBIOMICS Microbial Community DNA Standards II, Zymo Research). These standards were used to evaluate the quality of our sample processing pipeline and to identify potential contaminants from the reagents. Extraction wells were loaded with one spatula of homogenized gut content (∼0.25 g). DNA extraction protocol followed the manufacturer instructions. It included a chemical lysis consisting of 60 µL of manufacturer provided SL solution heated at 60 • C and added before a bead beating step, which corresponded to 10 min at 20 Hz using a Tissue Lyzer, then the plate was flipped and shake again for 10 min at 20 Hz. DNA recovery was based on magnetic beads and automated with a Kingfisher Flex robot. DNA was eluted in 200 µL of elution buffer before quantification of DNA quantity and quality using a NanoDrop 8000 spectrometer. Extracted

PCR Amplification
PCR amplification was done using primers selected for their compliance with the Earth Microbiome Project (Parada et al., 2016): 515F-Y (5 -GTGYCAGCMGCCGCGGTAA) and 926R (5 -CCGYCAATTYMTTTRAGTTT). The targeted sequence was 411 bp and corresponded to the V3-V4 regions of the prokaryotic 16S rRNA gene. PCR amplification was carried out in 96 well plates in triplicate for each DNA extract and was done in a 25 µL reaction volume. The PCR mix consisted of 9.75 µL of water, 0.75 µL of DMSO, 0.5 µL of each primer (concentration), 12.5 µL of Phusion ready-to-use Taq mix (Phusion High-Fidelity PCR Master Mix with GC Buffer), and 1 µL of DNA. The PCR cycle consisted of 35 cycles of 10 s denaturation at 98 • C, 1 min annealing at 58 • C, and 1 min 30 s of extension at 72 • C. Final extension was held for 10 min at 72 • C before keeping the reaction at 4 • C. The success of PCR amplification was checked on 2% agarose gel in TAE buffer and using a 100 bp DNA ladder. The wells left empty during DNA extraction served as negative controls for contamination of the PCR reactions. PCR triplicates were pooled and stored at −20 • C before sequencing. Amplicons library was constructed by the Genotoul platform 3 and sequencing was carried out using an Illumina MiSeq (2 × 250 bp) sequencer.

Amplicon Sequencing and Sequences Processing
The demultiplexed reads were processed using the R software environment and the package dada2 (Callahan et al., 2016). Briefly, the quality of the reads for each sample was inspected using graphic representations of their quality scores and reads were filtered based on their length and quality. Amplicon sequence variants (ASVs) were inferred using the dada algorithm (Divisive Amplicon Denoising Algorithm) after pooling dereplicated reads from all samples. Then, forward and reverse reads were merged and chimeric sequences were removed. After these steps, a third of the original sequences were kept for further analyses, on average. The taxonomic classification of ASVs was performed with the naive Bayesian RDP classifier implemented in dada2 and using the SILVA reference database nr_V132 (10.5281/zenodo.1172783). As many ASVs were not affiliated at the genus level, ASVs sequences were blasted against the NCBI 16S rRNA bacterial and archaeal database, and the best hit was used to assign ASVs at the genus level.
Several data cleaning steps were performed to remove poorly characterized and rare ASVs. First, the mock communities and blank samples were used to identify contamination from the reagents (e.g., extraction kit, polymerase), and these ASVs were removed from the data set (e.g., Ralstonia, Rhizobium). Second, ASVs not assigned to the bacterial domain, unclassified at the Class level, assigned to chloroplasts or mitochondria were filtered out. Third, ASVs present in less than two samples were removed.
Amplicon sequence variants sequences were aligned using mafft implemented in Qiime2 before being inserted in the Greengene reference phylogenetic tree (Janssen et al., 2018). The tree was then ultrametricized using pathd8 4 .

Identification of the Core Microbiome
The core microbiome of Sparidae fishes, i.e., the most frequent and/or abundant bacterial members of their microbiome, was estimated for the 12 Sparidae species altogether using the core identification algorithm from Magurran and Henderson (2003). This method accounts for both ASVs occurrence and abundance across communities and is based on the comparison of the observed abundance-occurrence distribution of ASVs with a random distribution under a stochastic Poisson model. We selected two individuals from each species in order to give equal weight to each species and the procedure was repeated 1000 times using different combinations of individuals for the species represented by more than two individuals. Then, we defined the size of the core as the median number of core ASVs (n = 842 ASVs) observed across all 1000 combinations and thus the core microbiome of Sparidae was defined as the 842 ASVs that were considered as core the highest number of times across the 1000 replicates. The abundance-occurrence distribution of

Relationships Between Microbiome Dissimilarity and Host Ecological, Morphological and Phylogenetic Dissimilarities
In order to test for phylosymbiosis signal and identify the determinants of the core microbiome in Sparidae fishes, we computed microbiome dissimilarity using complementary diversity indices at several ranks of the bacterial taxonomy. Interspecies dissimilarities were estimated after taking the average taxa abundance across all the individuals from each species.
We used six dissimilarity indices from the Hill numbers framework (Chao et al., 2014), as these combine taxonomic and phylogenetic dissimilarities within a common mathematical framework. In addition, Hill numbers allow the estimation of composition and structure dissimilarity by giving more or less weight to taxa abundances, depending on the value of a parameter, the order q. These six indices were computed for six ranks across the bacterial taxonomy (Phylum, Class, Order, Family, Genus, and ASVs). For that, all the ASVs were aggregated at a given taxonomic rank and we computed dissimilarity on aggregated data. In a second approach, we tested whether intrataxa dissimilarity was related with host dissimilarity. Here, the dataset was split into individual bacterial taxa for the highest taxonomic ranks (Phylum, Class, and Order) and dissimilarity was computed using only the ASVs from the focus taxa. For this analysis we used only the taxa that were represented by at least 10 ASVs (N = 31 taxa corresponding to 7 phyla, 10 classes, and 14 orders).
We used Mantel tests implemented in the R package ecodist (Goslee and Urban, 2007) to test the relationships between microbiome and host (phylogenetic, dietary, and morphological) dissimilarity matrices.

Identification of Gut Microbiome Determinants in Sparidae Fishes
Differences between species and diet types in terms of gut microbiome structure were tested using multivariate Welch MANOVA (W d test, Hamidi et al., 2019) implemented in the R package MicEco (Russel, 2020). This approach has the advantage, over classically used PERMANOVA, of being robust to heteroscedasticity in the data, i.e., differences in multivariate groups dispersion. Association of host species and diet types with particular bacterial taxa was assessed using the LEfSe approach (Segata et al., 2011).
We used constrained analysis on principal coordinates (CAP analysis, Anderson and Willis, 2003) to determine which of the diet, gut anatomy and host phylogeny explained the better the differences in microbiome structure. Then, we identified which microbial taxa correlated with particular diet items or traits. Again, this was done for all the combinations of dissimilarity indices and bacterial taxonomic ranks. CAP models were defined to explain microbiome dissimilarity by (i) the proportions of the 12 dietary items, (ii) the values of the two morphological traits, and (iii) the phylogenetic relatedness of the hosts. This later was included as phylogenetic eigenvectors resulting from the decomposition of the cophenetic distance matrix between hosts performed using the R package PVR (Diniz-Filho et al., 1998;Santos, 2018). The individual effect of each explaining variable was tested using ANOVA and the anova.cca function in R.

Phylogenetic Relationships and Morphological Similarity of Mediterranean Sparidae
The 12 studied Sparidae species shared a common ancestor around 60 million years ago and separated in two main branches which do not contain recent (<12 Mya) speciation events (Figure 1). Two genera contained more than one species, Diplodus (n = 4) and Pagellus (n = 2), but while the former is monophyletic, the latter is spread in the two main branches of the family tree. From a dietary aspect, ten species are considered as macrofauna hunters according to FishBase (including all Diplodus and Pagellus species), but this rough classification blurs a high dissimilarity in terms of diet composition, with mostly zooplankton for Oblada melanura to omnivory on various proportion of crustacean, annelids, and/or echinoderms among other species (Figure 1). The two remaining species have a unique diet: Boops boops is a selective plankton feeder and Sarpa salpa is a strict herbivore feeding on benthic macrophytes. Set aside S. salpa, only B. boops and Diplodus puntazzo consume more than 15% of plant material. B. boops is the closest species to S. salpa (divergence ∼18.4 Mya) and those two species also have the highest values for all gut traits (Supplementary Table 2). Gut length and relative gut length showed a significant phylogenetic signal (Abouheif test, p-value < 0.05, Supplementary Table 3

Microbial Taxa Associated With Host Species and Diet Types
Most of the 189 key ASVs identified by the LEfSe analysis ( Table 2) were associated with S. salpa (n = 26), its closest relative B. boops (n = 17) and Pagellus erythrinus (n = 11), or with selective plankton feeding (n = 71) and aquatic plants grazing (n = 96). The other species (maximum 3 ASVs) and macrofauna hunting (4 ASVs) were associated with a much lower number of ASVs. S. salpa and its diet were mostly associated with members of the Clostridiales (Firmicutes) from the families Lachnospiraceae (notably the genus Tyzerella) and Ruminococcaceae (Angelakisella, Ruminiclostridium), the Bacteroidetes family Rikenellaceae (genera Alistipes and Rikenella) and several Desulfovibrionaceae ASVs. B. boops and its diet were associated with the Firmicutes family Ruminococcaceae (genera UCG-008) and the Lachnospiraceae Tyzerella, the Rikenellaceae (genera Alistipes and Rikenella) and the genus Desulfovibrio. P. erythrinus was mostly associated with Planctomycetes from the Pirellulaceae family (Blastopirellula and Rhodopirellula) but also Alphaproteobacteria (Ahrensia, Roseobacter).

Relationships Between Host and Gut Microbiome Dissimilarity
We did not observe any phylosymbiosis signal as phylogenetic distance between fish species was not significantly correlated with their microbiome dissimilarity for any of the 36 combinations of dissimilarity indices and bacterial taxonomy ranks tested (Mantel tests, p-values > 0.05, Figure 2). Diet dissimilarity was related to microbiome dissimilarity in only one out of 36 tests, i.e., using taxonomic microbiome dissimilarity computed on presence-absence data at the bacterial order level (q = 0). Morphological dissimilarity based on gut traits was significantly related to microbiome dissimilarity in five instances (out of 72 Mantel tests, p-values < 0.05), all with abundance-weighted phylogenetic dissimilarity (q = 2) Only the bacterial families with a relative abundance higher than 10% in at least one species are represented. Empty cells means relative abundance was lower than 1%.
computed at higher ranks of bacterial taxonomy (Phylum = 1, Class = 2, Order = 2). We found 11 significant relationships between intrataxa abundance-weighted (q = 2) phylogenetic dissimilarity (Supplementary Figure 5) and host dissimilarity (phylogenetic, diet and morphological). Half of these relationships were negative (6/13) and, overall, the strength of the relationships was weak. Only one significant relationship involved host phylogenetic distance, which was found to be negatively related with dissimilarity in the order Desulfovibrionales. The strongest association was found between phylogenetic dissimilarity within the Bacillales order and Bacilli class (Firmicutes phylum) and dissimilarity in gut traits.

Association Between Sparidae Microbiome, Gut Morphological Traits, Diet Items, and Host Phylogeny
More than half (n = 19, 53%) of the CAP models testing the relationship between gut traits, diet composition, host phylogeny, and microbiome structure were significant (i.e., p-value < 0.05 in 19 out of 36 models, Supplementary Table 4). The first two axes of the significant CAP models explained between 41 and 91% of the variation in gut microbiome (Supplementary Figure 6). Among these 19 significant CAP models, the most influential variable was by far gut length, as it exhibited the highest number of significant relationships with microbiome components (n = 19) and the highest F-values (Supplementary Table 4). The proportion of Echinoderms and Mollusks in the diet was the second and third most influential traits, followed by relative gut length. Phylogenetic eigenvectors did not explain significant proportion of the variance in any of the CAP models.
Hereafter, we focus on the results obtained using abundanceweighted taxonomic dissimilarity (q = 1) as this method yielded a significant CAP model five out of six taxonomic ranks (i.e., not at the Class) and explained a high proportion of variance in microbiome composition, on average (88 ± 8%). The first CAP axis explained 26.2% of the variation and mostly separated the species with a long gut that do not consume Mollusks, i.e., S. salpa, from the other species ( Figure 3A). The second CAP axis explained 18.6% of the variation and separated two species, P. erythrinus and B. boops, from the others.
We found 131 significant correlations between gut traits or the proportions of diet items and the abundance of bacterial taxa at different ranks. Only four were negative, which corresponded to the negative association between gut length and Proteobacteria, Betaproteobacteriales and Burkholderiaceae, and between Annelids and Firmicutes. Gut length was the variable that correlated with the highest number of bacterial taxa (Figure 3B), i.e., 60 including 48 ASVs, and was the only variable associated with taxa across all the taxonomic ranks (i.e., phylum to ASVs). The strongest correlations included those with the Clostridia class and the Clostridiales order (r = 0.83, p-value < 0.001), the Lachnospiraceae family (r = 0.81, p-value < 0.001) and one of its ASV (r = 0.87, p-value < 0.001). The associations between ASVs and groups of samples (species or diet categories) were estimated using LEfSe analysis and the differences between groups were tested using Kruskal-Wallis test (FDR corrected). SPF, selective plankton feeder; APG, aquatic plants grazer; MH, macrofauna hunter.
Frontiers in Marine Science | www.frontiersin.org FIGURE 2 | Gut traits explain differences in the phylogenetic structure of the microbiome better than host phylogeny or diet. Relationships between dissimilarity of the gut microbiome (Y axes) and phylogenetic, morphological, and dietary dissimilarity (columns) between host fishes (X axis) using Mantel tests (for each subplot: r = correlation between X and Y variables and p-value = p-value of the Mantel test, bolded r and p-values correspond to significant relationships). Gut microbiome dissimilarity, here estimated using abundance-weighted Hill numbers phylogenetic index with q = 2, was computed at six different levels of taxonomic resolution (row names on the right of the figure). Each point represents a pair of fish species. Pairs within the Diplodus genus are depicted in blue, pair within the Pagellus genus is depicted in pink, pair between Sarpa salpa and Boops boops is depicted in black and other pairs that included S. salpa or B. boops are depicted in green and orange, respectively.

Composition of the Core Microbiome of Sparidae From the Mediterranean Sea
The core microbiome of the 12 Sparidae species studied here was typical of what is reported for marine fishes FIGURE 3 | Association between microbiome, gut morphological traits and diet items in Sparidae. (A) Constrained analysis on principal coordinates (CAP) analysis. Fish individuals are represented by colored dots, gut traits and diet items are represented by italicized labels and arrows, bacterial ASVs are represented by gray squares. Dissimilarity between microbiomes was estimated using abundance-weighted taxonomic dissimilarity (Hill diversity with q = 1). Only the variables identified as explaining a significant proportion of total microbiome variation are depicted (ANOVA.CCA, p-value < 0.05). (B) Number of taxa whose abundance are significantly correlated with variables that explained a significant proportion of microbiome variation. ASVs correspond to panel (A) while other ranks correspond to plots in Supplementary Figure 4. (C) The ASVs associated with explaining variables (white bars in panel B) were taxonomically assigned at the Family level and their proportion were estimated. (Sullam et al., 2012;Egerton et al., 2018), i.e., a dominance of Proteobacteria, Bacteroidetes, and Firmicutes, with these three phyla representing more than 84% of the sequences. Our results also confirm that the genus Vibrio, albeit dominant on the mucus, skin and gills of teleost (Chiarello et al., 2018;Ruiz-Rodríguez et al., 2020), represents only a small proportion of their gut microbiome. Regarding Sparidae' microbiome, the literature is highly skewed toward one species, Sparus aurata, which is used as a model in aquaculture (Dimitroglou et al., 2010;Estruch et al., 2015;Nikouli et al., 2019;Rosado et al., 2019), while microbiomes of most other Sparidae have not been described in the wild, with the exception of a recent study (Ruiz-Rodríguez et al., 2020). To our knowledge, our study provides the first report of the gut microbiome composition for five species: B. boops, D. puntazzo, Diplodus sargus, Lithognathus mormyrus, and Pagellus acarne. The structure of the gut microbiome significantly differed between the 12 studied species, with no genus and only four families (Rikenellaceae, Lachnospiraceae, Desulfovibrionaceae, and Vibrionaceae) representing more than 10% of the sequences in more than two host species. Species-specific gut microbiomes have been reported for other marine fishes both within and between families (Miyake et al., 2015;Ruiz-Rodríguez et al., 2020), but some studies also reported the absence of inter-specific difference because of high inter-individual variability in the gut microbiome (Nikouli et al., 2020).
It is worth mentioning that the method used here to identify the members of the core microbiome was based on the abundance-occurrence distribution of ASVs (Fillol et al., 2016;Jeanbille et al., 2016) and not on arbitrarily defined occurrence or abundance thresholds. This could explain the high number of ASVs identified as core (n = 842) compared with other studies (e.g., n = 15 in Sullam et al., 2015;0 < n < 110 in Chiarello et al., 2018;n = 7 in Nikouli et al., 2020). Indeed, in our method, even ASVs observed in a low number of species (low occurrence) but that were consistently abundant within individuals from these species were identified as core.

Weak Support for Phylosymbiosis in the Gut Microbiome of Sparidae Fishes
Overall, we found little support for the existence of gut phylosymbiosis in Sparidae fishes, i.e., we did not observe significant relationships between host phylogenetic relatedness and microbiome dissimilarity (using both Mantel tests and CAP analysis). Such result is in agreement with the observation of weak phylosymbiosis signal in skin microbiome of teleost fishes (Chiarello et al., 2018;Doane et al., 2020;Sylvain et al., 2020). In our case, the absence of phylosymbiosis, despite inter-species differences, is driven by the fact that the 10 carnivorous species are spread on the two main branches of the Sparidae phylogeny while they have similar gut traits and microbiome. For instance, P. erythrinus and Pagrus pagrus are phylogenetically distant but have a microbiome similar to those of most of the other Mediterranean Sparidae (i.e., the Diplodus cluster, the other Pagellus species, S. aurata and L. mormyrus). This microbiome convergence in carnivorous Sparidae resulted in species pairs with small microbiome dissimilarity despite high phylogenetic distance (i.e., points in the bottom right of the plots in Figure 2) that blurred the phylosymbiosis signal. On the contrary, the outlier regarding diet, gut traits and microbiome (S. sarpa) is not a phylogenetic outlier, which resulted in species pairs with small phylogenetic distance and strong microbiome dissimilarity (i.e., green points in Figure 2).

Determinants of the Gut Microbiome in Mediterranean Sparidae
Gut morphology and diet were the strongest determinants of the gut microbiome in Mediterranean Sparidae. These relationships were mostly driven by the fact that S. salpa, and to a lesser extent B. boops, were simultaneously microbiome and gut outliers. The importance of gut features in shaping the microbiome has been proposed to explain the observation of similar microbiome composition in wild and domesticated populations of zebrafish that experienced very different dietary and environmental conditions (Roeselers et al., 2011). Here, gut length and relative gut length were found to be conserved in host phylogeny, with one clade composed of S. salpa an B. boops having significantly longer guts than other Sparidae. While the relative gut length of these two species are 2.5 and 2.1 times higher than those of carnivorous species, respectively, their absolute gut length are 3.3 times longer for S. salpa (591 mm against an average of 177 ± 56 mm for the ten carnivores) and only 1.8 times longer for B. boops (329 mm) which is a much smaller species. These two species also appeared as outliers regarding their trophic ecology, with uncommon diet classification and a low diversity diet compared with the other studied species, along with the consumption of algae in S. salpa.
One limitation of our study lies in the use of stomach content data collected in other studies from the same region and recovered from FishBase. Individual analysis of stomach content was not feasible as the studied species are mostly diurnal and the fishes were caught during the night so their stomach was empty. Another approach to get the identity of ingested preys at the individual level corresponds to the metabarcoding of the gut content and feces (Soininen et al., 2009;Kartzinel et al., 2015), but this approach does not provide quantitative estimates of the ingested preys, contrarily with gut content analyses (Parravicini et al., 2020).
Sarpa salpa and B. boops, with their unique diet specialization (plant grazer and plankton feeder, respectively) and their extremely long guts, host bacterial taxa that were scarce in the other ten species (classified as macrofauna hunters). These taxa belong to three families of the Clostridiales order in particular (Lachnospiraceae and Ruminococcaceae), which are obligate anaerobes well-known for their associations with several herbivorous hosts (Llewellyn et al., 2014;Moran et al., 2019) including fishes (Clements et al., 2009;Sullam et al., 2012;Jones et al., 2018). In addition, S. salpa and B. boops traits and microbiome were associated with two genera known for their anaerobic lifestyle, Akkermansia and Alistipes, and more generally with their respective phyla, Verrucomicrobia and Bacteroidetes, also known to be abundant in the gut of herbivorous fishes (Sullam et al., 2012;Miyake et al., 2015;Ruiz-Rodríguez et al., 2020).
Herbivory in general is associated with modifications of the gut to accommodate the microbes involved in the digestion of plant-based material (Moran et al., 2019), either with differentiated parts (e.g., rumen, hindgut) or more simply through an elongation of the digestive tract that favor more anaerobic conditions within the distal part of the gut. Longer guts also increase food transit time, a necessary step to allow the fermentation process to take place and extract nutrients from plant material (Clements and Choat, 1995). In these regards, S. salpa is truly unique among the studied Sparidae as its gut is twice longer than the one of its closest, and yet distant, relative B. boops (divergence ∼18.4 Mya), allowing a long transit time (∼18 h, Goldenberg and Erzini, 2014). Furthermore, contrarily with B. boops that only consumes a limited proportion of plants/algae, S. salpa is specialized on these resources and consequently provides the raw material required for fermentation of plant polysaccharides to SCFAs (butyrate, acetate) by gut anaerobes (Boutard et al., 2014). As a result, Clostridiales represented 56% of the microbiome in S. salpa and only 30% in B. boops, while Lachnospiraceae represented 40 and 14% of the total, respectively. Lachnospiraceae is the most abundant taxa in the cow rumen where it represents more than a third of the sequences (Seshadri et al., 2018) and participates to plant fermentation as they are particularly efficient in cellulose degradation (Schwarz, 2001;Biddle et al., 2013).
The multifaceted ecological uniqueness exhibited by S. salpa has been reported in other clades such as the giant and red pandas, that are bamboo-eating specialists within a group of carnivores. Similarly, recent studies revealed that their typical straight, short and non-differentiated carnivore guts limit their ability to host the typical diversified microbiome of herbivores (Guo et al., 2020), contrarily with the bamboo lemur (Hapalemur griseus), a bamboo specialist hosting a much more diversified microbiome (McKenney et al., 2018b) thanks to its omnivorous digestive tract and a very long transit time (>30 h).

CONCLUSION
In this study, we highlight the presence of a dietary, anatomical and microbiological outlier within a carnivorous family of teleost fishes, the Sparidae. We revealed that the diet originality of the strict herbivore S. salpa is supported by a longer gut which increases the food transit time as well as anaerobic conditions favoring bacteria clades able to degrade plant materials. A key challenge is now to assess whether S. salpa' microbiome is comparable with the one of typically herbivorous fish clades such as Kyphosidae, Acanthuridae or Siganidae, both in terms of composition and efficiency to help the host assimilate nutrients. This comparison is particularly urgent to address because two Lessepsian rabbitfishes are spreading in the Eastern part of the Mediterranean basin and led in many places to overgrazing of macrophytes (Vergés et al., 2014).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Data Dryad, with the following doi: 10.5061/dryad.kprr4xh35.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because the fish individuals used for this study were already dead when collected from local fishermen.

AUTHOR CONTRIBUTIONS
AE, J-CA, AA, RS, and SV collected the samples. RS, AG, and LB provided fish traits data. AA and AE performed the molecular lab work. AE analyzed the data, performed the statistical analyses, and wrote the manuscript with the help of J-CA and SV. SV, AE, and J-CA conceived the study. All authors gave final approval for publication and agreed to be held accountable for the work performed therein.