The stress challenge and collection of blood for measuring plasma cortisol concentration was conducted with the approval of the IUCAC of the USDA/ARS National Center for Cool and Cold Water Aquaculture, protocol #50.

### Mapping population

The QTL mapping population was identified from a broodstock population at the NCCCWA in Leetown, West Virginia, USA. This selective breeding program was initiated in 2002 and designed to select for growth in even year classes and disease resistance in odd year classes
[13]. Stress responsiveness was initially evaluated to determine whether or not it is associated with performance traits affecting aquaculture production efficiency. To this end a total of 584 fish representing 64 families and 9 replicated families (8 fish per family) from the 2002 year class were evaluated for stress response to crowding according to the protocol of Pottinger and Carrick
[9] to identify phenotypic variation of 11.6-93.9 ng cortisol/mL and a significant positive association between plasma cortisol response and growth performance
[12]. In the 2004 year class, 8 fish from each of the top 15% of the growth selected families were similarly evaluated for stress response to identify a heritability of 46%; mean family values were used to develop a high and low responding (P_{1}) generation
[20]. In 2006, high and low responding P_{1} fish were crossed to create a F1 generation consisting of 7 full-sib (FS) families including 222 offspring from 12 P_{1} parents which were evaluated for their potential for use in understanding the genetic control of this trait
[21].

### Stress challenge

The stress challenge method for crowding was modified from Pottinger and Carrick
[9] as described by Weber and Silverstein
[12]. The challenge for the parental generation is described in Weber *et al.*[20] and for the offspring in Vallejo *et al.*[17] Briefly, fish were reared as individual families and were exposed to an artificial ambient photoperiod and reared in continuous flow spring water with temperatures between ~11.5 and 13.5°C and dissolved oxygen near saturation. Fish were placed in 120 l blue polypropylene tanks at approximately one month of age. Seven weeks prior to the start of the experiment, fish were tagged with passive integrated transponder (PIT) tags (Avid Identification Systems Inc., Norco CA) and split by family into tanks containing eight fish each. Fish in the stress study were fed Zeigler Gold (Zeigler Bros. Inc., Gardners PA) at 2% body weight/day. The eight fish from each tank were sampled four times at 4-week intervals when the fish were approximately 160 g body weight. Fish were not fed the day of sampling or the afternoon preceding the sampling. For the crowding stress challenge, fish from a single family were netted and transferred from a 120 l tank to a 6 l or 15 l tank and left undisturbed for three hours. The 6 l or 15 l tanks were used in an effort to keep fish densities in the challenge consistent as fish continue to grow throughout the experiment. After three hours of crowding, the fish were then netted and transferred into an anesthesia bath followed by blood collection. Plasma cortisol was measured (ng/ml) by tritium radioimmunoassay following procedures described by Redding *et al.*[22].

### Stress response derived traits

The parents and offspring fish used in this study (sires = 5; dams = 7; offspring = 222) had four repeated measurements of plasma cortisol recorded at about 4-week intervals as described elsewhere
[17]. These plasma cortisol measurements are highly variable and heavily impacted by environmental effects with a polygenic model heritability of ~0.26
[17]. So, in order to minimize the plasma cortisol variation due to non-genetic factors, and consequently increase the statistical power of QTL detection, and reduce the false positive rate of detected QTL, we derived two stress response phenotypes to use in QTL analyses.

First, we estimated animal breeding value (EBV) using four repeated measurements of plasma cortisol and fitting a mixed inheritance linear model under a Bayesian framework with software iBay version 1.46
[23]. We decided to fit a permanent environmental effect in the repeated measures mixed model analysis to account for the covariance between the records of an individual (i.e., repeated measurements), and capture individual variation across measurements in the EBV computation. The covariates body weight, body length and sexual maturity which had significant effect on the predictive power of the response variable plasma cortisol
[17] were also included in the mixed model to minimize the variance in the sampled population.

