Temporal changes in skin and gill microbiomes of Atlantic salmon in a recirculating aquaculture system – Why do they matter?

Mucosal surfaces are key components of teleost health, providing defence against opportunistic pathogens and other insults. Maintaining the integrity of mucosal surfaces and their associated microbial communities, especially the gill and skin that have large surface areas exposed to the environment is essential. Production of Atlantic salmon in land-based recirculating aquaculture systems (RAS) has increased significantly in recent years as it allows greater control over stability of the environment in which fish are reared, reduces water demand and minimises environmental impacts. However, little is known about the impact of the RAS environment upon the temporal dynamics of skin and gill mucosal microbiomes. In this study we examined microbial communities in gill mucus, skin mucus and rearing water throughout freshwater (FW) RAS production, and at 1-week and 4-weeks following transfer to seawater (SW) in open cage production using 16S rRNA sequencing. Microbial diversity and richness in skin and gill mucus of fish reared in a RAS system were temporally dynamic. Dynamics in richness and diversity were similar in the two mucosal tissues, and to some extent also mirrored that of the surrounding water. Dysbiosis indicated by an abrupt decline in diversity during FW production coincided with an increase in the relative abundance of two taxa belonging to the RAS-biofilter-associated nitrogen-cycling genus Hydrogenophaga in RAS tank water and this was also observed in gill and skin mucus. Extensive overlap in core taxa was observed between gill and skin mucus, but host-specific cores were non-existent during the dysbiotic event with all cores present in the rearing water. Diversity remained stable during the transition from FW to SW, but distinct community composition and core taxa were observed in the two environments. Although RAS are closely controlled, significant temporal variation could be observed in temperature as well as levels of CO 2 and nitrogen compounds, reflecting the increasing biological load within the system over time. The results presented here suggest that, in terms of microbiomes, dysbiosis may occur in both the RAS environment and fish mucosal surfaces over time, but microbial communities have the capability to recover.


