“Seeder fish” were recovered from a field outbreak of CGD on the west coast mainland of Scotland on 15^th July 2022. These infected animals were bred from a stock of naive Atlantic salmon smolts (Benchmark Genetics Iceland AKA StofnFiskur, strain) and transported to an isolation unit at the Marine Environmental Research Laboratory (MERL), Campbeltown Scotland which is part of the Institute of Aquaculture (IOA) and University of Stirling (UOS). On the 27^th September 2022, a subsample of this population (n=9 individuals) were sampled by technicians and diagnostic salmon gill material sent to Patogen AS (www.patogen.no). All samples were subjected to extraction of total RNA followed by an automated one-step standardized Real Time Polymerase Chain Reaction (RT-PCR) protocol using a relevant validated methodology (Norwegian Standard NS-EN ISO/IEC 17025). The samples (n = 9 individuals) were all confirmed positive for the following pathogens, Paramoeba perurans; Candidatus Branchiomonas cisticola, Paranucleospora theridion and Salmon Gill Poxvirus (SGPV). The remaining isolated animals (n=15 individuals) were used to propagate disease by sharing contaminated water from infected fish equally into two challenge tanks (n= 2 tanks) housing the Benchmark Icelandic year class 2021–2 family nucleus (naïve fish). The gill challenge was propagated, and gill health was meticulously monitored by the technicians under UK Home Office license. A total of two separate CGD challenges were conducted, with therapeutic freshwater treatment (4 hours duration) employed to separate them. Terminal samples were taken before any CGD-induced mortalities occurred after the second challenge. These samples included gill scores, fin clips for parental assignment and genotyping, fish body mass (grams) and maturation status.

The Atlantic salmon population selected for an CGD challenge test was derived from the 2021/2-year class and 13^th generation of Benchmark genetics Atlantic salmon breeding program in Iceland. A total of 195 full-sibling families were produced from a unique male-female pairing whereby, 190 families were produced using 1:2 paring (1 male: 2 female) and 5 families with 1:1 paring (1 male: 1 female), resulting in 195 full-sibling families and thus 95 paternal half-sib families. The families were selected from isolated family tanks and pooled into a single population after first feeding at Kollafjörður in Iceland and then transported by air freight. On May 5^th, 2022, a total of 3120 individual salmon parr (16 individuals per family) were shipped to the University of Stirling, (mean weight 26.0 g) and were kept in freshwater from May 22 to August 22, 2022, until the time of vaccination (Alpha Jet Micro 6, Pharmaq). On June 28, 2022, it was recorded that 0.7% of the fish had died, leaving 3097 fish weighing 97g. A total of 2952 smolts (approximately 15 individuals per family) were transferred to seawater on August 23, 2022, and were randomly placed into 2 x 14 m3 fiber glass seawater tanks at the Marine Environmental Research Laboratory at the University of Stirling in Scotland.

To detect the presence of the disease, and to mimic the multiple CGD infections which Atlantic salmon are exposed to in the commercial production environment (ocean net pens) two independent disease challenges were undertaken separated by a freshwater treatment. The first phase of the CGD challenge test began on September 14, 2022, in tanks with a salinity of 34 ppt and a temperature of 15 °C. In this phase, 2060 individuals (approximately 10 individuals per family) were challenged with CGD. Gill disease symptoms were observed in all 40 sampled fish 13 days post-infection, and the average CGD gill score was 1.4 (range 0.6 - 2.4). After completing the first CGD challenge, a freshwater treatment was conducted for four hours on September 30, 2022. Thereafter, a second CGD challenge began on October 1, 2022, with a total of 1946 fish remaining alive from the first challenge (9–10 individuals per family) entered a second CGD challenge. Ten days post-infection of the second CGD challenge, disease symptoms were observed in all 24 randomly sampled fish. The average gill score was 1.5 (range 0.6 - 2.6). A terminal sample was collected on October 25, 2022, whereby all remaining live fish were humanely dispatched post-CGD infection and sampled for gill score, DNA fin clip, body weight (grams) and sexual maturation recordings.