Second, we also estimated an animal effect using best linear unbiased prediction (BLUP3) and fitting three repeated measurements of plasma cortisol in a multivariate mixed model with the software ASReml version 2.0
[24]. The mixed model also included covariates indicated above that had significant effect on the predictive power of the response variable plasma cortisol. Here, we estimated heritability of each repeated measurement and genetic correlations between repeated measurements of plasma cortisol with ASReml version 2.0
[24], and noticed that measurement at time 1 was different than measurements at time 2, 3 and 4 in heritability and genetic correlations. So, we decided to calculate an index BLUP3 that was weighted by their relative heritability using measurements at time points 2, 3 and 4, respectively.

In order to determine the effect of using adjusted animal effects for mid-parent genetic effect in the QTL analysis, we performed QTL analysis using HS regression interval mapping as outlined above using offspring animal effects that were adjusted for mid-parent genetic effects. The adjustment of offspring animal effects for mid-parent genetic effect allowed accounting for the effects of relatives, and the use of a pure measure of offspring individual’s genetic value in the QTL analysis
[25, 26].

As expected, we determined that the estimated animal effects EBV and BLUP3 had a Pearson’s correlation of *r* = 0.84 (*P* <0.0001) with SAS Procedure REG
[27]. In this study, we used the estimated animal effects EBV and BLUP3 in the QTL analysis. The use of estimated animal effects such as EBV and daughter-yield deviation (DYD) in whole genome QTL scans has been well documented in livestock species
[25, 28–31].

### Genotyping and linkage analysis

A panel of 412 microsatellite loci identified from the NCCCWA Genetic Map
[19] were used to genotype the parents and offspring of seven FS QTL mapping families (sires = 5; dams = 7; offspring = 222; total marker genotyped fish = 234). Markers were either genotyped using the tailed protocol
[32] or by direct fluorescent labelling (with FAM, HEX, or NED) of the forward primer. Primer pairs were obtained from commercial sources (forward primers labelled with FAM or HEX from Alpha DNA, Montreal, Quebec, Canada, or NED from ABI, Foster City, CA, USA). PCR reactions consisted of 12 μl reaction volumes containing 12.5 ng DNA, 1.5-2.5 mM MgCl_{2}, 1.0 μM of each primer, 200 μM of dNTPs, 1X manufacturer’s reaction buffer and 0.5 units Taq DNA polymerase. Thermal cycling consisted of an initial denaturation at 95°C for 15 min followed by 30 cycles of 95°C for 1 min, annealing temperature for 45 s, 72°C extension for 45 s and a final extension at 72°C for 10 min. PCR products were visualized on agarose gels after staining with ethidium bromide. Markers were grouped in combinations of two or three markers based on differences in fluorescent dye color and amplicon size. Three μl of each PCR product was diluted with 20 μl of water, 1 μl of the diluted sample was added to 12.5 μl of loading mixture made up with 12 μl of HiDi formamide and 0.5 of Genscan 400 ROX internal size standard. Samples were denatured at 95°C for 5 min and kept on ice until loading on an automated DNA sequencer ABI 3730 DNA Analyzer (ABI, Foster City, CA, USA). Output files were analyzed using GeneMapper version 3.7 (ABI, Foster City, CA, USA), formatted using Microsoft Excel and stored in Microsoft Access database. Linkage maps were constructed for each chromosome using marker genotypes from all seven mapping families with the software MULTIMAP version 2.0
[33] using the approach of Rexroad *et al*[19].

### Testing loci for Mendelian segregation distortion

Before QTL analysis, all STR loci were tested for Mendelian segregation distortion (MSD). In outbreed populations, the progeny of an informative QTL mapping family can have any of these marker genotype proportions: 1:1; 1:2:1 and 1:1:1:1. Within individual FS families, the marker genotype counts were performed using a Perl script (Written by G. GAO, unpublished). Then, Chi-Square goodness-of-fit test of marker genotype counts to expected proportions under Mendelian segregation was performed with SAS Procedure FREQ (SAS, 2007) using a default significance level of α = 0.01. Loci with significant MSD were not used in the QTL analysis unless these loci had high quality marker genotypes (i.e., minimum genotyping errors). The rationale to use STR loci that had significant MSD in QTL genome scan is that loci with significant MSD can be linked to viability/survivability and stress response QTL.

### QTL analysis using HS regression interval mapping

