Genome-wide Association Studies Reveal Similar Genetic Architecture with Shared and Unique QTL for Bacterial Cold Water Disease Resistance in Two Rainbow Trout Breeding Populations

Bacterial cold water disease (BCWD) causes significant mortality and economic losses in salmonid aquaculture. In previous studies, we identified moderate-large effect QTL for BCWD resistance in rainbow trout (Oncorhynchus mykiss). However, the recent availability of a 57K SNP array and a genome physical map have enabled us to conduct genome-wide association studies (GWAS) that overcome several experimental limitations from our previous work. In the current study, we conducted GWAS for BCWD resistance in two rainbow trout breeding populations using two genotyping platforms, the 57K Affymetrix SNP array and restriction-associated DNA (RAD) sequencing. Overall, we identified 14 moderate-large effect QTL that explained up to 60.8% of the genetic variance in one of the two populations and 27.7% in the other. Four of these QTL were found in both populations explaining a substantial proportion of the variance, although major differences were also detected between the two populations. Our results confirm that BCWD resistance is controlled by the oligogenic inheritance of few moderate-large effect loci and a large-unknown number of loci each having a small effect on BCWD resistance. We detected differences in QTL number and genome location between two GWAS models (weighted single-step GBLUP and Bayes B), which highlights the utility of using different models to uncover QTL. The RAD-SNPs detected a greater number of QTL than the 57K SNP array in one population, suggesting that the RAD-SNPs may uncover polymorphisms that are more unique and informative for the specific population in which they were discovered.