Introduction
As the global population continues to rise, as does the need to produce sufficient quality protein to sustain healthy diets.Aquaculture production outputs have risen in line with this need, however such upscaling of production raises environmental concerns.Such concerns have prompted a drive for more sustainable and efficient systems, diversifying away from traditional open water net-based production to land-based recirculating aquaculture systems (RAS) which greatly reduce water usage and minimise risks of escape events and eutrophication (d'Orbcastel et al., 2009;Attramadal et al., 2014).The use of RAS systems for smolt production is now practiced in many countries involved in commercial Atlantic salmon (Salmo salar) production.Atlantic salmon are anadromous, spending periods of the life cycle in both fresh-and saltwater.The preparatory transition from a freshwater adapated juvenile (parr) to a saltwater adapted smolt, called the parrsmolt transformation, is a key stage in the life cycle of Atlantic salmon both in the wild and in production, and is accompanied by many physiological, morphological and behavioural changes (Björnsson et al., 2011;McCormick, 2012).However, little is known about how production in such tightly controlled systems impacts on the development and stability of the mucosal microbiome, which is an important component of the salmonid immune system.
In teleost fish, the skin, gill, gut and nares have a mucosal layer which contains a wide variety of host derived antimicrobial compounds and immunological components acting as the first line of defence against infectious pathogens (Salinas, 2015;Salinas et al., 2021).Coupled with the mucosa there is a complex and continually adapting community of symbionts and opportunistic microorganisms which are integral in maintaining mucosal health and has constant interaction with the host immune system (Kelly and Salinas, 2017).The threat of infectious disease is constant in the aquatic environment.Maintaining homeostasis in these microbial communities via competitive exclusion (Coyte et al., 2015) and the parallel monitoring of microbial composition by both innate and adaptive immune systems is key in preventing invasion by opportunistic pathogens and influencing community composition.In addition, fish are in constant contact with external environments which are rarely stable.Microbiome composition is sensitive to a number of environmental factors including temperature (Rosado et al., 2021), salinity (Lokesh and Kiron, 2016;Dehler et al., 2017), rearing system (Minich et al., 2020a), stress (Boutin et al., 2013;Minniti et al., 2017), as well as factors associated with host biology including developmental stage (Lokesh et al., 2019;Heys et al., 2020) and body site (Uren Webster et al., 2018;Minich et al., 2020a).
Despite the importance of the external mucosa, studies relating to skin and particuarly gill microbiomes in salmonids are few in number.A recent review of fish skin microbiomes (Gomez and Primm, 2021) identified just 6 papers related to salmonids analysing inter-individual variation (Boutin et al., 2014), tissue comparisons (Lowrey et al., 2015), fresh-and salt-water communities (Lokesh and Kiron, 2016;Hamilton et al., 2019), response to handling stress (Minniti et al., 2017) and relative contribution of environment and genetics (Uren Webster et al., 2018).Additionally, effect of culture system and mucosal proteome interactions (Minniti et al., 2019) have been analysed.A number of these aforementioned skin studies also included gill (Lowrey et al., 2015;Minich et al., 2020a) with additional studies in gill addressing the influence of fishmeal-free diets (Schmidt et al., 2016), comparing Flavobacterium psychrophilum susceptible and resistant lines (Brown et al., 2019) and the effect of ploidy and salmonid alphavirus infection (Brown et al., 2021).
Many of the previous studies have focussed on a single time point, which would assume a stable microbial community.Year-long temporal analyses in pond-reared European seabass (Rosado et al., 2021) and wild marine Pacific chub mackerel (Minich et al., 2020b) identified highly dynamic external microbiomes impacted strongly by water temperature.In this study we analyse the temporal dynamics of the gill and skin microbiomes in Atlantic salmon in a commerical RAS, considered a stable system, and post-seawater transfer (SWT).To our knowledge no studies to date have examined the temporal dynamics of the microbial communities in external mucosal layers of Atlantic salmon during FW production.Results showed mucosal microbiomes were temporally dynamic and that these changes mirrored those of the environmental communities in the water.A potentially dysbiotic event was observed during freswater RAS production with a dramatic increase in the denitrifying genus Hydrogenophaga, constituing up to 70% of the total relative abundance.This suggested fluctuations in nitrogen metabolism in the system, reflected in both host and water microbial communities.However, microbial communities recovered and no adverse impacts were observed on fish health or growth, highlighting redundancy or resilience in the mucosal microbiomes.

Fish maintenance and sampling schedule
Mixed sex juvenile Atlantic salmon were followed from parr to smolt stage in a single stream of a commercial recirculating aquaculture system and sampled at multiple timepoints as previously described (Lorgen-Ritchie et al., 2021).The sampling schedule is summarised here utilising degree days (dd), cumulative temperature over a number of days, a measure used to determine smolt windows where dd = 0 represents the onset of artificial spring photoperiod.Fish were first sampled at − 802 dd at the parr stage (FW1-14.05.2019).Fish were exposed to a 'winter' photoperiod (12 L:12D) for 714 dd (6-7 weeks) and sampled at − 276 dd (FW2-18.06.2019), just prior to the onset of spring photoperiod at − 48 dd (FW3-02.07.2019), and as smolts just prior to SWT (261 dd) following application of a "spring" signal in the form of constant light (24L:0D) ( FW4-18.07.2019FW4-18.07. to 25.07.2019)).Fish were transferred 334 DD at the end of July to a sea water cage site.At each FW sampling point, six fish were sampled from triplicate tanks (n = 3 tanks, 18 fish in total) with the exception of FW1 (n = 2 tanks, 12 fish in total).Fish (n = 2 cages, 12 fish in total) were sampled approximately 1 week and 4 weeks post-SWT.The one-week sampling point was chosen to represent the early phase of SW transfer and the four-week as a later phase where host physiology has had a chance to adapt to the new environment.Similar timepoints have been utilised in other studies investigating the impact of SW transfer (Lokesh and Kiron, 2016;Dehler et al., 2017).Sampling was only carried out at two SW timepoints to minimise interruption to the farm operation at a critical point in the production cycle.All dd represent mean values across replicate tanks.
Prior to sampling, fish were culled by lethal anaesthesia (MS-222, 1000 ppm, PHARMAQ, Norway) and blood sampled using heparinised syringes.Blood was centrifuged (2500RPM at 4 • C for 15 min) to separate the plasma which was stored at − 20 • C for later chloride analyses.Weight and length of all fish was recorded.Gill (4-6 gill filaments) were collected using fine scissors (primary gill arch) in SEI buffer (150 mM sucrose, 10 mM EDTA, 50 mM imidazole, pH 7.3) and flash frozen on liquid nitrogen for later Na + , K + − ATPase (NKA) enzyme activity analysis.
Gill mucus samples were taken by gently rolling a cotton swab across all gill arches on the left side of each fish and adding to a 2 mL sample collection tube containing 1.5 mL RNAlater™ (Ambion Inc., USA) (Clinton et al., 2021).Skin mucus samples were taken by gently rolling a cotton swab along the lateral line of the left side of each fish and again adding to a 2 mL sample collection tube containing 1.5 mL RNAlater™ (Ambion Inc., USA).Samples were stored at 4 • C for 24 h followed by longer term storage at − 80 • C. As previously described (Lorgen-Ritchie et al., 2021), 4 × 50 mL water samples of tank or cage water were collected at each sample point.Water samples were transported at room temperature and stored at -20 • C prior to filtration through 0.2 μM Whatman Cyclopore polycarbonate membrane filters (Sigma-Aldrich; WHA70634702) using a vacuum pump.Filters were stored at − 80 • C until extraction.
Plasma chloride (Cl − ) levels (resting levels, no seawater challenge) were determined by colorimetric titration according to manufacturer's protocol using silver nitrate as the reagent (926S Chloride Analyser, Sherwood Scientific, Cambridge).
Na + , K + − ATPase activity was determined with a kinetic assay run in 96-well microplates at 26 • C and read at a wavelength of 340 nm for 10 min according to the method of (McCormick, 1993).Protein concentrations were determined thereafter using a BCA (Bicinchoninic acid) Protein assay kit (SIGMA, Aldrich, UK).

DNA extraction
Mucus swab samples were thawed on ice and sliced in half lengthwise with a scalpel.Excess RNAlater™ was removed by gently squeezing swab heads between tissue to remove residual salt from the storage solution before transferring to a 2 mL Eppendorf tube for extraction.DNA was then extracted as described previously using the QIAamp Fast DNA Stool Mini Kit (Qiagen) according to the manufacturers protocol with modifications described in (Dehler et al., 2017).DNA was extracted from water filters using this same protocol as described in (Lorgen-Ritchie et al., 2021).