For each stress response trait, we performed combined sire-family (five HS families) and dam-family (seven HS families) HS regression analysis, separately, using the web-based software GridQTL
[34]. This software implements a multi-marker approach of interval mapping in HS families as described by Knott *et al.*[35]. This method of QTL analysis does not assume the parents had fixed QTL alleles, instead it relaxes the assumption of fixed QTL allele
[35]. A QTL with a gene substitution effect is fitted at 1-cM intervals along the chromosome using this one-QTL model, *y*
_{
ij
} = *a*
_{
i
} + *b*
_{
i
}
*x*
_{
ij
} + *e*
_{
ij
} where *y*
_{
ij
} is trait score of individual *j* from sire or dam *i; a*
_{
i
} is average effect for HS family *i; b*
_{
i
} is regression coefficient within HS family *i* (substitution effect of QTL); *x*
_{
ij
} is conditional probability for individual *j* within HS family *i* of inheriting allele 1 (or 2); and *e*
_{
ij
} is the residual error.

In the HS regression analysis, the likelihood ratio (LR) test statistic was defined as
where *ln* stands for natural logarithm, Î(*z*) is the likelihood function evaluated at the maximum likelihood estimate (MLE) for the full model that includes polygenic and QTL effects, and Î_{
r
}(*z*) is the MLE for the restricted model under which *r* parameters of the full model are assigned fixed values
[36]. The *P*-value was calculated assuming an *F*-value distributed with numerator DF equal to the number of sires or dams, and denominator DF equal to the total number of offspring minus twice the number of sires or dams
[35]. The chromosome-wide *F*-value (*F*
_{
ChromWide P=0.05
}) and experiment-wide *F*-value (*F*
_{
ExperWide P=0.05}) were estimated using 10,000 permutations with software GridQTL
[34]. The genome-wide significance level (*P*
_{
GenomeWide
}) for detected QTL was estimated as *P*
_{
GenomeWide
} = 1 − (1 − *P*)^{
g
}[37] where *P* is the nominal *P*-value, and *g* = 380 STR loci used in the HS regression analysis. The QTL with *F*-value ≥ *F*
_{
ChromWide P=0.05
} was defined as suggestive QTL (*); and QTL with *F*-value ≥ *F*
_{
ExperWide P=0.05
} or *P*
_{
GenomeWide
} ≤ 0.05 was defined as significant QTL (**). The proportion of phenotypic variance explained by the QTL was calculated as
where *MSE*
_{
full
} and *MSE*
_{
reduced
} are the mean squared error of the full and reduced model, respectively
[35]. The 95% QTL confidence intervals were estimated using 10,000 bootstraps with re-sampling with software GridQTL
[34].

Briefly, the HS regression analysis was performed following these steps: First, all chromosomes were genome scanned for EBV/BLUP3 QTL performing sire-family and dam-family regression analysis, separately, using one-QTL model. At this stage: (a) the chromosome-wide significance threshold level (*F*
_{
chromWideP=0.05}) to declare suggestive QTL was determined using 10,000 permutations; and (b) the 95% confidence interval (CI95) for each suggestive QTL was determined using 10,000 bootstraps with re-sampling. Second, chromosomes with suggestive QTL were re-scanned for additional QTL using one-QTL model that accounted for the effect of already detected QTL until not detecting more QTL in a chromosome. Third, the chromosomes with two detected QTL were re-scanned using two-QTL models. In the chromosomes with two detected QTL: (a) the (*F*
_{
chromWideP=0.05}) significance level for one QTL at a time was estimated while accounting for the effect of the other detected QTL using 10,000 permutations; and (b) the CI95 for each QTL was determined by fixing alternatively the effect of each detected QTL using 10,000 bootstraps with re-sampling
[38]. Fourth, the experiment-wide significance threshold level (*F*
_{
ExperWideP=0.05}) to declare significant QTL was determined by analyzing altogether marker genotype data from 29 chromosomes using one-QTL models and 10,000 permutations
[25, 26].

### QTL analysis using variance components approach

