Single-Step Genome-Wide Association Study for Resistance to Piscirickettsia salmonis in Rainbow Trout (Oncorhynchus mykiss).

One of the main pathogens affecting rainbow trout (Oncorhynchus mykiss) farming is the facultative intracellular bacteria Piscirickettsia salmonis Current treatments, such as antibiotics and vaccines, have not had the expected effectiveness in field conditions. Genetic improvement by means of selection for resistance is proposed as a viable alternative for control. Genomic information can be used to identify the genomic regions associated with resistance and enhance the genetic evaluation methods to speed up the genetic improvement for the trait. The objectives of this study were to i) identify the genomic regions associated with resistance to P. salmonis; and ii) identify candidate genes associated with the trait in rainbow trout. We experimentally challenged 2,130 rainbow trout with P. salmonis and genotyped them with a 57 K single nucleotide polymorphism (SNP) array. Resistance to P. salmonis was defined as time to death (TD) and as binary survival (BS). Significant heritabilities were estimated for TD and BS (0.48 ± 0.04 and 0.34 ± 0.04, respectively). A total of 2,047 fish and 26,068 SNPs passed quality control for samples and genotypes. Using a single-step genome wide association analysis (ssGWAS) we identified four genomic regions explaining over 1% of the genetic variance for TD and three for BS. Interestingly, the same genomic region located on Omy27 was found to explain the highest proportion of genetic variance for both traits (2.4 and 1.5% for TD and BS, respectively). The identified SNP in this region is located within an exon of a gene related with actin cytoskeletal organization, a protein exploited by P. salmonis during infection. Other important candidate genes identified are related with innate immune response and oxidative stress. The moderate heritability values estimated in the present study show it is possible to improve resistance to P. salmonis through artificial selection in the rainbow trout population studied here. Furthermore, our results suggest a polygenic genetic architecture for the trait and provide novel insights into the candidate genes underpinning resistance to P. salmonis in O. mykiss.

As in any intensive animal production system, infectious diseases are one of the main threats affecting the success and sustainability of aquaculture (Yáñez ). In the case of salmonid production, one of the major pathogens affecting productivity is the facultative intracellular bacteria Pisciricketssia salmonis, etiological agent of salmonid rickettsial syndrome (SRS). This bacterium was first identified in 1989 in Chile, in a farmed coho salmon (Oncorhynchus kisutch) population (Cvitanich ). Since then, mortalities resulting from SRS have been also identified in Atlantic salmon (Salmo salar) and rainbow trout (Oncorhynchus mykiss) in several countries, such as Scotland, Ireland, Norway and Chile (Fryer and Hedrick 2003). In Chile, SRS was responsible for 20.7, 67.9 and 92.6% of the mortalities associated with infectious diseases in S. salar, O. kisutch and O. mykiss, species respectively (Sernapesca 2018). To date, strategies for P. salmonis control and treatment are mainly based on vaccines and antibiotics. The effectiveness of both approaches has not been adequate (Rozas and Enríquez 2014). Therefore, it has been estimated that economic losses due to SRS mortalities, reached up to US$450 million in Chile in 2012 (Camussetti ). However, variables such as laboratory diagnosis screening expenses or loss of quality of the harvested fish and products were not considered, implying that the economic impact could be even higher.

Therefore, selective breeding could be a feasible alternative to enhance disease resistance; reducing mortality rates caused by P. salmonis, as well as improving animal health and productivity (Bishop and Woolliams 2014: Yáñez and Martínez 2010). However, the main requisite to include a trait into a genetic program is the presence of significant additive genetic variance for that particular trait within the population (Falconer and Mackay 1996). Previous studies estimated heritability values ranging from 0.11 to 0.41 for P. salmonis resistance in Atlantic salmon and coho salmon (Yáñez ; Yáñez ; Barría ). In the case of rainbow trout, Yoshida estimated heritabilities ranging from 0.39 to 0.57 for resistance to P. salmonis using day of death and 0.54 to 0.62 for binary survival as trait definitions. Altogether, these results demonstrate the possibility of improving this trait by means of artificial selection in different salmonid species.