Upon completion of the second phase of the CGD challenge test, fish survivors underwent a meticulous recording process. For gill phenotyping, the industry standard protocol and methodology utilized was adopted with some modifications to assess gill damage and severity (Taylor et al., 2009a). Each fish survivor was humanely dispatched in accordance with local laws and inspected under consistent bright light whereby gill structures were fully exposed to allow thorough examination of each gill arch, including the internal corners between the arches. This method allowed the generation of eight different scores (0-5) ranging from clear (0), with no infection present to heavily infected (5). In total, 16 gill scores were generated per survivor which were averaged per fish. The gill score description is presented in Table 1.

A tissue sample was collected from all the 1691 survivors at the termination of the second challenge test were sent to IdentiGEN (https://identigen.com/) in Dublin, Ireland for DNA extraction and genotyping using 67K Axiom Affymetrix SNP genotyping array BERGEN04 (www.thermofisher.com). The raw genotype data were quality filtered using the Affymetrix axiom analysis suite software using a sample call rate (CR) of ninety and a SNP CR-cut off of 95. Following this, a secondary filtering of samples and genotypes was conducted using Plink v 1.9 (Purcell & Chang, 2015) with the following specified filtering thresholds. Single nucleotide polymorphism markers with a minor allele frequency (MAF) lower than 5%, missing rate higher than 10%, and those not passing the Hardy–Weinberg equilibrium exact test (p< 1.0 × 10−6), were excluded from the data. Filtering was also performed at the individual level with animals removed from the dataset based on the missing genotype rate (>10%) and abnormal levels of heterozygosity (higher than ± 6 SD of population average). The filtering step retained a total of 62087 SNPs and after filtering quality control of the SNPs, a total of 1663 animals were available for analysis, (833 males and 830 female), which had genotype and phenotype records, with an average of eight records per family group.

Despite conducting two phases of the CGD challenge test infections, separated by a freshwater treatment, terminal sampling for gill scores was performed solely during the second infection. In the first infection, 12 random samples per tank were collected and recorded to detect the presence of the disease. Consequently, the subsequent data analysis was based on the gill scores obtained from the second phase of challenge test infection, comprising 1663 phenotypes.

A multivariate linear animal mixed model was used to estimate the variance components of gill scores for 1663 genotyped fish. The model, which included fixed effects related to the corresponding scoring person, tank, and random effect such as genomic breeding values (GEBVs), can be expressed as follows:

where y_i is the vector of observations for the three traits (y₁ average left and right gill score, y₂ is left gill score and y₃ right gill score). The X_i is the design matrix for fixed effect related to scoring person, and b(_i) contains the fixed effects for scoring person (three persons scoring gills), covariate body weight (g) at a time of gill scoring and experimental challenge tank (two tanks). Z_i is the design matrix for random effects associated with genomic breeding values, and g_i is the vector for genomic breeding values for each trait. The random effects are assumed to follow multivariate normal distributions:

g_i ~ MVN(0, V⨂G), where V is the additive genetic 3 × 3 (co) variance matrix between traits, ⨂ is the direct product operator, G is the genomic relationship matrix.

e_i ~ MVN(0, R⨂I_i) where R is the 3 × 3 residual (co) variance between the traits. and I_i is the identity matrix for the 1663 fish. The environment effect common to full-sibs has not included in any model because all families were pooled at a very early stage during production in Iceland prior to shipping to University of Stirling.

The genomic relationship matrix (G) was computed according to (VanRaden, 2008) method as follows:

where z is a centered (-2 p_i)  matrix of SNPs marker genotypes that were coded as 0 for homozygotes AA, 1 for heterozygote AB/BA and 2 for homozygote BB assuming B as the reference allele ( pi) and Nsnp is the total number of SNPs marker genotypes used in the analysis.

Variance component and breeding value for the 1663 animals recorded for gill scores were estimated using the Asreml package (Gilmour et al., 2002). The heritability (h²) was estimated separately for each trait using the following formula:

where σg2 is the additive genetic variance from G matrix for each trait; σe2 is the residual variance derived from the diagonal residual variance matrix (R) for each trait; and σp2=σg2+σe2.

Additionally, a multivariate linear animal mixed model, similar to the one described above, was employed to estimate genetic parameters and correlations for gill scores across all 16 gill arches (R1–R8 and L1–L8). Due to convergence issues, the model was run with 8 traits at a time until all heritability and covariance estimates were obtained.

To estimate the genetic correlation between gill score and weight of fish at scoring, a bivariate animal model was used. In this simplified version, fish weight was excluded as a covariate. The model was run separately for average left + right gill score, left gill score and right gill score.

The genetic correlation between all the traits were estimated using multivariate models described above and calculated as:

rg=σg1g2σg12σg22, where σg1g2 was the estimated covariance between both traits, σg12 is the additive genetic variance of trait 1, and σg22 is the genetic variance of trait 2. To test significance of genetic correlation between traits, the additive genetic correlation among the traits was assessed by likelihood ratio test LR= −2*(LogLR−LogLF) (Lynch and Walsh, 1998) using similar bivariate model as described above, where LogLR is the log of restricted likelihood of the model, where the genetic correlation among the traits is fixed to zero and LogLF is the log of the restricted likelihood of the model where the genetic correlation between the traits is not fixed. This statistic follows a mixture of distribution between 0 and 1 degrees of freedom.

To estimate the association of gill scores with the SNPs marker genotypes, genome-wide association analysis was performed using the linear mixed animal model implemented in the GCTA package with the option of “—mlma-loco” and the total variance explained by all SNPs was estimated using the options of “—reml” functions (Yang et al., 2011).

where y, is a vector of n (n=1663) gill scores, b_(i) contains the fixed effects for scoring person (three persons scoring gills), covariate body weight (g) at a time of gill scoring and experimental challenge tank (two tanks); W is the incidence matrix for the SNPs containing maker genotypes coded as 0 =AA, 1 =AB/BA, 2 =BB and α is the allele substitution effect for each SNPS; Z is the incidence matrix of genotyped individuals, u is the vector of genomic breeding values with u ~ N(0,Gσu2), e is the vector of random residual effects with e  ~ N(0,Iσe2).  The G matrix was computed as described above in the variance component section. Markers were considered genome significant when they crossed the Bonferroni thresholds for multiple testing (alpha =0.01) of 0.01/tg, where tg = 62510 markers were considered chromosome wide significant when Bonferroni threshold for multiple testing surpassed (alpha =0.01) of 0.01/tc, where tc = 2156 (average number of SNPs per chromosome). To evaluate the distribution of the test statistics and to assess potential genomic inflation arising from population structure or technical artifacts, quantile-quantile (QQ) plots were generated by comparing the distribution of observed versus expected -log₁₀ (p-values). The genomic inflation factor (Lambda, λ) was calculated based on the median of the observed chi-squared statistics and provided a quantitative measure of deviation from the expected distribution.

Cross validation was performed to estimate the accuracy genomic prediction of gill-scores using the same model as detailed under the variance component estimation section. A random sample of the 1333 animals were used as a training set, while the remaining 300 animals were used as validation sets. This cross-validation procedure was replicated 50 times with the prediction accuracy computed in each replicate. The prediction accuracy for each replicate was computed as follows:

where ρ is the Pearson moment correlation coefficient, GEBV is the estimated genomic breeding values, y_adj is the adjusted phenotype (y_adj = y−Xb), and h² is the heritability of the trait.

To determine parentage and family connections, we conducted a parentage assessment based on a run of homozygosity pairwise counts between offspring and parent genotypes. The assignment was considered successful if the runoff homozygosity counts are less than the defined threshold, representing the number of SNPs due to error, and if the sire and dam combination matched the recorded crossings. For a fast run of homozygosity pairwise counts between parent and offspring, the recommendations of (Ferdosi and Boerner, 2014) were employed. This analysis relied on genomic data from all parents and offspring that met the initial filtering criteria. After parentage analysis, there was an 85% success of assigning parents to the survivors of CGD infection that passed the quality filtering (1663). For this reason, it was not possible to correctly identify all the family of CGD-challenged fishes. Therefore, we phased the parent together with the offspring to classify them into sib-groups using the hsphase R package (Ferdosi et al., 2014). The CGD-challenged fishes were then classified into a total of 206 sib groups.

The descriptive statistics of 1,663 mean gill scores are shown in Table 2. The minimum and maximum possible average gill scores for the 16 areas of the gill-arches (8 parts on each sides of the left and right gill arches) are 0 and 7. For the total area of the gill (left gill + right gill), the gill scores had a mean score of 3.64 and the total gill scores had a range of 1.2 – 6.5. For the left gill, the mean score was 3.62 and ranged from 1.12 - 6.5. The average score for the right gill was 3.66 and ranged from 1.0 - 6.8. These results were far away from the minimum expected value and close to the expected maximum value. The SD and CV for the total scores were 0.89 and 0.24, respectively. For scores from the left gill only, the SD and CV were 0.92 and 0.25, while for the right gill only, the SD and CV were 0.91 and 0.25. All three scientists who performed the gill scoring had comparable gill-score results to the total scoring. The distribution of average gill scores and the number of individuals scored by each person is given in Figure 1.

The distribution of the average CGD gill scores for the three traits analyzed: Left gill + Right gill, Left gill, and Right gill is shown in Figure 1. The average gill scores showed a normal distribution for all three traits analyzed. For all three traits analyzed, the average from the total CGD gill-score had a large variation across the animals challenged.

The mean gill scores across the different anatomical partitions for left gill, right gill and left and right gill and mean weight per sib groups are represented in Figure 2. The primary Y-axis (left) represents the distribution of the family mean CGD gill scores, while the secondary Y-axis (right) shows the number of fish per sib group (blue line). The X-axis displays sib groups arranged in ascending order of their averaged CGD gill score. The figure illustrates substantial variation in mean CGD gill score, with some families exhibiting lower average scores while others show higher susceptibility and increased gill damage. Of note for the Y-axis the grey bars correspond to average score per sib group and mean weight respectively while the blue line represents total individuals represented per sib group.

The estimated genetic parameters of gill scores from the multi-variate analysis using the genomic relationship matrix are presented in Table 3. The heritability estimate was slightly higher when all scorings from the total areas of the gill were used for the analysis. Moreover, the heritability estimates from the left and right sides of the gill were not significantly different from each other. Genetic correlation among the three traits analyzed was close to 1, indicating all three areas of the gill were the same trait. All the three gill score traits had a negative genetic correlation with the weight of the fish. Moreover, the heritability and genetic correlation estimate from a multi-variate analysis of the gill scores from the 16 areas of the gill are presented in Table 4. The heritability estimate had a range between 0.16 and 0.26. The genetic correlation among the 16 scorings of the gill was greater than 0.9.

In the matrix, the diagonal elements represent heritability estimates along with their standard errors (h² ± SE), while the off-diagonal elements display genetic correlations between trait pairs, also accompanied by their respective standard errors (r_g ± SE). The additive genetic variance (σ^2g) for the mean gill score traits was approximately 0.20 ± 0.03 and for environmental variance (σ^2e) was 0.51 ± 0.01. For body weight, the additive genetic variance was estimated at 2985 ± 321.

The genomic based prediction accuracy for the three traits analyzed are given in Table 5.

The accuracy of prediction was comparable for the three traits analyzed. The accuracy of prediction ranged from 0.6 to 0.62 for all the CGD traits analyzed, while the bias was 1.26-1.33 for the traits measured.

Figure 3 illustrates the genome-wide association analysis of SNPs related to robustness against CGD infection. The horizontal red line represents the threshold for genome-wide significance (P< 1.6×10^-5), while SNPs highlighted in green indicate those exceeding the suggestive significance threshold (P< 1 × 10^-5). For total gill scores (Figure 3A), 2 SNPs on chromosome 24 surpassed the suggestive threshold. In addition, a single SNP on another chromosome 16 also reached this level. In the analysis of left gill scores (Figure 3B), one SNP on chromosome 24, three SNPs on chromosome 9, and one SNP on chromosome 5 exceeded the suggestive threshold. For the right gill scores (Figure 3C), only one SNP on chromosome 16 passed the suggestive significance threshold. In the genome wide association analysis conducted for all three traits, no SNPs surpassed the genome-wide significance threshold.

Quantile-quantile (QQ) plots were generated to visualize the distribution of the observed versus expected -log₁₀(p-values) and to assess potential genomic inflation factor attributable to the population architecture or technical artifacts. The plots demonstrated close alignment with the null expectations across the majority of the distribution, with modest deviations suggestive of true associations. Genomic inflation factors (λ) ranged from 1.18 to 1.20, suggesting marginal systematic bias for each of the traits assessed (Left gill + right gill = 1.20, Left gill = 1.18 Right gill =1.20 respectively Figures 3A–C).

Parentage assignment was conducted using genomic data from both parents and offspring used in the challenge. Due to the presence of ungenotyped parents and potential genotyping errors, complete family identification was not feasible for all individuals. To address this limitation, all experimental fish were grouped into 206 sib groups after phasing the SNPs information of all the fish and their parents using hsphase R package (Ferdosi et al., 2014). While the dataset contained more sib groups than families, this imbalance was not expected to bias the results, given that family-level identification was required to evaluate family effects and demonstrate family differences in gill severity to CGD. Moreover, statistical modelling indicated that the sib group effect was not significant in either the univariate and/or multivariate models. The presence of tank effect was minimized due to common rearing prior to the CGD challenge test, albeit other common environmental effects including, maternal and non-additive genetics effects may have contributed to the observed variation. Consequently, the use of sib groups in place of family assignments did not significantly impact the robustness of the study’s outcomes.

Based on a multivariate genetic analysis of CGD gill scores and genotypic information, it has been determined that the gill score is a trait with moderate heritability (h² = 0.28 - 0.29). This finding is consistent with previous genomic heritability estimates of h² = 0.28 and h² = 0.24, reported in Norwegian and Scottish Atlantic salmon challenge test experiments for AGD (Kube et al., 2012; Robledo et al., 2018; Boison et al., 2019). Interestingly, our study revealed that analyzing CGD gill scores separately for the left and right gill regions resulted in the same outcome as total scores. Furthermore, by subjecting each of the 16 CGD gill scores to a multivariate analysis as individual traits, we obtained heritability estimates ranging from 0.16 to 0.26. Notably, the genetic correlation across all 16 scores exceeded 0.90. These findings suggest that the subdivision of the gill into 16 distinct areas, as done in the current study, may not be necessary for phenotyping purposes. Rather, a division into eight or fewer areas would produce comparable outcomes. This result has significant importance for future gill challenge tests and could reduce the associated burden of gill phenotyping to either left or right side only.

In addition, the presence of substantial genetic variance (σ^2g = 0.2) as compared to environmental variance (σ^2e = 0.5) indicates that a considerable proportion of phenotypic variation in gill score is under genetic control. This is reflected by the moderate heritability estimates of gill damage to CGD in the Icelandic nucleus breeding program and demonstrates the potential for selection response to CGD within this breeding program. Due to the inability to reliably score gill condition post-mortem owing to decomposition, only surviving fish were phenotyped and genotyped. However, this approach may lead to an underestimation of genetic variance and heritability (h²) if the excluded mortalities had systematically different scores. While this limitation was necessary to ensure data quality, it introduces uncertainty in the variance component estimates and should be considered when interpreting the results. Furthermore, future selection candidates identified from the current study should undergo performance tests within a commercially relevant environment to further validate selection and gill damage to CGD.

The genome-wide estimated breeding value for gill scores presented revealed individual and family differences in gill severity to CGD infection within the Icelandic salmon breeding nucleus. This highlights the potential benefits of incorporating CGD gill scores into the breeding program index for fast-tracking this trait for improving disease resistance via selective breeding (Yáñez et al., 2014). Notably, these findings were based on the challenge test albeit the material used to propagate the disease challenge was isolated from a relevant field outbreak of gill disease in Scottish net pens. Previous studies indicate low genetic correlation (r_g = 0.02 - 0.34) between challenge tests and field testing for AGD resistance in the commercial environment (Gjerde et al., 2019; Lillehammer et al., 2019). This may be due to the multiple gill pathogens being present during field trials vs. single pathogen challenge tests and/or the complexity of the natural environment which may encompass all factors and their interactions that may contribute to increased pathogenicity (Norris et al., 2008). This suggests the outcomes in the event of an AGD field outbreak may differ depending on the disease type and composition. During bacterial disease challenge studies, favorable genetic correlations between challenge and field tests have been found (Gjøen et al., 1997). In addition, the individual and average sib groups’ breeding value estimate shows there is a large variation in gill severity to CGD gill infection among individuals and sib groups. Consequently, validation of these results using a complex gill disease field outbreak is essential in light of previous studies which low genetic correlations between challenge tests and field trials for AGD.

In the current study, we conducted a genome-wide association study (GWAS) analysis to identify the specific genomic regions associated with AGD gill scores in fish that survived the challenge test. Apart from a few SNPs on chromosomes 9, 16, 24 and 5 surpassing the suggestive level threshold, our analysis did not show any significant SNPs that passed the genome-wide threshold. In contrast to the findings of this study, previous studies have reported that significant SNPs were found on chromosomes 9 and 16 for the trait of mean gill score in an Atlantic salmon challenge test, while a significant SNP was reported on chromosome 9 and others for total AGD gill scores in a challenge test experiment (Boison et al., 2019). However, the current study did not observe genome-wide associations between gill score due to CGD and SNPs. The absence of genome-wide significant SNPs in our analysis may be attributed to the small effect sizes of the genetic variants and the limited size of our dataset. Our analysis was conducted using data from a multiple challenge test derived from a single year-class, specifically focusing on phenotyping after the second gill challenge. These constraints may have limited the statistical power to detect significant associations, or alternatively, may suggest that gill score due to CGD is polygenic in nature which is consistent with previous findings on AGD resistance (Ajasa et al., 2024a). Significant associations between AGD-related gill scores and specific genomic regions have been identified in both controlled challenge experiments and natural field outbreaks (Aslam et al., 2020). To better detect significant associations between gill severity to CGD infection and specific genomic regions, this may require larger datasets from multiple experiments and/or field outbreaks.

The genomic prediction accuracies for the three CGD-related traits were moderately high, with relatively low standard errors (Table 5). The highest accuracy was observed for the combined gill score, likely due to increased information content and reduced measurement noise compared to single-side assessments. These results align with previously reported within-population prediction accuracy for AGD in Atlantic salmon which ranged from 0.56 to 0.76 across different year classes (Ajasa et al., 2024b). Key indicators of model performance were also within the expected range. The reliability calculated as the ratio of GEBV variance to additive genetic variance was 0.55, with a corresponding accuracy of 0.74. These values are consistent with the cross-validation results and support the conclusion that the model’s predictive ability is appropriate. Overall, the results suggest that CGD-related traits have moderate heritability and are suitable for genomic selection, similar to other complex disease traits in salmon. As expected, genomic prediction accuracy is influenced by factors such as the extent of linkage disequilibrium (LD) between markers and QTL, the size of the reference population, trait heritability, and the distribution of QTL effects (Hayes et al., 2009). While the relationship between the training and validation sets does not exactly reflect the common practice where the validation set is constructed to mimic sib-based genomic selection by including at least one individual per family (Ajasa et al., 2024a), our approach of randomly masking 300 animals and repeating the process 50 times appeared to have created appropriate relationships between the training and validation sets based on accuracies obtained.

The results of our validation indicated a surprisingly high bias, with a regression coefficient slope of 6.8 between the adjusted phenotype and genomic estimated breeding values (GEBV) when using the full datasets, and a consistent value of 6.6 in the 50-iteration cross-validation sets. A bias of this magnitude suggests a considerable underestimation of true breeding values (TBV), which is highly unusual in quantitative genetics. We hypothesize that this inflated bias may stem from an overestimation of the covariance between adjusted phenotypes (y_adj) and GEBV, potentially due to properties of the CGD gill score data or significant differences between adjusted and observed phenotypes.

It is important to note that the slope derived from a scatter plot of GEBVs from a cross-validation against GEBVs from the full dataset is not a valid measure of bias. The regression coefficient for bias is specifically defined as the relationship between GEBV and TBV, not between two sets of GEBVs. Therefore, any slope from such a plot should not be interpreted as a measure of bias. Given the challenges in accurately estimating bias with the conventional method using adjusted phenotypes, it may be recommended to adopt the semi-parametric approach proposed by Legarra and Reverter (2018). This method offers a robust way to estimate population accuracy and prediction bias without relying on adjusted phenotypes. Instead, it evaluates a dispersion metric that captures how closely on average GEBVs from full and partial data sets align. This approach has recently been applied in aquaculture by Sae-Lim et al. (2025) and applying this method would provide a more reliable evaluation of the model’s predictive performance.