INTRODUCTION
Bacterial cold water disease (BCWD) causes significant mortality and economic losses in salmonid aquaculture (NEMATOLLAHI et al. 2003;BARNES AND BROWN 2011). The etiological agent of BCWD is a gram-negative bacterium, Flavobacterium psychrophilum (Fp), and current methods to control BCWD outbreaks are limited. At the National Center for Cool and Cold Water Aquaculture (NCCCWA), we have developed a selective breeding program to increase rainbow trout genetic resistance against BCWD and have shown that BCWD resistance is a moderately heritable trait that responds to selection (LEEDS et al. 2010). Furthermore, we revealed complex genetic architecture of BCWD resistance (VALLEJO et al. 2010) and identified several moderate-large effect quantitative trait loci (QTL) for this trait in the NCCCWA odd-and even-year rainbow trout selective-breeding populations (WIENS et al. 2013;VALLEJO et al. 2014a;LIU et al. 2015b;PALTI et al. 2015b). While those loci can be fine mapped to identify positional candidate genes, the complex genetic architecture of BCWD resistance and high genetic variability discovered in past studies (VALLEJO et al. 2014a) suggest that whole genomeenabled selection is more effective for improving genetic resistance against BCWD in rainbow trout aquaculture, and we were able to empirically demonstrate that whole genome-enabled selection can double the accuracy of predicted genetic merit of potential breeders compared to traditional family-based selection (VALLEJO et al. 2017).
For agricultural livestock species, single nucleotide polymorphism (SNP) chips have been the platform of choice for whole genome genotyping of at least 50K SNPs (MATUKUMALLI et al. 2009;RAMOS et al. 2009;GROENEN et al. 2011); including the recently developed 57K 5 SNP chip for rainbow trout (PALTI et al. 2015a). However, sequencing-by-genotyping methods that do not require a priori marker discovery or a reference genome sequence and are capable of simultaneous marker discovery and genotyping in many individuals were developed for genetic analyses (DAVEY et al. 2011). One such technique is restriction-site-associated DNA (RAD) sequencing (MILLER et al. 2007;BAIRD et al. 2008). In recent years, the method of RAD genotyping by sequencing has been widely used in salmonid species for SNP discovery, generating linkage maps, QTL mapping, genome-wide association studies (GWAS) and for evaluating genome-enabled selection (HECHT et al. 2012;HOUSTON et al. 2012;MILLER et al. 2012;HALE et al. 2013;HECHT et al. 2013;NARUM et al. 2013;BRIEUC et al. 2014;CAMPBELL et al. 2014;GONEN et al. 2014;HOUSTON et al. 2014;PALTI et al. 2014;LIU et al. 2015a;LIU et al. 2015b;PALTI et al. 2015b;VALLEJO et al. 2016).
Several moderate-large effect QTL associated with BCWD resistance have been identified on 24 of the 29 rainbow trout chromosomes using linkage analysis mapping (JOHNSON et al. 2008;WIENS et al. 2013;PALTI et al. 2014;QUILLET et al. 2014;VALLEJO et al. 2014a) and GWAS methods (CAMPBELL et al. 2014;LIU et al. 2015b;PALTI et al. 2015b). However, these previous studies had limitations on QTL detection. First, most of the reported genome-wide association analyses have tested one SNP at a time using single-regression or mixed linear models with a fixed SNP effect along with a random polygenic effect to capture the effects of all other genes. Although those studies have been successful in detecting associations, those associations typically explain only a small fraction of the trait genetic variance (VISSCHER et al. 2010). Conversely, in GWAS analysis using whole-genome selection models that simultaneously fit all SNPs as random effects, the SNPs jointly explain a larger proportion of the genetic 6 variance which highlights the utility of using multiple-regression GWAS models with the SNPs joined into genomic windows for accurate QTL mapping (HAYES et al. 2010;FAN et al. 2011;ONTERU et al. 2011). Second, most of the reported QTL for BCWD resistance were identified using linkage-based methods and GWAS was performed within individual segregating families with relatively small sample size and low detection power. Consequently about half of the reported findings from such QTL mapping studies are expected to represent false positives (HIGGINSON AND MUNAFO 2016). Third, the previous studies used lower density genotyping platforms and did not have a reference genome physical map for accurate prediction of the order and physical proximity of the genetic markers. Furthermore, to our knowledge, in the current study we are using for the first time SNP genotype data with the genome physical map coordinates for GWAS in rainbow trout.
There is ambiguity on the best computational algorithm when using multiple-regression based models in genomic selection (GS) and GWAS experiments. The genetic architecture of the trait and the population structure can have a significant impact on the accuracy of the genomic predictions and estimated marker effects. Therefore, it is important to compare the performance of the best competing algorithms on GS and GWAS when evaluating a trait with complex inheritance for the first time in a population. This will allow effective discovery of QTL underlying the genetic basis of the complex trait and control the type I error rate, which is often high in genome-wide discovery experiments.
In GWAS with models that fit all SNPs simultaneously, the genomic BLUP (GBLUP) method assumes a polygenic architecture of the trait and uses all SNP data to estimate the genomic relationship (G) matrix. In contrast, the Bayesian variable selection method assumes that the genetic variance is explained by a reduced number of markers with small-moderate or large effects (HABIER et al. 2007;HAYES et al. 2009; DE LOS CAMPOS et al. 2013;FERNANDO AND GARRICK 2013;HOWARD et al. 2015). Based on this assumption, GBLUP is expected not to perform as well as Bayesian variable selection models when the trait is controlled by few moderate-to-large effect QTL. The GBLUP method was modified into the single-step GBLUP (ssGBLUP) method which allows combining the pedigree (A) and genomic-derived relationships into an H relationship matrix (AGUILAR et al. 2010;LEGARRA et al. 2014), and to the weighted single-step GBLUP (wssGBLUP) method which emulates the Bayesian variable selection models by fitting in the multiple regression model selected SNPs that explain moderate-large fraction of the trait genetic variation (WANG et al. 2012).
The recent development of the 57K SNP array (PALTI et al. 2015a), a dense genetic linkage map with 47,939 SNP markers (GONZALEZ-PENA et al. 2016), and the release of the improved rainbow trout reference genome (GenBank assembly Accession GCA_002163495) have provided the needed tools for performing GWAS to identify genomic regions associated with BCWD resistance in rainbow trout. The main objectives of this study were to (1) identify and validate QTL associated with BCWD resistance in two commercially-relevant rainbow trout breeding populations; (2) characterize the genetic architecture of rainbow trout resistance to BCWD; (3) compare the QTL mapping efficiency and determine whether the Chip-SNP and RAD-SNP genotyping platforms detect the same QTL; and (4) compare the QTL mapping efficiency of two widely-used multiple-regression GWAS models.

