Amoebic Gill Disease (AGD) in Atlantic Salmon Investigated Through a Holo‐Omic Lens

TL;DR

This study explores how Atlantic salmon's gill microbiota and genetics interact to influence resistance to amoebic gill disease, using a combined genetic and microbial approach.

Contribution

The study introduces the gill microbiota as an extended resistance trait in AGD research, integrating 16S rRNA profiling with quantitative genetics.

Findings

The gill microbiota of AGD-affected salmon is dominated by Simkaniaceae and Arcobacteraceae.

Microbial diversity and Simkaniaceae abundance show moderate associations with resistance indicators like gill score and amoebic load.

Genome-wide associations suggest several genomic regions linked to resistance traits and microbiota characteristics.

Abstract

Interactions between host genetics and the resident microbiota are complex. Understanding these interactions offers interesting alternatives for addressing gill health and disease resistance in salmonids. Amoebic gill disease (AGD), caused by Neoparamoeba perurans, remains a threat to Atlantic salmon, particularly in aquaculture where prevention and treatment options are scant. Selective breeding or genetic engineering towards AGD resilience presents viable prevention strategies. While several studies have addressed AGD resistance in Atlantic salmon using transcriptomic and quantitative genetic approaches, the influence of the gill microbiota as a genotypic encoded phenotype for AGD resilience remains underexplored. Addressing this, we leveraged a holo‐omic approach using 16S rRNA profiling and quantitative genetics, treating the microbiota as an extended resistance trait. In this small‐scale exploratory work, we investigated the microbiota of AGD‐challenged Atlantic salmon in terms of two common resistance indicator traits: gill score and amoebic load. We then performed a GWAS using the traditional indicator traits and traits of the microbiota. We found that the gill microbiota of the AGD‐affected salmon was dominated by two bacterial families, Simkaniaceae and Arcobacteraceae. Additionally microbial diversity and the relative abundance of Simkaniaceae, potentially derived from the amoeba, showed moderate variation with indicator traits. We identified several genomic regions that showed suggestive association with gill score and traits of the microbiota, and explored genes in these regions in relation to AGD resistance. Although the underlying mechanisms shaping gill microbiota dynamics in gill disease remain unknown, this study highlights the potential of addressing AGD through an integrative approach that considers the interplay between host genetics, the microbiota, and their roles in disease resistance.

In this exploratory study, we integrated 16S rRNA profiling with quantitative genetics to examine the gill microbiota of Atlantic salmon during infection with Neoparamoeba perurans, the causative agent of amoebic gill disease (AGD), treating the microbiota as an extended resistance trait. Microbial communities were dominated by Simkaniaceae and Arcobacteraceae, and microbial diversity and the relative abundance of key taxa showed moderate associations with resistance indicators, gill score and amoebic load. Genome‐wide association analyses using both disease traits and microbial traits indicated several genomic regions with suggestive associations, providing a preliminary basis for the future implementation of microbiota data in AGD resistance studies.

Figures

Left panel shows the comparison of alpha diversity estimates based on hill numbers at orders of diversity (q = 0, 1, 2) which apply increasing weight to the relative abundance of taxa, among samples grouped by (A–C) amoebic load Ct value where the colours represent Ct < 30 (grey), 30 < Ct > 35 (pink) and Ct > 35 (purple), and (D–F) gill damage (light–moderate, gill score < 3.5 = yellow; moderate–advanced, gill score ≥ 3.5 = red). The top right panel (G) shows a Sankey diagram of the bacterial families constituting over 1% of the overall relative abundance in the dataset where the thickness of the lines represents the relative abundance. The bottom right panels (H, I) Show beta‐diversity PCoAs using Bray–Curtis distance of samples grouped by amoebic load Ct value and gill damage, respectively.

(A) −log2‐transformed relative abundances of the bacterial families constituting over 1% of the total relative abundance across samples on a gradient from black (low) to dark yellow (high). The dendrograms represent the hierarchical clustering of samples and taxa abundances using Euclidean distance and complete linkage. The bottom legends show the samples based on groupings of amoebic load Ct value and gill damage respectively. (B) Correlation between the −log2 transformed relative abundance of Simkaniaceae and amoebic load Ct value. (C) Differential abundance of Simkaniaceae in samples grouped by gill damage. The colours represent the gill damage groups, Light–moderate (red) and moderate–advanced (yellow).

(A) Merged GWAS for AGD resistance for the phenotypes with suggestive SNPs. The coloured dots represent the SNPs crossing the suggestive threshold and the colours represent the phenotype. The horizontal red line represents the genome wide significance threshold and the blue line represent the suggestive threshold. (B) q–q plot of the observed/expected p‐values for the phenotypes with suggestive SNPs, the colours represent the phenotypes. (C) GWAS plot for the peak observed on chromosome 10, the colours represent the different phenotypes where the peak was observed.

Violin plots and a PCoA depicting the differences in phenotypes among genotypes on the suggestive SNPS for (A) Gill score on SNP AX‐87955988, (B) 2nd axis of a beta diversity PCoA on SNP AX‐88127064, (C, D) Alpha diversity (q = 0) on SNPs AX‐88304640 and AX‐87891546 respectively.

Tables

TABLE 1. Variance component estimates for AGD resistance for phenotypes where suggestive SNPs were observed.

Phenotype	σa2±se	σp2±se	σe2±se	h2±se
Gill score	0.07 ± 0.1	0.32 ± 0.06	0.26 ± 0.11	0.2 ± 0.32
Beta‐diversity axis.2	0.01 ± 0.01	0.04 ± 0.01	0.03 ± 0.01	0.17 ± 0.35
Alpha‐diversity(q = 0)	13.21 ± 31.44	90.99 ± 17.12	77.78 ± 33.09	0.15 ± 0.34