PCR amplification and sequencing
PCR reactions and subsequent library preparations were carried out as previously described (Lorgen-Ritchie et al., 2021) amplifying variable regions 3 and 4 of the 16S rRNA gene with the 341F/785R primer pair (Klindworth et al., 2013).Illumina adapter overhang sequences were added to the 5 ′ end of each primer.The forward primer (341F) had the sequence 5 ′ TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACG GGNGGGCWGCAG, and the reverse primer (785R) 5 ′ GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTA CHVGGGTATCTAATCC with the bold underlined sequence being the locus-specific V3-V4 primers.In summary, triplicate PCR reactions were performed in a 10 μL reaction including 2 μL of each forward and reverse primer (1 μM stock, Sigma), 5 μL of 2× KAPA HiFi HotStart ReadyMix including high-fidelity polymerase (KAPA Biosystems Ltd., UK) and 1 μL of DNA.PCR conditions included an initial denaturation at 95 • C for 3 min, followed by 27 cycles of 30 s at 98 • C, 30 s at 57 • C and 30 s at 72 • C after which a final extension of 72 • C for 5 min was applied.A subset of resulting PCR products were separated on an Agilent 2200 Tapestation (Agilent Technologies, Italy) for quality control and to verify amplification of correct sized product.The Nextera XT Index Kit (Illumina, San Diego, CA) was used to attached dual indices and Illumina sequencing adapters (P5 and P7) by PCR to the amplicons to produce the final libraries.Libraries were quantified using a Quant-iT High-Sensitivity dsDNA Assay (Thermofisher, USA) using an Omega FLUOstar plate reader (BMG Labtech, UK).Final libraries were pooled equimolarly.The final library was denatured and diluted to 1.2 nM prior to loading onto a MiSeq flow cell and sequencing on the Illumina MiSeq platform (Illumina, San Diego, CA).Two flow cells were required to obtain sufficient depth and samples were randomised between the libraries loaded onto the cells.Sequencing was performed on an Illumina MiSeq platform using a 2 × 300 bp paired end protocol with a targeted sequencing depth of 80,000 reads per sample.

