### Ethics statement

Experimental procedures were approved by the All University Committee on Animal Use and Care at Michigan State University (AUF# 09/03-114-00).

### Animals

Animals from a three-generation resource pig population developed at Michigan State University were used for this study. This population is an F_{2} cross originated from 4 F_{0} Duroc sires and 15 F_{0} Pietrain dams. The full pedigree consists in a single large family of 19 F_{0}, 56 F_{1} (including 50 females and 6 males), and 954 F_{2} animals. Further details of population development and animal management are found in Edwards et al. [16, 17].

### Phenotypic data

Over 60 different phenotypes related to growth, body composition, carcass merit and meat quality were collected on the Michigan State University F_{2} Duroc x Pietrain resource population. In this study, we focused on carcass and meat quality phenotypes that were measured on or were directly related to longissimus dorsi (loin) muscle. Details of carcass and meat quality phenotype collection were published in Edwards et al. [17]. Briefly, carcass traits collected included loin muscle weight, and loin muscle pH and temperature at 45 min and 24 h postmortem. During carcass fabrication, measurements of loin muscle area and off-midline 10th-rib backfat thickness were also recorded. In addition, a section of the loin was further evaluated for meat quality traits. Traits included subjective and objective color, marbling and firmness. Samples were also evaluated for proximate composition, including moisture, intramuscular fat and protein. A trained sensory panel evaluated samples for juiciness, tenderness, connective tissue and off-flavor.

### Genotypic data

Animals from the Duroc x Pietrain resource population, including F_{0}, F_{1}, and F_{2} individuals, were genotyped for 124 dinucleotide microsatellites genetic markers (3-9 markers per chromosome) at a commercial laboratory (GeneSeek Inc., Lincoln, NE). This genotype information was used to derive breed of origin probabilities across the genome of F_{2} animals. In particular, probabilities of each F_{2} individual being homozygous for Duroc alleles (*P*
_{11}), homozygous for Pietrain alleles (*P*
_{22}), or heterozygous (*P*
_{12} *or P*
_{21}) were estimated at each microsatellite marker and at 11 equidistant inter-marker positions, yielding in total 1,279 putative QTL positions spanning the whole pig genome. Breed of origin probabilities were derived assuming that the parental breeds (i.e., Duroc and Pietrain) were fixed for alternative QTL alleles [18].

### Transcriptomic data

Longissimus dorsi (loin) muscle tissue was sampled from a total of 176 F_{2} individuals during slaughter. The transcriptome of this tissue was measured for each of the 176 F_{2} animals using a pig whole-genome 70-mer oligonucleotide microarray. This microarray includes 20,400 annotated oligonucleotides spanning the whole swine genome. Details regarding tissue sample collection, sample preparation, microarray hybridization and pre-processing data were reported in Steibel et al. [19]. The resulting normalized gene expression data (intensity values) were expressed in the log2 scale.

### Genome-wide linkage analysis

The dataset for analysis included several phenotypes, genotype information, and gene expression data for a total of 171 F_{2} individuals. Two complementary whole-genome scans were performed: first, we carried out a classical phenotypic QTL mapping (pQTL) integrating phenotypic and genotypic data, and second we performed an expression QTL mapping (eQTL) integrating transcriptional profiling with genotypic data.

For the pQTL mapping, the following linear model was fitted separately to each phenotype directly related to loin muscle (e.g., loin muscle weight, loin muscle area):

$$ {y}_{ijk}=\mu +se{x}_i+ grou{p}_j+ carcw{t}_k\cdot \beta +{c}_k\cdot \alpha +{e}_{ijk} $$

where *y*
_{
ijk
} is the phenotypic trait under study of the *k*
^{th} F_{2} animal within the combination of *sex*
_{
i
} and⋅ *group*
_{
j
}, ⋅ *μ* is the general mean, *sex*
_{
i
} represents the fixed effect of the sex of the *k*
^{th} animal, *group*
_{
j
} represents the fixed effect of the slaughter group of the *k*
^{th} animal, and *carcwt*
_{
k
} is the carcass weight of the *k*
^{th} animal as a linear covariate. As mentioned before, the additive QTL coefficient \( c \) was derived assuming that the parental breeds were fixed for alternative alleles. In particular, *c*
_{
k
} = *P*
_{11} − *P*
_{22} is the conditional expectation of the number of Duroc alleles carried by the *k*
^{th} animal. The significance of the additive pQTL effect *α* at each of the 1,279 putative pQTL positions for each phenotypic trait was tested using an *F*-test by comparing the full model to the reduced model without the QTL effect. Significance thresholds of 5 % at genome-wise level were determined through the use of permutation tests [20].

For the eQTL mapping, the following linear mixed model was fit to normalized log-intensity data:

$$ {w}_{ijkl}=\mu + dy{e}_i+ arra{y}_j+\operatorname{se}{x}_k+{c}_l\cdot \alpha +{e}_{ijkl} $$

where *w*
_{
ijkl
} is the normalized log-intensity for each oligonucleotide measured in the loin muscle of the *l*
^{th} animal, *μ* is the general mean, *dye*
_{
i
}, *array*
_{
j
}, and *sex*
_{
k
} are effects accounting for systematic variation in the microarray experiment of the *l*
^{th}animal; *dye* and *sex* were fitted as fixed effects, while *array* was fitted as a random effect. As described above, *c*
_{
l
} is the additive QTL coefficient of the *l*
^{th} animal calculated as *P*
_{11} − *P*
_{22}. The significance of the additive eQTL effect *α* at each of the 1,279 putative eQTL positions and for each expression trait was tested using a likelihood ratio test by comparing the aforementioned model to a reduced model without QTL effect. The *p*-values were corrected for multiple testing across all expression traits and positions using Benjamini and Hochberg procedure [21].

### Causal structural learning

Causal structures are represented here using graphical models; these models combine the rigor of a probabilistic approach with the intuitive representation of relationships given by graphs. Graphical models are composed of two parts: a set *V* of random variables describing the quantities of interest, and a graph *G* = (**V**, *E*) in which each vertex *ν* ∈ *V* is called node, and each edge *e* ∈ *E*, also called arc or link, is used to express the dependence structure of the data, i.e., the set of dependence relationships among the variables in *V* [22].

There are several structure learning algorithms that can be used to infer the network structure underlying a given set of correlated variables, assuming that conditional independencies in the joint probability distribution of these variables mirror *d*-separations in the causal structure (for more details, see [6, 23]). One of such algorithms is the Inductive Causation (IC) algorithm, which is able to search for a class of minimal causal structures that are compatible with the conditional independencies implied by the joint distribution of the data [24]. The IC algorithm, when applied to a set *V* of variables, can be described as follows:

Step 1. For each pair of variables *A* and *B* in *V*, search for set of variables *S*
_{
AB
} ⊂ *V* such that *A* and *B* are independent given **S**
_{
AB
}. If there is no such set, i.e., if *A* and *B* are dependent for every possible **S**
_{
AB
}, then place an undirected edge between *A* and *B*.

Step 2. For each pair of non-adjacent variables *A* and *B* with a common adjacent variable *C*, search for a possible set **S**
_{AB} containing *C* such that *A* and *B* are independent given **S**
_{AB}. If there is no such set, then assign the direction of the edges A ‐ C and C ‐ B as *A* → *C* and *C* ← *B*.

Step 3. In the partially directed graph returned by the previous two steps, orient as many of the undirected edges as possible in such a way that it does not result in (i) new v-structures (i.e. new unshielded colliders) or (ii) directed cycles.

Even though the IC algorithm provides the theoretical framework for causal structural learning using conditional independent tests, its application to practical problems with several variables is hampered due to the exponential number of possible conditional independence relationships to be tested. This has led to the development of more efficient algorithms. Here, we have used one of such algorithms, the Incremental Association Markov Blanket (IAMB) algorithm [25]. The IAMB algorithm first learns the Markov Blanket of each variable in the dataset; the Markov Blanket of a given variable *Y* is defined as the minimal set of variables conditioned on which all other variables are probabilistically independent of the target *Y*. This preliminary step reduces the number and the size of the subsets considered in the conditional tests, and hence results in a lower computational complexity without compromising the accuracy of the resulting causal network [25].

Practical application of the IAMB algorithm involves performing a set of statistical decisions using conditional independence tests. In the context of normally distributed variables, these tests are functions of the partial correlation coefficients *ρ*
_{
XY|W
} between *X* and *Y* given **W**. Here, we used the Fisher’s Z test, which involves a transformation of the linear correlation coefficient and is defined as:

$$ Z\left(X,Y\left|\mathbf{W}\right.\right)=\frac{1}{2}\cdot \sqrt{n-\left|\mathbf{W}\right|-3\cdot } \log \frac{1+{\widehat{\rho}}_{XY\left|\mathbf{W}\right.}}{1-{\widehat{\rho}}_{XY\left|\mathbf{W}\right.}} $$

which has an approximate normal distribution with mean zero and variance 1, i.e., *Z*(*X*, *Y*|**W**) ∼ *N*(0, 1). After the structure of the network was learned, the estimation of the parameters of the local distributions was performed using maximum likelihood. Since the variables under study are continuous, the causal parameters take the form of regression coefficients. Furthermore, the stability of the structure of the causal networks was evaluated using Jackknife resampling. By leaving out one observation per time from the dataset, we could evaluate the stability of each edge in the original network in terms of presence (binary variable; presence or absence in the resampled network) and direction (three possible outcomes; same direction as the original arrow, opposite direction, or undirected arc). All these analyses were performed using the bnlearn package [26] implemented in the R language/environment [27].