Interpopulation Variation in the Atlantic Salmon Microbiome Reflects Environmental and Genetic Diversity

Variation in the microbiome has a fundamental influence on host health, ecology, and evolution, but the scope and basis of this variation are not fully understood. We identified considerable variation in skin and gut microbial communities between seven wild and captive populations of Atlantic salmon, reflecting divergent environmental conditions and fish genetic diversity. In particular, we found very pronounced differences in the intestinal microbiomes of wild and hatchery-reared fish, likely reflecting differences in diet. Our results offer an insight into how the microbiome potentially contributes to the generation of local adaptations in this species and how domestication alters intestinal microbial communities, highlighting future research directions in these areas.


DNA extraction and 16S rRNA sequencing
DNA extraction from all gut and skin swab samples was performed using MoBio PowerSoil® DNA Isolation Kit (Cambio, Cambridge, UK) according to the manufacturer's instructions, with an additional incubation step of 10 minutes at 65 °C prior to bead beating (Llewellyn et al. 2016). Water samples were centrifuged at 5000xg for 1 hour at 4°C before DNA was extracted from the pellet. DNA concentration and purity was assessed using a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies, Wilmington, USA). 16S library preparation using Nextera XT Index kit was performed according the Illumina Metagenomic Sequencing Library Preparation guide using 12.5 ng total genomic DNA. The V4 hypervariable region of the bacterial 16S gene was amplified using the primers selected as the best candidates for bacterial and archaeal representation (Klindworth et al. 2013); 519F (5'-CAGCMGCCGCGGTAA) and 785R (5'-TACNVGGGTATCTAATCC), each with 5' tags for barcode attachment. Reaction conditions for the first PCR amplification consisted of an initial denaturation at 95°C for 3 min, followed by 25 cycles of 95°C for 30s, 55°C for 30s Fisher Scientific, UK) and pooled in equal concentrations, before sequencing using an Illumina MiSeq (300 bp PE reads).

16S rRNA bioinformatics analysis
Adaptor contamination and poor quality bases from the 3' end were removed from the raw sequence reads using a sliding window of 4 bp and a minimum quality score of Q=20 in Trimmomatic (Bolger et al. 2014). Microbial community analysis was then conducted using mother v1.37 (Kozich et al. 2013), Qiiime v1.9 (Caporaso et al. 2010) and R v3.3.2 (R_Core_Team 2014). Forward and reverse reads were then fully overlapped and merged, and then filtered to retain amplicons of the target size range (260-300 bp) and remove those containing ambiguous bases. Contigs were aligned to the Silva seed reference database (version 123)  which was first trimmed to include only the target V4 region. For contigs with high quality alignments, further noise reduction was performed using mothur's pre-clustering algorithm. Potential chimeras were then removed through implementation of UCHIME (Edgar et al. 2011) in mothur before taxonomic classification using the Silva reference taxonomy and removal of mitochondrial, eukaryote and chloroplast sequences.
Analysis of microbial community alpha diversity was performed at the operational taxonomic unit (OTU) level, following clustering of sequences in mothur based on 97% sequence similarity. Several samples had low numbers of good quality microbial sequence reads. In order to maximise sample inclusion whilst ensuring high Good's coverage (≥ 94%) for all included samples, reads were subsampled to a depth of 4012/sample. Seventy six gut and 81 skin samples (min of 10 samples per population), together with a water sample from each site, were retained and used for analysis of alpha diversity. We calculated two measures of alpha diversity (Chao1 richness and Shannon diversity) to examine variation in microbial diversity using mothur. We analysed variation in alpha diversity by linear mixed modelling with the lme4 package in R using group (hatchery vs. wild), fork length, condition factor, sex, individual heterozygosity and individual MHC heterozygosity as fixed factors, and population as a random factor to account for spatial autocorrelation. Given the influential role of water on the fish microbiome, we included the microbiome diversity of the water samples at each site as an offset to statistically control for the effects of the surrounding water on fish microbial diversity. Model simplification was achieved by single deletion tests using the drop1 command via maximum likelihood; the minimal adequate model was then refitted by restricted maximum likelihood (Zuur et al. 2009), and the Satterthwaite approximation was used to obtain approximate significance levels using the lmerTest package. Variance components were calculated using the VarCorr function.
Analysis of microbial community beta diversity was performed based on weighted Unifrac distances and the Bray-Curtis similarity index. Phylogenetic trees were constructed using all retained sequences, using the relaxed neighbour joining method implemented in Clearcut (Evans et al. 2006), within mothur, then used to calculate weighted Unifrac (Lozupone et al. 2011) distances between samples. Bray-Curtis distances between samples were calculated using mothur based on OTU assignment. Non-metric multidimensional scaling analysis was performed separately for gut and skin samples in mothur, for both measures of structural diversity. After checking that data met assumptions of homogeneity of variance using Betadisper, statistical analysis of Unifrac and Bray-Curtis distances was performed using Adonis (both within the Vegan package for R; (Oksanen et al. 2017)). We assessed effects of origin (wild/hatchery), population, fork length, condition factor, sex, individual heterozygosity and individual MHC heterozygosity on gut and skin structural community variation, and used the strata function to specify a nested model of population within in group origin. Additionally, HOMOVA (Stewart & Excoffier 1996), within mothur, was used to specifically quantify the degree of intra-population variation in community structure, and test for statistical differences in this degree of variance between populations. To investigate a hypothesised influence of fish genetic background on microbial community structure a Mantel test, using the Pearson method, was employed in mothur to test for correlation between individual level genetic distances and weighted Unifrac/ Bray-Curtis distances for both the gut and the skin.
OTUs that were present in at least 80% of all individuals, as well as separately in 80% of wild and in 80% of hatchery fish, were identified for the gut and for the skin samples using the compute_core_microbiome function in Qiime. Following filtering of singleton OTUs from the dataset, a Kruskal-Wallace test, incorporating FDR calculation for multiple testing, was implemented using the group_significance function in Qiime to identify skin and gut OTUs with significantly differential abundance between fish from a wild and hatchery origin. In addition to these OTU-level analyses further community composition analysis was performed at the Phylum level using the summarize_taxa function in Qiime, followed by cluster analysis of all gut, skin and water samples using the Bray-Curtis similarity index and visualisation using FigTree v1.4.3 (Rambaut 2007).