TABLE 2. SNP markers reported to be significant at the suggestive threshold in the GWAS.

SNP	AX‐87955988	AX‐88127064	AX‐88304640	AX‐87891546
Chr.	4	18	3	10
Position	15,608,124	69,701,288	78,377,066	12,153,141
Phenotype	Gill score	Beta diversity Axis.2	Alpha diversity (q = 0)	Alpha diversity (q = 0)
Freq.	0.336	0.188	0.155	0.289
SNP eff.	0.496	0.196	10.94	9.693
SE	0.119	0.048	2.712	2.279
p	2.94 × 10−5	3.80 × 10−5	5.47 × 10−5	2.11 × 10−5

Full text

Introduction

1

Salmonids are one of the most commercially important aquatic food products accounting for 20% of the total value of exported aquatic products in 2022 (FAO 2024). Diseases that occur in salmonid aquaculture are therefore a major threat to production, continued growth and sustainability of the industry. One important factor frequently mentioned in the context of salmonid aquaculture disease is gill health (Boerlage et al. 2020; Herrero et al. 2018; Mitchell and Rodger 2011; Oldham et al. 2016). The gills are not only the primary organ of respiration in salmonids, but also serve other important functions such as pH regulation, excretion of nitrogenous waste, hormone regulation and ionic regulation (Evans et al. 2005; Koppang et al. 2015). Being exposed to the surrounding water, the gills present one of the first lines of defence against invading pathogens and thereby also a potential point of entry for infectious agents (Cabillon and Lazado 2019; Koppang et al. 2015; Salinas 2015).

Amoebic gill disease (AGD) is one of the major threats to the gill health of marine salmonids and hence of major importance in salmon aquaculture (Oldham et al. 2016). Outbreaks of AGD can lead to severe economic losses and high mortality rates (Oldham et al. 2016) and in extreme scenarios mortality rates of up to 80% have been recorded in Atlantic salmon sea cages (Steinum et al. 2008). The primary causative agent of AGD is the amoeboid Neoparamoeba perurans (Crosbie et al. 2012). AGD pathology presents itself with white lesions, increased mucus production, epithelial hyperplasia and respiratory disruptions (Adams and Nowak 2003; Marcos‐López and Rodger 2020; Zilberg and Munday 2000). Moreover, AGD is often accompanied by other infectious agents resulting in a more complex disease phenotype often termed complex gill disease (CGD) (Boerlage et al. 2020; Herrero et al. 2018). AGD prevention and treatment options are scant, and currently limited to freshwater or hydrogen peroxide bathing (Marcos‐López and Rodger 2020; Rodger 2014). Therefore, other options such as selective breeding or genetic engineering to increase the resilience of Atlantic salmon stock to AGD may present additional prevention strategies.

A number of studies have addressed AGD resistance in Atlantic salmon using transcriptomic and quantitative genetic approaches. Robledo et al. (2018) identified two regions on chromosome ssa18 which indicated a suggestive association with AGD resistance, Aslam et al. (2020) identified SNPs on chromosomes ssa01, ssa02 and ssa05 which were significantly associated with AGD resistance and Boison et al. (2019) identified regions on chromosome ssa04, ssa09 and ssa13. Previous studies suggest that AGD resistance has a polygenic nature, with heritabilities estimated in the range of 0.09–0.48—typically characterised as moderately heritable—where the estimates differ based on the source population of Atlantic salmon and the number of re‐infection cycles (Aslam et al. 2020; Boison et al. 2019; Kube et al. 2012; Lillehammer et al. 2019; Robledo et al. 2018; Taylor et al. 2007, 2009). Studies in which transcriptomic analyses have been conducted have revealed the differential expression of genes involved in cell adhesion and proliferation, mucin secretion, cell cycle control, red blood cells, and immune response shedding further light on mechanisms of resistance (Boison et al. 2019; Botwright et al. 2021; Marcos‐López et al. 2018; Robledo et al. 2020; Wynne et al. 2008). Another layer may be added to the context of AGD resistance – that of the gill microbiota.

The microbiota is an often overlooked layer in the context of the host and has been exploited relatively little in the context of disease resistance. Recently Schaal et al. (2022) showed a potential link between host genetics, metabolic rate and gut microbiota in AGD, demonstrating that including the layer of the microbiota is a potentially valuable asset in terms of understanding disease resistance. Several studies have addressed the microbial communities of Atlantic salmon gills in the context of AGD (Birlanga et al. 2022; Botwright et al. 2021; Bowman and Nowak 2004; Clinton et al. 2024; Downes et al. 2018; Slinger et al. 2020, 2021). Collectively these studies indicate that the gill microbiota is generally dynamic – not unexpected given its close proximity to the surrounding water – but altered during an AGD outbreak. Reported differences include the presence or altered relative abundance of known pathogenic taxa and differences in diversity. It has further been suggested that some bacteria (Winogradskyella sp. in this case) may be able to enhance the severity of an AGD infection (Embar‐Gopinath et al. 2005, 2006). Moreover bacteria have been reported to live in close association or as endosymbionts of the amoeba itself (Benedicenti et al. 2019; MacPhail et al. 2021; Nylund et al. 2018), and even differences in microbial composition among virulent and attenuated amoeba cultures have been described (Ní Dhufaigh, Botwright, et al. 2021; Ní Dhufaigh, Dillon, et al. 2021). Taken together, the body of work regarding AGD and gill microbiota suggests that while AGD had been associated with alterations of the microbiota, the causal role and functional significance of these changes remain unresolved. Thus, these complex dynamics of the host, parasite and microbiota are challenging to address, especially in the context of host genetics and disease resistance.

