Effects of aging on the skin and gill microbiota of farmed seabass and seabream

Background Important changes in microbial composition related to sexual maturation have been already reported in the gut of several vertebrates including mammals, amphibians and fish. Such changes in fish 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 microbiota 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 microbiota in shaping the fish skin and gill microbiota. Results Microbial diversity, composition and predicted metabolic functions varied across fish maturity stages. Alpha-diversity in the seabass microbiota varied significantly between age groups and was higher in older fish. Conversely, in the seabream, no significant differences were found in alpha-diversity between age groups. Microbial structure varied significantly across age groups; moreover, high structural variation was also observed within groups. Different bacterial metabolic pathways were predicted to be enriched in the microbiota of both species. Finally, we found that the water microbiota was significantly distinct from the fish microbiota across all the studied age groups, although a high percentage of ASVs was shared with the skin and gill microbiotas. Conclusions We report important microbial differences in composition and potential functionality across 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 fish coping mechanisms during stressful conditions. Our results provide new evidence suggesting that growth and sexual maturation have an important role in shaping the microbiota of the fish external mucosae and highlight the importance of considering different life stages in microbiota studies. Supplementary Information The online version contains supplementary material available at 10.1186/s42523-020-00072-2.


Background
Research on animal microbial communities (microbiota) is growing exponentially as the link between microbiota and host health is strongly validated by emerging evidence [1][2][3][4][5][6][7][8]. Age-related uctuations in the microbiota 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 microbiota 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.
Importantly, most studies testing the role of age on sh microbiota 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 microbiota to environmental changes and the high interindividual microbiota 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 microbiota. 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] microbiota of juvenile and mature adult sh from several species, showed a general pattern of differentiation between life-stages with differences attributed to intraspeci c niche partitioning [14,35]. Additionally, increases in body weight were seen to be associated with an increase in the microbial structure (i.e. beta-diversity) of the skin and gill microbiota of wild rabbit sh [36].
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 the skin and gill bacterial microbiota 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 in their microbiota. Additionally, we investigated the impact of the microbial communities present in the water column in the skin and gill microbiota.

Results
Skin, gill and water microbial samples from the different age groups of both species were collected simultaneously (same day and time) from separate ponds (Additional le 1). Three age cohorts were sampled for the seabass, which included sh on their 1 st , 2 nd and 3 rd year of age, and two for the seabream, which included sh in their 2 nd and 3 rd 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 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 age groups.

Microbial diversity across age groups
Alpha-diversity. Microbial alpha-diversity was calculated using Shannon, Faith's phylogenetic diversity (PD), ACE and Simpson 1-D indices. In general, the skin microbiota showed higher mean values for the alphadiversity indices than the gill microbiota across all age groups in both sh species, except in the Simpson 1-D index of the late juveniles and adults of the seabass (Additional le 2). In seabass, the skin and gill microbiota of late juveniles and mature adult sh presented higher mean values of alpha-diversity than the microbiota of the early juveniles ( Figure 1A, Additional le 3). In seabream, skin and gill microbiota showed similar alpha-diversity in both cohorts ( Figure 1B Table 1). Moreover, high inter-individual variability within age groups was also observed. However, differences in beta-diversity dispersion within age groups for the skin and gill microbiotas of both species were small (Bray-Curtis distance, Figure 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) microbiotas of all studied age groups (Additional le 4). 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 (Additional le 4). The most abundant microbial phyla and genera found in both sh species varied between age groups and tissues ( Figure 3, Additional le 4). LME models showed that the relative abundance of all those phyla was signi cantly different between age groups, except in the gill microbiota of the seabream, where the relative abundance of Cyanobacteria did not vary (Additional le 5). LME analyses also revealed that 100% and 63% of the genera varied between age groups 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 5). 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 5).
Microbial predicted functional diversity across age groups About 462±18 KEGG pathways were inferred in the skin and gill microbiota of the seabass, while 455±4 pathways were inferred in the skin and gill microbiota of the seabream. Linear discriminant analysis of the metagenomic predictions performed in LEfSe showed that different pathways were signi cantly enriched for each age group in both species (Figure 4, Additional le 6). Overall, there were more enriched pathways in the older age groups (seabass late juveniles and mature adults of both species) than in the younger age groups of both species (Figure 4, Additional le 6).
While there were no signi cantly enriched pathways in the skin microbiota of early juvenile seabass, enriched pathways in the gill microbiota of this age group were related to metabolic regulator biosynthesis, purine nucleotide degradation, sugar degradation and fermentation of pyruvate. In the skin microbiota of late juveniles seabass, enriched pathways were related to thiamine biosynthesis, aldehyde degradation and L-arabinose degradation; while in the gill microbiota they 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 gill microbiota was also enriched by pathways related to the biosynthesis of chlorophyll, folate, hemiterpene, L-alanine, Ltyrosine, NAD, secondary metabolite and ubiquinol, chloroaromatic compound degradation, fermentation to lactate and glycolysis ( Figure 4, Additional le 6).
In the skin microbiota 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 microbiota 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 microbiota; whereas pyrimidine and purine deoxyribonucleotide de novo biosynthesis, autotrophic CO 2 xation and fermentation of pyruvate were enriched in the gill microbiota ( Figure 4, Additional le 6).

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

