Effects of disease, antibiotic treatment and recovery trajectory on the microbiome of farmed seabass (Dicentrarchus labrax)

The mucosal surfaces of fish harbour microbial communities that can act as the first-line of defense against pathogens. Infectious diseases are one of the main constraints to aquaculture growth leading to huge economic losses. Despite their negative impacts on microbial diversity and overall fish health, antibiotics are still the method of choice to treat many such diseases. Here, we use 16 rRNA V4 metataxonomics to study over a 6 week period the dynamics of the gill and skin microbiomes of farmed seabass before, during and after a natural disease outbreak and subsequent antibiotic treatment with oxytetracycline. Photobacterium damselae was identified as the most probable causative agent of disease. Both infection and antibiotic treatment caused significant, although asymmetrical, changes in the microbiome composition of the gills and skin. The most dramatic changes in microbial taxonomic abundance occurred between healthy and diseased fish. Disease led to a decrease in the bacterial core diversity in the skin, whereas in the gills there was both an increase and a shift in core diversity. Oxytetracycline caused a decrease in core diversity in the gill and an increase in the skin. Severe loss of core diversity in fish mucosae demonstrates the disruptive impact of disease and antibiotic treatment on the microbial communities of healthy fish.


Results
Approximately, a total of 3.6 million raw reads were generated, while the number of sequences per sample ranged from 2,354 to 50,564 (Supplementary Table S1). A total of 6,485 unique ASVs were detected, but after normalization and depletion of Archaea and Algae ASVs, a total of 3,827 ASVs (1,560,279 sequences) and 3,741 ASVs (1,904,115 sequences) were analyzed for the gill and skin microbiomes, respectively (Supplementary Table S1). Taxa showing a mean proportion ≥4% in any state were considered the most abundant taxa. Analyses of alphaand beta-diversity showed no significant differences (P > 0.05; data not shown) between samples from the two healthy time points as well as between the samples from the three recovery time points (Fig. 1) for both gill and skin microbiomes. Therefore, in all our subsequent analyses, samples from Aug 21 (Healthy 1) and Aug 29 (Healthy 2) were combined into the "healthy" state; and samples from Sep 19 (Recovery 1), Sep 26 (Recovery 2) and Oct 3 (Recovery 3) were combined into the "recovery" state in order to increase sample size (Fig. 1).
Gill bacterial composition and diversity. No significant differences were detected in alpha-diversity across all states (RRPP, P > 0.5), with the exception of the Shannon index (RRPP, P = 0.04) ( Table 1, Fig. 2A). There were, however, significant differences between healthy and recovery states for all alpha-diversity indices (RRPP, P ≤ 0.03; Table 1). Beta-diversity also varied greatly between states (PCoA, Fig. 3A), with significant differences both in overall and pairwise comparisons in almost all the tests (Adonis, P < 0.05; Table 1).
Bacteroidetes, Proteobacteria and Verrucomicrobia were the most abundant phyla retrieved from the seabass gill microbiome across states, accounting for 83% to 93% of the sequences altogether ( Table 2). The most abundant genera were NS3a marine group, Polaribacter 4, Pseudomonas and Rubritalea, which were present in all four states ( Table 2, Supplementary Fig. S1A). Other relatively abundant genera were: (a) Polynucleobacter, which accounted for 4-7% of the sequences in the diseased, treatment and recovery states, but only 0.2% in the healthy state; (b) Stenotrophomonas, represented by 5% of the sequences in the healthy state, and 2-3% in the remainder states; and (c) Photobacterium, which accounted for 5% of the sequences in the diseased state but ≤1 in all other states ( Table 2, Supplementary Fig. S1A).
Taxa mean proportions varied between states: (i) in healthy versus diseased states 7 taxa increased and 3 decreased; (ii) in diseased versus treatment states 3 taxa increased, 2 decreased and 6 remained almost constant; and (iii) finally, in treatment versus recovery states 3 taxa decreased, 6 increased and 1 remained constant ( Table 2, Fig. 4A). The 3 most abundant bacterial phyla and the 8 most abundant genera all varied significantly (P ≤ 0.04) in their mean proportions across the four studied states (Fig. 4A, Table 1). In addition, 9 of these taxa varied significantly (P ≤ 0.04) between healthy and disease states, whereas 5 varied significantly (P ≤ 0.04) between treatment and recovery states (Table 1). Only 2 genera varied significantly (P ≤ 0.03) between disease and treatment states (Table 1). Of the 55 ASVs recovered from the gill core microbiome across the four states, 21 were present in the healthy state, 26 in the diseased state, 5 in the treatment state and 10 in the recovery state (Fig. 5A). Four of these ASVs were unique to the healthy state, 11 unique to the diseased state and one to the recovery state (Fig. 5A). Of the 11 unique ASVs recovered from the gill core microbiome of diseased fish, one was identified as Photobacterium damselae. There were 8 other ASVs belonging to the Photobacterium genus, of which 7 were unique to the diseased state and one was found in all four states (Supplementary Table S1). This suggests that P. damselae is the most likely causative agent of the disease in the diseased fishes.
Skin bacterial composition and diversity. Alpha-diversity estimates varied significantly across all states (RRPP, P ≤ 0.002; Table 1, Fig. 2B) and between states in the skin microbiome. They decreased significantly between healthy and diseased fish and increased significantly between diseased and treatment states (RRPP, P ≤ 0.003; Table 1, Fig. 2B). Beta-diversity estimates show significant differences across all states and between states (Adonis, P < 0.05; Table 1).
As for the skin microbiome, Bacteroidetes, Proteobacteria and Verrucomicrobia were the most abundant phyla retrieved across states, accounting for 87% to 93% of the sequences altogether ( Table 2). The genera NS3a marine group, Polaribacter 4, Pseudomonas and Stenotrophomonas were the most abundant in all four states ( Table 2, Supplementary Fig. S1). Moreover, Pseudoalteromonas accounted for 5% of the sequences in the treatment state, but only 0.1-1% in the remaining states (Table 2, Supplementary Fig. S1).     Table 1). All taxa varied significantly between healthy and diseased states (P ≤ 0.01); all except one varied significantly between diseased and treatment states (P ≤ 0.02); and 2 genera varied significantly (P ≤ 0.02) between treatment and recovery states (Table 1).
A total of 43 ASVs formed the core microbiome of all four states, 17 were present in the healthy state, 8 were present in the diseased state, 33 in the treatment state and 8 in the recovery state (Fig. 5B). It is worth noticing that 2 ASVs were unique to the healthy state and 16 ASVs were unique to the treatment state.