The development of next generation sequencing technologies has facilitated the identification of thousands of single nucleotide polymorphisms (SNPs) segregating along the genome of several animals, including aquaculture species (Yáñez ). Thus, using a genotyping by sequencing (GBS) approach in conjunction with genome-wide association studies, some authors evaluated genomic regions associated with resistance to bacterial infections in aquaculture species (Liu ; Palti ; Palaiokostas ; Barría ). However, in salmonid species, the use of SNP panels has been the most used alternative for genotyping a high number of individuals with thousands of genetic variants simultaneously. This has been made simpler by the development of high density SNP arrays for Atlantic salmon (Houston ; Yáñez ) and rainbow trout (Palti ). The use of these SNP panels have also allowed the comparison of the accuracy of estimated breeding values (EBV) using genomic selection to pedigree-based genetic evaluations for resistance to infectious diseases in Atlantic salmon (Ødegård ; Tsai ; Bangera ; Correa ), coho salmon (Barría ) and rainbow trout (Vallejo ; Vallejo ; Yoshida ; 2018b). SNP arrays have also enabled the dissection of the genetic architecture of resistance to bacterial diseases in salmonids. For instance, genomic regions and candidate genes associated with resistance to P. salmonis in Atlantic and coho salmon (Correa ; Barría ), and bacterial cold water disease (BCWD) in rainbow trout (Vallejo ) have been identified.

To date there are no studies aimed at identifying genomic regions or candidate genes associated with resistance to P. salmonis in rainbow trout populations. Therefore, the main objectives of the current study were to i) identify genomic regions associated with resistance to P. salmonis in a farmed rainbow trout population, and ii) identify candidate genes associated with the trait.

Materials and Methods

Population and experimental challenge

The population used in this study was comprised by a rainbow trout (Oncorhynchus mykiss) broodstock (year-class 2011), owned by Aguas Claras (Puerto Montt, Chile) and belonging to a genetic improvement program run by Aquainnovo (Puerto Montt, Chile). This population was artificially selected for growth, appearance-related traits and carcass quality for three generations. For a detailed description about rearing conditions and population management please see Flores-Mara , Rodríguez and Neto .

Fish from 105 full-sib families (48 half-sib families) with an average weight of 7.0 ± 1.5 g, were PIT-tagged for individual traceability of families. After tagging, fish were maintained in a single tank until they were transferred to Aquainnovo’s Aquaculture Technology Center Patagonia in August 2012. Fish were acclimated for 20 days in a 15m3 tank, prior to experimental challenge. A random sample of fish were selected to evaluate the sanitary status of the population by means of qRT-PCR for Infectious Salmon Anemia virus (ISAV), Infectious Pancreatic Necrosis virus (IPNV), and Renibacterium salmoninarum, and culture for Flavobacterium spp. Later, a total of 2,130 juveniles (with an average of 23 individuals per family and ranging from 17 to 27 fish per family), were intraperitoneally (IP) injected with 0.2ml of a lethal dose (LD50) of a LF-89 strain of P. salmonis inoculum. Post injection, fish were equally distributed into three different tanks, considering similar family distribution into each replicate (with 5 to 9 fish per family in each tank). Environmental parameters were measured throughout the experiment and the challenge continued until the mortality curve showed a plateau. Daily mortality was recorded, and body weight was measured for each fish at time of death or at the end of the experiment (FW). Surviving fish were killed and body weight was also recorded. Fin clips from all fish were sampled and stored in 95% ethanol at -80° until they were genotyped.

Genotyping