Looking at this complex system in a more holistic way – namely the host and the associated microbiota in the same framework – may have potential in shedding further light on mechanisms of AGD resistance. Holo‐omics which integrates combined ‐omics layers from the host and its associated microbiota (Nyholm et al. 2020; Odriozola et al. 2024) provides a unique framework for investigating host–microbiota interactions in the holobiont – the host and its microbiome together (Baedke et al. 2020; Bordenstein and Theis 2015; Zilber‐Rosenberg and Rosenberg 2008). Holo‐omics is promising for offering new insights into ecology and evolution (Alberdi et al. 2022; Foster et al. 2017; Nyholm et al. 2020; Theis et al. 2016) and has great potential to be applied in industries such as aquaculture (Limborg et al. 2018).

We intended to investigate whether we could leverage a holo‐omic framework to shed more light on AGD resistance/susceptibility. Here, we utilise a combined approach of 16S metabarcoding and quantitative genetics to address AGD resistance. We characterise the gill microbiota of Atlantic salmon subjected to an AGD disease challenge and investigate microbiota composition in relation to traditional AGD phenotype scoring methods (gill score and amoebic load). We then performed a genome‐wide association study (GWAS) to investigate SNPs associated with the phenotypic scoring methods and the microbiota as two distinct, but interlinked, phenotypes.

Methods and Materials

2

The AGD challenge data used in this study has been published and described by Robledo et al. (2020). Briefly, a cohabitation AGD challenge was performed using 797 Atlantic post‐smolt salmon (~18 months, mean weight after challenge 463.5 g) from 132 nuclear families from a commercial breeding programme (Landcatch, UK). Neoparamoeba perurans was sourced from naturally infected farmed salmon in the West of Scotland. Fish with a similar level of AGD infection (assessed by gill damage) were used as seeder fish at a ratio of 15% seeder to naïve fish. The pathogen was maintained exclusively in vivo and was not cultured in the laboratory prior to the challenge to mimic realistic production settings. The experimental challenge was conducted at the experimental facilities of University of Stirling's Marine Environmental Research Laboratory, Machrihanish (Scotland, UK) in a 4 m^3^ seawater tank (water temperature between 13°C and 14°C, at 33–35 ppt of salinity), supplied with ambient flow‐through seawater filtered to approximately 90 μm for the duration of the trial. The challenge consisted of three separate cycles of infection with a recovery period after two infections, based on Taylor et al. (2009), in which the gill scores for the third infection round were shown to best predict survival outcomes. A freshwater treatment was applied 21 days after the start of the challenge followed by a week of recovery and by the addition of seeder fish. During the challenge the fish were checked visually four times daily and in the third cycle of infection the disease was allowed to progress until the sampling point.

At sampling fish were terminated by an overdose of anaesthetic (Phenoxyethanol, 0.5 mg/L) and gill damage was recorded for both left and right gills, using all gill arches, by a single operator. Gill damage, as used throughout this manuscript, refers to gross, macroscopic scoring of gill lesions and does not imply histological assessment. Gill damage was scored from 0 for clear gills (healthy red gills, no gross sign of infection) to 5 for heavy AGD infection level (extensive lesions covering most of the gill surface) according to Taylor et al. (2016), the mean score of the two gills was used as the resistance phenotype. The second left anterior gill was dissected out and stored in ethanol for amoebic load estimation by qPCR using N. perurans specific primers and 16S bacterial profiling.

V3‐V4 16S rRNA Bacterial Profiling

2.1

Out of 797 a total of 77 samples were available for 16S bacterial profiling. This subset was determined by availability and sequencing costs, but covers the full spectrum of gill scores, providing representative biological coverage. Equal sizes of gill tissue were cut off each sampled salmon and transferred to a lysing E‐matrix tube containing 0.9 mL of DNA/RNA shield (Zymo Research). Prior to DNA extraction samples were lysed for 5 min in a Tissuelyzer (Qiagen) at 30 GHz after which they were spun down at 16.000 g. In a randomised order, 400 μL of lysate from samples was transferred to a 96 deep‐well plate. DNA extraction was then performed using the Quick‐DNA MagBead Plus kit (Zymo Research) following the recommendations of the manufacturer. Negative extraction controls containing 400 μL DNA/RNA were included in the process. The DNA concentration was measured using a Qubit fluorometer (Invitrogen). Prior to amplification a qPCR was performed on a subset of the samples to check for inhibitor presence and optimal cycle number. Initial amplification was performed using the forward and reverse primers 341F (5′‐CCTAYGGGRBGCASCAG‐3′) and 806R (5′‐GGACTACNNGGGTATCTAAT‐3′; Yu et al. 2005) with Illumina adapter overhangs, targeting the V3‐V4 region of the bacterial 16SrRNA. Each amplification reaction consisted of 7.5 μL ddH_2_O, 9.5 μL AccuPrime SuperMix II (Invitrogen), 1 μL 10 μM forward 16S primer, 1 μL 10 μM 16S reverse 16S primer, and 1 μL sample DNA, totaling 20 μL. The amplification PCR settings were as follows: incubation at 95°C for 10 min, followed by 30 cycles of denaturation at 95°C for 15 s, annealing at 53°C for 20 s and extension at 68°C for 40 s; after the completed cycles this was followed by a 10 min final extension at 68°C. The Nextera XT index kit v2 Set A (Illumina) was used to uniquely index the PCR products. The indexed PCR products were then visualised on a 1% agarose gel and subsequently pooled in an equimolar fashion based on band strength and purified using SPRI beads. Libraries were quantified using the Agilent BioAnalyzer 2.100. Paired‐end sequencing was performed using the Illumina MiSeq platform, reagent kit v3 at 600 cycles by the GeoGenetics Sequencing Core, University of Copenhagen, Globe Institute. Negative extraction controls and PCR controls were included in the library construction and sequenced for downstream quality control.