Ethics statement
Protocols for this study were reviewed and approved by the NCCCWA Institutional Animal Care and Use Committee (Kearneysville, WV).

Rainbow trout rearing and BCWD challenge
Details of the 21-day BCWD challenge have been reported elsewhere (SILVERSTEIN et al. 2009;LEEDS et al. 2010). Mortalities were removed and recorded daily and fin clipped. Fish that survived the challenge were euthanized in a lethal dose of Tricaine methane sulfonate (Tricaine-S, Western Chemical, Inc., Ferndale, WA) and fin clipped. Fin clips from all fish (mortalities and survivors) were individually kept in 95% ethanol until DNA was extracted using established protocols (PALTI et al. 2006).

Rainbow trout populations used in GWAS
Fish used in this study were sampled from two populations, and all analyses were performed separately by population. The first sample included fish with genotypes and phenotypes from 10 full-sib (FS) families sampled from a total of 71 pedigreed FS families with phenotype data from year-class (YC) 2005 of the NCCCWA BCWD resistant line (NCCCWA) and it was described in our previous GS study (Vallejo et al., 2016). Briefly, the YC 2005 families represented the base generation of the breeding line, and thus had not previously been selected for BCWD resistance. Each family had n = 39−80 fish evaluated in the laboratory BCWD challenge in one or two tanks per family. The phenotypic dataset included disease resistance phenotypes from n = 1 0 phenotype STATUS had two categories: 1 = the fish died during the 21 days post challenge evaluation period; and 2 = the fish survived for the duration of the challenge. The DAYS and STATUS records were analyzed separately using univariate GWAS models described below.

SNP genotyping platforms
The fish sampled from TLUM and NCCCWA populations were genotyped using the Rainbow Trout Affymetrix 57K SNP array (Chip) following previously described procedures (PALTI et al. 2015a) and the samples were genotyped by a commercial service provider (Geneseek, Inc., Lincoln, NE) following the Axiom genotyping procedures described by the array manufacturer (Affymetrix). The quality control (QC) bioinformatics pipeline applied to the Chip-SNP genotype data collected in the TLUM (VALLEJO et al. 2017) andNCCCWA (VALLEJO et al. 2016) populations were already described. After genotype data QC, a total of 41,868 and 49,468 SNPs were included in the TLUM and NCCCWA raw Chip genotype datasets, respectively.
The fish sampled from the NCCCWA population were also genotyped by sequencing of restriction-site-associated DNA (RAD) tag libraries as we have previously described (VALLEJO et al. 2016). After genotype data QC, a total of 24,465 RAD-SNPs were included in the raw RAD genotype dataset. The raw sequence data from the RAD libraries were deposited in the NCBI SRA database (Accession SRP063932).
Before performing GWAS analyses, the raw marker genotype datasets were further QC filtered using algorithms implemented in the software BLUPF90 (MISZTAL et al. 2015). For the 1 1 frequency higher than 0.05, and departures from Hardy-Weinberg equilibrium less than 0.15 (difference between expected and observed frequency of heterozygotes). Parent-progeny pairs were tested for discrepant homozygous SNPs, and SNPs with a conflict rate of more than 1% were discarded from further analysis. For the RAD data, the QC retained SNPs with a genotype calling rate higher than 0.70. Following this final QC step,33,838 SNPs,39,112 SNPs and 9,534 SNPs were retained for analyses of the TLUM, NCCCWA (Chip) and NCCCWA (RAD) datasets, respectively. Next, we determined the physical map position (GenBank assembly Accession GCA_002163495) of each of the QC filtered markers and found that a small fraction did not have a physical map location. The numbers of effective genotyped markers and effective genotyped fish that were used with each specific GWAS model and genotyping platform in the evaluated populations are shown in Table 1.