Discussion
In this study, we investigated the dynamics of the gill and skin microbiomes in 140 samples of the farmed seabass Dicentrarchus labrax during a natural disease outbreak and subsequent antibiotic treatment with oxytetracycline. We used high-throughput sequencing technology to generate 16S rRNA bacterial ASVs and examine changes in microbial composition and diversity over six weeks. We identified Photobacterium damselae as the most probable causative agent of disease.
The most abundant taxa found in the gill and skin microbiomes of healthy farmed seabass (Dicentrarchus labrax) belonged to the Proteobacteria, Bacteroidetes and Verrucomicrobia phyla. These phyla have been previously described as the most abundant in the gill and skin microbiomes of several teleosts 14,50-52 , including the seabass 16,53,54 . At the genus level, the most abundant taxa were the NS3a marine group, Polaribacter 4, Pseudomonas, and Stenotrophomonas in the gills and skin (Table 2, Supplementary Fig. S1, Fig. 4), and Rubritalea in the gills ( Table 2, Supplementary Fig. S1, Fig. 4). These results are mainly in accordance with previously described microbiomes of healthy seabass 16,54 , including fish retrieved from the same farmed population during winter months 16 . However, one of the most abundant genera in the healthy seabass gill microbiome was Polynucleobacter 16 , which  www.nature.com/scientificreports www.nature.com/scientificreports/ in the present study only accounted for 0.2% of the sequences in the healthy state, and 4-7% in the other three studied states (Table 2). Another compositional difference was the high abundance of Stenotrophomonas found in both tissues in apparently healthy individuals (Table 2) in this study, but not in Rosado et al. 16 . Several environmental factors known to impact microbiome composition, such as seasonality 55,56 and water temperature 57 , could be driving these differences between our two studies.
The composition and diversity of the gill and skin seabass microbiomes varied differently during infection. Whereas in the skin there was a significant decrease in alpha-diversity between healthy and diseased fish, there were no significant differences in the gill microbiome. An overall decrease in microbial richness was also reported for the skin of Atlantic salmon as a result of infection with salmonid alphavirus 6 and sea lice 15 ; but interestingly, as in the present study, Legrand et al. 14 reported significant differences in microbial richness between the skin of healthy and enteritis-infected yellowtail kingfish, but not in the gills.
Significant changes in beta-diversity occurred in both gills and skin, showing clear signs of dysbiosis in both tissues. In the skin microbiome of diseased fish, the abundance of taxa from the non-pathogenic NS3a marine group and Polaribacter 4 decreased, whereas the pathogenic Pseudomonas and Stenotrophomonas significantly increased. Pseudomonas spp. almost doubled their abundance and largely dominated the skin microbiome of diseased fish. While the genus Stenotrophomonas contains important globally emergent fish pathogens (e.g. Stenotrophomonas maltophilia 58,59 ), Pseudomonas harbors both opportunist fish pathogens (e.g. P. baetica, P. chlororaphis 60,61 ; amongst others 8,62 ) and taxa with known antimicrobial activity against fish pathogens (e.g. Flavobacterium psychrophilum 63 ). For example, P. fluorescens is an important pathogen of carp and salmon 64,65 , but is also known to inhibit the growth of Saprolegnia, an oomycete that causes huge losses in aquaculture 66,67 . Importantly, a ten-fold increase of Pseudoalteromonas, which was not amongst the most abundant taxa in healthy fish, occurred in the skin of diseased fish. Species from this genus can inhibit the growth of both Vibrio spp. and Photobacterium damselae 6,68-71 , hence an increase of Pseudoalteromonas could lead to a decrease of the other two genera, as we have seen in the skin microbiomes of seabass transitioning from healthy to diseased states (from 2% to 0.7% and from 0.3% to 0.2%, respectively). In the gills of diseased fish, the majority of the most abundant bacterial genera in the healthy state (NS3a marine group, Polaribacter 4, Pseudomonas, Rubritalea and Stenotrophomonas) decreased significantly in abundance during infection, with the exception of Polynucleobacter. Amongst the most abundant taxa in the gill, only Photobacterium spp. was exclusively associated with diseased fish, where it showed a 25-fold increase. Similarly, all studies addressing the effects of parasitic infection on fish microbiomes reported significant changes in microbial composition 6,14,15 . Importantly, all of these studies reported an increase of potentially pathogenic taxa, which highlights the opportunistic nature of such pathogens 6,14,15 . Although Photobacterium damselae was only highly abundant in the diseased gill microbiome, the  www.nature.com/scientificreports www.nature.com/scientificreports/ tissue with more significant shifts in overall bacterial composition (alpha-diversity) between healthy and diseased states was the skin. This is not totally unexpected, since it has been shown that this pathogen can unequally affect the microbiome of distinctive mucosal surfaces such as the skin and gill 14 .
The effects of the disease in the core microbiomes were also significant and again different between tissues, with a shift of core species in the gill and a decrease of core diversity in the skin from healthy to diseased states. A shift of the microbial assemblages with enrichment of specific groups was also described for the gill microbiome of the yellowtail killifish as a result of enteritis 14 .
Antibiotics administration can negatively impact host physiology in different ways (e.g., inhibiting mitochondrial gene expression 72 ; decreasing enzymatic activity 46 ), leading to dysbiosis and the emergence of antibiotic resistant bacteria 35,[42][43][44]73 . Specifically, the reported effects of oxytetracycline in the gut microbiome of the Atlantic salmon showed a clear reduction in taxonomic diversity, becoming almost exclusively composed of the oxytetracycline resistant Aeromonas spp., which include the salmon pathogens Aeromonas sobria and A. salmonicida 74 . Similarly, in zebrafish, long-term exposure (6 weeks) to environmental concentrations of oxytetracycline, prompted both a decrease in gut microbial diversity and higher mortality when fish were challenged with the pathogen A. hydrophila 46 . The impact of broad-spectrum antibiotics in the skin microbiome of Gambusia affinis have also been assessed 41,42 . In this case, the use of rifampicin led to a decrease of diversity in the skin microbiome after 2.6 days of antibiotic administration. Additionally, as reported for zebrafish and Atlantic salmon, fish subjected to rifampicin antibiotic administration were more susceptible to infection due to osmotic stress and exhibited less growth compared to the control group, an effect that lasted one month after treatment 41,42 . A key difference with the present study is that the fish used by Carlson et al. 41,42 were healthy before antibiotic administration. Importantly, our results showed that skin core diversity was higher in healthy than in recovery individuals, indicating a negative effect of disease and antibiotic use.
In the present study, administration of oxytetracycline resulted in a dramatic reduction of Photobacterium abundance in the gill microbiome, with this genus no longer being one of the most abundant taxa in the treatment and recovery states. This was expected given the reported sensitivity of P. damselae to several antibiotics, including oxytetracycline 38 . Pseudoalteromonas, however, remained one of the most abundant taxa in the skin microbiome during treatment perhaps due to the host innate immune response mediated by the skin microbiome, given the ability of this genus to produce antimicrobial metabolites that are correlated with host homeostasis 68 .
Previous studies on the impact of antibiotics on fish skin microbiomes showed that, even though stabilization of bacterial communities during recovery occurs, neither diversity nor composition returns to healthy-like values in the short term (after 1 week) 41,42 . Here the relative frequency of the most abundant taxa found in the skin microbiome of the seabass during the recovery period, which corresponded to 3 weeks, was similar to that in healthy individuals (P ≥ 0.1 for all taxa except the NS3a marine group). In the gill microbiome, however, www.nature.com/scientificreports www.nature.com/scientificreports/ differences in taxa proportions between the healthy and recovery states were significant for almost all of the most abundant taxa. Hence, although dysbiosis due to infection was more noticeable in the skin than in the gill, the microbial communities present in the skin seem to be more resilient than those of the gill. Importantly, although the abundance of Photobacterium damselae in the gill seemed to have been controlled through antibiotic administration, it increased significantly in the recovery state, surpassing its initial proportion in the healthy state.
In summary, the mucosal surfaces of fish, such as the gill and the skin, are constantly exposed to several pathogens in the aquatic environment and are crucial to prevent and/or control disease 1 . It has been shown that both infectious diseases and antibiotic treatment lead to a decrease in microbial diversity, which translates into a decrease in host immunity 14,42 . Here we described microbial changes in the gill and skin of adult seabass in response to a natural disease outbreak followed by a succeeding treatment with oxytetracycline. We showed that the gill and skin microbiomes are highly disturbed by both infection and antibiotic treatment, ultimately decreasing their diversity.