The genomic DNA from the sampled fin clips was extracted using a commercial kit (DNeasy Blood & Tissue Kit, Qiagen), following the manufacturer’s instructions. Genotyping was performed using a commercial 57K SNP array (Affymetrix AxiomTM myDesignTM) developed by the National Center for Cool and Cold water Aquaculture at USDA and Aquagen (Palti ).

Quality control (QC) was assessed through Affymetrix’s Axiom Analysis Software, using default settings. Then, a second QC using Plink software (Purcell ) was applied to remove SNPs with a genotype call rate lower than 0.90, minor allele frequency (MAF) < 0.01 and deviated from Hardy-Weinberg Equilibrium (P < 1x10−6). Individuals with a call rate lower than 0.90 were also removed from further analyses.

Trait definition

Resistance to P. salmonis was defined as time to death (TD), measured in days, with values ranging from 1 until the end of challenge test. Additionally, resistance to P. salmonis was also defined as binary survival (BS), with a value of 1 or 0 based on if the fish died or survived until the end of the challenge.

Genome-Wide association study

A single-step GWAS (ssGWAS) analysis was performed to identify genomic regions associated with resistance to P. salmonis. This approach considered fish with both phenotypes and genotypes and also individuals with phenotypes but no genotypes in the analysis (Wang ). The pedigree and genotypic data in ssGWAS are connected through the H matrix. Thus, the H matrix combines both the pedigree and the genomic relationship matrices (Aguilar ). Thus, the inverse of the H matrix is:where is the inverse of the numerator relationship matrix, considering all the phenotyped animals, is the inverse of the pedigree-based relationship matrix considering only the genotyped animals, and is the inverse genomic relationship matrix. The following model was used for GWAS analysis:where y is the vector of phenotypes (for TS and BS), is the vector of fixed effects (tank as factor and final body weight as a covariate), a is the vector of random animal effects (including the genotype information), e is the vector of residuals, and and are the incidence matrices for fixed and random effects, respectively.

The variance of and are estimated as follow:where is the identity matrix and and are the genetic additive and residual variances, respectively.

A linear model (AIREML) and a threshold model (THRGIBBS1F90) were used for TD and BS, respectively. Both models were fitted using BLUPF90 programs (Misztal ). For the latter, a total of 200,000 Markov Chain Monte Carlo (MCMC) iterations were used, the first 20,000 were discarded as burn-in iterations and from the remaining 180,000 samples, we saved one from every 50. Therefore, the analyses included 3,600 independent samples.

For TS and BS, the heritability was estimated as follows:where is the additive genetic variance estimated using the H matrix, and is the residual variance.

The genetic correlation between both resistance definitions was calculated using the following formula indicated by Falconer and MacKay (1996):where is the additive genetic covariance between resistance as TD and as BS, while and correspond to the additive genetic variance of TD and BS, respectively.

To identify genomic regions associated with each trait, we estimated the percentage of the genetic variance (PGV) explained by windows of 20 adjacent SNPs. Then, if a 20 SNP window explained more than 1% of the PGV, we considered that region as associated with resistance to P. salmonis, as proposed by Gonzalez-Pena ).

Candidate genes

The candidate genes were identified by searching 500kb up and downstream from the SNP explaining the highest PGV within each associated window. For this purpose, we used the last version of the Oncorhynchus mykiss reference genome (GCA_002163495.1). The criteria for selecting candidate genes lies in the function of the protein that encodes each gene found, mainly related to immune response, DNA repair, stress response and similar pathways.

Linkage disequilibrium

Whole genome linkage disequilibrium (LD) was calculated using Plink v 1.09 and measured as Pearson’s squared correlation (r2) among SNP genotypes. The pair-wise LD was estimated based on the formula proposed by Hill and Robertson (1968) and estimated for each SNP pair. Based on the physical distance among these SNP pairs, we created bins of 100kb. The decay and extent of the LD was visualized by plotting the mean r2 for each bin, from 0 to 10Mb.