Discussion
We characterized the skin and gill microbiota 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 indicate that sh age, in particular sexual maturation and growth, in uence skin and gill microbial diversity (Table 1; Figure 1), composition (Additional le 4; Figure 3, 5) and predicted functions (Additional le 6; Figure 4).

Microbial diversity across age groups
Fish growth and sexual maturation are usually accompanied by extreme morphological and physiological changes [e.g. 44,45]. Importantly, some of these 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 microbiota [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 showed that increases in body weight are accompanied by increased 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 microbial diversity, composition and predicted functionality observed in the present study.
The skin and gill microbiota of older age groups of seabass showed signi cantly higher alpha-diversity than early juveniles. Additionally, a higher percentage of signi cant differences in the relative abundance of the most abundant phyla and genera occurred between early juveniles and older age groups (67±27% and 55±38% in the skin and gill, respectively; Additional le 5). This suggests that the skin and gill microbiota of the seabass were highly dynamic, diversifying with age. Conversely, the skin and gill microbiotas of seabream juveniles and adults showed similar alpha-diversity means, even though a high percentage of the most abundant phyla and genera varied between age groups (70±42% and 63±18% in the skin and gill, respectively; Additional le 5). Variation in microbiota alpha-diversity between different age groups has been previously reported for many sh species. For example, studies on the zebra sh and salmon gut microbiota, 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]. Moreover, the relative abundance of bacterial groups that are predominant in the microbiota changed with aging in other sh species as well [e.g. 14, 17, 20].
The differences found in the present study in microbial 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]), being particularly evident in longitudinal studies encompassing several months [13,17,20]. High inter-individual variability within age groups was also previously reported for other sh species [e.g. 14, 28,30]. However, our results showed that sh age only explained a low percentage of the variation in the bacterial community structure (R2 < 0.1). This suggests that microbial differences between age groups are small at the community level, but clearly noticeable at the species level, with a high proportion of the predominant bacterial taxa (61±39% and 46±32% in the skin and gill of seabass, respectively; 70±42% and 63±18% in the skin and gill of the seabream, respectively) changing their abundance with sexual maturation. Although our statistical models accounted for sampling dates as a random factor, other biotic and abiotic factors (e.g., variation in the environment and individual weights) could be responsible for most of the variation observed in community structure [e.g. 12,19,28].
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 microbiotas of both sh species across age groups. Although metabolic information is particularly limited for sh microbiotas, studies on other vertebrates, mainly in humans and their gut microbiotas, are starting to shed light on the bene cial outcomes speci c microbial metabolic functions have on the host health and physiology [e.g. 50].
Notwithstanding that present results should be interpreted with caution since PICRUSt2 results are limited by the currently available genomes and biased towards human health microorganisms [51], one could suggest that some of the enriched metabolic pathways found in our analysis could also improve the seabass and seabream health and physiology. The protective role of the microbiota is often related to the production of secondary metabolites that provide chemical defense and mediate bacterial diversity [4].
Secondary metabolites with antimicrobial activity have been previously isolated from microbial species inhabiting the gut microbiome of sh [52]. Here, the biosynthesis of secondary metabolites that have been associated with antimicrobial activity, including hemiterpene [e.g. 53], were enriched in mature adults of both species. Additionally, the biosynthesis of chlorophyll and several amino acids, herein enriched in older age groups of both sh species, have also been found to be expressed in the skin and gut of healthy humans [e.g. 54; 55]. Importantly, amino acid biosynthesis was reported in the gut microbiota of grass carp when fed a protein-de cient diet, suggesting a metabolic role of the gut microbiota towards sh nutrition [56]. The biosynthesis of vitamins, here enriched in older age groups of seabass, has been found bene cial for human skin (e.g Saxena et al., 2018) and gut mucosa, including folate and thiamine [57]. In addition, polyamines are bacterial metabolites known to have several bene ts towards gut mucosa recovery [58].
Here, these pathways were enriched in the juveniles of seabream.
It is also 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. 59]. Biosynthesis of fatty acids and unsaturated fatty acids were two of the predicted metabolic pathways enriched in the microbiota of mature seabream. These same pathways have also been enriched in previous analyses of the skin and gut microbiota in the atlantic salmon [60,61] and in the skin microbiota of the common snook [62] 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 chloroaromatic compounds by bacteria is essential to remove them from the environment and to prevent absorption through the skin and gills in aquatic animals [63][64][65].
Following alpha-diversity patterns, sh from older age groups, particularly in the seabass, had greater enrichment of predicted functions related mainly to the biosynthesis and degradation of compounds; as well as, to a lesser extent, metabolism and energy cycles. We hypothesize that the increase in microbial diversity observed as sh ages leads to wider functional diversity. This could prove bene cial to those shes, given the key physiological modi cations older sh groups are experiencing during sexual maturation and growth.

