Characterization of the skin and gill microbiomes in farmed seabass and seabream across different age groups

Background Important changes in microbiome composition related to sexual maturation have been already reported in the gut of several vertebrates including mammals, amphibians and sh. Such changes in sh are linked to reproduction and growth during developmental stages, diet transitions and critical life events. We used amplicon (16S rRNA) high-throughput sequencing to characterize the skin and gill bacterial microbiomes of farmed seabass and seabream belonging to three different developmental age groups: early and late juveniles and mature adults. We also assessed the impact of the surrounding estuarine water microbiome in shaping the sh skin and gill microbiomes. Results Microbiome diversity, composition and potential metabolic functions varied across sh maturity stages. Alpha-diversity in the seabass microbiome varied signicantly between age groups and was higher in older sh. Conversely, in the seabream, no signicant differences were found in alpha diversity between age groups, although it was higher in the skin of juveniles. Microbiome structure varied signicantly across age groups. Different bacterial metabolic pathways were predicted to be enriched in the microbiomes of both species. Finally, we found that the water microbiome is signicantly distinct from all the sh microbiomes across the studied age groups, although a high percentage of ASVs is shared with the skin and gill microbiomes. Conclusions We report important microbial differences in composition and potential functionality across the different ages of farmed seabass and seabream. These differences may be related to somatic growth and the onset of sexual maturation. Importantly, some of the inferred metabolic pathways could enhance the host coping mechanisms during stressful conditions. Our results provide new evidence suggesting that growth and sexual maturation have an important role in shaping the external mucosa microbiomes of sh and highlight the importance of considering different life stages in microbiome studies. and seabream, respectively). We, thus, conclude that growth and sexual maturation are likely the main drivers of the differences found herein. Overall, our results were in line with what has been previously found in the skin [35, 36] and gill [14, 36] microbiomes of several wild reef sh, suggesting this could be a general pattern across sh. Our results also highlight the importance of considering sexual maturation as a key factor shaping external sh mucosa microbiomes, especially in studies focusing on farmed sh, where the microbiome and disease dynamics can be very important.


Abstract
Background Important changes in microbiome composition related to sexual maturation have been already reported in the gut of several vertebrates including mammals, amphibians and sh. Such changes in sh are linked to reproduction and growth during developmental stages, diet transitions and critical life events. We used amplicon (16S rRNA) high-throughput sequencing to characterize the skin and gill bacterial microbiomes of farmed seabass and seabream belonging to three different developmental age groups: early and late juveniles and mature adults. We also assessed the impact of the surrounding estuarine water microbiome in shaping the sh skin and gill microbiomes.
Results Microbiome diversity, composition and potential metabolic functions varied across sh maturity stages. Alpha-diversity in the seabass microbiome varied signi cantly between age groups and was higher in older sh. Conversely, in the seabream, no signi cant differences were found in alpha diversity between age groups, although it was higher in the skin of juveniles. Microbiome structure varied signi cantly across age groups. Different bacterial metabolic pathways were predicted to be enriched in the microbiomes of both species.
Finally, we found that the water microbiome is signi cantly distinct from all the sh microbiomes across the studied age groups, although a high percentage of ASVs is shared with the skin and gill microbiomes.
Conclusions We report important microbial differences in composition and potential functionality across the different ages of farmed seabass and seabream. These differences may be related to somatic growth and the onset of sexual maturation. Importantly, some of the inferred metabolic pathways could enhance the host coping mechanisms during stressful conditions. Our results provide new evidence suggesting that growth and sexual maturation have an important role in shaping the external mucosa microbiomes of sh and highlight the importance of considering different life stages in microbiome studies.