We also performed a genome scan for genetic loci linked to stress response phenotypes (animal BLUP3 index and EBV for plasma cortisol) using a variance components (VARCOMP) method of QTL analysis with software SOLAR version 4.0
[39]. The aim on using a second linkage-based method of QTL analysis was to reduce the rate of false positive claims of QTL. This robust VARCOMP method of QTL analysis was performed using the seven FS families altogether (total parents and offspring *n* =234). This combined family QTL analysis enables accounting for the impact of common environmental effects on QTL mapping which increases the resemblance between full-sibs as they share the same environment during early stages of rainbow trout rearing.

In order to perform multipoint QTL analysis, the multipoint identical by descent (MIBD) relationship matrix for each pair of individuals was estimated with the program LOKI version 2.4.5
[40]. We used LOKI because SOLAR cannot estimate MIBDs in complex pedigrees with double grandparent-grandchild relationships as those used in this study. The VARCOMP approach performs multipoint genome scan at 1-cM intervals. At each testing interval, SOLAR calculates the residual genetic variance or proportion of the total variance due to the polygenic component (*h*
_{
u
}
^{2}); the heritability associated with the QTL or proportion of the total variance due to the QTL (*h*
_{
q
}
^{2}); and the logarithm of odds (LOD) score as *LOD* = *log*
_{10}
*L*(*QTL*)/*L*(*polygenic*)], where *L(QTL)* stands for the likelihood of the full model that includes the polygenic and QTL effect. The genome-wide significance level of detected QTL was estimated using this expression *P*
_{
GenomeWide
} = 1 − (1 − *P*)^{
g
}[37] which assumes the use of sparse genetic maps; here *P* is the nominal *P*-value, and g = 365 STR loci used with VARCOMP analysis (average marker distance = 6.6 cM). The QTL with LOD ≥ 2 was defined as suggestive QTL (*); and the QTL with LOD ≥ 3 or *P*
_{
GenomeWide
} ≤ 0.05 was defined as significant QTL (**).

### QTL analysis using LDLA method

We performed a genome scan for stress response QTL using an LDLA method with the module from the web-based software GridQTL
[34, 41]. We reasoned that the rate of false positive QTL may be reduced by following up the linkage-based QTL analyses with an LDLA-based method of QTL analysis because the QTL signals have to conform with both LD and LA assumptions
[42].

The LR test statistic is defined as indicated above in the section of QTL mapping using HS regression analysis. The LDLA method performs LR tests at 1-cM intervals along the chromosome. Here, the LR test statistic is defined as *LR* = *L*(*QTL*)/*L*(*polygenic*) where *L(QTL)* is the likelihood of the full model that includes the polygenic and QTL effects. The nominal *P*-value was estimated assuming the LR test statistic follows a chi-square distribution with two degrees of freedom
[43]. The genome-wide significance level of detected QTL was also estimated as *P*
_{
GenomeWide
} = 1 − (1 − *P*)^{
g
}[37] where *P* is the nominal *P*-value, and g = 365 STR loci used with LDLA analysis. The QTL with LR ≥ 13.82 or nominal *P* ≤ 0.001 was declared as suggestive QTL (*), and the QTL with *P*
_{
GenomeWide
} ≤ 0.05 was declared as significant QTL (**). The variance components due to polygenic (*σ*
_{
u
}
^{2}), additive QTL (*σ*
_{
a
}
^{2}), dominant QTL (*σ*
_{
d
}
^{2}) and residual error (*σ*
_{
e
}
^{2}) were used to estimate the proportion of total variance due to the polygenic component *h*
_{
u
}
^{2} = *σ*
_{
u
}
^{2}/(*σ*
_{
u
}
^{2} + *σ*
_{
a
}
^{2} + *σ*
_{
d
}
^{2} + *σ*
_{
e
}
^{2})], the proportion of total variance due to additive effect QTL *h*
_{
addQTL
}
^{2} = *σ*
_{
a
}
^{2}/(*σ*
_{
u
}
^{2} + *σ*
_{
a
}
^{2} + *σ*
_{
d
}
^{2} + *σ*
_{
e
}
^{2})], and the proportion of total variance due to dominance effect QTL *h*
_{
domQTL
}
^{2} = *σ*
_{
d
}
^{2}/(*σ*
_{
u
}
^{2} + *σ*
_{
a
}
^{2} + *σ*
_{
d
}
^{2} + *σ*
_{
e
}
^{2})], respectively.