Bioinformatic Processing and Analysis of 16S Amplicon Sequencing Data

2.2

Approximately 12.9 million read pairs were generated for the 83 samples and controls. Initial quality checks were performed using FastQC (Andrews 2010) and MultiQC (Ewels et al. 2016). Adapter and initial quality trimming were subsequently performed using TrimGalore (Krueger et al. 2021). The DADA2 pipeline (Callahan et al. 2016) was used for further quality filtering, trimming, ASV inference, and removal of chimeric sequences. Taxonomy was assigned in the DADA2 framework using the Silva/v138 database training set (Quast et al. 2013). LULU (Frøslev et al. 2017) was then applied for post‐clustering curation of the ASVs to minimise errors after which Decontam (Davis et al. 2018) was applied to remove contaminant ASVs from the dataset (Table S3 lists the contaminant taxa). The ASVs were further processed using the R package Phyloseq (McMurdie and Holmes 2013). ASVs assigned to Eukaryota, Mitochondria or Chloroplast were removed and ASVs were then agglomerated on the family level based on lack of resolution on the genus level and to minimise noise. ASVs were normalised using simple proportions. A Hill numbers framework using the R package hilldiv (Alberdi and Gilbert 2019) was used to carry out alpha diversity analyses, applying different orders of diversity (q), where q = 0 treats all taxa equally, q = 1 gives proportional weight to taxa, and q = 2 emphasises dominant taxa. Differential abundance analysis was performed using the Wilcoxon rank‐sum test with an FDR adjusted p‐value. Beta‐diversity estimations were performed using a PCoA with Bray–Curtis distance and statistical differences were estimated using PERMANOVA with Bray‐Curtis distance from the R package Vegan (Oksanen et al. 2020). Statistical analyses were performed in R/4.3.1 (R Core Team 2023). Additional packages used for visualisation were ggplot2 (Wickham 2016), ggpubr (Kassambara 2020), and patchwork (Pedersen 2024).

Genotyping and Estimation of Genetic Parameters

2.3

The fish were genotyped using a SNP array with 47K SNPs (Houston et al. 2014), after DNA extraction from fin clip tissue samples using the DNeasy 96 tissue DNA extraction kit (Qiagen, UK) as described in Robledo et al. (2018). Standard quality control was performed using PLINK v1.90b6.21 (Purcell et al. 2007) where individuals with less than 95% of their SNPs genotyped were removed, SNPs with a minor allele frequency (MAF) less than 0.05, a call rate of less than 95% and SNPs that deviated from the Hardy–Weinberg equilibrium (p < 10^−6^) were removed. This resulted in 26,459 SNPs for 58 individuals. Next the software GCTA (Yang et al. 2011) was used to construct a genomic relationship matrix (GRM) with:

$$[eqn]$$

where gjk denotes the relationship between individual j and k. N refers to the total number of SNPs (26,459), the copy number of reference allele for a ith SNP for the jth and kth fish is denoted by xij and xik and pi is the reference allele frequency.

The Average Information Restricted Maximum Likelihood (AI‐REML) (Gilmour et al. 1995) algorithm applied in GCTA was then used for the estimation of genetic parameters based on the GRM using the following linear mixed model approach:

$$[eqn]$$

where yi denotes a phenotype for the ith fish, μ is the mean in the population, ui is the random additive genetic value for individual i, following a normal distribution uN0Gσg2, where G is the GRM made with GCTA (Equation 1) and σg2 the estimated genetic variance. ei is the residual effect following a normal independent distribution eN0Iσe2 where σe2 denotes the residual variance.

A GWAS was then performed in GCTA using a mixed linear model association with the leave‐one‐chromosome‐out (loco) approach.

$$[eqn]$$

where yi denotes a phenotype for the ith fish, μ is the mean in the population, for the jth SNP aj denotes the additive genetic effect of the reference allele with the genotype for individual i (gij) represented with 0,1 or 2. ei is the residual effect following a normal independent distribution eN0Iσe2 where σe2 denotes the residual variance. ui is then a random vector of polygenic effects following a normal distribution uN0Gσg2, where G (because of the mlma‐loco approach) is a partial GRM constructed with 28 chromosomes after leaving out the chromosome containing the jth SNP, and σg2 is the estimated genetic variance. A Bonferroni correction (α = 0.05) was used to determine the genome wide significance threshold (−log10α/n) where n is the number of SNPs, and the chromosome‐wide significance threshold (−log10α/n/29). The phenotypes used were both the traditional AGD scoring methods (amoebic load and gill score) and traits of the microbiota which were; the unweighted and weighted alpha diversity measures, the first two axis of a beta diversity PCoA using Bray–Curtis distance, and inverse normally transformed (INT) relative abundances of families that contributed over 3% of the total relative abundance in the dataset.

Results

3

Sequencing and Sample Overview

3.1