Fish and water microbiota comparisons
The water microbiota of shponds were signi cantly distinct and more diverse than the skin and gill microbiota 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,[66][67][68], gills [14,36], gut [7,15,18,21,69], 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 microbiota [17,21], others have also shown that water microbiota do not in uence directly the microbiota of sh mucosa [7, 8, 13-15, 18, 19, 22, 28, 30, 34, 36, 66-68, 70, 71]. Importantly, a previous study of the skin microbiota of the seabass and seabream also showed signi cant differences with plaktonic communities [67]. However, in that study only a low number of OTUs (3%) was shared between skin and water microbiota, 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.
Microbial dissimilarities depicted by PCoAs showed that, although signi cantly different, the skin microbiota of both species clustered more closely to the water microbiota than the gill microbiota. 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 microbiota (p<0.03), but not the skin microbiota. 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 microbiota.

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 microbiota 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, we demonstrate that, to some extent, changes occurring later in life can also be correlated with aging factors such as growth and sexual maturation. We uncovered important differences in the diversity, composition and predicted function of the skin and gill microbiota across age groups of farmed seabass and seabream.
Fish included in this study were exposed to a heterogeneous estuarine environment, that varies seasonally as well as daily (e.g., salinity and temperature uctuations); hence, although we observed signi cant differences in microbial community structure due to age, other unmeasured biotic and abiotic factors may have more deeply structured the sh microbiota. 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 attributed to age found herein. Overall, our results were in line with what has been previously found in the skin [35,36] and gill [14,36] microbiotas 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 mucosal microbiota, especially in studies focusing on farmed sh, where the microbiota and disease dynamics can be very important.