The LD within 1Mb window of specific genomic regions (500kb upstream and 500kB downstream from the SNP of interest) was estimated as the average r2 among all the SNPs located in the window.

Data availability

Raw genotypes, phenotype data and the location of each SNP used in this study in the Oncorhynchus mykiss genome are available from the online repository figshare (https://doi.org/10.6084/m9.figshare.8146076). Table S1 containing all the genes located within 1Mb windows surrounding the SNPs explaining the highest proportion of genetic variance is available at the same link mentioned above.

Results

Descriptive statistics and heritabilities

Summary statistics for resistance to P. salmonis measured as TD and BS, and for FW are shown in Table 1. The first death was recorded on day 10 post intraperitoneal injection; the last on day 32. Average TD was 23.26 (SD = 7.86) days. At the end of the experimental challenge the proportion of non-survivor fish, taking into account the 105 challenged families, was 0.59 (SD = 0.49). Furthermore, the cumulative mortality ranged from 7.7 to 100%, indicating considerable phenotypic variation for resistance to P. salmonis among families in the rainbow trout population. Cumulative mortality within each replicate tank was 59.4, 65.1 and 64.7%. Mortality peaked on days 12, 15 and 19 post injection. Average final body weight was 173.80 (SD = 52.27) g. This trait ranged considerably among challenged fish, with a minimum of 46.10g and maximum 448g.

Table 1

Summary statistics for time to death (TD), binary survival (BS) and final weight (FW) measured in 2,130 rainbow trout individuals after an experimental challenge against Piscirickettsia salmonis

Trait	Mean	SD	CV(%)	Min	Max
TD (days)	23.26	7.86	33.27	10	32
BS	0.59	0.49	0.83	0	1
FW (g)	173.80	52.27	30.07	46.10	448

Variance components for TD and BS are shown in Table 2. Significant heritability values were estimated for both trait definitions. Thus, 0.48 0.04 and 0.34 0.04 were estimated for TD and BS, respectively. Furthermore, a high genetic correlation was found between both traits (-0.96 0.01).

Table 2

Genetic parameters and heritabilities for resistance to Piscirickettsia salmonis defined as time to death (TD) and binary survival (BS)

Trait	σa2a	σe2b	h^2(SE)c
TD	25.95	28.92	0.48(0.04)
BS	6.27x10−2	1.21x10−1	0.34(0.04)

Additive genetic variance.

Residual variance.

Heritability and standard error.

Genome-wide association study

From all genotyped animals, 2,047 passed quality control (representing 97.10% of the total). A total of 26,068 SNPs remained in the set for further analyses (∼64.68%). Figure 1 shows the Manhattan plot for the proportion of variation explained for resistance to P. salmonis measured as TD and BS. We identified four genomic regions associated with resistance as TD. These regions were located on Omy03, Omy14, Omy24 and Omy27. For BS, we identified three genomic regions associated with the trait. These were found on Omy05, Omy27 and Omy30. The genomic region located on Omy27 was found to be associated with resistance to P. salmonis for both TD and BS. In both cases, this common genomic region explains the highest proportion of genetic variance for each trait, with 2.4 and 1.5% for TD and BS, respectively. The SNP explaining the highest proportion of the genetic variance (Affx-88923370) is the same for both TD and BS (Table 3).

Figure 1

Genome-wide association analysis for resistance to Piscirickettsia salmonis in rainbow trout (Oncorhynchus mykiss). Resistance was defined as time to death (A) and as binary survival (B).

Table 3

Top markers associated with Piscirickettsia salmonis resistance defined as TD and BS in rainbow trout

Ranking	Name	Chra	Pos (Bp)	Alleleb	-log(pval)	PGVc	Genesd
Time to death
1	Affx-88923370	27	9998276	A / G	5.65	2.43	usp2, nlrc3, tap, pitpna
2	Affx-88916453	3	14818380	T / C	1.56	1.41	stl2, aicda, il11, gstk1
3	Affx-88922612	14	10975036	T / G	3.12	1.21	tlr4, tax1bp1, satb1
4	Affx-88927397	24	11828385	C / A	1.55	1.02	a2m, pou2af1, nfkbiz
Binary survival
1	Affx-88923370	27	9998276	A / G	NAe	1.50	usp2, nlrc3, tap, pitpna
2	Affx-88951679	5	68055053	T / C	NA	1.12	faslg, prdx6, plpp6
3	Affx-88908715	29	32519588	C / T	NA	1.01	mapk12, gls

Chromosome.

Resistant/Susceptible allele.

Percentage of genetic variance.

Summary of the genes located within 1Mb window.

Not assessed.

Using the O. mykiss reference genome (GCA_002163495) we identified candidate genes associated with resistance to P. salmonis. Table 3 shows a summary of the genes located proximate to the SNPs explaining the highest proportion of the genetic variance within each genomic region.

Among the candidate genes flanking the most important SNP on Omy03 for TD, we found Gluthatione S-transferease kappa 1 (gstk1) and Interleukin-11 (il11). These genes are involved in the response to oxidative stress and the immune response to bacterial infections, respectively (Oruc et al. 2004; Wang et al. 2005). On Omy14, we found the Toll-like receptor 4 (tlr4) gene, which has been suggested to act as a bacteria sensor (Palti 2011). On Omy24, we found alpha-2-macroglobulin-like (a2m), which is part of a broad-spectrum protease inhibitor, and it has been suggested to play a role in the defense against Cryptobia salmositica in rainbow trout (Zuo and Woo 1997).

Also, we found POU class 2 associating factor 1 (pou2af1), which has been described as a coactivator of transcription factors that regulate immunoglobulin expression of B cells (Teitell 2003), and NF-kappa-B inhibitor zeta-like (nfkbiz), a regulator of pathogen recognition, phagocytosis and production of cytokines by dendritic cells (Rozas-Serri et al. 2017).

For BS, on Omy05 we found fas ligand (faslg), whose protein has been suggested as an important mediator of anti-bacterial innate immune response, by inducing apoptosis of target cells and recruiting phagocytic cells (Kaur et al. 2004). On the same chromosome we found Peroxiredoxin-6-like (prdx6), one of the six different isoforms that conforms the peroxiredoxins group, which are antioxidants proteins that protect cells from oxidative damage and is likely to be involved in protective response against a bacterial infection in Scophthalmus maximus (Zheng et al. 2010).

On Omy29, MAPK12 was found; previous studies described that MAPK12 is involved on the signaling pathways responsible for TNF-α secretion from rainbow trout macrophages (Roher et al. 2011). Glutaminase kidney isoform, mitochondrial-like (gls) was also found on Omy29, which is part of a family of enzymes that play a role in nucleotide, amino acid and urea biosynthesis (Kumada et al. 1993).

On Omy27 we found genes related with innate immune response regulation and some molecules related with metabolic processes and apoptosis. However, the SNP explaining the highest proportion of genetic variance is located within an exon of the gene Smoothelin protein 2 (Smtnl2) which remains poorly characterized both in humans and fishes, but it is believed to have a role in actin cytoskeleton organization.

The complete list of genes located within the 1Mb window flanking the SNPs explaining the highest proportion of genetic variance for resistance to P. salmonis, is shown in Table S1.

The LD was estimated for each window containing SNPs associated with P. salmonis resistance. Table 4 shows the number of SNPs along the 1Mb windows associated with the trait, which varied from 7 (Omy03) to 32 (Omy29). Regarding the average r2, estimations varied from 0.21 to 0.69 for almost all the associated windows. However, the genomic regions located on Omy03 and Omy14 showed a value of 0.14 and 0.18, respectively, below the 0.2 threshold suggested by (Meuwissen ).

Table 4

Estimated average linkage disequilibrium (measured as r2) for the regions flanking the top markers associated with P. salmonis resistance as time to death (TD) and binary survival (BS)

Ranking	Name	Chra	Window size (Bp)	Number of SNPsb	r2
Time to death
1	Affx-88923370	27	1189027	30	0.21
2	Affx-88916453	03	1037096	07	0.14
3	Affx-88922612	14	972244	29	0.18
4	Affx-88927397	24	1283717	24	0.32
Binary survival
1	Affx-88923370	27	1189027	30	0.21
2	Affx-88951679	05	985641	22	0.69
3	Affx-88908715	29	1020166	32	0.24

Chromosome.

Number of SNPs markers within the window.

Discussion

In the current study we show significant genetic variation for resistance to P. salmonis in a farmed rainbow trout population. A moderate to high heritability was estimated for resistance as TD (0.48) and BS (0.34). These estimates are higher than those reported in previous studies carried out for resistance to other bacterial diseases in aquaculture species, with heritabilities ranging from 0.22 to 0.38 (Ødegård ; Palaiokostas ; Vallejo ). In the case of P. salmonis resistance, several studies have evaluated the presence of genetic variation in different salmonid species. Thus, similar estimates have been shown for Atlantic salmon, when using pedigree or genomic data, with values ranging from 0.19 to 0.39 (Yáñez ; Yáñez ; Correa ; Bangera ). In the case of coho salmon, heritability estimates range from 0.16 to 0.27 when resistance is defined as a linear or binary trait (Yáñez ; Barría ).

Recent studies in rainbow trout, using pedigree and genome-based genetic evaluation approaches, estimated heritabilities ranging from 0.39 to 0.57 for TD and from 0.54 to 0.62 for BS (Yoshida ); values which are within the range of our estimations. Moreover, our results suggest a higher effect of the additive genetic component on the phenotypic variance for resistance to P. salmonis in rainbow trout when compared to S. salar and O. kisutch, which would imply potentially faster genetic progress for the improvement of resistance to P. salmonis by means of artificial selection in the rainbow trout population used in the present study.

The current rainbow trout population was founded by using 18 sires and 48 dams in 2005 (Bassini ; Flores-Mara ), representing more animals than those that were used for establishing a coho salmon breeding population in Chile in 1997 (Barría ; Dufflocq ; Yáñez ). Thus, the higher heritability estimates for P. salmonis resistance in rainbow trout than in coho salmon, could be explained by i) a major genetic variability due the higher number of animals and ii) a fixation of the alleles on the coho salmon population as result of a higher number generations under artificial selection.