Background
Research on animal microbial communities (microbiomes) is growing exponentially as the link between microbiome and host health is strongly validated by emerging evidence [1][2][3][4][5][6][7][8]. Age-related uctuations in microbiomes are well studied in humans and are considered as "natural, inevitable and benign" [9]. Critical microbial changes occur during infancy and old age, coinciding with stages when the immune system is also more fragile [9]. Results linking changes in the gut microbiome to reproduction and growth (e.g. monkeys, [10]) or disease resistance in early life stages (e.g. amphibians, [11]) have been also found in other vertebrates.
Ecological factors, such as diet transitions [e.g. 21] or critical life events (e.g. habitat transition, [20]), which in turn are intrinsically linked to sexual maturation, also play a major role in shaping the sh gut microbiome. Importantly, most studies testing the role of age on sh microbiomes were cross-sectional and based on a single time point or short time window [e.g. 16,18,19,21,22]. Thus, given the high susceptibility of the sh microbiome to environmental changes and the high inter-individual microbiome variability [e.g. [28][29][30], the compound effect of all these factors could be hard to interpret [31].
Fish skin and gills and their associated mucous and microbes form a natural physical and chemical barrier to pathogens [4,32,33]. Despite this protective role, little is known about potential host developmental effects on skin and gill microbiomes. Filling out this knowledge gap is especially important in sh farming, where diseases are a main concern causing large mortality rates [e.g. 34]. Two previous studies in wild reef sh comparing the gill [14] and skin [35] microbiomes of juvenile and mature adult sh from several species, showed a general pattern of differentiation between life-stages with the differences being attributed to intraspeci c niche partitioning [14,35]. Additionally, increases in body weight were seen to be associated with an increase in the microbiome structure (i.e. beta-diversity) of the skin and gill microbiomes of wild rabbit sh [36].
The European seabass (Dicentrarchus labrax) and the gilthead seabream (Sparus aurata) are two of the most important farmed sh in Europe (global production of 191,003 and 185,980 tns, respectively, in 2016, [37]). The gilthead seabream is a protandric hermaphrodite, maturing rst as males between years 1 and 2, with sex reversal occurring in the following 2-3 years [38][39][40]. The European seabass reaches sexual maturity between years 2 and 3 in males, and after year 3 in females [41][42][43]. Typically, in semi-extensive production systems, both sh are reared until they reach their rst commercial size (18-24 months). However, demand for larger sh sizes has been increasing, meaning that both species can reach sexual maturation before harvest.
Here we used amplicon (16S rRNA) high-throughput sequencing to characterize over six months the skin and gill bacterial microbiomes of farmed seabass and seabream from different ages (juvenile stages and mature adults). Our main aim was to describe differences in composition, structure and potential metabolic functions. Additionally, we investigated the impact of the microbial communities present in the water column in the skin and gill microbiomes.

Results
Skin, gill and water microbial samples from the different age groups of both species were collected simultaneously (same day) from separate ponds. Three age cohorts were sampled for the seabass, which included sh on their 1st, 2nd and 3rd year of age, and two for the seabream, which included sh in their 2nd and 3rd year of age. Due to non-invasive sampling, we have coupled available information from the literature [38][39][40][41][42][43] with data provided by the sh farmers about the weight and age of maturation of both species, to classify samples into age groups.
The three seabass age cohorts were then classi ed as early juveniles, late juveniles and mature adults, respectively; while the two seabream age cohorts were classi ed as juveniles and mature adults, respectively -see the Materials and Methods section for more details. Differences in the average weight estimated for each age group at the beginning and end of our sampling indicate a 245% growth for the seabass early juvenile group, an 83% growth for the late juvenile group and a 43% growth for the mature adult group. For the seabream a 143% growth was estimated for juveniles and a 16% growth for the mature adults. Descriptive analyses were performed for each age group separately and comparative statistical analyses were performed between groups.
Microbiome diversity across age groups Alpha-diversity. Microbial alpha-diversity was calculated using Shannon, Faith's phylogenetic diversity (PD), ACE and Fisher indices. The skin microbiome showed higher alpha-diversity than the gill microbiome across all age groups in both sh species (Additional le 1). In seabass, the skin and gill microbiomes of late juveniles and mature adults sh presented higher alpha-diversity than the microbiomes of the early juveniles ( Fig. 1A, Additional le 2). In seabream, the skin microbiome of juveniles presented higher alpha-diversity than the microbiome of mature sh, while the gill microbiome showed similar diversity in both cohorts (Fig. 1B, Additional le 2). Linear Mixed Effects (LME) model analysis (diversity ~ age group + (1|sampling date) showed most alpha-diversity estimates varied signi cantly between seabass age groups in both tissues. Pairwise comparisons between age groups in seabass showed signi cant differences in alpha-diversity for almost all of the early vs late juvenile and early juvenile vs mature adult comparisons (p < 0.05, Table 1), while late juvenile vs mature adult comparisons were never signi cant (p > 0.05, Table 1) in both tissues. In the seabream only the Shannon and PD indices of the gill microbiomes varied signi cantly between juveniles and mature adults (p < 0.04, Table 1). Table 1 Mean alpha-diversity values, and alpha-and beta-diversity comparisons for the skin and gill microbiomes of the different age groups of seabass Dicentrarchus labrax and seabream Sparus aurata. Variation in alpha-diversity was assessed using Linear Mixed Effect models, with age groups as a xed factor and sampling time as a random factor. Differences in beta-diversity were assessed using PERMANOVA. For each linear model effect test (alpha-diversity) we report the F statistic and signi cance (P value) and for each PERMANOVA test (beta-diversity) we report the R2 statistics and signi cance (P value Beta-diversity. Microbial structure was estimated using phylogenetic UniFrac (unweighted and weighted) and Bray-Curtis distances. The PERMANOVA analyses of dissimilarities (diversity ~ age group, strata = sampling date) showed signi cant differences between the age groups of both species (p < 0.02, Table 1), except for the UniFrac Weighted distance between the gills of early and late seabass juveniles (p = 0.1, Table 1), seabass late juveniles and mature adults (p = 0.2, Table 1), and the skin of juveniles and seabream adults (p = 0.3, Table 1). Principal Coordinate Analyses (PCoAs) were used to visualize microbial structure (dissimilarity) and depicted the differences between early and late juvenile/mature seabass groups, in both tissues (Bray-Curtis distance, Fig. 2). For the seabream, however, differences between age groups were not evident (Fig. 2).
Bacterial taxa. Proteobacteria and Bacteroidetes were the most abundant (≥ 5%) phyla in the skin (averaging 41 ± 4% and 39 ± 2% of the sequences in seabass and 55 ± 4% and 31 ± 4% in seabream) and gill (averaging 52 ± 7% and 25 ± 5% in seabass and 69 ± 4% and 12 ± 1% in seabream) microbiomes of all studied age groups ( Table 2). The NS3a marine group and a genus belonging to the Flavobacteriacea family were the most abundant (≥ 5%) genera in the skin (10 ± 1 and 11 ± 2, respectively) and gill (6 ± 1 both) of all the age groups in seabass; while Burkholderia-Caballeronia-Paraburkholderia was the most abundant genus in the skin (17 ± 1) and gill (25 ± 0) of both age groups in seabream ( Table 2). The most abundant microbial phyla and genera found in both sh species varied between age groups and tissues (Fig. 3, Table 2). LME models showed that the relative abundance of all those phyla was signi cantly different between age groups, except in the gill microbiome of the seabream, where the relative abundance of Cyanobacteria did not vary (Additional le 3). LME analyses also revealed that 100% and 63% of the genera varied in the skin and gill of the seabass, respectively, while 40% and 50% varied in the skin and gill of the seabream, respectively (Additional le 3). Pairwise comparisons of taxa across age groups in seabass yielded a higher percentage of signi cant differences between early juveniles and mature adults in both tissues (100% in the skin and 38% in the gill) than between early and late juveniles (67% in the skin and 13% in the gill), or between late juveniles and mature adults (0% in the skin and 25% in the gill) (Additional le 3). Table 2 Relative mean proportions (%) of the most abundant phyla and genera (≥ 5%) in the skin and gill microbiomes of the different age groups of the seabass Dicentrarchus labrax and the seabream Sparus aurata, and in the water column. Taxa with a ≥ 5% relative mean proportion in a group are indicated in bold. Unknown genera are identi ed as u.g.

Seabass Seabream
Skin Microbial predicted functional diversity across age groups About 462 ± 18 KEGG pathways were inferred in the skin and gill microbiomes of the seabass, while 455 ± 4 pathways were inferred in the skin and gill microbiomes of the seabream. Linear discriminant analysis of the metagenomic predictions performed in LEfSe, showed that different pathways were signi cantly enriched for each species and for each age group in both species (Fig. 4, Additional le 4).
While there were no signi cantly enriched pathways in the skin of early juvenile seabass, enriched pathways in the gills of this age group were related to metabolic regulator biosynthesis, purine nucleotide degradation, sugar degradation and fermentation of pyruvate. In the skin of late juveniles seabass enriched pathways were related to thiamine biosynthesis, aldehyde degradation and L-arabinose degradation; while in the gills were related to denitri cation, galactose degradation and nitrogen compound metabolism. In mature seabass, pyrimidine and purine deoxyribonucleotide de novo biosynthesis were enriched in both tissues. Additionally, the gills were also enriched by pathways related to the biosynthesis of chlorophyll, folate, hemiterpene, L-alanine, L-tyrosine, NAD, secondary metabolite and ubiquinol, chloroaromatic compound degradation, fermentation to lactate and glycolysis (Fig. 4, Additional le 4).
In the skin of seabream juveniles, enriched pathways were related to amine and polyamine biosynthesis and degradation, choline biosynthesis, and sugar acid and toluene degradation; whereas in the gill only pyrimidine and purine deoxyribonucleotide de novo biosynthesis were identi ed. The enriched pathways of the seabream mature adults were related to fatty acid, L-methionine, NAD, palmitate, palmitoleate, siderophore, stearate and unsaturated fatty acid biosynthesis, pyrimidine and purine nucleotide salvage, aspartate superpathway and TCA cycle in the skin; whereas pyrimidine and purine deoxyribonucleotide de novo biosynthesis, autotrophic CO2 xation and fermentation of pyruvate were enriched in the gill (Fig. 4, Additional le 4).

Fish and water microbiome comparisons
The microbiome of shpond water showed higher alpha-diversity than the skin and gill microbiomes of seabass and seabream, except when compared to the Shannon index estimated for the seabass late juveniles (Additional le 1). The analyses of dissimilarities between the skin and gill microbiomes and the water microbiome were statistically signi cant for all pairwise comparisons (PERMANOVA, p < 0.001, Table 3). Moreover, results from Mantel tests revealed a correlation between gill and water microbiomes of seabass and seabream across age groups (p < 0.03, Table 3), except in the case of late juvenile seabass (p > 0.05, Table 3). PCoAs showed that the water microbiome clustered more closely to the skin microbiome than to the gill microbiomes in both shes (Additional le 5). In both species, the percentage of ASVs shared between skin and water microbiomes, and between gill and water microbiomes was very similar (14%±1 and 15%±1 of ASVs (amplicon sequence variants), respectively) (Fig. 5).

Discussion
We characterized the skin and gill microbiomes of different age groups of farmed European seabass and gilthead seabream using 16S rRNA amplicon high-throughput sequencing. By taking into account potential environmental and seasonal effects, the results of the present study show that sh age in uences skin and gill microbiome diversity and structure (  Fig. 3, 5) and predicted functions (Additional le 4; Fig. 4).

Microbiome diversity across age groups
Fish growth and sexual maturation is usually accompanied by extreme morphological and physiological changes [e.g. 44,45]. Importantly, some of the changes reported for the skin and gills have been suggested to also affect their microbiota. For example, changes in epidermal structure derived from sexual maturation (e.g. increases in the number, size and activity of the mucous cells) have been reported in several sh species [e.g. 44,46], and suggested to contribute to a higher infection with Saprolegnia fungus in the cases of the sea trout and brown trout [47]. Likewise, changes in the hormones expressed in the skin alter the biochemistry of the skin mucous and also potentially affect its microbiome [48]. Fish growth and sexual maturation also impact gill morphology and function in some sh species. For example, the ability to osmoregulate at different salinities was seen to increase throughout the developmental stages of the seabass (between larva and juvenile individuals, [45]. Additionally, body size was also identi ed as the main factor affecting morphological variation in gill rakes and the size of their pores in the Silver Carp and Gizzard Shad, suggesting that the overall ltering ability of these species is related to size and maturation [49]. Importantly, a recent study in rabbit sh showed that increases in body weight are accompanied by increases in the microbial community structure of the skin and gill of rabbit sh [36]. We thus hypothesize that such physiological and morphological changes occurring during sh growth have led to the changes in microbiome diversity, composition and predicted functionality observed in the present study.
The skin and gill microbiomes of older age groups of seabass showed signi cantly higher alpha-diversity than early juveniles. Although all of the most abundant phyla were maintained between age groups, the skin and gill microbiomes of the seabass were highly dynamic, diversifying with age. Conversely, the skin microbiome of seabream juveniles showed a tendency to exhibit higher alpha-diversity than adults, though these differences were not signi cant. Variation in microbiome alpha-diversity between different age groups has been previously reported for many sh species. For example, studies on the zebra sh and salmon gut microbiome, have reported differences between mature and immature life stages; however those differences also coincide with other major ecological changes in the sh, such as diet [17] or environment transitions [20].
The differences found in the present study in microbiome structure across age groups, which were consistently signi cant in both species, have been already reported in other sh (e.g., several reef sh [14,35]; Salmo salar [19]), mainly in longitudinal studies several months long [13,17,20].

Microbial predicted functional diversity across age groups
The predicted functional analysis suggests that distinct signi cantly enriched metabolic pathways are expressed in skin and gill microbiomes of both sh species across age groups. Following alpha-diversity patterns, the number and diversity of enriched pathways was higher in mature seabass adults when compared to juveniles, especially in the gill. In seabream, on the other hand, there were essentially no differences in both microbial diversity and number of enriched pathways between age groups. However, one must interpret these results with caution, since PICRUSt2 results are limited by the currently available genomes and biased towards human health microorganisms [50]. However, it is worth noticing that some of the enriched metabolic pathways detected in the present study could be driven by the high environmental variability of the Alvor estuary where these sh are reared. In estuaries, salinity variations occur on a daily basis due to tides and pollutants can be prevalent [e.g. 51]. Biosynthesis of fatty acids and unsaturated fatty acids were two of the predicted metabolic pathways enriched in the microbiome of mature seabream skin. These same pathways have also been enriched in previous analyses of the skin and gut microbiomes in the atlantic salmon [52,53] and in the skin microbiome of the common snook [54] when transitioning between freshwater and seawater. Additionally, two of the predicted metabolic pathways identi ed in both sh species were related to degradation of toxic compounds. Speci cally, biodegradation of the highly prevalent toxic pollutants toluene and chloroaromaric compounds by bacteria is essential to remove them from the environment and to prevent absorption through the skin and gills in aquatic animals [55][56][57].

Fish and water microbiome comparisons
The water microbiome of shponds were signi cantly distinct and more diverse than the skin and gill microbiomes of both sh, regardless of their age. It is known that free-living microbial communities retain higher richness than host-associated communities [31], with many studies showing a higher bacterial diversity in water relative to sh skin [28,30,36,[58][59][60], gills [14,36], gut [7,15,18,21,61], stomach [36], hindgut [36] and whole larvae [22]. Although some studies in sh have shown that the microbial communities found in the water tend to be recovered in the larval gut microbiome [17,21], others have also shown that water microbiomes do not in uence directly the microbiomes of sh mucosa [7, 8, 13-15, 18, 19, 22, 28, 30, 34, 36, 58-60, 62, 63]. Importantly, a previous study of the skin microbiome of the seabass and seabream [59] also showed signi cant differences with plaktonic communities. However, in that study only a low number of OTUs (3%) was shared between skin and water microbiomes, whereas in the present study higher percentages of ASVs were shared between the skin (14%±1) and the gill (15%±1) of both sh species and the surrounding water.
Microbiome dissimilarities depicted by PCoAs showed that, although signi cantly different, the skin microbiome of both species clustered more closely to the water microbiome than the gill microbiome. However, only a small percentage of the variation (PC 1 -average 18%±2; PC 2 -average 10%±1) was explained by this analysis. On the other hand, the results from the Mantel tests showed a correlation between the water and gill microbiomes (p < 0.03), but not the skin microbiomes. This suggests that although both skin and gill are permanently in contact with water, the gill environment may be more susceptible to variations in the water microbiome.

Conclusions
Skin and gill are important mucosal barriers that protect the sh from the external environment. They are in permanent contact with the water column and thus prone to pathogenic bacterioplankton colonization. However, most studies so far investigating microbiome changes related to sh age have either strictly focused on early life stages (i.e., larvae development) or on the gut microbiota. In the present study important differences were uncovered in the diversity, composition, and predicted function of the skin and gill microbiomes across age groups of farmed seabass and seabream. Besides the increments in biomass recorded at the end of our sampling and the onset of sexual maturity, the estimated growth rate of each cohort also changed. Growth rate decreased drastically with age, being much higher in juveniles (243% and 83% for early and late seabass juveniles, and 143% in seabream) relative to adults (43% and 16% in adult seabass and seabream, respectively). We, thus, conclude that growth and sexual maturation are likely the main drivers of the differences found herein. Overall, our results were in line with what has been previously found in the skin [35,36] and gill [14,36] microbiomes of several wild reef sh, suggesting this could be a general pattern across sh. Our results also highlight the importance of considering sexual maturation as a key factor shaping external sh mucosa microbiomes, especially in studies focusing on farmed sh, where the microbiome and disease dynamics can be very important.

Fish species, sampling and preparation
Fish were sampled at a semi-intensive open-water farm in the Alvor Estuary (Ria Formosa, Portimão, Portugal). In this sh farm, seabass and seabream production can take up to 36 months, so having a healthy mucosa during this time is of utter importance. The gilthead seabream is a protandric hermaphrodite, maturing rst as males between years 1 and 2 in the wild, with sex reversal occurring in the following 2-3 years [38][39][40]. The European seabass reaches sexual maturity between years 2 and 3 in males, and after year 3 in females [41][42][43]. In this particular sh farm, seabass typically reaches sexual maturity at approximately 275 g, whereas for seabream maturity is usually attained at 300 g. We monitored the skin and gill microbiomes of seabass and seabream of different age cohorts, including juveniles and adults. Due to sampling restrictions within the sh farm, sampling was strictly non-invasive and sh could not be dissected to con rm sexual maturation. The categorization of the age group cohorts was based on previous studies [e.g. 38,41] and the weight at maturity records available at this farm.
Samples were collected every other week (12 sampling time points) between August 2017 and January 2018 (6 months). We simultaneously sampled three seabass age groups cohorts with approximately one year old difference. Fish were categorized as early juveniles (9 months and an average weight of 22 g at the beginning of the study and 15 months and an average weight of 76 g at the last sampling date), late juveniles (18 months and an average weight of 151 g at the beginning of the study and 24 months and an average weight of 277 g at the last sampling date), and mature adults (32 months and an average weight of 467 g at the beginning of the study and 38 months and an average weight of 669 g at the last sampling date). We also simultaneously sampled two seabream cohorts categorized as juveniles (15 months and an average weight of 103 g initially and 21 months and an average weight of 250 g at the last sampling date), and mature adults (37 months and an average weight of 411 g at the beginning of the study and 37 months and an average weight of 476 g at the last sampling date). Seabream of an intermediate age were not available.
Each age group and species was reared in separated but not distant open-water ponds (maximum 344 m and 380 m apart for seabass and seabream, respectively). In this sh farm, all ponds shared the same in ow of estuarine, which circulates between ponds and is naturally recycled. Hence, sh share roughly the same water quality and environment. Additionally, sh of each species were bought from commercial hatcheries where genetic background is limited.
Fish were caught from each tank using a sh line, and gill and skin samples were non-invasively taken using sterile swabs (Medical Wire & Equipment, UK). The right laments between the rst and second arches of the gill and the right upper lateral part of the sh skin from head to tail were swabbed. Afterwards sh were released unharmed. Water samples (1 L) were collected from the ve different culture ponds at the same time as sh swabbing was performed, except during the month of December, when no water samples could be collected. Water samples were ltered through 0.2 µm lters on collection day. Swabs and lters were immediately frozen at -20ºC and then transported in dry ice to the CIBIO-InBIO laboratory where they were kept at -80ºC until processing.
Five sh were sampled per week per age group, totaling 60 individuals per species and age group. A total of 360 seabass samples (60 skin and 60 gills x 3 age groups) plus 29 water samples from their corresponding shponds and a total of 240 seabream samples (60 skin and 60 gills x 2 age groups) plus 16 water samples from their corresponding shponds were processed. The seabass and their corresponding water samples were processed using the PowerSoil DNA Isolation Kit (QIAGEN, Netherlands), while seabream and their corresponding water samples were processed using the PureLink Microbiome DNA Puri cation Kit (ThermoFisher Scienti c, UK). We used two different DNA extraction kits due to supply shortage at the time of extraction. This technical difference did not impact the goals of our study since we studied each sh species separately (i.e., microbiomes are not compared between sh species). DNA concentration and quality were measured in a NanoDropTM 2000 Spectrophotometer (ThermoFisher Scienti c, USA). DNA extractions were shipped on dry ice to the University of Michigan Medical School (USA) for ampli cation and sequencing according to the protocol of Kozich et al. [64]. Each sample was ampli ed for the V4 hyper-variable region of the 16S rRNA gene (~ 250 bp). All amplicon libraries were pooled and sequenced in a single run of the Illumina MiSeq sequencing platform.
Approximately 8,313,608 and 6,943,265 16S rRNA sequences were retrieved for seabass and seabream, respectively. The number of sequences per sample ranged from 726 to 46,001 in seabass and from 5,145 to 151,713 in seabream. After normalization and removal of non-bacterial reads, 8,724 and 5,754 ASVs were assigned to the skin and gill, respectively, of seabass; while 5,308 and 3,423 ASVs were assigned to the skin and gill, respectively, of seabream. A total of 2,543 ASVs were retrieved from the water samples collected in seabass shponds, while 1,440 ASVs were retrieved from the waters of seabream shponds. Taxa showing a mean relative proportion ≥ 5% in any group were considered the most abundant in that group.

Data processing and statistical analysis
Raw FASTQ les were denoised using the DADA2 pipeline in R with the paramenters for ltering and trimming being trimLeft = 20, truncLen = c(220,200), maxN = 0, maxEE = c(2,2), truncQ = 2 [65]. A midpoint rooted tree of ASVs was estimated using the Quantitative Insights Into Microbial Ecology 2 package (QIIME2; release 2019.7). A table containing amplicon sequence variants (ASVs) was constructed and taxonomic inferences made against the SILVA (138 release) reference database [66]. ASV abundances were normalized using the negative binomial distribution [67], which accounts for library size differences and biological variability.
Microbial taxonomic alpha-diversity (intra-sample) was calculated using Shannon, Faith's phylogenetic diversity (PD), ACE and Fisher indices as implemented in the R package phyloseq [68]. Variation in microbial composition (alpha-diversity) and the mean proportions of the most abundant taxa (≥ 5% of all reads) were assessed using Linear Mixed Effects models (LME) with the lmer R package [69]. Since we were interested in assessing whether microbial diversity varied across sh age groups (predictor), we used age groups as a xed factor and sampling date (with 12 sampling time points) as a random factor. The nal general LME formula was expressed as: microbial diversity ~ sh age group + (1|sampling time point). Microbial structure (beta-diversity) was estimated using phylogenetic Unifrac (unweighted and weighted) and Bray-Curtis distances. Dissimilarity in microbial structure between samples was visualized using principal coordinates analysis (PCoA). Additionally, differences in community structure driven by sh age group were further tested using permutational multivariate analysis of variance (PERMANOVA) as implemented in the adonis function of the vegan R package [70]. We used the strata argument to permutate sampling dates and ran 1,000 permutations.
Previous sh studies of skin and gill microbiomes [e.g. 7,8,36,61], including seabass and seabream [71], have shown remarkable differences in microbial composition and structure across host species and tissues. Additionally, a previous study by our group [72] showed that disease and antibiotic treatment in seabass leads to asymmetrical shifts in skin and gill microbial communities. Therefore, all our statistical analyses were carried out separately for each sh species and tissue.
To assess to what extent water microbial communities shaped skin and gill microbiomes across sh age groups, we estimated the number of shared ASVs between sh and water microbiomes and constructed Venn diagrams in R. PERMANOVA and mantel testes [73] were used to assess differences in community structure and correlations between tissues and water microbiomes, respectively, in both species.
Finally, microbial potential metabolic functions were predicted using the metagenomic Phylogenetic Investigation of Communities by Reconstruction of Unobserved States software (PICRUSt2) embedded in QIIME2 [74], applying a weighted nearest sequenced taxon index (NSTI) cutoff of 0.03. Predicted metagenomes were collapsed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway metadata [75]. Differentially abundant metabolic pathways in the skin and gill microbiomes of seabass and seabream across age groups were identi ed using linear discriminant analysis (LDA) in LEfSe, using age groups as classes [76]. As suggested by the authors, we used a Pvalue cut-off of 0.05 and a LDA effect size cut-off of 2 [76].

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and material
The datasets generated and/or analyzed during the current study are available in the NCBI Sequence Read Archive (SRA) database within the BioProject ID XXXXX (will be added upon acceptance).

Competing interests
The authors declare that they have no competing interests. were performed with a P value and LDA score cut-offs of 0.05 and of 2, respectively. Additional le 5: PCoA plot computed using Bray-Curtis distances for water, skin and gills microbiomes of the seabass Dicentrarchus labrax (A) and the seabream Sparus aurata (B). Each dot represents a microbiome sample and is coloured by tissue/origin (skin, gill and water). Figure 1 Mean values and standard deviations of Shannon alpha-diversity estimates plotted for the early juveniles/juveniles (green), late juveniles (yellow) and mature adults (orange) of the seabass Dicentrarchus labrax (A) and the seabream Sparus aurata (B). Pairwise comparisons of alpha-diversity were assessed using Linear Mixed Effect models with age groups as a xed factor and sampling time as a random factor. Statistically signi cant differences are denoted with an asterisk and non statistically signi cant differences with "ns".