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 fish 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 first as males between years 1 and 2 in the wild, with sex reversal occurring in the following 2–3 years [38–40]. The European seabass reaches sexual maturity between years 2 and 3 in males, and after year 3 in females [41–43]. In this particular fish 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 fish farm, sampling was strictly non-invasive and fish could not be dissected to confirm 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 fish farm, all ponds shared the same inflow of estuarine, which circulates between ponds and is naturally recycled. Hence, fish share roughly the same water quality and environment. Additionally, fish of each species were bought from commercial hatcheries where genetic background is limited.
Fish were caught from each tank using a fish line, and gill and skin samples were non-invasively taken using sterile swabs (Medical Wire & Equipment, UK). The right filaments between the first and second arches of the gill and the right upper lateral part of the fish skin from head to tail were swabbed. Afterwards fish were released unharmed. Water samples (1 L) were collected from the five different culture ponds at the same time as fish swabbing was performed, except during the month of December, when no water samples could be collected. Water samples were filtered through 0.2 µm filters on collection day. Swabs and filters 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 fish 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 fishponds and a total of 240 seabream samples (60 skin and 60 gills x 2 age groups) plus 16 water samples from their corresponding fishponds 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 Purification Kit (ThermoFisher Scientific, 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 fish species separately (i.e., microbiomes are not compared between fish species). DNA concentration and quality were measured in a NanoDropTM 2000 Spectrophotometer (ThermoFisher Scientific, USA). DNA extractions were shipped on dry ice to the University of Michigan Medical School (USA) for amplification and sequencing according to the protocol of Kozich et al. [64]. Each sample was amplified 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 fishponds, while 1,440 ASVs were retrieved from the waters of seabream fishponds. 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 files were denoised using the DADA2 pipeline in R with the paramenters for filtering 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 fish age groups (predictor), we used age groups as a fixed factor and sampling date (with 12 sampling time points) as a random factor. The final general LME formula was expressed as: microbial diversity ~ fish 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 fish 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 fish 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 fish species and tissue.
To assess to what extent water microbial communities shaped skin and gill microbiomes across fish age groups, we estimated the number of shared ASVs between fish 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 identified using linear discriminant analysis (LDA) in LEfSe, using age groups as classes [76]. As suggested by the authors, we used a P-value cut-off of 0.05 and a LDA effect size cut-off of 2 [76].