The effect of the genetic architecture of a trait (among other variables) on the accuracy of breeding values obtained through genomic selection (GS) is widely known (Daetwyler ; Goddard 2009). Previous studies in salmonid species (Atlantic salmon and coho salmon), suggest that resistance to P. salmonis is a polygenic trait (Correa ; Barría ). Based on the 26K SNPs which passed QC, our study similarly suggests a polygenic nature for resistance to P. salmonis resistance in rainbow trout (i.e., several loci involved with the trait, with a small effect each). Thus, it is expected that, when comparing genome-based Best Linear Unbiased Predictor (GBLUP) method against a Bayesian approach (e.g. Bayes C), the former would show increased accuracy of the estimated breeding values over the latter (Habier ; Hayes ) for the current rainbow trout population. Nonetheless, as predicted by Yoshida this was true only at low SNP densities (i.e., 0.5 to 10 K). When 20K and 27K were used, Bayes C outperformed GBLUP accuracies. The authors suggested that this could be due to an oligogenic architecture of the resistance trait, or that Bayes C had higher effectiveness in capturing the linkage disequilibrium between the SNPs and a QTL when more SNPs were used.

Although a high genetic correlation between TD and BS was found, only one common genomic region surpassed the 1% threshold for both traits (Omy27). The relatively low correspondence between the two traits in terms of the identified QTL might be partially explained based on the different approaches used for estimation of QTL. While for TD we fitted a linear model, for BS we fitted a threshold model. Thus, QTL detected in the observed scale in the case of TD might not be detected on the underlying scale for BS and vice versa. However, there are several genomic regions in common (e.g., Omy13, Omy19, Omy26) below the 1% threshold for TD and BS, and other in which the threshold was surpassed only in one trait (e.g., Omy03, Omy05, Omy14). Thus, increasing the statistical power might allow increasing the number of common regions associated to P. salmonis resistance.