Material And Methods
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 to reach minimum commercial size, 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 275g, whereas for seabream maturity is usually attained at 300g. We monitored the skin and gill microbiota of seabass and seabream of different age cohorts, including juveniles and adults. Due to sampling restrictions within the sh farm, sampling was strictly noninvasive 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.
We collected samples 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 from this species were categorized as early juveniles (9 months and an average weight of 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 from an intermediate age-cohort 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; Additional le 1). In this sh farm, all ponds share the same water in ow, which is taken from a single point in the estuary. Water in each pond is naturally recycled at each high tide (twice a day) and never shared between ponds. Hence, sh share roughly the same water quality and environment. Additionally, fry were bought from commercial hatcheries where genetic variation is usually limited [72].
Fish were caught from each tank using a sh line, and skin and gill samples were non-invasively taken using sterile swabs (Medical Wire & Equipment, UK). We swabbed 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. Afterwards sh were released unharmed. We collected water samples (1 L) 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.
We sampled ve sh per week per age group, totaling 60 individuals per species and age group. We processed 360 seabass samples (60 skin and 60 gills x 3 age groups) plus 29 water samples from their corresponding shponds. We also processed 240 seabream samples (60 skin and 60 gills x 2 age groups) plus 16 water samples from their corresponding shponds. We processed seabass and corresponding water samples using the PowerSoil DNA Isolation Kit (QIAGEN, Netherlands), while seabream and 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., microbiota is not compared between sh species). We measured DNA concentration 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 [74]. We estimated a midpoint rooted tree of ASVs using the Quantitative Insights Into Microbial Ecology 2 package (QIIME2; release 2019.7). We constructed a table containing amplicon sequence variants (ASVs) and made taxonomic inferences against the SILVA (138 release) reference database [75]. We normalized ASV abundances using the negative binomial distribution [76], 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 Simpson 1-D indices as implemented in the R package phyloseq [77]. We assessed variation in microbial composition (alpha-diversity) and the mean proportions of the most abundant taxa (≥5% of all sequences) using Linear Mixed Effects models (LME) with the lmer R package [78]. 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 [79]. We used the strata argument to permutate sampling dates and ran 1,000 permutations.
Previous sh studies of skin and gill microbiota [e.g. 7,8,36,69], including seabass and seabream [80], have shown remarkable differences in microbial composition and structure across host species and tissues. Although the two species studied here are farmed in the same location, provided feeds were different throughout the sampling period. For this reason, and because sh samples were processed using two different DNA extraction kits, we did not compare microbiotas between shes [e.g. 21,81]. Additionally, a previous study by our group [82] showed that disease and antibiotic treatment in seabass leads to asymmetrical shifts in skin and gill microbial communities. Therefore, we carried out all our statistical analyses separately for each sh species and tissue.
Microbial potential metabolic functions were predicted using the metagenomic Phylogenetic Investigation of Communities by Reconstruction of Unobserved States software (PICRUSt2) embedded in QIIME2 [83], 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 [84]. We identi ed differentially abundant metabolic pathways in the skin and gill microbiota of seabass and seabream across age groups using linear discriminant analysis (LDA) in LEfSe, using age groups as classes [85]. As suggested by the authors, we used a P-value cut-off of 0.05 and a LDA effect size cut-off of 2 [85].
Finally, to assess to what extent water microbial communities shaped skin and gill microbiota across sh age groups, we estimated the number of shared ASVs between sh and water microbiota and constructed Venn diagrams in R. We used PERMANOVA and mantel testes [86] to assess differences in community structure and correlations between tissues and water microbiota, respectively, in both species. 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. 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 are denoted with "ns".   Table 1: Mean alpha-diversity values, and alpha-and beta-diversity comparisons for the skin and gill microbiota of the different age groups of seabass Dicentrarchus labrax and seabream Sparus aurata (n=60 per species x age group x tissue). Variation in alpha-diversity was assessed using Linear Mixed Effect models, with age groups as a fixed 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 significance (P value) and for each PERMANOVA test (betadiversity) we report the R2 statistics and significance (P value). Significant differences are indicated in bold. EJ: early juveniles; LJ: late juveniles; MA: mature adults; J: juveniles.