Genome-wide association analyses
This study was conducted to identify chromosomal regions that have the greatest impact on the variation of BCWD resistance because it is difficult to infer individual SNP effects from a multiple regression model that fits markers simultaneously at a genome-wide scale (GARRICK AND FERNANDO 2013). Therefore, instead of using the markers effect to make an inference on a particular locus, we used the markers effect to make an inference about a particular genomic region that encompasses a number of contiguous loci in association with the trait (FAN et al. 2011;FERNANDO et al. 2014). Thus the GWAS analysis was performed separately for each population and BCWD resistance phenotype using multiple-regression GWAS models that use simultaneously all the SNPs in the association test. Two multiple-regression GWAS methods 1 2 were used: Bayesian variable selection BayesB (BayesB) (FERNANDO AND GARRICK 2009) and weighted single-step GBLUP (wssGBLUP) (WANG et al. 2012;MISZTAL et al. 2015).
Before starting the search for genomic regions associated with BCWD resistance, we tested 0.5 and 1Mb exclusive-windows in GWAS with BayesB model and found that 1Mb windows provided Manhattan plots with less noisy baseline which was in agreement with other GWAS reports (KIZILKAYA et al. 2013;SAATCHI et al. 2013); so we decided to use 1Mb exclusive-windows in the GWAS performed with BayesB. We also used 1Mb sliding-windows in GWAS performed with wssGBLUP which are more informative than exclusive-windows.

Bayesian variable selection model BayesB
The BayesB model uses only fish that had both genotype and phenotype records. The TLUM population sample included n = 1473 fish from 50 YC 2013 families with phenotype and genotype records, and the NCCWA sample included n = 577 fish from 10 YC 2005 families with phenotype and genotype records ( Table 1) In GWAS with wssGBLUP, the weights for each SNP are 1's for the 1 st iteration which means that all SNPs have the same weight (i.e., standard single-step GBLUP). For the next iterations (2 nd , 3 rd , etc.), the weights are SNP specific variances calculated using both the SNP allele-substitution effect estimated in the previous iteration and their corresponding allele frequencies (WANG et al. 2012). In this study, we decided to use results from the 2 nd iteration because they provide the highest accuracy genomic predictions (VALLEJO et al. 2016) and marker effects (WANG et al. 2012;IRANO et al. 2016;MELO et al. 2016).
In GWAS with TLUM population sample, the linear and threshold models for DAYS and STATUS, respectively, included the effects of population mean, random animal, random common environment and random error. The full-sib fish progeny from each family was allocated into two tanks for BCWD challenge evaluation, so the variable tank nested within family was used to model the common environment effect. The linear and threshold models used with the NCCCWA sample were similar to those used with the TLUM sample. The linear model for DAYS and threshold model for STATUS were fitted using computer applications 1 5 implemented in the software BLUPF90 (MISZTAL et al. 2015). The binary phenotype disease STATUS was analyzed with a threshold model using a Bayesian approach which included a single chain with a total of 270,000 iterations; the first 20,000 iterations were discarded as burnin iterations; then from the remaining 250,000 samples, one from every 50 samples were saved.
Thus, 5000 independent samples were used in the analysis. The proper mixing and convergence of these MCMC iterations were also assessed using the R package CODA (PLUMMER et al. 2006).

Criteria to declare QTL associated with BCWD resistance
The results from the GWAS performed with BayesB and wssGBLUP were used to identify genomic windows and QTL associated with BCWD resistance. A two stage approach was used to identify a QTL associated with BCWD resistance. First, the genomic windows with explained genetic variance (EGV) greater than 1% and 2% in the TLUM and NCCCWA populations, respectively, were declared as genomic regions associated with BCWD resistance. The threshold to declare a QTL was raised to EGV≥2% in the relatively small NCCCWA sample (n = 577) to control the type I error rate. Second, to determine if neighboring or overlapping windows on the same chromosome belong to the same QTL region we used the following criteria: all windows associated with BCWD resistance that were bounded within a region smaller than 20Mb and were less than 10Mb apart from another associated window were grouped to a single QTL region. The QTL nomenclature we used was based on the chromosome number and QTL region within the chromosome, where the region with the lowest genome assembly position numbers