Sequencing data bioinformatics
Analysis of sequence data were carried out using DADA2 (Callahan et al., 2016) in RStudio v1.1.456using R v3.6.1 as described previously (Lorgen-Ritchie et al., 2021).Briefly, adapters and primers were removed using TrimGalore! (https://github.com/FelixKrueger/TrimGalore) and reads with an overall Phred quality score less than 30 were discarded.Forward reads were truncated to 250 bp and reverse reads to 200 bp.Remaining reads were denoised, merged, screened for chimeric sequences which were subsequently removed, and assigned as distinct amplicon sequence variants (ASVs) using DADA2.Samples with less than 1000 reads were excluded from further analysis.Taxonomic classification of ASVs was carried out using the assignTaxonomy function of DADA2 using the Silva reference taxonomy v132 (Quast et al., 2013) and default parameters.Assignment of species was also conducted using the Silva species assignment v132, allowing for assignment of multiple species.Known non-relevant co-amplified products including mitochondrial, eukaryotic, cyanobacteria and chloroplast sequences were excluded from analysis as were those ASVs that were only sequenced as singletons.Samples with less than 500 reads following filtering were excluded from further analysis.Taxonomic composition of triplicate positive controls was in agreement with the mock community reference (Fig. S1).

Statistical analysis
Statistical analysis was carried out in RStudio v1.1.456using R v3.6.1 and the package phyloseq (McMurdie and Holmes, 2013) as described previously (Lorgen-Ritchie et al., 2021).Growth and water quality parameters were analysed by one-way ANOVA with Tukey's HSD post-hoc test.All samples were subsampled to an equal depth of 2622 reads before calculation of alpha and beta diversity.Differences in alpha diversity across sampling points was determined by Kruskal-Wallis comparisons of Shannon (Shannon, 1948) and Chao1 measurements (Chao, 1984) measurements followed by pairwise testing using the Wilcoxon rank sum test.
Community structure (beta diversity), determined by Bray-Curtis dissimilarity distance (Bray and Curtis, 1957) was visualized using non-metric multidimensional scaling (NMDS) ordination plots, implemented using the Vegan package (Oksanen et al., 2017) and plotted using ggplot2 (Wickham, 2016).Data ellipses based upon an assumed multivariate t-distribution were drawn at a level of 0.95 with stat-ellipse in ggplot2 to provide a visual summary.Permutational multivariate statistical analysis of community separation (PERMANOVA) was carried out using the Adonis function in the Vegan package and pairwise comparisons computed using adonis.pair in the EcolUtils package (Salazar, 2020).Core microbiota were identified using the microbiome R package (Lahti and Shetty, 2017) with a prevalence cut-off of 90% and a lower relative abundance limit of 0.1%.Log2 relative abundances of core ASVs across samples were presented in heatmaps drawn with Pheatmap (Kolde, 2012) within R, using Euclidean distance clustering of ASVs.
In order to identify functional pathways based upon 16S communities, Piphillin was used to normalize the non-rarefied amplicon data by 16S rRNA gene copy number and to infer metagenomic contents (Iwai et al., 2016;Narayan et al., 2020).A sequence identity cut-off of 99% was implemented.The inferred metagenomic functions were assigned using the Kyoto Encyclopedia of Genes and Genomes database (KEGG; May2020 Release) and KEGGREST (Tenenbaum, 2019) was utilised to obtain KEGG pathway names and BRITE hierarchies from pathway identifiers.STAMP v2.1.3(Parks et al., 2014) was used to test for statistically significant differences in pathway contributions to parent terms using Welch's t-test corrected for multiple-testing by Benjamini-Hochberg false discovery rate (FDR).Differences were considered significant at q < 0.05.

Sequencing outputs
In gill mucus samples 3,318,899 raw reads were obtained, for both forward and reverse reads with a mean read depth of 36,877 ± 3520 (SE).After quality filtering, denoising and chimera removal in DADA2, 1,260,333 reads with a mean of 14,004 ± 1479 (SE) per sample were retained.For skin samples, 4,398,912 raw reads were obtained at a mean read depth of 48,877 ± 7453 (SE) with 1,796,959 reads at a mean read depth of 19,966 ± 3550 (SE) remaining following refinement.For water samples, 873,217 reads with a mean read depth of 72,768 ± 11,805 (SE) per sample were obtained and 516,367 reads at a mean read depth of 43,031 ± 7037 (SE) per sample were retained following refinement.Samples yielding less than 1000 reads were excluded from further analysis.Of the 90 gill samples sequenced, 86 were retained for downstream analysis while 87 out of 90 skin samples and all 12 water samples were retained.

Fish growth and condition factor
Length and weight increased significantly throughout FW (p < 0.001), but no further significant increases were observed in the first four weeks post-SWT, considered to be due to initial loss of appetite upon transfer to seawater (Table 1).In line with smoltification, condition factor showed a decline throughout the study period (p < 0.001).At the time of seawater transfer, fish had a mean smolt index score of 3.86 ± 0.25 and gill NKA activity of 7.76 ± 0.50 μmol mg prot − 1 h − 1 respectively.Specific growth rate was positive throughout FW sampling points, but significantly declined upon initial SWT (p < 0.05).SGR across the entire study equated to 1.12% body weight per day.

RAS water quality
A number of water quality parameters were measured in header tanks immediately upstream of rearing tanks in the FW RAS and these varied over time throughout the duration of the study (Table 2).

Alpha diversity
Alpha diversity was temporally dynamic in both gill (Shannon F = 45.8, p < 0.001; Fig. 1A1) and skin (Shannon F = 62.2, p < 0.001; Fig. 1A2) mucus.Initial Shannon diversity was higher in skin than in gill, with levels in skin comparable to that of the water in the RAS.Diversity showed a sharp and significant decline between FW1 and FW2 in skin and gill mucus (both p < 0.001).Following this decline, diversity showed an increasing temporal trend from FW2 to FW4 in skin and gill mucus (both p < 0.001), recovering to FW1 levels by FW3 in gill, and FW4 in skin.Diversity remained constant post SWT (FW4 to SW1) and over time in SW in skin and gill.
Chao1 species richness was also temporally dynamic in gill (Chao1 F = 28.5, p < 0.001; Fig. 1B1) and skin (Chao1 F = 42.1, p < 0.001; Fig. 1B2).Richness in skin mucus showed a significant decline between FW1 to FW2 (p < 0.001) followed by a subsequent increase from FW2 to FW4 (p < 0.001).In gill however, no significant decline in richness was observed between FW1 and FW2, but an increasing trend was observed from FW2 to FW4 (p < 0.05) and significantly higher richness was observed at FW3 compared to FW2 (p < 0.01).Consistent with observations in diversity, no changes in richness were observed following SWT (FW4 to SW1) and richness remained constant at both sampling points in SW in gill mucus, but a significant decline was observed in skin mucus from SW1 to SW2 (p < 0.05).No significant impact of sampling point was observed on alpha metrics in all water samples (Figs.1A3 and  B3), likely due to lack of statistical power with small sample numbers, but temporal patterns mirrored those observed in gill and skin mucus with a decline in diversity from FW1 to FW2 followed by recovery to FW4.Diversity and richness in water samples declined from FW4 to SW1 and remained consistent from SW1 to SW2.

Beta diversity
PERMANOVA analysis revealed significant differences in the microbiome structure in gill (F 5,71 = 31.126,R 2 = 0.687, p < 0.001; Fig. 2A) and skin (F 5,78 = 28.452,R 2 = 0.646, p < 0.001; Fig. 2B) mucus over time.The largest separation was evident between the FW and SW samples, and all pairwise comparisons of sampling points were significantly different (p < 0.001).Less separation was observed between sources (i.e.gill, skin or water) than between sampling points within each sample type.PERMANOVA analysis revealed significant differences in the microbiome structure between gill mucus, skin mucus and water samples (F 2,169 = 3.892, R 2 = 0.044, p < 0.001; Fig. 2C).Pairwise comparisons identified significant differences in beta diversity between all source contrasts: gill vs skin (p < 0.01); gill vs water (p < 0.05); skin vs water (p < 0.01).SW water samples showed greater separation from gill and skin samples than FW water samples.

Community composition in skin and gill mucus, water and diet
Microbial community composition at the phylum (Fig. S2) and ASV (Table S2) level were determined at each sampling point for gill mucus, skin mucus and water samples.Relative abundances of dominant taxa at the phylum level are presented in Fig. 3A.Phyla present at <1% were grouped as 'Other'.Seven phyla were present at >1% in gill mucus, six in skin mucus and seven in water.Proteobacteria and Bacteroidetes were the dominating phyla throughout in all sample types.Proteobacteria alone constituted 85.8% of the total abundance in gill mucus, 82.3% in skin mucus and 73.8% in tank water at FW2.This proportion further increased to 81.9% in tank water at FW3. Proteobacteria and Bacteroidetes constituted >90% of total relative abundance in skin mucus at all sampling points, but in both gill and skin mucus, a temporal increase in the proportion of Bacteroidetes and a decrease in Proteobacteria was observed from FW2 to FW4, which was more pronounced in gill.In contrast, Proteobacteria increased in relative abundance from FW1 to FW3 in tank water before declining at FW4 while Bacteroidetes showed the inverse pattern.Post-SWT, Verrucomicrobia increased in prevalence in gill mucus (SW: 21.3%; SW2: 13.4%) and constituted a higher relative abundance than Bacteroidetes (SW: 16.4%; SW2: 11.8%).Verrucomicrobia was observed at >1% in water samples at all sampling points, but peaked at FW1 (12.7%).Relative abundance of Actinobacteria increased in water post-SWT, becoming more dominant than Bacteroidetes by SW2 (26.5%).Plactomycetes relative abundance was also found at a higher level in SW at SW2 (9.0%).
Relative abundances of dominant taxa at the ASV level are presented in Fig. 3B.When considering composition at the ASV level, it became apparent that the clear increase in abundance in Proteobacteria observed in gill mucus at FW2 was a result of domination by two ASVs belonging to the genus Hydrogenophaga (ASV7 and ASV18).These two ASVs constituted 71.7% of total relative abundance in gill mucus at FW2 and 42.3% at FW3.These two ASVs also dominated skin mucus (FW2: 62.3%; FW3: 40.8%) and tank water communities (FW2: 58.8%; FW3: 72.9%).The top 40 ASVs generally constituted a greater proportion of the total community in gill mucus compared to skin mucus, particularly at FW4 and SW1.

'Core' microbiota
Taxa were considered as 'core' if present in at least 90% of samples at a minimum of 0.1% relative abundance.Cores were considered separately for each mucosal tissue across all temporal samples, by salinity, and finally by individual sampling point (Table S3 and Fig. S3).A total of 108 cores were identified in in gill mucosal samples and 140 in skin.Of these taxa 83 were identified as core in both tissues while gill harboured 25 unique core taxa, and skin 57.No overall core taxa were identified when considering all sampling points together.Considering FW only, two core taxa were identified in gill mucus; ASV7 -Hydrogenophaga and ASV24 -Flavobacterium.ASV7 was also identified as a core FW taxon in   skin mucus alongside ASV37 -Flavobacterium.In SW, five core ASVs (belonging to Gammaproteobacteria) were identified in both mucosal surfaces: ASV99 and ASV229 -Acinetobacter; ASV22 -Psychrobacter; ASV74 -Citrobacter; ASV201 -Pseudomonas.Two additional SW cores were identified in gill mucus (ASV78 -Rubritalea; ASV89 -Shewanella) and one additional skin mucus SW core was observed (ASV227 -Stenotrophomonas).Upon removal of the dysbiotic sampling points (FW2 and FW3) from the core analysis, ASV99 (Acinetobacter) was identified as an overall core in skin mucus.ASV99 was also identified as a FW core in skin alongside ASV167, also of the genus Acinetobacter.ASV7 was no longer considered a FW core in either gill or skin mucus.Considering sampling points individually a total of 4, 5, 10, 14, 30 and 60 core taxa were identified in gill mucus at samplings from FW1 to SW2 respectively, while in skin 44, 6, 14, 15, 31 and 62 were identified.
We next considered taxa unique to or shared between individual sampling points in each mucosal surface, and the proportions of which were unique from the bacterial community in the surrounding water (Fig. 4).At FW1, 61.4% of taxa identified in gill mucus were unique to FW1 compared to 25% in gill mucus.Of those unique to FW1, 18.5% in skin mucus were not shared with water, while in gill mucus, all were identified in water.Of those shared with other sampling points, 66.7% in gill mucus were unique from the water microbiome compared to 23.5% in skin mucus.Proportions of unique taxa fell to zero at FW2 indicating no taxa unique to this sampling point, and all identified core taxa were shared with the surrounding tank water, as was also the case at FW3.From FW2 onward, the number of core taxa along with the percentage of taxa unique to individual sampling points rose steadily to a maximum of 91.7% in gill at SW2 and 83.9% in skin at SW1.At FW4, the proportion of taxa unique from the surrounding water increased to 33.3% of taxa shared with multiple timepoints in skin mucus, and to 11.1% of those unique to FW4, but remained at zero in gill mucus.By SW1, 100% of taxa shared between timepoints were unique to the surrounding water in both mucosal surfaces, while of those unique to SW1, 84.6% in skin and 44.0% in gill were absent from water.These proportions rose to 94.5 and 98.0% in gill and skin respectively by SW2.

Microbiota associated with nitrogen cycling in RAS
Bacteria known to be associated with nitrogen cycling in RAS biofilters were identified in water, gill and skin samples (Fig. 5).Ammonia oxidising bacteria (AOB; Nitrosomonas and Nitrosospira) and nitrite oxidising bacteria (NOB; Nitrospira and Nitrobacter) were identified at similarly low levels in all sample types.Despite low levels throughout, a peak in Nitrosomonas was observed at FW4 in water and also in Nitrospira in water, gill and skin samples.The heterotrophic de-nitrifier Pseudomonas peaked in all sample types at the initial sampling point, but another de-nitrifier, Comamonas, was notably absent, or at very low levels, throughout the study.A dramatic increase in the relative abundance of the autotrophic de-nitrifier Hydrogenophaga was observed at FW2 and FW3 in water, reaching over 70% of total relative abundance at FW3.This surge was also observed at comparable magnitude in gill and skin, but peak relative abundance was reached at FW2 rather than FW3 as found in the water.

Functional enrichment for processes in skin and gill microbial communities
Temporal fluctuations in microbial community composition may be indicative of a change in microbiome functionality over time which can be inferred from 16S sequences using programmes which match 16S data to bacterial genomes and thus a suite of potentially activated genes.Although such methods have inherent pitfalls including a lack of sequenced aquatic genomes and no incorporation of transcriptomic data, they can be useful in presenting an overview of potential functionality.One such programme, Piphillin, inferred 375 pathways from 4554 ASVs present in Atlantic salmon gill mucus at an identity cut-off of 99%.Removing human diseases and top-level terms, 294 pathways remained and 155 (52.7%) of these pathways were related to metabolism.In skin, 376 pathways were inferred from 5332 ASVs, with 156 of 295 pathways (52.9%) related to metabolism.Contributions of level 2 terms to Metabolism were temporally dynamic in both gill and skin mucus and temporal dynamics were highly similar between the two sample types (Fig. S4).
Statistical analysis of Piphillin outputs using STAMP revealed reciprocal temporal dynamics in the contributions of Metabolism level 2 terms 'Xenobiotic degradation and metabolism' and 'Carbohydrate metabolism' (p < 0.001) over time, coincident with the peak in Hydrogenophaga at FW2 and FW3 in both gill and skin (Fig. 6).Contributions of 'Metabolism of terpenoids and polyketides' and 'Lipid metabolism' were also significantly increased in both mucosal surfaces during the increase in Hydrogenophaga abundance.Notably, in skin, 'Amino acid metabolism' constituted a large percentage of all pathways and showed a significant increase in contribution at FW2.Amino acid metabolism will have a relationship with nitrogen and ammonia levels, relating back to the observed increases in nitrogen cycle products in tank water.
Analysing functional data at the individual pathways related to 'Xenobiotic metabolism and degradation', the contribution of 'Polycyclic aromatic hydrocarbon degradation' was the term showing the largest increase in contribution coincident with Hydrogenophaga dominance at FW2 in both gill and skin mucus, followed by 'Aminobenzoate degradation', 'Benzoate degradation', 'Chlorocyclohexane and chlorobenzene degradation' and 'Atrazine degradation' (Fig. S5).For amino acid pathways in skin, the largest change in contribution was observed for "Valine, leucine and isoleucine degradation", followed by Phenylalanine metabolism" and the same pattern was observed in gill (Fig. S6).

Discussion
Microbial diversity, richness and core taxa in skin and gill mucus of juvenile Atlantic salmon reared in a RAS system were temporally dynamic.Dynamics in richness and diversity were largely similar in the two mucosal tissues, and to some extent mirrored those of the surrounding water.The proportion of sampling point-specific cores became increasingly distinct from water microbial communities over time.An abrupt decline in diversity during FW production coincided with a clear increase in abundance of two taxa belonging to the RAS-biofilter-related genus Hydrogenophaga in RAS tank water as well as in gill and skin mucus.Although this dominance resulted in a temporary loss of diversity in mucosa, there were no obvious physical implications and specific growth rate was comparable with a second RAS stream (unpublished data).Functional analysis showed an increased contribution of metabolic pathways relating to degradation and metabolism of xenobiotic compounds and amino acid metabolism in both skin and gill coincident with the domination of Hydrogenophaga in FW.To our knowledge this is the first temporal study of external mucosal surfaces of Atlantic salmon during FW smoltification in a RAS and results show that the microbial communities in these systems can be subject to abrupt and dramatic changes over a short time period, likely reflecting water quality/chemistry and biofilter fluctuations.

Alpha and beta diversity are temporally dynamic in gill and skin mucus in RAS
In an enclosed RAS, incremental increases in particulate organic matter as a consequences of increasing organic load may provide additional substrate, increasing the carrying capacity and thus enabling bacterial communities to expand in number (Attramadal et al., 2014;Fossmark et al., 2020).Microbial diversity in skin and gill declined significantly between FW1 and FW2 before increasing throughout FW4 to a level maintained post-SWT and similar dynamics were observed for richness.These findings differ from those of a previous study where the skin microbiome of Atlantic salmon reared in flow-through tanks in FW exhibited lower evenness, and higher richness and diversity posttransfer to a SW flow-through system (Lokesh and Kiron, 2016).Similarly, in a wild population of Arctic charr (Salvelinus alpinus), Shannon diversity was also lower post-SW migration (Hamilton et al., 2019).Such differing observations may stem from factors including differences in FW rearing system (flow-through vs RAS), post-SWT on-growing system (flow-through vs open sea loch), FW rearing temperature (12 • C vs 15 • C) or SW temperature (12 • C vs ambient).Indeed, a recent study showed significantly higher (2-fold) microbial richness in skin mucus from Atlantic salmon parr reared in two FW RAS hatcheries when compared to those reared in a flow-through hatchery (Minich et al., 2020a).
Seawater transfer has been shown to re-shape the skin microbiome in Atlantic salmon reared in FW flow-through tanks (Lokesh and Kiron, 2016) and distinct communities were observed in Arctic charr in FW and SW (Hamilton et al., 2019).Such shifts may in effect 're-set' the external microbial communities, resulting in a period of tolerance in mucosal immune function to allow a new community to establish as a consequence of introduction to novel microbes encountered in the saline environment.Distinct microbial community composition pre-and post-SWT were also observed in fish reared in a FW RAS as well as between sampling points during both FW and SW phases.Although statistical analysis identified distinct community composition by both mucosal tissue and sampling point, greater distinction was observed temporally than by tissue.Fish body site was the strongest driver of community composition in samples taken at a single sampling point from both RAS and FT systems (Minich et al., 2020a) and a study of five mucosal sites (skin, gill, olfactory rosettes, anterior and posterior gut) in rainbow trout (Oncorhynchus mykiss) also identified a significant effect of body site on community composition with skin hosting the highest diversity (Lowrey et al., 2015).The stability of gill and skin microbial communities showed an increase in richness over time (38 sampling points across a period of one year) in wild marine Pacific chub mackerel (Scomber japonicus) with gill and skin showing very similar dynamics (Minich et al., 2020b).In free-living marine microbes, a study conducted over a 6-year period identified day length as the primary driver of diversity in microbial communities and identified richness to be highest during the winter in the North Atlantic (Gilbert et al., 2012).
Interestingly, in the current study, the identified dysbiosis occurred during a time when fish were exposed to a 'winter' photoperiod as a signal to induce smoltification (McCormick, 2012).The occurrence of dysbiosis could potentially be a direct consequence of the short photoperiod, associated with changes in osmoregulatory physiology known to accompany smoltification or the result of an interaction between environmental and host factors.For example, water temperature is also a strong predictor of microbiome dynamics with dysbiotic events and associated proliferation of potentially pathogenic taxa associated with warmer temperatures or shifts from cold to warm temperatures annually (Rosado et al., 2021).Indeed, in the current study, a period of increased FW rearing water temperature coincided with a reduction in microbial diversity in both water and external mucosa.Such dynamics are also important to consider in aquaculture in terms of year-round production where different cohorts of fish are transferred to sea at different times of the year and thus experience different temperature dynamics.The fish in this study were transferred to sea in summer and as such microbiome dynamics during the seawater transfer period would likely not be applicable to a cohort of fish transferred in the winter months.

Shared prevalent taxa in gill and skin mucus of Atlantic salmon
Proteobacteria and Bacteroidetes were the dominating phyla observed throughout the study in both gill and skin mucus, in keeping with previous studies (Lowrey et al., 2015;Lokesh and Kiron, 2016;Brown et al., 2019;Minich et al., 2020a;Brown et al., 2021).The relative abundance of Proteobacteria in the skin microbiome was reported to increase post-SWT in Atlantic salmon reared in a FW flow-through system (Lokesh and Kiron, 2016).Contrasting results were obtained in the present study where FW production was in a RAS, and Proteobacteria made up 67.4-82.3% of the total community throughout both FW and SW sampling points.In gill mucus, Proteobacteria abundances were more temporally dynamic and the relative abundance almost doubled between FW1 and FW2, driven by two ASVs belonging to the genus Hydrogenophaga.In SW gill mucosal communities, a third phylum, Verrucomicrobia, became prevalent, driven by the genus Rubritalea.
Core taxa were defined as those present in 90% of individuals at a minimum of 0.1% relative abundance.A core microbiome of 11 OTUs was previously found in Atlantic salmon skin across different populations of wild and hatchery fish which consisted predominantly of Proteobacteria (Uren Webster et al., 2018).In another study, 75 core OTUs were identified in FW, and 139, 97 and 100 at one, two and four weeks post-SWT with 19 shared core OTUs across all sampling points (Lokesh and Kiron, 2016).In the present study, no core taxa were identified across either gill or skin when considering all sampling points together, but it is important to consider that previous studies utilised OTUs as opposed to more specific ASVs.Single ASVs belonging to the genera Hydrogenophaga (Gammaproteobacteria) and Flavobacterium (Bacteroidia) were identified as core taxa across FW in both gill and skin, although each mucosal tissue harboured a different core Flavobacterium ASV and overall relative abundance of the genus was higher in gill than in skin.This is consistent with a study in rainbow trout which found Flavobacterium to comprise up to 61.7% of the gill and 24.0% of the skin microbiome (Lowrey et al., 2015).A lower relative abundance of Flavobacterium on salmon skin in SW compared to FW has also been reported (Lokesh and Kiron, 2016) and indeed in the current study, Flavobacterium was not identified as a core in SW in either skin or gill.The presence of Hydrogenophaga as a core FW taxon is potentially related to the association of this genus with RAS biofilters (Rurangwa and Verdegem, 2015).

Dysbiosis results in loss of host-specific core taxa
Considering sampling points individually, the number of cores were not stable over time, neither were the proportions of core taxa unique to each sampling point, suggesting dynamic communities in gill and skin mucus with few sustained cores at the ASV level.Higher numbers of core taxa at SW sampling points than those in FW could be a result of free water flow between sea cages, however, as much as 98% of core taxa were not found in the corresponding SW samples indicating that surrounding water alone could not be responsible for the expansion of cores at SW sampling points.Coincident with temporary domination by two Hydrogenophaga taxa and a loss of alpha diversity and richness at FW2, core taxa in gill and skin mucus which were host-specific (i.e.not identified in the surrounding water) disappeared, suggesting a dysbiosis in mucosal membranes.Alteration of the microbiome does not necessarily have a negative impact on the host but does provide the opportunity for opportunistic taxa to establish themselves if competitive exclusion in the host microbiome is disturbed.
Analysis of shared and unique taxa in gill and skin mucus also indicated the existence of a common microbiome in external mucosal surfaces and supports the suggestion of microbial exchange between the two external niches (Rosado et al., 2021).The proportion of cores shared by gill and skin mucus showed an increasing trend as time proceeded in FW, which initially declined post-SWT before re-establishing.This is in contrast to previous work which observed a lack of core taxa across mucosal layers including gill and skin in rainbow trout (Lowrey et al., 2015), however this study analysed just 6 fish which will have an impact on the assignment to core taxa.Interestingly, in this study, skin mucus shared no unique core taxa with seawater while gill shared 22-53% of unique cores, suggesting a stronger influence of the seawater microbial community on microbial communities in gill mucus than those of skin, perhaps due to active movement of water across the gills during osmotic regulation (Hoar, 1988).Thus, the observed higher proportions of core taxa shared with water in FW compared to SW could reflect changes in osmotic gradients with salinity transfer.Care must be taken when interpreting these results as a single pooled water sample was sequenced at each sampling point in SW.
Water microbial communities are known to influence richness and diversity of the fish skin microbiome to a greater extent than the gut microbiome where dietary diversity is thought to have a more pronounced effect (Boutin et al., 2013;Giatsis et al., 2014;Uren Webster et al., 2018).However, many studies to date indicated distinct communities in fish skin and water (Minniti et al., 2017;Uren Webster et al., 2018), despite the direct contact of skin with water.The current study suggests that shifts in surrounding water microbial communities can impact the microbiomes of skin and gill mucus and a clear overlap between temporal dynamics was observed in a FW RAS.However, a high proportion of core ASVs present in skin and gill mucus were not identified in corresponding water samples, particularly in SW, supporting selection of the mucosal microbiota from that of the environment and the potential for maintenance of essential components of the microbial community.

Nitrogen-cycling taxa identified in external microbiomes
Diets in aquaculture are protein-rich in order to enhance growth, and this can result in an increase in metabolic nitrogen waste (Ip et al., 2004).Ammonia is toxic to fish (Schram et al., 2010) and concentrations in RAS should be kept low to prevent damage to tissues.Temporal accumulation of organic matter and nitrogenous compounds are a consequence of the closed-nature of RAS, and microbes play a key role in maintaining water quality in RAS.Previous studies revealed distinct microbial communities in a RAS biofilter, tank water and mucosal surfaces, but also observed that various microbes responsible for nitrogen cycling, present in biofilters were also detectable in fish mucus, suggesting that some microbes involved in nitrogen cycling are not simply restricted to the biofilter, but can also circulate through RAS systems, potentially colonizing fish mucosal surfaces (Schmidt et al., 2016;Minich et al., 2020a).When organic matter accumulates in a RAS, heterotrophic blooms can occur, outcompeting nitrifying microbes as heterotrophs obtain carbon and energy from organic matter (Leonard et al., 2000;Rurangwa and Verdegem, 2015).Primary heterotrophic microbes associated with denitrification in RAS include Pseudomonas and Paracoccus, and primary autotrophs include Rhodobacter and Hydrogenophaga (Rurangwa and Verdegem, 2015).Dominance by two ASVs assigned to the facultative autotrophic de-nitrifier Hydrogenophaga (Xing et al., 2018) in RAS tank water samples was observed at FW2 and FW3 and these were also dominat in gill and skin microbiomes, supporting that the water microbiome does have the potential to impact external mucosal microbial communities (Uren Webster et al., 2018;Lokesh et al., 2019;Gupta et al., 2019).However, diversity recovered more rapidly in the mucosal layers than in the tank water, supporting the role of the host in regulating mucosal microbiomes.
In the current study we observed rising trends in ammonia, nitrate and CO 2 levels in water in the RAS over production time, with the largest increase in nitrate coinciding with increased Hydrogenophaga abundance at FW2, suggesting that RAS tank water chemistry may directly impact external microbiomes.This is in line with findings in an array of fish species residing in a heterogeneous inland water system which identified a role for sporadic nutrient pollution events in driving dysbiosis in external microbiomes (Krotman et al., 2020) and increased presence of Hydrogenophaga coincided with an algal bloom associated with pH and nitrogen/phosphorous levels (Li et al., 2012).Coincident with the temporary dysbiosis and dominance of Hydrogenophaga, functional analysis identified an increased contribution of metabolic pathways relating to 'Xenobiotics metabolism and degradation' in gill and skin mucus at the expense of 'Carbohydrate metabolism'.Amino acid metabolism pathways also increased in contribution, particularly "Valine, leucine and isoleucine degradation" and phenylalanine metabolism.Xenobiotic compounds are those non-natural to the environment and this suggests that dysbiosis could be related to the temporal dynamics of organic and inorganic compounds in the RAS (Vadstein et al., 2012).At the pathway level, 'Polycyclic aromatic hydrocarbon (PAH) degradation' showed an increased contribution coincident with dysbiosis.

Conclusions
Atlantic salmon gill and skin mucus-associated microbial communities were temporally dynamic in a RAS.An overdominance in RAS biofilter-associated Hydrogenophaga in RAS tank water significantly impacted diversity, richness and the presence of host-associated core taxa in the external mucosal surfaces, but diversity was quickly reestablished in mucosal tissues.Distinct shifts in microbial communities were also observed in both skin and gill mucus following transfer to a completely novel seawater environment where fish require the plasticity to re-establish a new community suited to the SW niche.The results of this study highlight the importance of considering temporal dynamics, often not considered in earlier microbiome research, as well as the potential for deviations in RAS water microbial composition to temporarily destabilise mucosa-associated communities.Although RAS are considered stable in terms of environment, changes in water quality and chemistry over time can have significant and acute impacts on the microbial community not just in water, but also in external mucosal tissues.Longer-term impacts of distinct community shifts on health and immune function are yet to be determined.

Fig. 1 .Fig. 2 .
Fig. 1.Alpha diversity (A: Shannon) and richness (B: Chao1) comparisons.(A1) Shannon: gill mucus, (A2) Shannon: skin mucus and (A3) Shannon water across sampling points.(B1) Chao1: gill mucus, (B2) Chao1: skin mucus and (B3) Chao1: water.Red circles indicate individual samples for FW RAS and blue triangles are individual samples from SW cages.Superscripts indicate significant differences between sampling points derived from pairwise testing.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. (A) Total relative abundance of each phyla across sampling points in gill mucus, skin mucus and water samples.Only phyla constituting >1% of total relative abundance at each timepoint are presented individually, with those at lower abundances grouped together in 'Other'.(B) Relative abundance of the 40 most abundant taxa in each sample type, coloured by genus.

Fig. 4 .
Fig. 4. Shared and unique core ASVs.Number of core ASVs either shared between sampling points (hatched columns) or unique (filled columns) to each sampling point in (A) gill and (B) skin plotted on y-axis 1.The percentage of these shared (dashed line) or unique (solid line) cores not identified in corresponding water samples is plotted on y-axis 2. Taxa were designated as core if present in >90% of individuals at a minimum relative abundance of 0.1%.

Fig. 6 .
Fig. 6.Functional annotation inferred by Piphillin.Significant differences in 'Metabolism' level 2 terms as proportion of parent terms prior to and during Hydrogenophaga dominance (FW1 and FW2) in gill and skin mucus.Benjamini-Hochberg adjusted p-values are presented (q).

Table 1
Length, weight, condition factor and specific growth rate in Atlantic salmon smolts during smolting (n = 2 or 3 tanks, 12 to 18 fish total, means ± SD).Fork length was measured in centimetres and weight in grams.CF: condition factor, SGR: specific growth rate, NKA: ATPase activity.Degree days expressed relative to onset of spring photoperiod (24 L:0D).Significant changes in metrics over the study period were determined by ANOVA and superscript letters denote pairwise significance determined using Tukey's HSD test.Modified from Lorgen-Ritchie et al. (2021) with the addition of smolt score, plasma chloride and gill NKA.

Table 2
Water quality parameters within the FW RAS stream.Measurements made daily in header tanks immediately upstream of rearing tanks.Data presented for sampling days; May 14th (FW1), June 18th (FW2), July 2nd (FW3) and July 25th (FW4).Mean, minimum (min) and maximum (max) values determined from May 14th to July 25th 2019 are also presented.