A 16S amplicon library was sequenced consisting of DNA extracted from 77 gill samples including negative extraction and PCR controls and a positive control totalling 12.9 million read pairs related to 83 samples. After quality filtering and trimming of sequences, ASV inference, taxonomic identification, removal of contaminant sequences and post‐clustering curation, we were left with a dataset comprised of 844 ASVs. 35% of the ASVs were not assigned at the genus level; hence the ASVs were merged at the family level and analysed accordingly. After filtering out samples with a low ASV count and samples not assigned both Ct value for amoebic load and a gill score we were left with 59 samples containing a total of 180 bacterial families. Most of the samples had an assigned gill score around 3–4 while amoebic load Ct values were more evenly distributed among the 59 samples. For analysis of microbiota composition, we grouped the samples by gill score into light–moderate (gill score < 3.5) and moderate–advanced (gill score ⋝ 3.5) by degree of gill damage as described in (Taylor et al. 2016) and by amoebic load into Ct < 30, 30 < Ct < 35, Ct > 35 to appropriately cover the range of amoebic load Ct values.

Microbiota Profiles Vary by AGD Scoring Methods

3.2

We first investigated the general composition of the microbiota on the family level, since 35.5% of the total abundance of ASVs was unassigned on the genus level. After filtering, we found that 61.44% of the overall relative abundance is attributable to two families namely, Simkaniaceae (31.3%) and Arcobacteraceae (29.79%). After Simkaniaceae and Arcobacteraceae the families Vibrionaceae, Marinomonadaceae, Pseudoalteromonadaceae, Rhodobacteraceae, Flavobacteriaceae each contributed over 3% of the overall relative abundance. The gill microbiota of the AGD challenged salmon in this study can therefore be described as dominated by two Families namely Simkaniaceae and Arcobacteraceae and a few families contributing between 1% and 3% of relative abundance (Figure 1G).

Left panel shows the comparison of alpha diversity estimates based on hill numbers at orders of diversity (q = 0, 1, 2) which apply increasing weight to the relative abundance of taxa, among samples grouped by (A–C) amoebic load Ct value where the colours represent Ct < 30 (grey), 30 < Ct > 35 (pink) and Ct > 35 (purple), and (D–F) gill damage (light–moderate, gill score < 3.5 = yellow; moderate–advanced, gill score ≥ 3.5 = red). The top right panel (G) shows a Sankey diagram of the bacterial families constituting over 1% of the overall relative abundance in the dataset where the thickness of the lines represents the relative abundance. The bottom right panels (H, I) Show beta‐diversity PCoAs using Bray–Curtis distance of samples grouped by amoebic load Ct value and gill damage, respectively.

We then looked at the diversity among the samples grouped by amoebic load and gill score. Using alpha diversity metrics with hill numbers and orders of diversity (Alberdi and Gilbert 2019) we found significant differences in diversity between samples grouped by amoebic load on the Family level (Figure 1A–C). These significant differences were observed for the two orders of diversity q = 1 (p‐adjusted = 0.024) and q = 2 (p‐adjusted = 0.017), where more weight is attributed to relative abundance. No significant difference was found on the ASV level (Figure S1). When we grouped samples according to the severity of gill damage no significant difference was seen in alpha diversity on the family level (Figure 1D–F) or ASV level (Figure S1). Overall, the weighted microbial diversity measures (q = 1 and q = 2) trend to increase with a lower amoebic load (Figure 1A–C), indicating that microbial communities in fish with a lower amoebic load may be dominated by a broader set of common taxa, rather than a few highly abundant ones Moreover, the alpha‐diversity measures indicate that the samples were largely dominated by a few bacterial families confirming previous observations.

The beta‐diversity did not differ among samples grouped by amoebic load (Figure 1H,I) but differed slightly among samples grouped by gill damage (PERMANOVA, R ^2^ = 0.057, p = 0.003). Considering the little difference in ordination space and a significant result in a test for homogeneity of dispersion (p = 0.006), the effect observed may be due to differences in group dispersions.

Simkaniaceae Bacterium Correlates with Disease Status

3.3

Clustering of taxa abundances across samples revealed no general patterns in terms of amoebic load or gill damage (Figures 2A and S2), but microbiota‐specific correlations. Differential abundance analysis of the bacterial families constituting more than 1% of the overall relative abundance shows that one bacterial family, Simkaniaceae, is differentially abundant (p.adjust = 0.037) among samples grouped by gill damage (Figure 2C). No differentially abundant bacterial families were detected in samples grouped by amoebic load. Notably the representative ASV of the Simkaniaceae has a 99.78% percentage identity (using BLAST) with Candidatus Syngnamydia Salmonis, which may be an intracellular bacterium of the amoeba as it has been found in co‐culture with the amoeba and identified in AGD before (Nylund et al. 2018, 2015). Interestingly there is a significant correlation (R ^2^ = 0.29, p = 0.024) between the amoebic load Ct value and the relative abundance of Simkaniaceae (Figure 2B) suggesting a co‐occurrence pattern with this bacterium and the amoeba.

(A) −log2‐transformed relative abundances of the bacterial families constituting over 1% of the total relative abundance across samples on a gradient from black (low) to dark yellow (high). The dendrograms represent the hierarchical clustering of samples and taxa abundances using Euclidean distance and complete linkage. The bottom legends show the samples based on groupings of amoebic load Ct value and gill damage respectively. (B) Correlation between the −log2 transformed relative abundance of Simkaniaceae and amoebic load Ct value. (C) Differential abundance of Simkaniaceae in samples grouped by gill damage. The colours represent the gill damage groups, Light–moderate (red) and moderate–advanced (yellow).

GWAS Identifies Suggestive Association With Gill Score and Bacterial Traits

3.4

We then performed a GWAS using both the traditional AGD scoring methods (amoebic load and gill score) as a phenotype together with chosen phenotypic metrics from the microbiota (see Methods). One sample was discarded after filtering leaving 58 samples for the GWAS.