QTL segregating in both NCCCWA and TLUM populations
In order to identify QTL that might be overlapping between the two populations and between the two genotyping platforms in the NCCCWA population, we assigned the genome physical map positions to all flanking markers of the genomic windows associated with BCWD resistance using the rainbow trout reference genome sequence (GenBank assembly Accession GCA_002163495) and searched for overlapping QTL regions within each chromosome using the flanking markers physical map genome coordinates.

Data availability
1 7 The authors state that all data necessary for supporting the conclusions of this research article are included within the article and it's Supplemental Material. Figure S1 shows Manhattan plot with GWAS results for survival DAYS in TLUM sample genotyped with 57K Chip-SNP. Figure S2 shows Manhattan plot with GWAS results for survival DAYS in NCCCWA sample genotyped with 57K Chip-SNP. Figure S3 shows Manhattan plot with GWAS results for survival DAYS in NCCCWA sample genotyped with RAD-SNPs. Table S1 contains summary of all genomic windows and QTL detected in this GWAS. Table S2 contains list of QTL associated with survival DAYS found in TLUM population using the 57K Chip-SNP. Table S3 contains list of QTL associated with survival DAYS found in NCCCWA population using the 57K Chip-SNP.

Heritability of BCWD resistance
The heritability or proportion of phenotypic variance explained by the markers for survival DAYS and the binary survival STATUS were previously reported (VALLEJO et al. 2016;VALLEJO et al. 2017). Briefly, here they were moderate and relatively constant with a range of 0.24−0.34 (DAYS) and 0.23−0.35 (STATUS) ( Table 1); and the mean heritability for DAYS and 1 8 STATUS were similar (0.29). The heritability of STATUS estimated on the underlying scale of liability using a threshold model was transformed to the observed scale of disease survival STATUS using already described procedures (VALLEJO et al. 2017). Overall, for both BCWD phenotypes and across genotyping platforms and populations, the mean heritability estimated with wssGBLUP (0.32) was slightly higher than that estimated with BayesB (0.27).

QTL associated with BCWD resistance in the TLUM population (57K SNP Chip)
We have previously shown that the BCWD resistance phenotypes DAYS and STATUS yielded similar results in QTL mapping (VALLEJO et al. 2014a;LIU et al. 2015b;PALTI et al. 2015b) and genomic selection experiments (VALLEJO et al. 2016;VALLEJO et al. 2017). In this study, we observed also that STATUS and DAYS were affected by similar QTL regions with few exceptions (Overall, DAYS detected three more QTL than STATUS) (Table S1) In the TLUM population, a total of 45 windows with EGV ≥ 1% were detected on chromosomes Omy3, 5, 8, 13 and 25 (Table 2; Figure 1; Table S1 and Figure S1). Fourteen windows were detected with BayesB with EGV up to 57.6% and 31 windows were detected with 1 9 Four QTL (3.2, 8.1, 13.2 and 25.1) were detected by both GWAS models (Table S1) Table 2).

QTL for BCWD resistance in the NCCCWA population (57K SNP Chip)
In the NCCWA population using the 57K SNP chip, we detected a total of 11 windows associated with BCWD resistance on chromosomes Omy3, 5, 10, 22 and 25 (Table 3; Figure 2; Table S1 and Figure S2). Two and nine QTL windows were detected with BayesB and wssGBLUP, respectively, and each GWAS model explained up to 5.6% and 16.1% of the genetic variance, respectively.
Only one QTL (3.2) was detected by both GWAS models (Table S1). We did not detect any BayesB model specific QTL, and four QTL (5.1, 10.1, 22.1 and 25.1) were detected only with the wssGBLUP method.

0
Overall, we detected five QTL (3.2, 5.1, 10.1, 22.1 and 25.1) associated with BCWD resistance that explained up to 18.2% of the genetic variance in the NCCCWA population when accounting only for the largest EGV window in each QTL (Table 3 and Table S3).

QTL for BCWD resistance in the NCCCWA population detected with RAD SNPs
In the NCCCWA population using the RAD SNP genotypes, we detected a total of 18 windows associated with BCWD resistance on chromosomes Omy3, 5, 10, 11, 13, 15 and 25 (Table 4; Figure 3; Table S1 and Figure S3). Four windows were detected by BayesB with EGV up to 17.3% and 14 windows were detected by wssGBLUP with EGV up to 26.4%.

