Genic and non-genic SNP contributions to additive and dominance genetic effects in purebred and crossbred pig traits

Ethics approval

Data used for this study was originated in the UtOpIGe project ANR-10-GENOM_BTV-015. All experimental protocols of the project were approved on the 01/23/2013 (R-2012-NM-01) by the local ethic committee (Comité Rennais d’Ethique en matière d’Expérimentation Animale). The Comité Rennais d’Ethique en matière d’Expérimentation Animale is registered in the French National Committee of ethical reflexion on animal experimentation number 7 (Comité National de Refléxion Ethique sur l’Expérimentation Animale). All methods were carried out in accordance with guidelines and the French regulation in Animal research (articles R214-87 to R214-137 of the French rural code https://www.recherche-animale.org/sites/default/files/c_rural_2013.pdf) updated by the 2013–118 decree and five orders from February 1st 2013, published on February 7th, according to the 2010/63 directive from the EU. This regulation is under the responsibility of the French Ministry of Agriculture.

The study was carried out in compliance with the ARRIVE guidelines (http://www.nc3rs.org.uk/page.asp?id=1357).

Animal material

Animals were provided by the three French breeding companies of the Alliance R&D association (composed of Axiom, Choice Genetics France, Nucléus and IFIP) involved in the UtOpIGe project ANR-10-GENOM_BTV-015. We used 636 purebred Piétrain (PB) and 720 crossbred Piétrain × Large White (CB) entire “intact” male piglets produced on selection and multiplication farms and tested at a single test station INRA UE3P France Génétique Porc phenotyping station (UE3P, INRA, 2018. Unité expérimentale Physiologie ET Phénotypage des Porcs, France, https://doi.org/10.15454/1.5573932732039927E12). Both, PB and CB animals were descendants of 90 Piétrain boars. They entered the test station facilities in Le Rheu (France) at approximately 5 weeks of age and were slaughtered at a fixed weight of 110 kg (at 5–6 months of age).

Phenotypes

Animals were misspelled weighed at the beginning (when the animals reached approximately 35 kg) and at the end of the test period (when they reached 110 kg). Average daily gain (ADG, kg/day) was calculated as the body weight gained during the test period divided by the duration of the period. Rib back fat thickness (BFT, mm) was measured on carcass with the Capteur Gras Maigre method¹⁹. At the slaughterhouse, carcasses were chilled in a cooling room at 4 °C for 24 h. Ultimate pH of the semimembranosus dorsi muscle (PHS, pH units) was measured using a Xerolyt electrode (Mettler-Toledo, Australia) and a Sydel pH meter (Sydel, France) at 24 h post- mortem. Further details regarding complete data collection and experimental design can be found in Tusell et al.¹⁰.

Genotypes

Animals were genotyped using the Illumina Porcine SNP60 Bead Chip (Illumina, Inc., San Diego). Single nucleotide polymorphisms (SNPs) with a call rate lower than 0.90 and a minor allele frequency lower than 0.05 were removed. For the remaining SNPs. Missing genotypes were imputed using a naïve method that sampled genotypes with probability weights of the allele frequencies at each locus. The missing genotypes were then replaced with these sampled genotypes. Animals that presented Mendelian inconsistencies with their parents were discarded. After quality control, 46,816 SNPs were included in the analyses. Due to separate data edition, the number of animals with records and number of SNPs retained for the analyses slightly differed for each trait. Summary statistics of the three analyzed phenotypes in the two populations are presented in Table 1.

SNP annotation

Chromosome information of SNPs (i.e. the map file containing SNP-ID and RS-Numbers) was downloaded from the Animal Genome Database (https://www.animalgenome.org/repository/pig/). The SNPs that did not have RS number or a non-unique RS number were discarded leading to 39,727 SNPs available for the analysis. The physical positions of the SNPs obtained through Ensembl database (https://useast.ensembl.org/Sus_scrofa/Tools/VEP) for pig (Sus Scrofa, Assembly: Sscrofa11.1, accession date: December 2018) allowed to locate each SNP into a genic or an intergenic region.

SNPs were classified in three categories. Genic region category (Genic, 19,672 SNPs) encompassed the SNPs annotated into the following categories: introns (15,824 SNPs), synonymous (321 SNPs), upstream the gene (1,595 SNPs), downstream the gene (1,338 SNPs), 5′ (77 SNPs) and 3′ untranslated regions (UTR) (382 SNPs), missense (108 SNPs), other exon mutations (27 SNPs). Intergenic region category (Intergenic, 20,055 SNPs) encompassed all SNPs annotated outside these genic regions. All region category (All) included the 39,727 SNPs.

Statistical analysis

Genic and intergenic additive and dominance variances in purebred and crossbred populations

To explore the additive and dominance variance explained by the SNPs of the categories (i.e., Intergenic, Genic and All SNPs) in ADG, BF and PHS in the two populations, the following general univariate model was fitted separately for each trait, population and genomic region^9,14:

where (mathbf{y}) was the phenotypic value of individuals, ({varvec{upbeta}}) is a vector of systematic effects and (mathbf{p}) a vector of pen nested within batch random effects. Term (mathbf{e}) is the vector of residual effects. Terms ({mathbf{X}}_{1}) and ({mathbf{X}}_{2}) are incidence matrices that assign systematic and nested within batch effects to the phenotypes, respectively. Model for ADG included the effects of weight at the beginning of test (covariate), and pen nested within batch effect (random effect, 65 levels). Model for BFT included the effects of hot carcass weight (covariate) and pen nested within batch effect (random effect, 62 levels). Model for PHS included the effects of hot carcass weight (covariate) and date of slaughter (systematic effect, 43 levels). Term (mathbf{f}) is a vector of inbreeding coefficients calculated as the average homozygosity per individual and (b) is the inbreeding depression coefficient¹⁴. Terms (mathbf{a}) and (mathbf{d}) are the vectors of animal additive and dominance genotypic effects, respectively. Terms (mathbf{Z}) and (mathbf{W}) are incidence matrices relating additive and dominance genotypic effects to either PB or CB animals with −1, 0, 1 (additive) and 0, 1, 0 (dominance) values for the AA, Aa and aa genotypes, respectively. Additive and dominance genotypic (co)variances were modelled as (mathbf{G}=frac{mathbf{Z}{mathbf{Z}}^{boldsymbol{^{prime}}}}{left{tr([mathbf{Z}{mathbf{Z}}^{boldsymbol{^{prime}}}])/nright}}{sigma }_{A*}^{2}) and (mathbf{D}=frac{mathbf{W}mathbf{W}}{left{tr([mathbf{W}{mathbf{W}}^{boldsymbol{^{prime}}}])/nright}}{sigma }_{D*}^{2})¹², where ({sigma }_{A*}^{2}) and ({sigma }_{D*}^{2}) are the estimated variance components, and (n) the number of animals.

A Bayesian framework was adopted for inference. The prior distributions for the parameters of the model were (P({varvec{beta}},b)sim k), (Pleft({varvec{p}}|{sigma }_{p}^{2} right)sim N(0, {mathbf{I}sigma }_{p}^{2})), (P(mathbf{a}|{mathbf{G}, sigma }_{A*}^{2})sim N(0, mathbf{G}{sigma }_{A*}^{2})), (P(mathbf{d}|{mathbf{D}, sigma }_{D*}^{2})sim N(0, mathbf{D}{sigma }_{D*}^{2})) and (Pleft({varvec{e}}|{sigma }_{e}^{2}right)sim N(0, {mathbf{I}sigma }_{e}^{2})) where (k) is a constant, (mathbf{I}) is an identity matrix, and ({sigma }_{p}^{2}) and ({sigma }_{e}^{2}) are the nested within batch and residual variances, respectively.

The variance components were estimated either using all the 39,727 available SNPs (All), the 19,672 SNPs located in the genic regions (Genic) or the 20,055 SNPs located in the intergenic regions (Intergenic). Thus, for each genomic region, the SNPs included in (mathbf{Z}) and (mathbf{W}) differed according to the genomic region used (i.e., All, Genic or intergenic SNPs).

Following¹², the variance components ({sigma }_{Aboldsymbol{*}}^{2}) and ({sigma }_{D*}^{2}) estimated in the genotypic models were then used to retrieve the additive and dominance SNP variances that, together with the allelic frequencies of each population, allowed to obtain the additive and dominance deviation variances for the two populations across the three analysed traits. Hence, with the SNPs of the three regions, three different models (using either Genic, Intergenic or All SNPs) were implemented per trait (ADG, BFT and PHS) and per population (i.e. PB or CB) leading to 18 different models.

Predictive ability

Predictive ability of a model including only additive genetic effects and an inbreeding coefficient was compared to a model including additive and dominance effects and an inbreeding coefficient. These two models (additive model or additive and dominance model) were run separately for each combination of population (PB or CB), trait (ADG, BFT and PHS) and genomic region (Genic, Intergenic and All markers) leading to 36 different models. All models included the same systematic and non-genetic random effects described above for each of the three traits. Predictive ability of the models was evaluated by cross-validation (CV). Specifically, a four-fold CV scheme was used by attributing animals randomly to one of four separate subsets. From these four subsets, three folds were combined to create a training set and the remaining fold was used as testing set. Each of the four subsets was applied as a testing set only once. Because of the small size of the sample, the four-fold CV was replicated 10 times at random, and results were averaged over replications²⁰. Predictive abilities were assessed via Pearson’s correlation between pre-adjusted phenotypes and predicted phenotypes in the testing sets.

Additive and dominance genotypic correlations between PB and CB populations

Additive and dominance genotypic correlations between PB and CB performances were estimated using all SNPs available and separately for each trait. In this case, PB and CB performances of the trait (i.e., ADG, BFT or PHS) were considered to be different traits in PB and CB populations and were jointly analyzed with a bivariate genotypic model accounting for additive and dominance effects¹² and a genomic inbreeding coefficient¹⁴. The following general two-trait model was applied:

Index (k) (for (k=PB, CB)) is used to denote either the PB ((k=PB)) or the CB ((k=CB)) populations. Term ({mathbf{y}}_{k}) is a vector of phenotypes, ({{varvec{upbeta}}}_{k}) is a vector of systematic effects, ({mathbf{p}}_{k}) is a vector of pen nested within batch effects (only included in ADG and BFT models), ({mathbf{f}}_{k}) is a vector of inbreeding coefficients and its corresponding inbreeding depression coefficient (({b}_{k})), and ({mathbf{e}}_{k}) is a vector of residual effects. Terms ({mathbf{X}}_{1,k}), ({mathbf{X}}_{2,k}), ({mathbf{Z}}_{A,k}) and ({mathbf{Z}}_{D,k}) are incidence matrices that assign systematic, pen nested within batch effects, additive genotypic effects, and dominance genotypic effects to the phenotypes, respectively. There were no correlations assumed between pen nested within batch effects, between residual effects or between these effects and other random effects. The (co)variance matrix for the pen nested within batch effects is (varleft[begin{array}{l}{mathbf{p}}_{PB}\ {mathbf{p}}_{CB}end{array}right]=mathbf{I}otimes mathbf{P}=mathbf{I}otimes left[begin{array}{ll}{sigma }_{p,PB}^{2}& 0\ 0& {sigma }_{p,CB}^{2}end{array}right]) and the (co)variance matrix for the residuals is (varleft[begin{array}{l}{mathbf{e}}_{PB}\ {mathbf{e}}_{CB}end{array}right]=mathbf{I}otimes mathbf{R}=mathbf{I}otimes left[begin{array}{ll}{sigma }_{e,PB}^{2}& 0\ 0& {sigma }_{e,CB}^{2}end{array}right]), where ({sigma }_{p,k}^{2} and {sigma }_{e,k}^{2}) are the pen nested within batch and residual variances in the PB and CB populations, respectively. Vector ({mathbf{a}}_{k}) is the vector of additive genetic effects and vector ({mathbf{d}}_{k}) is the vector of dominance genotypic effects. The (co)variance matrix for the genotypic additive effects is:

(varleft[begin{array}{l}{mathbf{a}}_{PB}\ {mathbf{a}}_{CB}end{array}right])= ({mathbf{G}}_{0}otimes mathbf{G}=left[begin{array}{ll}{sigma }_{A*PB}^{2}& {sigma }_{A*PB,CB}\ {sigma }_{A*PB,CB}& {sigma }_{A*CB}^{2}end{array}right]otimes frac{mathbf{Z}{mathbf{Z}}^{mathbf{^{prime}}}}{left{tr([mathbf{Z}{mathbf{Z}}^{mathbf{^{prime}}}])/nright}})where ({sigma }_{A*k}^{2}) is the additive genotypic variance in either PB or CB, and ({sigma }_{A*PB,CB}) is the additive genotypic covariance between PB and CB populations. Similarly, the (co)variance matrix for the genotypic dominance effects is:

where ({sigma }_{D*k}^{2}) and ({sigma }_{D*PB,CB}) are the dominance genotypic variances in PB and CB, and the additive genotypic covariance between PB and CB populations. Terms (mathbf{Z}) and (mathbf{W}) are incidence matrices relating additive and dominance genotypic effects to the PB and CB animals coded as described in the models above. In this bivariate genotypic model, the (mathbf{Z}) and (mathbf{W}) matrices include all PB and CB animals and the model estimates additive and dominance genotypic covariances which cannot be interpreted in the same way as the statistical covariance of breeding values and dominance deviations¹². The prior distributions for the parameters of the model were (P({varvec{beta}},b)sim k), (Pleft({varvec{p}}left|mathbf{P}right.right)sim N(0, mathbf{I}otimes mathbf{P})), (P(mathbf{a}left|{mathbf{G}}_{0}right.)sim {varvec{M}}{varvec{V}}N(0, {mathbf{G}}_{0}otimes mathbf{G})), (P(mathbf{d}left|{mathbf{D}}_{0}right.)sim {varvec{M}}{varvec{V}}N(0, {mathbf{D}}_{0}otimes mathbf{D})) and (Pleft({varvec{e}}left|mathbf{R}right.right)sim {varvec{M}}{varvec{V}}N(0, mathbf{I}otimes mathbf{R})).

Additive and dominance genotypic correlations between BFT, ADG and PHS in PB and CB populations

A tri-trait genotypic model including additive and dominance effects and a genomic inbreeding coefficient was used to estimate additive and dominance genotypic correlations between BFT, ADG and PHS within each PB and CB population. This approach allowed estimating (co)variances of additive and dominance genotypic effects between these three traits within purebred and crossbred populations. To achieve that, the following tri-trait model was applied separately in each population:

Index (j) is used to denote either ADG, BFT or PHS trait (i.e.,(k=ADG, BFT,PHS)). Term ({mathbf{y}}_{j}) is the phenotypic value of individuals, ({{varvec{upbeta}}}_{{varvec{j}}}) is a vector of systematic effects, ({mathbf{p}}_{ADG}) and ({mathbf{p}}_{BFT}) are the vector of pen nested within batch effects for ADG and BFT. Terms ({mathbf{a}}_{j}) and ({mathbf{d}}_{j}) are the vectors of additive and dominance genotypic effects, respectively. Term ({mathbf{e}}_{j}) is the vector of residual effects. Matrices ({mathbf{X}}_{1,j}), ({mathbf{X}}_{2,ADG}), ({mathbf{X}}_{2,BFT}), ({mathbf{Z}}_{A,j}) and ({mathbf{Z}}_{D,j}) are incidence matrices that assign the corresponding systematic effects and pen nested within batch effects, additive and dominance genotypic random effects to the phenotypes. There was no correlation between pen nested within batch effects, between residual effects or between these effects and other random effects. The (co)variance matrix for the pen nested within batch effects is (varleft[begin{array}{l}{mathbf{p}}_{ADG}\ {mathbf{p}}_{BFT}\ {mathbf{p}}_{PHS}end{array}right]=mathbf{I}otimes mathbf{P}=mathbf{I}otimes left[ begin{array}{lll}{sigma }_{p,ADG}^{2}& 0& 0\ 0& {sigma }_{p,BFT}^{2}& 0\ 0& 0& 0end{array}right],) where ({sigma }_{p,ADG}^{2}) and ({sigma }_{p,BFT}^{2}) are the pen nested within batch random variances for ADG and BFT, respectively. The (co)variance matrix for the residuals is (varleft[begin{array}{l}{mathbf{e}}_{ADG}\ {mathbf{e}}_{BFT}\ {mathbf{e}}_{PHS}end{array}right]=mathbf{I}otimes mathbf{R}=mathbf{I}otimes left[begin{array}{lll}{sigma }_{e,ADG}^{2}& 0& 0\ 0& {sigma }_{e,BFT}^{2}& 0\ 0& 0& {sigma }_{e,PHS}^{2}end{array}right],) where ({sigma }_{e,j}^{2}) is the residual variance for each trait. Pen nested within batch effect and the residuals were assumed to be uncorrelated due to overparameterization of the models.

The (co)variance matrix for the additive genotypic effects is

whereas the (co)variance matrix for the dominance genotypic effects is

Terms ({sigma }_{A*j}^{2}) and ({sigma }_{D*j}^{2}) are the additive and dominance genotypic variances for each trait, respectively, and ({sigma }_{A*i,j}) and ({sigma }_{D*i,j}) are the additive and dominance genotypic covariances between the traits.

Terms (mathbf{Z}) and (mathbf{W}) are incidence matrices relating additive and dominance genotypic effects of either the PB or the CB animals coded as described in the models above, respectively. The prior distributions for the parameters of the model were (P({varvec{beta}},b)sim k), (Pleft({varvec{p}}left|mathbf{P}right.right)sim {varvec{M}}{varvec{V}}N(0, mathbf{I}otimes mathbf{P})), (P(mathbf{a}left|{mathbf{G}}_{0}right.)sim {varvec{M}}{varvec{V}}N(0, {mathbf{G}}_{0}otimes mathbf{G})), (P(mathbf{d}left|{mathbf{D}}_{0}right.)sim {varvec{M}}{varvec{V}}N(0, {mathbf{D}}_{0}otimes mathbf{D})) and (Pleft({varvec{e}}left|mathbf{R}right.right)sim {varvec{M}}{varvec{V}}N(0, mathbf{I}otimes mathbf{R})).

Parameter inference

A Bayesian framework was adopted for inference. Flat prior distributions were assumed for the elements of the matrices of the variance components of all the models. Marginal posterior distributions of the parameters of interest were estimated via Gibbs’s sampling algorithm. Mean, highest density interval at 95% (HPD95%) and probability of the parameter of being higher or lesser that certain value, were obtained from these marginal posterior distributions. Univariate models were initially implemented using BGLR software (https://github.com/gdlc/BGLR-R)²¹, and the multiple-trait models were later run with GIBBS2f90 software (http://nce.ads.uga.edu/wiki/doku.php?id=start)²². For each proposed model, single chains of 200,000 iterations were run by discarding the first 20,000 iterations. The burn-in was determined by visual inspection of the chains and by the procedures of Raftery and Lewis²³ and Geweke²⁴. Samples of the parameters of interest were saved every 10 rounds.