Methods
Ethical statement. This study monitored a natural infection and subsequent antibiotic treatment as part of routine procedures in a commercial fish farm. All animals were handled by the fish farm employees, our sampling through swabbing was non-invasive and fish were released unharmed with no mortalities observed. According to the Portuguese legislation DL N° 113/2013, our work does not involve animal experimentation and therefore is exempted from the need of ethical approval.
Experimental design, sample collection and preparation. Ten individuals of seabass were collected once a week between August 21 and October 3, 2016, from the same rearing tank in a commercial fish farm located in the estuarine environment of the Ria Formosa (Portimão), southern Portugal. Fish were hatched at September 26, 2014 and entered the growth facility at March 6, 2015. Fish were kept in an open water circulation system in a semi-intensive farming facility, where water is supplied to each tank from the estuary. Fish were kept at a density of ca. 3 kg/m 3 corresponding to roughly 100 fish/tank with fish weighting on average 281 g. Given that it was not possible to tag individual fish and the unlikelihood of re-sampling the same individuals every week, a subset of samples believed to be representative of the population was chosen, i.e. 10 individuals (~ 10%), and for statistical purposes individuals were considered as pseudo-replicates. All fish were fed with the same commercial feed and they shared the same clinical history. Individuals were randomly caught using a fishing pole and skin and gill swabs were collected immediately using tubed sterile dry swabs (Medical Wire & Equipment, UK). Skin samples were taken by swabbing several times along the right upper lateral part of the fish from head to tail, while gill swabs were taken from the right filaments between the first and second arch. Due to the non-invasive nature of our sampling procedure, it was not possible to ascertain the sex of the individuals sampled; however, we do not expect this to impact our conclusions since, to the best of our knowledge, no gender bias in microbiome composition has ever been reported for skin or gill of piscine hosts. Swabs were immediately stored at −20 °C until www.nature.com/scientificreports www.nature.com/scientificreports/ transported on dry ice to the CIBIO-InBIO laboratory by airmail where they were kept at −80 °C until further processing.
To assess gill and skin microbiome dynamics in the seabass during a disease outbreak and under oxytetracycline treatment after infection, fish were sampled in 4 different states: healthy, diseased, treatment and recovery (Fig. 1). During the healthy state (August 21 and 29), all fish specimens were considered healthy due to a lack of visible disease symptoms, such as external lesions or behavioural alterations. On September 8 fish began to die in the farming tanks, showing symptoms of disease, and treatment with oxytetracycline antibiotic (a broad-spectrum tetracycline) was initiated, being administrated at 35 g/Kg through commercial feed for at lasted 8 days. On the same day, smears from spleen and kidney were collected for culture using Bionor kits DE020, MONO-VA-50 for Vibrio anguillarum and DL020, MONO-Pp-50 for Photobacterium piscicida. Agglutination essays were not conclusive and, at that stage, the causative agent of the disease was unknown. We do not have samples from September 8, hence we used the samples from our closest time point, September 5, which we classified as potentially diseased (i.e., diseased state). Antibiotic treatment lasted until September 16 and fish were sampled on September 12; this sample point corresponded to the treatment state. Then, three additional time points were sampled (September 19 and 26, and October 3), when fish were no longer dying or presented signs of infection; these three time points corresponded to the recovery state.
Total DNA from 140 fish samples (70 skin and 70 gills) was extracted using the PowerSoil DNA Isolation Kit (QIAGEN, Netherlands), following the manufacturer's protocol. DNA extractions were shipped in dry ice to the University of Michigan Medical School (USA) for amplification and sequencing on a single run of the Illumina MiSeq platform according to the protocol of Kozich et al. 75 Each sample was amplified for the V4 (~250 bp) hypervariable region of the 16S rRNA gene, using the primers in Caporaso et al. 76 This region has been widely used to characterize microbiomes from vertebrates (Earth Microbiome Project 77 ), including fish 42,78-80 . Data and statistical analysis. Raw FASTQ files were analyzed using the Quantitative Insights Into Microbial Ecology 2 (QIIME2; release 2018.4) platform. Clean sequences were aligned against the SILVA (132 release) reference database 81 using the DADA2 pipeline 82 . A feature table containing amplicon sequence variants (ASVs) was constructed and normalized using the negative binomial distribution 83 . The core microbiome was assessed at the ASV level for the gill and skin of seabass for each state (healthy, diseased, treatment and recovery) separately. An ASV was considered as part of the core microbiome if present in 100% of the samples from each state. Core diversity is here defined as number of ASVs represented in a given group.
Microbial alpha-diversity (intra-sample) was calculated using Shannon, ACE, Fisher and Faith's phylogenetic diversity (PD) indices as implemented in the R package phyloseq 84 . Microbial beta-diversity (inter-sample) was estimated using phylogenetic Unifrac (unweighted and weighted) and Bray-Curtis distances. Dissimilarity between samples was assessed by principal coordinates analysis (PCoA). Variation in microbial alpha-diversity and the mean proportions of the most abundant taxa (with more than 4% of all reads) were assessed using linear models with randomized residuals in a permutation procedure (RRPP). Differences in community composition (beta-diversity) were tested using permutational multivariate analysis of variance (PERMANOVA) with 1,000 permutations as implemented in the adonis function of the R vegan package 85 . All statistical analyses were carried out separately for the gills and skin. All statistical analyses were performed in R-studio v1.0.143 86 .

Data availability
The raw sequences are available at NCBI Sequence Read Archive (SRA) database within the BioProject ID PRJNA575053.