Based on the underlying genetic architecture of the resistance trait in the current population, here we used a single-step GWAS approach which converts the genomic breeding values of genotyped animals obtained by single step genomic evaluation into SNP effects (Wang ). This approach uses available pedigree data and besides including animals with both phenotype and genotype data, it also allows including animals with only phenotype records (those animals with missing genotype). This method is also adapted to estimate SNP effects and variances with fast computing times, generating robust estimates with practical simplicity. The use of 20 SNPs sliding windows may improve the accuracy and precision on identifying QTL when compared with a SNP-by-SNP approach (Wang ; Zhang ).

The use of a uniform-sized sliding-window approach, instead of using adjacent windows, allowed us to accounting for the linkage disequilibrium among the markers. The selection of an optimal window size is complicated and it depends on the population under study (Tang ). A 20 SNPs per window has been frequently used for the identification of genomic regions associated with different traits in different rainbow trout breeding populations (Gonzalez-Pena ; Neto ). Although Gonzalez-Pena used adjacent windows, the authors found that using a 20-SNPs sized window showed the lowest signal-to-noise-ratio than when using higher SNPs number. Thus we used a similar windows size for comparative purposes with previous studies in rainbow trout populations.

As showed by Wang , power and precision of ssGWAS, evaluated by the correlation between true simulated QTL effects and the sum of m adjacent SNPs, is high (up to 0.81 ± 0.02) when using 1500 genotyped animals, indicating a high correspondence between true and estimated QTL. Here, we used 2047 fish with genotypes, which may be generating a similar or somewhat higher statistical power and precision than the one presented by Wang . This accuracy was maximum when m = 8 SNPs and decreased sharply when m = 40 SNPs. When m = 16 SNPs the estimated accuracy reached up to 0.80 ± 0.03. Thus, we assume that the likelihood of the identified QTL being true positives can be also considered high (∼0.80). Furthermore, Zhang suggested that using SNP windows is a better approach than using single SNP, given that the true number of QTL is unknown. Thus, we believe that using a 20 SNPs windows-size represents an appropriate approach for identifying genomic regions and genes likely associated with resistance to P. salmonis.