No meaningful inference could be made from heritability estimates due to large standard errors, likely due to the modest sample size (Table 1 for the phenotypes with suggestive SNPs, and Table S1 and Figure S3 for all phenotypes tested). This is also reflected in the expected/observed p‐value q–q plots where there is a slight deflation towards the tail of the q–q plots (Figure 3B). While no genome‐wide or study‐wide significant SNPs were found for any of the traits tested, we did detect some outlier SNPs at the suggestive significance threshold (Table 2).

(A) Merged GWAS for AGD resistance for the phenotypes with suggestive SNPs. The coloured dots represent the SNPs crossing the suggestive threshold and the colours represent the phenotype. The horizontal red line represents the genome wide significance threshold and the blue line represent the suggestive threshold. (B) q–q plot of the observed/expected p‐values for the phenotypes with suggestive SNPs, the colours represent the phenotypes. (C) GWAS plot for the peak observed on chromosome 10, the colours represent the different phenotypes where the peak was observed.

A single suggestive SNP (AX‐87955988) was found on chromosome 4 when gill score was used as a phenotype (Figure 3A). Looking at the genotypes of that suggestive SNP we can see that individuals with the AA genotype tend to have a lower gill score compared to the individuals carrying the B allele (Figure 4A). Applying the microbiota diversity measures as phenotypes revealed a suggestive SNP (AX‐88127064) on chromosome 18 for the 2nd axis of the beta diversity PCoA (Figure 3A). No obvious differences were observed on the 2nd axis on the PCoA other than individuals with the AA genotype seeming to cluster more towards the lower end of the 2nd axis compared with individuals with the B allele (Figure 4B).

Violin plots and a PCoA depicting the differences in phenotypes among genotypes on the suggestive SNPS for (A) Gill score on SNP AX‐87955988, (B) 2nd axis of a beta diversity PCoA on SNP AX‐88127064, (C, D) Alpha diversity (q = 0) on SNPs AX‐88304640 and AX‐87891546 respectively.

Two SNPs (AX‐88304640 and AX‐87891546) on chromosome 3 and 10 were significant at the suggestive threshold for richness estimates (Figure 3A). In both instances individuals with the BB genotype have lower diversity (at q = 0) than do individuals carrying the A allele (Figure 4C,D).

No other suggestive SNPs were found for the relative abundances of the other bacterial families tested. Interestingly, a peak is present on chromosome 10 for both the 2nd axis of the PCoA and the richness estimates (Figure 3C). Furthermore, when looking at the INT transformed relative abundances as phenotypes (Figure 4C) the peak on chromosome 10 is also present for the abundance of the family Arcobacteraceae although no suggestive SNP is present (top SNP, AX‐88038371 at position 15,555,246 bp). No other suggestive SNPs were present.

A published list of differentially expressed genes in gill and head kidney between AGD‐susceptible and more resistant Atlantic salmon (Robledo et al. 2020) from the same population was used to guide candidate gene exploration. Looking in a 2 Mb range of the suggestive SNPs and the peak on chromosome 10 where the peak is present, we identified several genes (Table S2).

First, the gene cullin‐4b, a part of the cullin family of proteins, essential for protein degradation through ubiquitination (Fan et al. 2022; Sarikas et al. 2011) was found on chromosome 4 close to the SNP associated with mean gill score and reported to be more expressed in the head–kidney of susceptible fish (Robledo et al. 2020). Second, the gene slc4a2 encodes an anion exchange protein involved in the exchange of CL‐ and HCO3‐ across cell membranes and likely contributes to pH regulation in fish (Esbaugh et al. 2012; Romero et al. 2004; Shmukler et al. 2008) is located in proximity to the SNP on chromosome 4 associated with alpha diversity (q = 0) and reported to be more expressed in resistant fish gills (Robledo et al. 2020). Third, slc34a2 encodes a phosphate carrier protein involved in cellular homeostasis of inorganic phosphate (Hilfiker et al. 1998; Werner et al. 2016) and was found in proximity to the peak on chromosome 10 and was reported to have higher expression in the head‐kidney of AGD susceptible fish (Robledo et al. 2020). Finally, eps15 coding for epidermal growth factor substrate 15, is involved in the internalisation of the epidermal growth factor receptor (EGFR) and secretion and endocytosis in general (van Henegouwen and Paul 2009), was found close to the peak on chromosome 10 and was reported to be more expressed in the gills of susceptible fish. Literature from murine models and humans suggests that EGFR is involved in processes surrounding key indicators of AGD including increased mucus production and hyperplasia (Burgel and Nadel 2004; Takeyama et al. 2008).

Discussion

4

In this study we applied a holo‐omic approach to investigate AGD resistance in Atlantic salmon, integrating host genetic variation, gill microbiota composition, and potential amoebic symbionts. We described the composition of the gill microbiota of AGD‐affected salmon and characterised the differences in microbiota composition based on the severity of the disease. We then conducted a GWAS with both traditional resistance indicator traits and traits of the gill microbiota composition. Here we discuss the perspectives of the results and address the limitations of the study, with an eye towards inspiring future implementation of microbiota data in similar disease resistance studies.

Influence of Environmental Factors, Sampling Procedure and Simkaniaceae Presence

4.1

The Atlantic salmon gill microbiota has been suggested to be altered in response to AGD (Birlanga et al. 2022; Clinton et al. 2024; Slinger et al. 2020). While our observations are in agreement with microbiota differences in relation to AGD, the results differ somewhat from the microbiota composition described in other AGD‐related gill microbiota studies. This difference may be due to more transient environmental factors such as differences in water temperature, pH, salinity, geographic location or aquaculture environment (Lokesh and Kiron 2016; Sylvain et al. 2016; Minich et al. 2020; Amill et al. 2024). However, the differences may also be due to the procedure used to sample the gill microbiota.