DISCUSSION
In GWAS studies, the use of correct statistical models and computer algorithms is paramount to successfully identify the underlying genetic basis of resistance to complex diseases in livestock and aquaculture species. To date there have been several reported GWAS using single-marker association tests in fin fish species (CAMPBELL et al. 2014;AYLLON et al. 2015;GENG et al. 2015;GONEN et al. 2015;LIU et al. 2015b;PALTI et al. 2015b;TSAI et al. 2015;TSAI et al. 2016) which generally are associated with high type I error rate because single-marker methods do not account for linkage disequilibrium (LD) between physically linked loci in the association test. In this study, we performed GWAS for loci associated with BCWD resistance using multiple- In the current GWAS, we detected 14 QTL associated with BCWD resistance in two commercially-relevant rainbow trout breeding populations from which 11 were validated QTL from previous studies and three were novel QTL. Here we confirmed that BCWD resistance is controlled by the oligogenic inheritance of few moderate-large effect loci and a large-unknown number of loci with small effects on BCWD resistance. However, despite the similar genetic architecture for this trait in both populations and the detection of overlapping QTL, we still detected major QTL differences between the two populations. We also found that the RAD and Chip genotyping platforms did not detect the same QTL in the NCCCWA population, and overall the RAD platform detected a greater number of QTL than the Chip platform. In addition, 2 2 the wssGBLUP and the BayesB multiple-regression GWAS models did not detect the same QTL, which highlights the utility of using different GWAS models to effectively optimize the discovery of QTL.

The genetic architecture of BCWD resistance in rainbow trout
Previously, we predicted that six to 10 QTL explaining 83% to 89% of phenotypic variance with either additive or dominant disease-resistant alleles plus polygenic effects may underlie the genetic architecture of BCWD resistance in the same NCCCWA population using Bayesian complex segregation analysis of phenotype and pedigree records (VALLEJO et al. 2010). In the current study we were able to confirm our prediction on the genetic architecture of the trait.
With 10 families from that original NCCCWA population, we uncovered 10 moderate-large effect QTL that explained up to 27.7% of the additive genetic variance for BCWD resistance (Table S1). Similarly, in the TLUM odd-year population, we detected 8 moderate-large effect QTL that explained up to 60.8% of the additive genetic variance for BCWD resistance.
Four QTL regions located on chromosomes Omy3, 5, 13 and 25 are segregating in both populations (Table S5). The shared QTL regions explain a substantial proportion of the additive genetic variance for BCWD resistance in the two populations (up to 18% and 38.6% of the genetic variance in NCCCWA and TLUM, respectively); suggesting a common underlying genetic architecture for BCWD resistance in the two populations. However, major differences were also detected between the two populations. Six QTL, which explained up to 9.7% of the genetic variance and are located on chromosomes Omy5, 10, 11, 15, 22 and 25 were found only in the NCCCWA population (Table S6). Conversely, four QTL which explained up to 22.2% of 2 3 the genetic variance and are located on chromosomes Omy3, 8 and 13 were only found in the TLUM population. Overall, our GWAS results confirmed the hypothesis that BCWD resistance is controlled by the oligogenic inheritance of several moderate-large effect QTL and many small effect polygenic loci (VALLEJO et al. 2010;VALLEJO et al. 2014a;LIU et al. 2015b;PALTI et al. 2015b).
Further fine-mapping of the BCWD-QTL position and eventual identification of putative candidate genes or disease-causal mutations would be advantageous for applying marker-based selection and advancing the understanding of the mechanisms of genetic resistance to BCWD in rainbow trout populations. This can be achieved by genotyping and disease testing a greater number of SNPs from positions within and near the major QTL regions, by re-sequencing highly characterized BCWD resistant and susceptible individuals as was successfully done in the search for the IPNV resistance gene in Atlantic salmon (MOEN et al. 2015). In addition, positional and functional candidate genes for the QTL can be generated by interrogating the newest version of the rainbow trout reference genome sequence (GenBank assembly Accession GCA_002163495).