Resistance to bacterial infections implies a modulation of the host immune response to inhibit or reduce the replication rate of the pathogen (Doeschl-Wilson and Kyriazakis 2012). The infection process caused by P. salmonis uses clathrin for internalization and then the actin cytoskeleton for vacuole generation (Ramírez ). Similar pathways have been observed in other mammalian intracellular gram-negative bacteria (Manon ; Valencia-Gallardo ). Within the region associated with TD on Omy03 we identified a gene coding for the receptor DC-SIGN related with the immune response and expressed on macrophage and dendritic-cell surfaces (Ahmed ). It has been previously described that Mycobacterium tuberculosis, interferes with the Toll-like receptor signaling by DC-SIGN, inhibiting interleukin-12 production (Gorvel ), a proinflammatory cytokine, which plays a key role in the performance of phagocytes in teleost fish (Álvarez ).

As mentioned before, endocytosis mediated by clathrin is the main pathway used by P. salmonis for cell invasion. Clathrin recruits, among other cell components, AP-2; which is regulated by NECAP-1 (Ritter ), a gene flanking the SNP explaining the highest proportion of genetic variance in Omy03 for resistance measured as TD. Similarly, on this chromosome we also found the gene glutathione S-transferase kappa 1 (gstk1), which is a member of the glutathione S-transferase family (GST), involved in cellular detoxification, and expressed in cells to reduce oxidative stress-related damage (Morel and Aninat 2011), a consequence of P. salmonis infection (Rozas and Enríquez 2014), and differentially expressed in Atlantic salmon after P. salmonis exposure (Rise ). A candidate gene related to resistance measured as BS found on Omy05 is, the fas ligand gene (faslg), a member of the TNF superfamily. The Fas/FasL pathway is essential for immune system regulation, including apoptosis induced by T cell activation and cytotoxic T lymphocytes (Siegel ).