In this study we used gill biopsies rather than relying on alternative strategies such as gill mucus scrapes or mucosal swabs. These strategies have been reported to be reliable alternatives for assessing gill microbiota in salmonids, but are also reported to have compositional differences when compared to each other (Birlanga et al. 2022; Clinton et al. 2021; Clokie et al. 2022). Although gill biopsies are undoubtedly an invasive sampling technique it may be appropriate in certain cases, especially in the case of endosymbiotic bacteria which may reside and propagate in the gill tissues or are derived directly from parasites and that would be missed from only sampling a mucosal swab. Perhaps including a combination of mucosal swabs along with gill biopsies will give a better representation of the entirety of the gill microbiota.

The most striking difference from other studies was the high relative abundance of Simkaniaceae. Candidatus Syngnamydia salmonis, a member of the Simkaniaceae has been found in association with Neoparamoeba perurans, the AGD‐causing parasite, and has even been grown in co‐culture with the parasite (Nylund et al. 2018, 2015). Moreover, Simkaniaceae have been reported in other AGD‐related studies (Birlanga et al. 2022; Clinton et al. 2024), although in lower abundances. Although the representative ASV of the Simkaniaceae in this study shows a high percentage identity with Candidatus Syngnamydia salmonis, we cannot conclusively confirm its identity. These observations support the previous notion that our results may represent both the internal and external parts of the Atlantic salmon gill microbiota along with a potential endosymbiont of the amoebic parasite itself. As suggested in Nylund et al. (2018), Simkaniaceae is likely not universal to all N. perurans and may thus be opportunists in this regard. The observed correlation between Simkaniaceae and amoebic load in this study suggests a potential relationship. However, whether the presence/absence or abundance of this Simkaniaceae has an influence on the virulence of the amoeba or disease progression of AGD – as has been shown with another bacterium (Embar‐Gopinath et al. 2005, 2006) – warrants further investigation.

Due to the quantitative study design, a control group of unaffected salmon was not included. While the absence of a control group does not affect the validity of the analyses, a baseline group could have strengthened certain interpretations. For example, a baseline group would likely give a better clue as to whether Simkaniaceae are truly derived from the N. perurans – although this is very likely given previous studies – and clarify whether the microbiota of AGD‐affected salmon is more perturbed compared to an uninfected state. Although the subset of fish analysed primarily exhibited moderate to high gill scores, the data still encompass a spectrum of microbial traits across individuals, allowing for comparisons and analysis of variation in relation to infection severity. Future work incorporating temporal or longitudinal sampling as AGD progresses would provide additional insight into relationships between microbial dynamics and infection severity.

Microbiota Diversity in Relation to AGD Scoring Methods

4.2

The results of the microbiota analysis differed based on the AGD scoring methods. Specifically, one diversity metric differed among samples grouped by amoebic load but not by gill damage and vice versa. This might indicate that the typical resistance phenotypes usually used for AGD may not fully capture the phenotypic nature of the disease and argues for including other layers like that of the microbiota in elucidating resistance mechanisms.

The differences observed among samples grouped by amoebic load in alpha diversity measures indicated that the diversity is reduced in samples with higher amoebic load. This is likely due to certain higher abundance taxa since the differences observed are clearer in the diversity metrics where the abundance of taxa is taken into account. This is supported by the observation that the bacteria belonging to the family Simkaniaceae– the most abundant bacteria in our dataset – correlated significantly with amoebic load value (Figure 3B).

Integrating Microbiota Traits in GWAS for AGD Resistance

4.3

A goal of this study was to use a holo‐omic approach by conducting a GWAS using both the traditional AGD scoring methods and to include traits of the microbiota composition as resistance‐informative phenotypes. It should be noted that the GWAS is underpowered due to the small sample size, which was reflected in the q–q plots of expected versus observed p‐values (Figure 4B). Due to unreliable heritability estimates we decided against estimating the percentage of explained variability of each suggestive SNP as the estimates would likely be under‐/over‐estimated due to the high SE of heritability and lead to spurious results. Despite a modest sample size and high standard error of heritabilities, our results conform with previous studies on AGD resistance as the heritability estimate of the gill score resistance trait was 0.20.

Interpreting microbiota traits in the context of a GWAS on a disease‐challenged population is complex. While changes in the microbiota between healthy Atlantic salmon gills and AGD‐infected Atlantic salmon gills are documented (Birlanga et al. 2022; Botwright et al. 2021; Clinton et al. 2024), it is unclear what constitutes a desirable or healthy gill microbiota in the context of disease resistance. Is it high or low microbial diversity, a specific microbiota community composition, or the presence of a few specific microbes? One possibility is that genetically conveyed resistance allows for the microbiota to stay in a more robust state compared to more susceptible genotypes, suggesting that host genetics influence microbiota composition in a way that enhances resistance. Thus, including traits of the microbiota as covariates with the more conventional resistance traits may add more holistic insights into host genetic‐microbiota‐parasite interactions compared to traditional GWAS approaches only linking host genotype and phenotype (see e.g., Brealey et al. (2022)). While we do not directly assess microbiota robustness, future studies could use longitudinal sampling during AGD infection to examine how host genetics shape temporal stability and shifts in microbiota composition, providing insight into host‐linked microbial resilience.

Opportunities and Challenges in Microbiota‐Based GWAS

4.4