Comparing the two SNP genotyping technologies in the NCCCWA population
Overall, the RAD genotyping technology (18 windows; EGV= 32.8%; SNP genotyping technologies (QTL 3.2, 5.1, 10.1 and 25.1). The overall better performance of RAD-SNPs than Chip-SNPs in detecting QTL associated with BCWD resistance in the NCCCWA population may be due to sample ascertainment bias effects considering that the 57K SNP Chip was developed using a collection of samples from different rainbow trout populations to maximize SNP polymorphism and discovery (PALTI et al. 2015a). Therefore, the polymorphic markers in the SNP Chip might be less informative than the SNPs we genotyped with the RAD technology, which were specifically discovered in the sampled families from the NCCCWA population, and are therefore more informative for characterizing genome loci polymorphisms in this dataset.

Comparing BayesB and wssGBLUP models
Overall, we noticed that the wssGBLUP detected higher number of windows (54) associated with BCWD resistance than the BayesB (20) across the three datasets (TLUM-SNP, NCCCWA-SNP and NCCCWA-RAD) we used here (Table S1). We also noticed that the QTL detection power of BayesB was more negatively impacted by sample size reduction than wssGBLUP: BayesB detected 14, 2 and 4 QTL windows and wssGBLUP detected 31, 9 and 14 QTL windows in TLUM-SNP, NCCCWA-SNP and NCCCWA-RAD datasets, respectively. Performing GS with the Chip genotyped SNPs, we have shown that BayesB predicts GEBVs with higher accuracy than wssGBLUP when using a training sample size of n = 1473 (VALLEJO et al. 2017); however, wssGBLUP outperforms BayesB when using a smaller training sample size of n = 583 (VALLEJO et al. 2016). So, in agreement with these previous GS results, it seems also that the QTL detection power of GWAS with BayesB is more sensitive to sample size reduction than wssGBLUP. We think that the difference in power robustness to sample size reduction of these 2 5 GWAS models is due to their algorithmic differences. The BayesB model does not explicitly use the available pedigree information and includes in the analysis only animals that had both genotype and phenotype records and as a consequence is less power robust to sample size reduction and needs a minimum sample size for optimal performance, i.e. n ≥ 1000 animals. In contrast, the power of GWAS analysis with wssGBLUP is more robust to sample size reduction because it uses the available pedigree information and includes in the analysis animals that have both phenotype and genotype records, and also animals that have only phenotype records (those with missing genotype data).
In Thus, these results highlight the importance of using at least two different GWAS algorithms to efficiently uncover the underlying genetic basis of resistance against BCWD in the studied populations.

Comparing QTL detected in TLUM and NCCCWA populations
In spite of the smaller sample size of the NCCCWA population in comparison to the TLUM population, we detected 10 QTL in the NCCCWA population and only eight QTL in the TLUM 2 6 population (Table S1). We hypothesize that because the TLUM sample size was much larger than the NCCCWA sample size, it is likely that the type I error rate was smaller in the TLUM sample than in the NCCCWA sample. Therefore, we predict that most of the QTL detected in the TLUM population are real, but some of the QTL detected in the NCCCWA population are false positives. Specifically, the NCCCWA population had two unique QTL (5.2 and 15.1) that were also not reported in past studies. Thus, those two NCCCWA-specific QTL may very be false positives.

Comparing our GWAS results with previous studies
Eleven of the 14 QTL we detected in this study were also reported in previous studies in the NCCCWA germplasm and in other populations (Table S7)  SNPs associated with BCWD resistance in another commercial rainbow trout population that were about 1Mb from our QTL Omy8.1 (Table S1 and Table S7); they also reported QTL that (2) they can also represent false positive results due to limitations and weaknesses of experimental-design and power of analysis as we describe here.

Conclusion
This GWAS is the most comprehensive genome-wide scan for QTL associated with BCWD resistance performed to date in two commercially-relevant rainbow trout breeding populations, using two whole-genome SNP genotyping platforms and two multiple-regression GWAS models. We identified a total of 14 moderate-large effect QTL associated with resistance to BCWD resistance, and four of those QTL were segregating in the two populations. These GWAS results confirmed that the genetic architecture of BCWD resistance is controlled by the oligogenic inheritance of few moderate-large effect genes and many small effect resistance loci.
Overall, the wssGBLUP detected higher number of QTL than the BayesB and both GWAS models did not detect the same QTL which highlights the utility of using two different GWAS algorithms to effectively discover QTL. The RAD genotyping platform detected higher number of QTL than the Chip technology and also both genotyping platforms did not detect the same QTL in the NCCCWA population. These GWAS results will advance the biological and 2 8 functional analysis of positional candidate genes using the annotation of the new rainbow trout reference genome (GenBank assembly Accession GCA_002163495). They are also likely to be used in the implementation of more efficient selective breeding strategies which will utilize the QTL-flanking SNPs in genome-enabled selection for BCWD resistance in rainbow trout aquaculture.  s  t  e  n  s  e  n  ,  I  .  A  g  u  i  l  a  r  a  n  d  I  .  M  i  s  z  t  a  l  ,  2  0  1  4  S  i  n  g  l  e  S  t  e  p  ,  a  g  e  n  e  r  a  l  a  p  p  r  o  a  c  h  f  o  r  g  e  n  o  m  i  c   s  e  l  e  c  t  i  o  n  .  L  i  v  e  s  t  o  c  k  S  c  i  e  n  c  e  1  6  6 : 5 4 -6 5 .  b  i  e  t  a  l  .  ,  2  0  1  5  G  e  n  o  m  e  w  i  d  e  a  s  s  o  c  i  a  t  i  o  n  a  n  d   g  e  n  o  m  i  c  p  r  e  d  i  c  t  i  o  n  f  o  r  g  r  o  w  t  h  t  r  a  i  t  s  i  n  j  u  v  e  n  i  l  e  f  a  r  m  e  d  A  t  l  a  n  t  i  c  s  a  l  m  o  n  u  s  i  n  g  a  h  i  g  h  d  e  n  s  i  t   b From each QTL, the window with the highest explained genetic variance is presented in this Table. c GWAS was conducted using Bayesian variable selection model BayesB (BayesB) and weighted single-step GBLUP (wssGBLUP) methods. BayesB used 1Mb exclusive-consecutive windows and wssGBLUP used 1Mb moving-sliding windows.  a From each QTL, the window with the highest explained genetic variance is presented in this Table. b GWAS was conducted using Bayesian variable selection model BayesB (BayesB) and weighted single-step GBLUP (wssGBLUP) methods. BayesB used 1Mb exclusive-consecutive windows and wssGBLUP used 1Mb moving-sliding windows.
c Explained genetic variance by tested window (%). d SNP positions in base pairs (bp) based on rainbow trout reference genome sequence (GenBank assembly Accession GCA_002163495). a From each QTL, the window with the highest explained genetic variance is presented in this Table. b GWAS was conducted using Bayesian variable selection model BayesB (BayesB) and weighted single-step GBLUP (wssGBLUP) methods. BayesB used 1Mb exclusive-consecutive windows and wssGBLUP used 1Mb moving-sliding windows.

2
c Explained genetic variance by tested window (%). d SNP positions in base pairs (bp) based on rainbow trout reference genome sequence (GenBank assembly Accession GCA_002163495).  Tables   Table S1 Summary of QTL associated1 with BCWD resistance detected using two GWAS methods and two SNP genotyping platforms in NCCCWA and TLUM populations. (XLSX file)

Table S2
Summary of QTL associated with BCWD survival DAYS in the Troutlodge US May (TLUM) population. (DOCX file)

Table S3
Summary of QTL associated with BCWD survival DAYS in NCCCWA population detected using the 57K Chip-SNP. (DOCX file)

Table S4
Summary of QTL associated with BCWD survival DAYS in NCCCWA population detected using RAD-SNPs genotyping. (DOCX file)

Table S5
QTL for BCWD resistance that are shared or segregating in both NCCCWA and TLUM populations. (XLSX file)

Table S6
QTL for BCWD resistance that are private to either NCCCWA or TLUM population. (XLSX file)

Table S7
Summary of QTL for BCWD resistance reported in previous studies in rainbow trout populations. (XLSX file)