For both resistance definitions, the same chromosome and identical SNP was identified as the marker explaining the highest genetic variation for the trait, which makes this QTL an interesting region. Within this region we found the gene phosphatidylinositol transfer protein alpha (pitpna), which belongs to the phosphatidylinositol family (ptdlns) (Piscatelli ), responsible for phospholipid transfer between cellular membranes (Thornbrough ), which in turn are regulators of cell signal transduction, membrane trafficking and cytoskeleton organization (Hilbi and Haas 2012). The latter process is affected by P. salmonis once inside the macrophages (Ramírez ). Similar to P. salmonis, Legionella pneumophila also replicates inside macrophages, and manipulates the vesicle generation inside the cell by joining with ptdlns 5 (Hilbi and Haas 2012).

Additionally, in this region we found the gene nlr family card domain containing 3 (nlrc3). Previously, Álvarez described differential expression of nlrc3 in rainbow trout in response to bacterial lipopolysaccharides (lps), specifically in the skin, liver and gills. This pattern has also been observed in Atlantic salmon during an infection with P. salmonis (Tacchi ), and is therefore a potential mechanism used by this bacteria to evade the immune response.

The gene tapsain (tap) is also involved in the immune response, transporting cytosolic peptides generated by the proteasome to load on MHC class I (Procko ). On Omy27, we found a gene that encodes a protein related to tapsain (TAPBPR), which negatively regulates tap; generating a reduction in immune response efficiency (Boyle ).

Similar to what has been found in a farmed Atlantic salmon population (Barría ), LD of 0.2 was estimated at minimum marker distance of ∼50Kb (Supplementary Figure 1). Based on O. mykiss genome length of 2.2Gb, we estimate that at least 44K SNPs would be necessary to perform a proper whole genome association scan. Despite the 26K SNPs used in the current study is below this number, it shed light into the genomic regions likely associated with P. salmonis resistance in a farmed rainbow trout population.

Among the genomic regions of interest, a high variation in average LD was estimated (from 0.14 to 0.69). However, it should be noticed that these values comprises a 1Mb window-size. Similarly, Vallejo found high levels of LD ranging between 0.21 to 0.44 across the rainbow trout genome. The candidate genes associated with resistance to P. salmonis are at a smaller distance with respect to the associated SNP, which in turn is located in the midpoint of the interval. Thus, a higher correlation between the SNP genotypes and the putative genes are expected within the evaluated genomic regions.

We expect that in the near future, the identification and validation of causative mutations affecting some of the candidate genes presented here, by means of functional studies, will provide a better understanding of resistance against this and other infectious diseases in rainbow trout and other salmonid species. These studies will be facilitated through international collaborative initiatives such as the Functional Annotation of All Salmonid Genomes, FAASG (Macqueen ).

Conclusions

To the best of our knowledge this is the first report identifying candidate genes related to resistance to P. salmonis in a farmed rainbow trout population. Genes likely related with resistance were identified close to SNPs explaining the highest proportion of genetic variance. Furthermore, we identified one common genomic region associated with resistance using both a linear and binary trait. Our results show that this trait is controlled by multiple genes each with a small effect. Therefore, a genomic selection approach is suggested as the best method to improve this trait by means of artificial selection.