When utilising the relative abundance of individual bacterial families – or any other taxonomic rank – in disease resistance‐based GWAS analyses, having a hypothesis‐driven rationale for testing a particular trait strengthens the context of analysis. Although exploratory data analysis can also provide valuable insights– as this study shows. For instance, in the case of the relative abundance of Arcobacteraceae, the second most abundant bacterial family, a suggestive SNP was identified on chromosome ssa10. This raises the question: is this genetic effect associated with AGD resistance, or does it reflect a different heritable bacterial population trait influenced by host genetics? Conversely, for Simkaniaceae, while no suggestive SNP was observed, it may be reasonable to hypothesise that its relative abundance could serve as a reasonably robust resistance indicator trait, given its documented presence in co‐culture with the AGD‐causing parasite (Nylund et al. 2018). However, the broader ecological context of Simkaniaceae, within the gill environment, in relation to its association with amoebae, and its global prevalence remains poorly understood, and its potential role therefore remains speculative. Nevertheless, our study presents a path for future work to clarify the function and relevance of gill microbes in the context of disease resistance. It has been noted that GWAS studies applying 16S microbiota data may be problematic to some extent, due to the hierarchical nature of taxonomic data and the relative abundance of microbial taxa (Awany et al. 2018; Bruijning et al. 2023). That is because the relative abundances of taxa can be represented across multiple levels (Phylum, Family, Genus, Species, etc.) and thus each host SNP may affect these levels differently. The interrelatedness across taxonomic levels may therefore obfuscate the potential to directly suggest resistance indicators from a taxon point of view. However, assuming that the composition of the microbiota is influenced by host genetic factors (Alberdi et al. 2025) we argue that using microbiota traits as a health indicator in a quantitative genetics framework better captures genetic influence on the microbiota that directly contributes to resistance mechanisms. Including the microbiota for scoring resistance phenotypes introduces a more holistic and integrated view of resistance mechanisms.

Towards Functional Insights

4.5

In a holo‐omic context the compositional taxonomy‐centric nature of 16S metabarcoding can be a limiting factor. While a 16S metabarcoding approach provides a birds‐eye view over microbiota dynamics, the approach fails to capture information regarding the gene content, functions and activity of the microbiota. Although methods have been developed to infer functional insights based on 16S metabarcoding data (Aßhauer et al. 2015; Douglas et al. 2020; Langille et al. 2013), these methods are constrained in their ability to accurately infer functional capacity (Sevigny et al. 2019), particularly in animal and environmental studies (Sun et al. 2020) or when considering the taxonomic resolution of the 16S amplicon data (Alberdi et al. 2022; Antony‐Babu et al. 2017; Welch et al. 2002). The limitations of marker‐gene based prediction methods are overcome by integrating long‐read metagenomic sequencing and metatranscriptomics layers. These layers include whole bacterial genomes, associated functional pathways, microbial genetic variants, metabolism predictions and gene expression, all of which may be applied as phenotypes in a GWAS (Doolittle and Booth 2017; Sanna et al. 2022). Thus, including these ‐omic layers, although costly, may give a clearer picture of the functional dynamics involved in coordinating host‐microbiota responses such as the gill microbiota of AGD‐affected salmon and the mechanisms underlying AGD resistance/susceptibility.

Highlighting the added value of including layers from both the host and microbiota, we used a set of differentially expressed genes from the same population (Robledo et al. 2020), based on their location in proximity to the suggestive SNPs. This host‐transcriptomic layer can assist in elucidating mechanisms of resistance and identifying genes of interest in a putative QTL identified by considering the layer of the microbiota as an extended host phenotype. This revealed interesting candidate genes involved in cellular processes and homeostasis that have not been documented before in relation to AGD (Boison et al. 2019; Robledo et al. 2020). However, given the limitations of this study regarding sample size and the inability to compare differences in gene expression with regard to SNP allelic frequencies, these results should be taken with caution.

Impact and Prospects for Future Research

4.6

The limited knowledge surrounding gill microbiota dynamics and its role in disease as well as the complex nature of AGD and its resistance mechanisms, highlights the need for novel approaches. The idea of incorporating the microbiota to study AGD resistance mechanisms through a holo‐omic lens simultaneously poses interesting challenges and questions, which our findings begin to address. The results of the study indicate differences in the microbiota traits based on traditional AGD scoring methods, showcasing the potential of using traits of the microbiota to study AGD disease resistance mechanisms. The techniques applied here can be of value in further understanding AGD resistance or susceptibility. However further larger‐scale studies are needed to corroborate the results and draw more confirmed conclusions about the underlying mechanisms governing concerted host‐microbiota disease responses.

Author Contributions

Eiríkur Andri Thormar: conceptualization (equal), data curation (lead), formal analysis (lead), investigation (lead), methodology (equal), visualization (lead), writing – original draft (lead), writing – review and editing (lead). Clémence Fraslin: data curation (equal), formal analysis (equal), investigation (equal), methodology (equal), supervision (equal), visualization (equal), writing – original draft (equal), writing – review and editing (equal). Morten T. Limborg: conceptualization (equal), data curation (equal), funding acquisition (equal), resources (equal), supervision (equal), writing – original draft (equal), writing – review and editing (equal). Diego Robledo: conceptualization (equal), data curation (equal), funding acquisition (equal), project administration (lead), resources (equal), supervision (lead), writing – review and editing (equal).

Ethics Statement

All animals were reared in accordance with relevant national and EU legislation concerning health and welfare. The challenge experiment was performed by the Marine Environmental Research Laboratory (Machrihanish, UK) under the approval of the ethics review committee of the University of Stirling (Stirling, UK) and according to Home Office licence requirements. Landcatch is an accredited participant in the RSPCA Freedom Foods standard, the Scottish Salmon Producers Organisation Code of Good Practice, and the EU Code‐EFABAR Code of Good Practice for Farm Animal Breeding and Reproduction Organisations.

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Data S1: ece372484‐sup‐0001‐Supinfo.docx.