## ABSTRACT

Structural equation models allow causal effects between 2 or more variables to be considered and can postulate unidirectional (recursive models; RM) or bidirectional (simultaneous models) causality between variables. This review evaluated the properties of RM in animal breeding and how to interpret the genetic parameters and the corresponding estimated breeding values. In many cases, RM and mixed multitrait models (MTM) are statistically equivalent, although subject to the assumption of variance-covariance matrices and restrictions imposed for achieving model identification. Inference under RM requires imposing some restrictions on the (co)variance matrix or on the location parameters. The estimates of the variance components and the breeding values can be transformed from RM to MTM, although the biological interpretation differs. In the MTM, the breeding values predict the full influence of the additive genetic effects on the traits and should be used for breeding purposes. In contrast, the RM breeding values express the additive genetic effect while holding the causal traits constant. The differences between the additive genetic effect in RM and MTM can be used to identify the genomic regions that affect the additive genetic variation of traits directly or causally mediated for another trait or traits. Furthermore, we presented some extensions of the RM that are useful for modeling quantitative traits with alternative assumptions. The equivalence of RM and MTM can be used to infer causal effects on sequentially expressed traits by manipulating the residual (co)variance matrix under the MTM. Further, RM can be implemented to analyze causality between traits that might differ among subgroups or within the parametric space of the independent traits. In addition, RM can be expanded to create models that introduce some degree of regularization in the recursive structure that aims to estimate a large number of recursive parameters. Finally, RM can be used in some cases for operational reasons, although there is no causality between traits.

## Key words

## INTRODUCTION

Mixed multitrait models (

**MTM**;Henderson, 1984

) have been the standard for the genetic evaluation of livestock for more than 40 yr. They consider genetic and environmental covariances between traits in predicting the breeding values. However, the covariance represents the association between variables, only, that can be produced by a causal relationship or by the presence of cofounders (Pearl, 2009

). Establishing causality between 2 variables poses a controversial difficulty in statistical inference. A causal relationship exists between y_{1}and y_{2}if and only if an intervention in the population to change y_{1}will change y_{2}(Pearl, 2009

). In that case, the association between y_{1}and y_{2}is a causal relationship; however, establishing causality from the association is challenging because forces other than causality, such as confounding variables, can produce an association from y_{1}to y_{2}(or vice versa). Even if a causal effect exists, the association is not directional; this is, an association does not identify whether y_{1}causes y_{2}or vice versa.Since the early 20th century, researchers have attempted to tease out causality based on statistical methods ranging from pathway analyses () to more complex structural equation models developed in many different fields (

Shipley, 2000

; Hershberger, 2003

; Van Dijk, 2003

). Structural equation models can postulate unidirectional causality (recursive models; **RM**) from y_{1}into y_{2}(or vice versa), or mutual causality between y_{1}and y_{2}(simultaneous models). In the fields of quantitative genetics and animal breeding,Gianola and Sorensen, 2004

were the first to suggest the use of structural equation models in the realm of mixed models and proposed a Bayesian implementation to sample values of the posterior distribution of the structural parameters with a Markov Chain Monte Carlo approach. Since then, structural equation models have been widely used in animal breeding to infer causal relationships between phenotypes and to estimate genetic parameters. Most applications have assumed unidirectional relationships and used RM; however, their statistical and biological interpretations remain debated. In most cases, RM inference can be achieved within an MTM framework, although it can present identification problems in certain circumstances, and some restrictions are needed to ensure proper parameter estimation. The objective of this review was to clarify some controversial aspects of recursive models that have not been fully addressed.## RECURSIVE MODELS IN ANIMAL BREEDING

No human or animal subjects were used, so this analysis did not require approval by an Institutional Animal Care and Use Committee or Institutional Review Board. Following the notation proposed by

where

Breeding values $\left({\text{u}}^{\prime}=\left\{{{\text{u}}^{\prime}}_{1},{{\text{u}}^{\prime}}_{2},\dots ,{{\text{u}}^{\prime}}_{n},{{\text{u}}^{\prime}}_{n+1},\dots ,{{\text{u}}^{\prime}}_{s}\right\}\right)$ and residual effects $\left({\text{e}}^{\prime}=\left\{{{\text{e}}^{\prime}}_{1},{{\text{e}}^{\prime}}_{2},\dots ,{{\text{e}}^{\prime}}_{n}\right\}\right)$ are distributed following multivariate Gaussian distributions:

where n is the number of multivariate records, s is the number of individuals including recorded (n) and unrecorded (s − n),

Gianola and Sorensen, 2004

, the statistical definition of RM is as follows:
$\mathrm{\Lambda}{\mathbf{y}}_{\mathrm{i}}={\mathbf{X}}_{\mathrm{i}}\mathbf{b}+{\mathbf{u}}_{\mathrm{i}}+{\mathbf{e}}_{\mathrm{i}},$

[1]

where

**b**is the vector of fixed effects,**y**_{i},**u**_{i}, and**e**_{i}are m × 1 vectors of phenotypic measurements, additive genetic effects and residuals of the m traits associated with the ith multivariate record, and**X**_{i}is its corresponding incidence matrix.**Λ**is an m × m matrix of the recursive parameters with ones on the diagonal and minus the recursive effects of the ith trait on the jth trait (−λ_{1→j}) in some or all the elements below the diagonal, as follows:$\mathrm{\Lambda}=\left[\begin{array}{cccc}1& 0& \dots & 0\\ -{\lambda}_{1\to 2}& 1& \dots & 0\\ \dots & \dots & \dots & \dots \\ -{\lambda}_{1\to \mathrm{k}}& -{\lambda}_{2\to \mathrm{k}}& \dots & 1\end{array}\right].$

Breeding values $\left({\text{u}}^{\prime}=\left\{{{\text{u}}^{\prime}}_{1},{{\text{u}}^{\prime}}_{2},\dots ,{{\text{u}}^{\prime}}_{n},{{\text{u}}^{\prime}}_{n+1},\dots ,{{\text{u}}^{\prime}}_{s}\right\}\right)$ and residual effects $\left({\text{e}}^{\prime}=\left\{{{\text{e}}^{\prime}}_{1},{{\text{e}}^{\prime}}_{2},\dots ,{{\text{e}}^{\prime}}_{n}\right\}\right)$ are distributed following multivariate Gaussian distributions:

$\text{u}~\text{N}\left(0,\text{A}\otimes \text{G}\right)\text{and}\text{e}~\text{N}\left(0,\text{I}\otimes \text{R}\right),$

where n is the number of multivariate records, s is the number of individuals including recorded (n) and unrecorded (s − n),

**G**and**R**are m × m genetic and residual (co)variances matrices,**I**is the identity matrix of respective order, and**A**is the numerator relationship matrix (which can be replaced by a genomic relationship matrix in a genomic selection approach).This RM can be reparametrized as a standard mixed model (MTM) by multiplying all terms of [1] by

with

where

**Λ**^{−1}, as follows:${\mathbf{y}}_{\mathrm{i}}={\mathrm{\Lambda}}^{-1}{\mathbf{X}}_{\mathrm{i}}\mathbf{b}+{\mathrm{\Lambda}}^{-1}{\mathbf{u}}_{\mathbf{i}}+{\mathrm{\Lambda}}^{-1}{\mathbf{e}}_{\mathbf{i}}={\mathrm{\Lambda}}^{-1}{\mathbf{X}}_{\mathrm{i}}\mathbf{b}+{\mathbf{u}}_{\mathbf{i}}^{\ast}+{\mathbf{e}}_{\mathbf{i}}^{\ast}$

with

${\mathbf{u}}^{\ast}\sim \mathrm{N}\left(0,\mathbf{A}\otimes {\mathbf{G}}^{\ast}\right)\phantom{\rule{thickmathspace}{0ex}}\mathrm{a}\mathrm{n}\mathrm{d}\phantom{\rule{thickmathspace}{0ex}}{\mathbf{e}}^{\ast}\sim \mathrm{N}\left(0,\mathbf{I}\otimes {\mathbf{R}}^{\ast}\right),$

where

**u***and**e***are the vector of breeding values and residuals under the MTM, and ${\mathbf{G}}^{\ast}={\mathrm{\Lambda}}^{-1}\mathbf{G}{\mathrm{\Lambda}}^{-{1}^{\prime}}\text{'}$ and ${\mathbf{R}}^{\ast}={\mathrm{\Lambda}}^{-1}\mathbf{R}{\mathrm{\Lambda}}^{-{1}^{\prime}}.$## IDENTIFICATION

The fully conditional density of data given the parameters of the model under an RM is as follows:

and, under the MTM, is as follows:

As ${\mathbf{R}}^{\ast -1}={\mathrm{\Lambda}}^{\prime}{\mathbf{R}}^{-1}\mathrm{\Lambda},$ $\phantom{\rule{thickmathspace}{0ex}}{\mathbf{u}}_{\mathrm{i}}^{\ast}={\mathrm{\Lambda}}^{-1}{\mathbf{u}}_{\mathbf{i}},$ and, if each factor affects all traits (

where f is the number of factors. It can be shown that fully conditional densities of the 2 models are completely equivalent, and the models can only differ due to the assumed prior distributions within a Bayesian framework. The MTM is completely identifiable, and

$\begin{array}{c}\mathrm{p}\left(\mathbf{y}|\mathbf{b},\mathbf{u},\mathrm{\Lambda},\mathbf{R}\right)\\ \propto \prod _{\mathrm{i}=1}^{\mathrm{n}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{1}{\text{2}}{\left(\mathrm{\Lambda}{\mathbf{y}}_{\mathbf{i}}-{\mathbf{X}}_{\mathbf{i}}\mathbf{b}-{\mathbf{u}}_{\mathbf{i}}\right)}^{\mathrm{\prime}}{\mathbf{R}}^{-1}\left(\mathrm{\Lambda}{\mathbf{y}}_{\mathbf{i}}-{\mathbf{X}}_{\mathbf{i}}\mathbf{b}-{\mathbf{u}}_{\mathbf{i}}\right)\right\}\\ =\prod _{\mathrm{i}=1}^{\mathrm{n}}\mathrm{e}\mathrm{x}\mathrm{p}\{-\frac{1}{\text{2}}{\left({\mathbf{y}}_{\mathbf{i}}-{\mathrm{\Lambda}}^{-1}{\mathbf{X}}_{\mathbf{i}}\mathbf{b}-{\mathrm{\Lambda}}^{-1}{\mathbf{u}}_{\mathbf{i}}\right)}^{\mathrm{\prime}}\left(\mathrm{\Lambda}{\mathbf{R}}^{-1}\mathrm{\Lambda}\right)\\ \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\left({\mathbf{y}}_{\mathbf{i}}-{\mathrm{\Lambda}}^{-1}{\mathbf{X}}_{\mathbf{i}}\mathbf{b}-{\mathrm{\Lambda}}^{-1}{\mathbf{u}}_{\mathbf{i}}\right)\}\end{array}\text{'}$

and, under the MTM, is as follows:

$\begin{array}{c}\mathrm{p}\left(\mathbf{y}|{\mathbf{b}}^{\ast},{\mathbf{u}}^{\ast},{\mathbf{R}}^{\ast}\right)\\ \propto \prod _{\mathrm{i}=1}^{\mathrm{n}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{1}{\text{2}}{\left({\mathbf{y}}_{\mathbf{i}}-{\mathbf{X}}_{\mathbf{i}}{\mathbf{b}}^{\ast}-{\mathbf{u}}_{\mathrm{i}}^{\ast}\right)}^{\mathrm{\prime}}{\mathbf{R}}^{\ast -1}\left({\mathbf{y}}_{\mathbf{i}}-{\mathbf{X}}_{\mathbf{i}}{\mathbf{b}}^{\ast}-{\mathbf{u}}_{\mathrm{i}}^{\ast}\right)\right\}\end{array}$

As ${\mathbf{R}}^{\ast -1}={\mathrm{\Lambda}}^{\prime}{\mathbf{R}}^{-1}\mathrm{\Lambda},$ $\phantom{\rule{thickmathspace}{0ex}}{\mathbf{u}}_{\mathrm{i}}^{\ast}={\mathrm{\Lambda}}^{-1}{\mathbf{u}}_{\mathbf{i}},$ and, if each factor affects all traits (

Wu et al., 2010

),
${\mathrm{\Lambda}}^{-1}{\mathbf{X}}_{\mathbf{i}}\mathbf{b}={\mathbf{X}}_{\mathbf{i}}{\mathrm{\Lambda}}^{-1}\left(\begin{array}{c}{\mathbf{b}}_{1}\\ \dots \\ {\mathbf{b}}_{\mathbf{f}}\end{array}\right)={\mathbf{X}}_{\mathbf{i}}{\mathbf{b}}^{\ast},$

where f is the number of factors. It can be shown that fully conditional densities of the 2 models are completely equivalent, and the models can only differ due to the assumed prior distributions within a Bayesian framework. The MTM is completely identifiable, and

**G***and**R***can be inferred from likelihood or Bayesian approaches. Instead, the RM has up to m^{2}/2 − m/2 additional recursive parameters, and there are an infinite number of combinations of**R**,**G,**and the recursive parameters that generate the same likelihood. Therefore, to achieve the inference of the parameters involved, some restrictions must be imposed. In animal breeding models, 2 strategies have been used to achieve statistical identification, (1) imposing restrictions on the elements of the (co)variance component matrices (Varona et al., 2007

) or (2) imposing constraints in the form of linear combinations of explanatory variables to dispose of instrumental variables (Gianola and Sorensen, 2004

). To illustrate these 2 types of constraints, let us assume a simple recursive model that has 2 traits (y_{1}and y_{2}) in 4 scenarios (Figure 1).For simplicity, assume that the means of the 2 traits are known. In the first scenario (Figure 1a), the only relationship between y

_{1}and y_{2}is defined by the recursive parameter (λ) because the residuals are not correlated. Therefore, there are 3 sources of information available, the variances of y_{1}and y_{2}and the covariance between them, which can be used to estimate ${\sigma}_{1}^{2},$ ${\sigma}_{2}^{2},$ and λ, and the model is completely identifiable. In the second scenario (Figure 1b), the residuals are correlated; therefore, the covariance between y_{1}and y_{2}depends on λ and σ_{12}, and we have only the same 3 sources of information to infer 4 parameters $\left({\sigma}_{1}^{2},{\sigma}_{2}^{2},{\sigma}_{12},\mathrm{a}\mathrm{n}\mathrm{d}\lambda \right)$ that make the model unidentifiable. One way to make the model identifiable is to include an auxiliary instrumental variable (z; Figure 1c). The instrumental variables must be (i) correlated with the explanatory variable (y_{1}) and (ii) not correlated with the residuals. In that scenario, 3 additional sources of information are available; namely, the covariance between z and y_{1}, the covariance between z and y_{2,}and the variance of z and 2 new unknowns $\left(\mathrm{w}\mathrm{a}\mathrm{n}\mathrm{d}{\sigma}_{z}^{2}\right).$ However, in the scenario where the variable z affects y_{2}, rather than y_{1}(Figure 1d), the covariance between z and y_{1}is null. Then, it does not provide any information, and the model is not identifiable.In short, the only ways to make recursive models identifiable are (1) imposing constraints on the (co)variance matrix (scenario a), or (2) having an instrumental auxiliary variable that influences the dependent traits only through the independent variable (scenario c). Identification of auxiliary variables requires finding a systematic effect or a covariate (in the jargon of mixed models) that explains a significant proportion of the variation in the independent traits, and influences the dependent traits through the recursive relationship, only. In the animal breeding scope, finding a variable that matches these 2 conditions can be difficult because traits of economic interest are often affected by common environmental factors. Nevertheless, there are some examples. For instance, milk yield is genetically correlated with days open, and have a recursive effect because cows that have high energy demands in lactation can incur a negative energy balance, which impairs reproductive function. In that case, the energy composition of feed can be an auxiliary variable. The amount of energy in the ration (z) influences the amount of milk (y

_{1}) a cow produces; however, in the absence of lactation, reproductive performance (y_{2}) is expected to be normal as long as the minimum feed requirements are met. Thus, the amount of energy in the ration (z) affects milk yield (y_{1}), not reproductive performance (y_{2}), but high levels of y_{1}would impair y_{2}.In more complex models, such as multiple-trait RM models that have several recursive dependencies, identifying auxiliary variables can be extremely cumbersome. However, forcing the residual covariance matrix to be diagonal guarantees the identification in multiple-trait RM, although we must bear in mind that this is a strict assumption in most cases. Phenotypic causal relationships between traits rarely is the sole cause of environmental covariance between traits, there are usually other shared external environmental effects. In addition, mixed models can have other restrictions such as setting to zero the genetic or the phenotypic covariance between traits, as proposed by

Varona et al., 2007

and Jamrozik and Schaeffer, 2011

.The equivalence between MTM and RM holds under the assumption of fully recursive relationships. If the number of constrains imposed is greater than the number required to be identifiable or if some of the recursive relationships are fixed, the model is over-identified. Over-identified models are often nested models with respect to the complete models (

Bentler and Satorra, 2010

), and the parameter space of **G***and**R***is restricted to a subspace. Therefore, they can be compared with frequentist or Bayesian models comparison tools; however, these tests compare model fit to the available data, only, and cannot differentiate causality from association. To illustrate that, one can test whether 2 traits (y_{1}and y_{2}) are associated but cannot determine whether y_{1}has an effect on y_{2}or vice versa, or whether cofounders affect both traits simultaneously.## RECURSIVE PARAMETERS: INTERPRETATION AND IMPLICATIONS

In RM, the recursive parameter λ

To evaluate the effect of the recursive parameter on the heritability and the genetic correlation, we set a simple bi-character scenario as a case of study, where ${\sigma}_{u1}^{2}={\sigma}_{u2}^{2}={\sigma}_{e1}^{2}={\sigma}_{e2}^{2}=1.$ The heritability of trait 2 and the genetic correlation between traits 1 and 2 can be calculated as follows (

For simplicity, λ

_{1→2}is the expected change in trait y_{2}for each unit of change in trait y_{1}, which affects other important parameters such as heritability and genetic correlation. Let us assume a simple recursive model that has 2 traits in which${\mathrm{y}}_{\mathrm{i}1}={{\mathbf{x}}^{\prime}}_{\mathrm{i}1}{\beta}_{1}+{\mathrm{u}}_{\mathrm{i}1}+{\mathrm{e}}_{\mathrm{i}1}$

[2]

${\mathrm{y}}_{\mathrm{i}2}={\lambda}_{1\to 2}{\mathrm{y}}_{\mathrm{i}1}+{{\mathbf{x}}^{\prime}}_{\mathrm{i}2}{\beta}_{2}+{\mathrm{u}}_{\mathrm{i}2}+{\mathrm{e}}_{\mathrm{i}2}$

[3]

To evaluate the effect of the recursive parameter on the heritability and the genetic correlation, we set a simple bi-character scenario as a case of study, where ${\sigma}_{u1}^{2}={\sigma}_{u2}^{2}={\sigma}_{e1}^{2}={\sigma}_{e2}^{2}=1.$ The heritability of trait 2 and the genetic correlation between traits 1 and 2 can be calculated as follows (

Gianola and Sorensen, 2004

):
${\mathrm{h}}_{{\mathrm{u}}_{2}^{\ast}}^{2}=\frac{{\sigma}_{\mathrm{u}2}^{2}+2{\lambda}_{1\to 2}^{}{\sigma}_{\mathrm{u}12}^{}+{\lambda}_{1\to 2}^{2}{\sigma}_{\mathrm{u}1}^{2}}{{\sigma}_{\mathrm{u}2}^{2}+{\sigma}_{\mathrm{e}2}^{2}+2{\lambda}_{1\to 2}^{}\left({\sigma}_{\mathrm{u}12}^{}+{\sigma}_{\mathrm{e}12}^{}\right)+{\lambda}_{1\to 2}^{2}\left({\sigma}_{\mathrm{u}1}^{2}+{\sigma}_{\mathrm{e}1}^{2}\right)},$

$\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\left({\mathrm{u}}_{1},{\mathrm{u}}_{2}^{\ast}\right)=\frac{{\sigma}_{\mathrm{u}12\phantom{\rule{thinmathspace}{0ex}}}+{\lambda}_{1\to 2}{\sigma}_{\mathrm{u}1}^{2}}{\sqrt{{\sigma}_{\mathrm{u}1}^{2}\left({\sigma}_{\mathrm{u}2}^{2}+2{\lambda}_{1\to 2}{\sigma}_{\mathrm{u}12}+{\lambda}_{1\to 2}^{2}{\sigma}_{\mathrm{u}1}^{2}\right)}}.$

For simplicity, λ

_{1→2}is expressed as units of the independent trait residual standard deviation $\left({\sigma}_{\mathrm{e}1}\right).$ Figure 2 shows the landscape of heritability of trait 2 and the recursive effect of trait 1 on trait 2 at various levels of the genetic covariance between the 2 traits. In the presence of recursiveness, the heritability of the dependent trait increases if λ_{1→2}has the same sign or direction as the genetic covariance. The stronger the genetic covariance, the larger the increase in h^{2}. If, however, the recursive parameter and the genetic covariance are in opposite directions, the heritability of the dependent trait decreases. In general, the stronger genetic correlation the larger effect is shown on the change in heritability.Genetic correlation under recursiveness is depicted in Figure 3. The relationship between λ

_{1→2}and the genetic correlation follows a sigmoidal curve. The presence of a recursive effect causes an increase (or decrease) in the genetic correlation if it has the same (or different) direction as the genetic covariance. The recursive parameter has the largest effect on h^{2}and the genetic correlation for values up to 1 genetic standard deviation. Values larger than that tend to have a lesser effect on the 2 genetic parameters because λ_{1→2}affects in a quadratic manner to ${\mathrm{h}}_{{\mathrm{u}}_{2}^{\ast}}^{2};$ therefore, if λ_{1→2}> 1, the phenotypic variability of trait 1 explained by nongenetic factors has a larger effect on the total phenotypic variability of trait 2. The generating R code is presented in Supplemental File S1 (https://figshare.com/articles/journal_contribution/APPENDIX_A_docx/22040597;Varona and González-Recio, 2023

).The example shows that, even though an RM is equivalent to a bi-character trait model, the presence of a causal effect implies that the heritability of the dependent trait and the genetic correlation have a restricted parametric space. Under recursiveness,
${\mathbf{u}}_{2}^{\ast}=({\mathbf{u}}_{2}+{\lambda}_{1\to 2}{\mathbf{u}}_{1}),$ the heritability of the dependent trait (u

_{2}) is affected by ${\lambda}_{1\to 2}^{2}{\sigma}_{\mathrm{u}1}^{2}$ in the numerator and by ${\lambda}_{1\to 2}^{2}\left({\sigma}_{\mathrm{u}1}^{2}+{\sigma}_{\mathrm{e}1}^{2}\right)$ in the denominator as additional terms in the heritability formula. Thus, for ${\lambda}_{1\to 2}^{2}>0,$ ${\mathrm{h}}_{{\mathrm{u}}_{2}^{\ast}}^{2}$ increases proportionally to ${\sigma}_{\mathrm{u}1}^{2}\phantom{\rule{thickmathspace}{0ex}}$ and ${\sigma}_{\mathrm{e}1}^{2}.$ Even if ${\sigma}_{\mathrm{u}2}^{2}=0,$ the heritability of ${\mathbf{u}}_{2}^{\ast}$ is greater than 0 as long as ${\lambda}_{1\to 2}$ is not 0 and ${\mathrm{h}}_{{\mathrm{u}}_{1}}^{2}$ is not null. Hence, even in the hypothetical case in which y_{2}is fully influenced by the environment (no genetic control), an MTM would estimate a nonzero heritability due to the recursive effect of y_{1}on y_{2}and the nonzero heritability of y_{1}.Similarly, the genetic correlation is restricted to a parameter space that depends on the magnitude of the recursive parameter and the genetic covariance in the MTM. The larger the recursive effect, the larger the genetic correlation between the 2 traits, independent of the genetic covariance between them. In this case, the genetic correlation does not come from pleiotropic effects of the genes that affect both traits, nor from linkage disequilibrium between different genes that affect each trait; it rather comes from a recursive nongenetic effect of trait 1 that affects trait 2. The residual correlation follows the same reasoning as the genetic correlation, and it is also restricted to a given parameter space that depends on the recursive parameter and the residual covariance. The causal relationship between both traits produces an inflation or deflation in the residual correlation if the SMM that is not due to the causal phenotype rather than to environmental factors.

## BREEDING VALUES: INTERPRETATION

The breeding values in the MTM for the ith individual (

**u**_{i}*) can be obtained from the breeding values in the RM (**u**_{i}) because ${{\mathbf{u}}_{\mathbf{i}}}^{\ast}={\mathrm{\Lambda}}^{-1}{\mathbf{u}}_{\mathbf{i}}$ and vice versa ${\mathbf{u}}_{\mathbf{i}}=\mathrm{\Lambda}\phantom{\rule{thickmathspace}{0ex}}{{\mathbf{u}}_{\mathbf{i}}}^{\ast}.$ Although MTM and RM may be statistically equivalent, their interpretations are very different. Mixed multitrait models assume that genetic and residual effects of traits are correlated and caused by pleiotropy or linkage disequilibrium between QTLs (genetic correlation) and common environmental effects (residual correlation). In contrast, the RM postulate a causal and unidirectional influence of some traits on others, which affects genetic and residual covariances between traits, but which are not completely determined by common genetic or environmental effects. Thus, the interpretation of the breeding values of traits that are affected by the phenotypic influence of other traits (dependent traits) differs between RM and MTM. In RM, the breeding values of the dependent traits must be understood as the effect of the genes that act directly on the traits, and not indirectly through the phenotypic influence of another trait. Therefore, breeding values in RM must be interpreted as the breeding values corrected for their causal influence or while holding physically the value of the independent traits (Valente et al., 2013

). In contrast, in MTM, breeding values reflect the entire additive effects of the genes on traits, even though they directly or indirectly affect the final phenotype. For example, dominant or more aggressive cows may produce more milk than compliant cows. Dominant cows may have access to more amount of fed, crowding out compliant cows who might eat less and produce less milk. Therefore, aggressiveness can have a positive recursive effect on milk production, even if there is no pleiotropy or LD between genes affecting both traits (null genetic correlation). In MTM, the breeding values for milk yield include the additive genetic effects related to aggressiveness and milk production, whereas in RM, they only include the additive genetic effects for milk yield without the influence of aggressiveness.For practical purposes, in a breeding scheme that has the objective of increasing (or decreasing) the phenotypic mean of a dependent trait, the breeding values generated by an MTM should be used. However, the statistical equivalence of MTM and RM can be invoked to calculate the breeding values from an RM
$({{\mathbf{u}}_{\mathrm{i}}}^{\ast}={\mathrm{\Lambda}}^{-1}{\mathbf{u}}_{\mathrm{i}}).$ Furthermore, an RM allows the split of the breeding values into direct genetic effects (

**u**_{i}) that are attributed to genes that directly affect the trait, and indirect genetic effects (**u**_{i}* −**u**_{i}), which are caused by genes that influence the independent trait and affect the dependent trait through its phenotypic influence. If selection is based on the predicted breeding values from the RM (**u**_{i}), the response to selection in the dependent traits will occur by acting on the direct genetic effects. Note, however, that it may cause a correlated genetic response in the independent trait because of the presence of additive genetic covariance. If the breeding plan aims to modify the dependent trait and keep the independent trait constant, restricted selection indexes (Kempthorne and Nordskog, 1959

) should be used.Similarly, estimates of the additive genetic (co)variances components in an RM (

**G**) and their equivalents in an MTM $({\mathbf{G}}^{\ast}={\mathrm{\Lambda}}^{-1}\mathbf{G}{\mathrm{\Lambda}}^{-{1}^{\prime}})$ allow the decomposition of the genetic variance between traits. The estimates of additive genetic variance from an RM reflect the additive genetic variance caused by genes acting directly on the trait. The difference between the estimates from an MTM and an RM is the additive genetic variance that the dependent traits incorporated through the phenotypic influence of other traits. For instance, assuming the same simple recursive models defined in equations [2] and [3], the additive variance in an MTM can be split into the additive variance in an RM $\left({\sigma}_{\mathrm{u}2}^{2}\right)$ and the additive variance generated by the causal influence of trait 2 $\left(2{\lambda}_{1\to 2}{\sigma}_{\mathrm{u}12}+{\lambda}_{1\to 2}^{2}{\sigma}_{\mathrm{u}1}^{2}\right).$ The same decomposition can be achieved for the additive genetic covariance between traits, which can be divided into the covariance generated by pleiotropic genes that affect the 2 traits $\left({\sigma}_{\mathrm{u}12}\right),$ simultaneously, and the additive genetic covariance caused by the phenotypic influence of the independent trait $\left({\lambda}_{1\to 2}{\sigma}_{\mathrm{u}1}^{2}\right).$Selection on independent traits produces a correlated response in the dependent trait because of both sources of genetic covariance. However, an environmental change or an external intervention (

Rosa and Valente, 2013

), that directly increases (or reduces) the independent traits, will modify the phenotypic performance of dependent traits. This phenotypic change is solely through the phenotypic influence between traits expressed by the recursive parameter, and will not imply any change in the genes associated with the genetic variation. For example, milk yield in dairy cattle is affected by the amount of energy in the diet and has a causal effect on the SCC in milk. An external intervention that modifies the diet of cows by increasing the energy in the ration will increase milk production. It will also indirectly modify the SCC, although there would be no direct influence of the energy of the diet on the SCC.External interventions can also modify the causal structure between traits. In the above example of the dominant and compliant cows, if they are raised under feed restriction or without ensuring full access to feeder, a breeding program selecting for larger production yields would also increase aggressiveness. Genetic correlation from an MTM may show positive estimates, because it does not differentiate between the causal phenotypic effect and the genetic correlation (association due to genes). However, if more feed bins and plenty of feed availability is supplied to cows, more compliant individuals can eat as much as the more aggressive individuals. Hence, the external intervention removes the phenotypic causal effect that favor dominant animals in detriment of other cows.

The ability to predict the expected outcome for those external interventions requires distinguishing between causality and association, which cannot be achieved using standard goodness-of-fit tests given the statistical equivalence of the RM and the MTM. The distinction between them can be extremely complex because of the possibility of confounding effects that affect several traits (

Rubin, 1974

). Some authors (Shipley, 2000

; Pearl, 2009

) have developed methods that can achieve that goal. An adaptation of the Inductive Causation (IC) algorithm (Pearl, 2009

) to the scope of mixed models by Valente et al., 2010

is the most widely used algorithm for inferring the causal structure within a set of phenotypic traits. Briefly, the IC algorithm has the following 3 steps:
- 1.For each pair of traits (y
_{1}and y_{2}), search for a set of traits such that y_{1}is independent of y_{2}, given this set. If y_{1}and y_{2}are dependent in each possible set, they are declared as adjacent and are connected by an undirected edge. - 2.If 2 nonadjacent variables (y
_{1}and y_{2}) are connected to an extra-adjacent variable (y_{3}) and are not independent for all sets of variables that contain y_{3}, connect y_{3}to y_{1}and y_{2}by a directed edge. The result is a partially oriented graph that has directed and undirected edges. - 3.Change as many undirected to directed edges as possible to avoid generating new colliders or cycles.

Valente et al., 2010

, within a Bayesian framework, those partial correlations can be obtained from the posterior distribution of the residual (co)variances matrix calculated in a standard multivariate mixed model. If the HPD (highest posterior density) of a partial correlation includes zero, the 2 traits are assumed to be conditionally independent. A more detailed description of the algorithm is described by Valente et al., 2010

, Valente et al., 2011

) and Rosa and Valente, 2013

. For each pair of traits, the algorithm chooses among independence, causality, or associated confounders. However, it cannot be ruled out that one trait has an effect on another, and, at the same time, there are confounding factors that affect both traits simultaneously. Therefore, neither the RM nor the MTM may reflect the real relationship between the 2 traits, and both models are subject to very strong assumptions. Again, the only way to measure the magnitude of the effects of causality effects or the influence of cofounders is to have an external intervention that modifies the independent trait and allows quantifying the effects on the dependent trait.## RECURSIVE MODELS IN ANIMAL BREEDING

After the introduction of recursive models in quantitative genetics and animal breeding (

Gianola and Sorensen, 2004

), they have been applied in numerous studies. For instance, recursive models have been applied to calving traits (López de Maturana et al., 2010

; Inoue et al., 2017

), milk traits (Rehbein et al., 2013

; Bouwman et al., 2014

; Santos et al., 2015

; Tiezzi et al., 2015

), growth traits (González-Rodríguez et al., 2014

; Leal-Gutiérrez et al., 2018

; Krugmann et al., 2020

), diseases (König et al., 2008

; Rehbein et al., 2013

; Quigley et al., 2018

), and microbiota traits (Saborío-Montero et al., 2020

) in different species. The methods behind recursive models and the wide range of applications have been reviewed elsewhere (- Saborío-Montero A.
- Gutiérrez-Rivas M.
- García-Rodríguez A.
- Atxaerandio R.
- Goiri I.
- López de Maturana E.
- Jiménez-Montero J.A.
- Alenda R.
- González-Recio O.

Structural equation models to disentangle the biological relationship between microbiota and complex traits: Methane production in dairy cattle as a case of study.

Wu et al., 2010

; Valente et al., 2013

; Cantet et al., 2017

; Inoue, 2020

).Recursive models were initially formulated based on the assumption that all traits have a continuous distribution; however, in animal breeding, it is quite common for the traits of interest to be categorical and have 2 or more possible outputs. For those types of traits, either threshold models (

Gianola, 1982

) or models that assume a non-Gaussian distribution for the residuals (Tempelman and Gianola, 1993

; Varona and Sorensen, 2010

) are used. Bivariate recursive models when the independent traits follow a Gaussian distribution has been described for Gaussian – Threshold Models (López de Maturana et al., 2007

; König et al., 2008

; Wu et al., 2008

), Gaussian – Binomial Models (Varona and Sorensen, 2014

), and Gaussian – Multiplicative Binomial models (Varona et al., 2020

). Furthermore, some authors (Heringstad et al., 2009

) have presented models that have complex causal structures and include several Gaussian and Categorical traits. In those studies, the most frequent assumption to ensure identification was to set the residual covariances between traits to zero.Traditionally, RM have proposed causal effects between phenotypes; however, in some cases, the number of traits can be many and clustered in latent pseudophenotypes or latent variables that summarize each group of phenotypes. For example,

Peñagaricano et al., 2015

analyzed 19 traits associated with postweaning growth, carcass quality, fat composition, or meat quality in pigs, which they summarized in 5 latent variables through factor analysis and were used to develop an RM. The use of latent variables might be useful to understand the causal structure of large sets of variables, such as those generated in massive phenotypic devices (Fernandes et al., 2020

) or those that result from the omics analysis (Lippolis et al., 2019

), as proposed by Saborío-Montero et al., 2021

.- Saborío-Montero A.
- López-García A.
- Gutiérrez-Rivas M.
- Atxaerandio R.
- Goiri I.
- García-Rodriguez A.
- Jiménez-Montero J.A.
- González C.
- Tamames J.
- Puente-Sánchez F.
- Varona L.
- Serrano M.
- Ovilo C.
- González-Recio O.

A dimensional reduction approach to modulate the core ruminal microbiome associated with methane emissions via selective breeding.

Recently, the use of RM in GWAS analysis has gained popularity. GWAS has become the main method for identifying genomic regions that harbor the genes that cause the genetic variation in quantitative traits. The influence of those genomic regions on the phenotypic traits can be caused by the direct influence of the genes within the genomic region or mediated by another trait that has a causal dependency on them. For that reason, some (

Momen et al., 2019

; Wang et al., 2020

) have developed RM for use in GWAS that provide several types of signals, direct, indirect, and overall. The direct signals point to the genomic regions that directly affect the dependent trait, and the indirect signals identify the regions that affect the trait through the phenotypic dependency between traits. The overall signals indicate the genomic regions that affect the genetic variability of the traits, either directly or indirectly.## EXTENSIONS OF RECURSIVE MODELS

### Recursive Inference from the Output of an MTM

The inference results on the (co)variance components and the recursive parameters of an RM can be transformed into the results of an MTM. Instead, in the opposite direction, as the number of parameters increases, there are an infinite number of combinations of the (co)variance components and recursive parameters, and the transformation cannot be performed directly. In some scenarios and under some specific restrictions, however, the statistical equivalence of MTM and RM can be used to obtain the recursive parameters and the (co)variance components (

**R**and**G**) in RM from the estimates of the**R*** and**G*** achieved in the MTM.#### Complete Recursiveness

In this scenario,

**R**is assumed to be diagonal, the traits are ordered by fixing the independent variables before the dependent variables, and every preceding trait has a causal influence on all the subsequent traits (i.e.**Λ**has not any zero element below the diagonal). Then,**R***can**LDL'**factorized, where**D**corresponds to the**R**matrix and**L**is**Λ**^{−1}. Subsequently, the**G**matrix can be obtained as $\mathbf{G}=\mathrm{\Lambda}{\mathbf{G}}^{\ast}{\mathrm{\Lambda}}^{\prime}.$ An important advantage of that approach is that it can be used with data sets that have missing information, i.e., when only some traits of a given individual are available and the rest are missing.#### Recursiveness in a Sequential Scenario

A more complex case that can be inferred from the output of an MTM is the scenario in which the traits are recorded sequentially. Suppose we evaluate n+m traits that are divided into 2 groups that are expressed at 2 different times (n traits at time A and m traits at time B). One possible assumption is that the n traits at a first time (A) have a causal effect on the m traits at time B (See Figure 4), and that there are no cofounders (i.e., residual covariances are equal to zero) between the traits of the different groups.

If an MTM is used, it is possible to estimate

**R***and**G***, which are (n+m) x (n+m) matrices. They must be transformed into**R**and**G**matrices and a set of recursive effects (λ) between the n traits recorded at time A and the m traits recorded at time B. The**R***matrix can be divided into 2 symmetric matrices of n x n (**A**) and m x m (**C**) elements and one rectangular matrix of m x n (**B**) and its transpose (**B'**). Subsequently, a block**LDL'**decomposition can be applied: ${\mathbf{R}}^{\ast}=\left[\begin{array}{cc}\mathbf{A}& {\mathbf{B}}^{\prime}\\ \mathbf{B}& \mathbf{C}\end{array}\right]=\left[\begin{array}{cc}\mathbf{I}& 0\\ \mathbf{B}{\mathbf{A}}^{-1}& \mathbf{I}\end{array}\right]\left[\begin{array}{cc}\mathbf{A}& 0\\ 0& \mathbf{C}-\mathbf{B}{\mathbf{A}}^{-1}{\mathbf{B}}^{\prime}\end{array}\right]\left[\begin{array}{cc}\mathbf{I}& {{\mathbf{A}}^{-1}}^{\mathrm{\prime}}{\mathbf{B}}^{\prime}\\ 0& \mathbf{I}\end{array}\right]={\mathrm{\Lambda}}^{-1}\mathbf{R}{\mathrm{\Lambda}}^{-{1}^{\prime}}$and, as before, $\mathbf{G}=\mathrm{\Lambda}{\mathbf{G}}^{\ast}{\mathrm{\Lambda}}^{\prime}.$ If the number of time points is >2, additional block**LDL'**decompositions can be used sequentially. A numerical example is presented in the Appendix.### Parsimonious Recursive Models

Recursive models posit that independent traits can be related to dependent traits by causal effects that are governed by a set of recursive parameters that must be inferred. In some cases, the number of independent traits that can have a causal effect can be very large. For instance, a plethora of new omics data (e.g., transcriptomics, proteomics, microbiomics) have dramatically increased the availability of new biological variables (

Another type of regularization for the recursive parameter inference was proposed by

Lippolis et al., 2019

). All of them can be understood as intermediate traits that affect the phenotypes of economic interest and can be treated as independent traits that have a causal relationship with the target traits. The number of recursive parameters (nr) will be vast, and the amount of statistical information to infer each one would be limited. Therefore, it is possible to define a model that introduces some degree of regularization in the recursive relationships. Christensen et al., 2021

suggested introducing a shrinkage or a prior distribution within the Bayesian framework by assuming a Gaussian prior distribution for the set of recursive relationships
$\lambda =\left\{{\lambda}_{1},{\lambda}_{2},\dots ,{\lambda}_{\mathrm{n}\mathrm{r}}\right\}$ between the intermediate traits and the target (**y**) as follows:_{T}${\mathbf{y}}_{\mathbf{T}}=\sum _{\mathrm{i}=1}^{\mathrm{n}\mathrm{r}}{\lambda}_{\mathrm{i}}{\mathbf{y}}_{\mathbf{i}}+\mathbf{X}{\mathbf{b}}_{\mathbf{T}}+\mathbf{Z}{\mathbf{u}}_{\mathbf{T}}+{\mathbf{e}}_{\mathbf{T}}$

${\mathbf{y}}_{\mathbf{i}}=\mathbf{X}{\mathbf{b}}_{\mathbf{i}}+\mathbf{Z}{\mathbf{u}}_{\mathbf{i}}+{\mathbf{e}}_{\mathbf{i}}$

$\lambda \sim \mathrm{N}\left(0,\mathbf{I}{\sigma}_{\lambda}^{2}\right).$

Another type of regularization for the recursive parameter inference was proposed by

Onogi et al., 2019

. They developed a recursive model that proposes that the weights throughout the growth process can be partially phenotypically determined by the initial weight. Because the influence of the initial weight is expected to depend on its temporal distance from each recorded weight, they defined a recursive function based on B-splines that modeled the evolution of the recursive effect over time. The model was applied to a cow growth data set, which showed that the influence of initial weight was high in the early growth period and decreased thereafter.### Heterogeneous Causal Structures

If there is population stratification, phenotypic causality may differ between subpopulations. In that case, RM can be reformulated to identify an alternative causal relationship between the traits. Strata might be known (e. g., male or female) or unknown (e. g., healthy or sick) categories of individuals. With known stratification, the use of an RM is straightforward with standard REML or Bayesian approaches, although equivalence with MTM no longer holds because the covariance structure in the MTM differs between subpopulations. If the distribution of individuals among subpopulations is unknown (

Jamrozik and Schaeffer, 2010

), an additional latent variable (Wu et al., 2007

) about the unknown strata of the individuals must be included. The inference can be achieved, under a Markov Chain Monte Carlo approach to Bayesian inference, by using a Metropolis-Hasting step (Hastings, 1970

) within a Gibbs sampler (Gelfand and Smith, 1990

) to sample the latent variables.### Nonlinear Relationship Between Traits

The traditional definitions of MTM and RM imply a linear relationship between traits; however, the causal or non-causal relationship between traits can be nonlinear (

Fuerst-Waltl et al., 1997

, Fuerst-Waltl et al., 1998

; Mulder et al., 2015

). A special case of heterogeneous causal structures occurs if the causal influence of the independent trait on the dependent trait varies along its range of values. For example, López de Maturana et al., 2009

, López de Maturana et al., 2010

) developed a series of recursive models that had a heterogeneous structural coefficient between gestation length, calving difficulty, and perinatal mortality in dairy cattle. They defined 4 segments along the parameter space of gestation length (261–267, 268–273, 274–279, 280–291 d), and showed that the recursive relationships differed among the segments, which confirmed that the relationship between traits was nonlinear. Similarly, Jamrozik et al., 2010

proposed a model that defines heterogeneous structural coefficients between milk yield and somatic cell score across (the first 3 lactations) and between (4 intervals) lactations. Likewise, Ibáñez-Escriche et al., 2010

described a change-point recursive model for the relationship between litter size and the number of stillbirths in pigs. The study did not fix the segments within the litter size parameter space *a-priori*, and used change-point techniques (Chib, 1998

) to infer the size of the regions of the parameter space associated with each recursive parameter. The results suggested that the causal effect was higher for large litters than for small litters. Another way to model nonlinear relationships between the dependent and the independent trait is to use a nonlinear function (González-Rodríguez et al., 2014

; Varona and Sorensen, 2014

; Varona et al., 2020

). However, in complex analyses that have more than 2 traits, the model can be affected by identification problems (Gianola and Sorensen, 2004

).A nonlinear causality between traits implies that the genetic parameters of the dependent trait vary along the parameter space of the independent trait. To illustrate this, let us assume a nonlinear recursive model with 2 traits:

where $\mathrm{f}\left({\mathrm{y}}_{\mathrm{i}1}\right)={\lambda}_{1}{\mathrm{y}}_{\mathrm{i}1}+{\lambda}_{2}{\mathrm{y}}_{\mathrm{i}1}^{2},$ the independent trait (y

The total genetic $\left({\mathbf{G}}^{\ast}={{\mathrm{\Lambda}}_{\mathrm{y}1}}^{-1}\mathbf{G}{{\mathrm{\Lambda}}_{\mathrm{y}1}}^{-{1}^{\prime}}\right)$ and residual $\left({\mathbf{R}}^{\ast}={{\mathrm{\Lambda}}_{\mathrm{y}1}}^{-1}\mathbf{R}{\mathrm{\Lambda}}_{\mathrm{y}1}{{1}^{-{1}^{2}}}^{}\right)$ (co)variance matrix can be calculated along the parametric space by using

The plots of the recursive function, the heritability of the dependent traits, and residual and genetic correlations along the parametric space of the independent trait (80 to 120 units) are presented in Figure 5.

${\mathrm{y}}_{\mathrm{i}1}={{\mathbf{x}}^{\prime}}_{\mathrm{i}1}{\beta}_{1}+{\mathrm{u}}_{\mathrm{i}1}+{\mathrm{e}}_{\mathrm{i}1}$

${\mathrm{y}}_{\mathrm{i}2}=\mathrm{f}\left({\mathrm{y}}_{\mathrm{i}1}\right)+{{\mathbf{x}}^{\prime}}_{\mathrm{i}2}{\beta}_{2}+{\mathrm{u}}_{\mathrm{i}2}+{\mathrm{e}}_{\mathrm{i}2},$

where $\mathrm{f}\left({\mathrm{y}}_{\mathrm{i}1}\right)={\lambda}_{1}{\mathrm{y}}_{\mathrm{i}1}+{\lambda}_{2}{\mathrm{y}}_{\mathrm{i}1}^{2},$ the independent trait (y

_{1}) takes values from 80 to 120 units, the recursive parameters are ${\lambda}_{1}=1.00,$ and ${\lambda}_{2}=-0.005,$ and the covariance matrices are as follows:$\mathbf{R}=\left(\begin{array}{cc}{\sigma}_{\mathrm{e}1}^{2}& {\sigma}_{e12}\\ {\sigma}_{e12}& {\sigma}_{\mathrm{e}2}^{2}\end{array}\right)=\left(\begin{array}{cc}100& 0\\ 0& 100\end{array}\right)$

$\mathbf{G}=\left(\begin{array}{cc}{\sigma}_{\mathrm{a}1}^{2}& {\sigma}_{\mathrm{a}12}\\ {\sigma}_{\mathrm{a}12}& {\sigma}_{\mathrm{a}2}^{2}\end{array}\right)=\left(\begin{array}{cc}30& 10\\ 10& 30\end{array}\right).$

The total genetic $\left({\mathbf{G}}^{\ast}={{\mathrm{\Lambda}}_{\mathrm{y}1}}^{-1}\mathbf{G}{{\mathrm{\Lambda}}_{\mathrm{y}1}}^{-{1}^{\prime}}\right)$ and residual $\left({\mathbf{R}}^{\ast}={{\mathrm{\Lambda}}_{\mathrm{y}1}}^{-1}\mathbf{R}{\mathrm{\Lambda}}_{\mathrm{y}1}{{1}^{-{1}^{2}}}^{}\right)$ (co)variance matrix can be calculated along the parametric space by using

$\begin{array}{c}{\mathrm{\Lambda}}_{\mathrm{y}1}=\left[\begin{array}{cc}1& 0\\ -\frac{\mathrm{\partial}\mathrm{f}\left({\mathrm{y}}_{\mathrm{i}1}\right)}{{\mathrm{y}}_{\mathrm{i}1}}& 1\end{array}\right]\phantom{\rule{thickmathspace}{0ex}}{{\mathrm{\Lambda}}_{\mathrm{y}1}}^{-1}=\left[\begin{array}{cc}1& 0\\ \frac{\mathrm{\partial}\mathrm{f}\left({\mathrm{y}}_{\mathrm{i}1}\right)}{{\mathrm{y}}_{\mathrm{i}1}}& 1\end{array}\right]\phantom{\rule{thickmathspace}{0ex}}\\ \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\frac{\mathrm{\partial}\mathrm{f}\left({\mathrm{y}}_{\mathrm{i}1}\right)}{{\mathrm{y}}_{\mathrm{i}1}}={\lambda}_{1}+2{\lambda}_{2}{\mathrm{y}}_{1}.\end{array}$

The plots of the recursive function, the heritability of the dependent traits, and residual and genetic correlations along the parametric space of the independent trait (80 to 120 units) are presented in Figure 5.

It can be seen that a nonlinear recursive function (Figure 5a) leads to a modification of the heritability along the parametric space of y

_{1}. In this example (Figure 5b), the heritability decreases from >0.25 at y_{1}= 80 to <0.21 at y_{1}= 120, and the genetic and the residual correlations decrease from 0.50 to 0.15 (Figure 5c) and from 0.20 to −0.20 (Figure 5d) when y_{1}changes from 80 to 120 units. It implies that the genetic and residual correlations between the 2 traits are also modified by an external intervention that modifies the phenotypic mean of y_{1}. Let us imagine that, in our population, the average independent trait is 80 and, at that point in the parametric space, it is strongly correlated with the dependent trait. Therefore, it can be used to achieve a correlated response under a breeding scheme. However, if there are some external interventions or selection processes that increase the independent trait to 120, the genetic correlation would be reduced to 0.15, and the correlated response in the dependent trait will be much lower.Note that these statements are conditioned on a model that has very strict assumptions, e.g., the nonlinear relationship is caused by the phenotypic effect of the independent trait on the dependent trait, only. Nevertheless, nonlinear recursive models might provide a new and interesting area of research about the relationships among the traits of livestock populations.

### Recursive Models as Operational Tools

Recursive (or structural equation) models have been developed to describe the causal relationships between variables (or traits); however, RM can be useful for operational reasons, even if there is no causality between traits, because they may perform better than MTM in terms of computational requirements, convergence properties, or stability of estimates (

Jamrozik and Schaeffer, 2011

). Moreover, they can be used to propose alternative models for complex problems. For instance, modeling non-genetic (or residual) sources of covariance in bivariate models that includes traits with Gaussian and non-Gaussian (Binomial, Poisson) residuals is challenging; however, some Gaussian and non-Gaussian traits have a strong phenotypic covariance caused by nongenetic factors. In fact, the assumption of statistical independence of the residuals in cases where there is nongenetic covariance can lead to undesirable bias in the predictions of breeding values. Varona and Sorensen, 2014

suggested using an RM and setting the non-Gaussian trait as dependent and the Gaussian trait as independent. Statistical equivalence of the MTM and the RM can be invoked to reconstruct the breeding values of the non-Gaussian trait. With that approach, the covariance between the Gaussian and non-Gaussian residuals is taken into account in the prediction of breeding values and bias is avoided.Another operational application of RM has been developed for the analysis of Residual feed intake (

Koch et al., 1963

), which is one of the most commonly used approaches to assess feed efficiency in livestock populations. The usual method for analyzing genetic variation in residual feed intake has a 2-stage approach (Berry and Crowley, 2013

). First, a linear regression model relating dry matter intake (DMI) to the energy sinks is used. Second, the estimated residuals (residual feed intake) are analyzed using an MTM. Recently, Jamrozik et al., 2021

and Wu et al., 2021

suggested an RM that jointly analyzed DMI and all energy sinks by postulating a recursive relationship between each energy variable and the observed DMI. As such, the systematic and breeding values for the residual feed intake are obtained directly, without the need for the 2-stage approach.A third scenario that can benefit from operational RM is that in which traits are defined as the ratio of 2 Gaussian traits (y

_{2}/y_{1}); e.g., the percentages of fat and protein in milk, the feed conversion ratio. Ratio traits have a Cauchy distribution only if the 2 traits are Gaussian and uncorrelated; otherwise, the distribution cannot be assigned to any standard distribution. Therefore, usually, breeding values for ratio traits are obtained from ad hoc models. To solve this problem,Jamrozik et al., 2017

proposed using an RM to model those kinds of traits as y_{2}= λ × y_{1}. Note that, in the RM, breeding values for the dependent trait (y_{2}) can be understood as the additive effect of the genes after holding the independent trait (y_{1};Valente et al., 2013

). Therefore, the breeding values provided by the model for y_{2}are proxies for the breeding values of the ratio (y_{2}/y_{1}). The main advantage of that approach is that it avoids the use of statistical approximations to the ratio trait, and it can be easily implemented in most mixed model software.## FINAL REMARKS AND CONCLUSIONS

In recent decades, RM have gained popularity in animal breeding. Recursiveness between 2 traits under quantitative genetic analyses implies that the genetic parameters such as heritability of the dependent trait or genetic correlation are influenced by both the genetic variability of the independent trait and the recursive effect between the 2 traits. Moreover, it must be noted that the heritability and correlation estimates under MTM may change if the environmental circumstances determining this nongenetic causal link are modified.

In most cases, RM are statistically equivalent to MTM but, because RM have additional parameters, their implementation requires imposing constraints that imply several strict assumptions. Essentially, there are 2 types of constraints, either to impose independence of the residuals between traits or to dispose of an auxiliary variable. In animal breeding, both of those assumptions are difficult to meet. The first requires that the only cause for the statistical association between the 2 traits is phenotypic causality and that there are no environmental effects acting on both traits, simultaneously. The second constrain requires the disposition of an auxiliary variable (systematic effect or covariate) that affects the independent trait without any direct effect of the dependent trait.

The estimates of variance components and breeding values can be transformed from an RM to an MTM, but the biological interpretation differs. In MTM, the breeding values predict the complete effects of genes on the traits and those values should be used for selection, whereas, in RM, breeding values reflect the additive genetic effects while holding constant the causal trait. If there is a recursive effect between 2 (or more) traits, it may be preferred to consider it in an RM framework for a more reliable prediction of the breeding values. The implementation of RM allows the splitting of the breeding values and the genetic variances and covariances into direct and indirect effects. It ensures a deeper interpretation of the direct and correlated response to selection, mainly under nonlinear recursive effects, or when the recursive effect goes in opposite direction to the additive genetic covariance. Furthermore, RM allow a more thorough understanding of the distribution of the genetic effects throughout the genome.

This review presented the particularities of the models, their interpretation, and some interesting and novel applications. Among these applications we mentioned (1) the statistical equivalence of RM and MTM can be invoked to obtain the estimates of the recursive relationships from the MTM estimates, (2) RM can be expanded to describe more complex relationships between traits including heterogeneous causal structures and nonlinear relationship between traits, (3) RM can be expanded to models that introduce some degree of regularization in the recursive structure that aims to estimate a very large number of recursive parameters, and (4) RM can be used for operational reasons because of their simplicity for implementation, even though there is no causality between traits. In conclusion, the use of RM in animal breeding is expected to increase due to their properties for causal inference and the possibility to use them as operational models. A plethora of applications of RM are expected with the advent of multiomic (including phenomics) data. We hope that this review will contribute to the development of future applications of RM in animal breeding.

## ACKNOWLEDGMENTS

This study received no external funding. The authors wish to acknowledge the valuable suggestions of the reviewers and the help of Bruce MacWhirter (Madrid, Spain) in the English edition. The authors have not stated any conflicts of interest.

## APPENDIX

Example of the calculation the additive and residual (co)variances matrices in a recursive model (RMM) from the estimates under a standard mixed models (SMM) in traits recorded sequentially at 3 time points.

Let us suppose that we dispose phenotypes for 6 traits recorded at time 1 (y

_{1},y_{2}), 2 (y_{3},y_{4}) and 3 (y_{5},y_{6}) and there is causal relationship between the traits recorded at time 1 with the traits recorded at times 2 and 3, and between traits recorded at time 2 with the recorded at time 3 as is described in Figure A1.Let us suppose that we dispose of the following estimates of the additive (

First, we use a first Block LDL decomposition as:

where

then

Now, we use a new Block LDL decomposition for ${\mathbf{R}}_{1}$ as:

where

Then, the residual (co)variance matrix with the RMM (

The ${\mathrm{\Lambda}}^{-1}$ matrix is

and

Therefore, the recursive effects between the traits recorded at time 1 and the traits recorded at time 2 are 0.143, the recursive effects between traits recorded at time 1 and traits recorded at time 3 are 0.033, and between traits recorded at time 2 with the recorded at time 3 are 0.133.

**G***) and the residual (**R***) (co)variances matrices achieved from a standard multitrait mixed model (SMM):${\mathbf{G}}^{\ast}=\left[\begin{array}{cccccc}2& 1& 1& 1& 1& 1\\ 1& 2& 1& 1& 1& 1\\ 1& 1& 2& 1& 1& 1\\ 1& 1& 1& 2& 1& 1\\ 1& 1& 1& 1& 2& 1\\ 1& 1& 1& 1& 1& 2\end{array}\right]$

${\mathbf{R}}^{\ast}=\left[\begin{array}{cccccc}5& 2& 1& 1& 0.5& 0.5\\ 2& 5& 1& 1& 0.5& 0.5\\ 1& 1& 5& 2& 1& 1\\ 1& 1& 2& 5& 1& 1\\ 0.5& 0.5& 1& 1& 5& 2\\ 0.5& 0.5& 1& 1& 2& 5\end{array}\right].$

First, we use a first Block LDL decomposition as:

${\mathbf{R}}^{\ast}=\left[\begin{array}{cc}\mathbf{A}& {\mathbf{B}}^{\prime}\\ \mathbf{B}& \mathbf{C}\end{array}\right]=\left[\begin{array}{cc}\mathbf{I}& 0\\ \mathbf{B}{\mathbf{A}}^{-1}& \mathbf{I}\end{array}\right]\left[\begin{array}{cc}\mathbf{A}& 0\\ 0& \mathbf{C}-\mathbf{B}{\mathbf{A}}^{-1}{\mathbf{B}}^{\prime}\end{array}\right]\left[\begin{array}{cc}\mathbf{I}& {\mathbf{A}}^{-1}{\mathbf{B}}^{\prime}\\ 0& \mathbf{I}\end{array}\right]={\mathrm{\Lambda}}_{1}^{-1}{\mathbf{R}}_{1}{\mathrm{\Lambda}}_{1}^{-{1}^{\prime}},$

where

$\mathbf{A}=\left[\begin{array}{cc}5& 2\\ 2& 5\end{array}\right]\phantom{\rule{thickmathspace}{0ex}}\mathbf{B}=\left[\begin{array}{cc}1& 1\\ 1& 1\\ 0.5& 0.5\\ 0.5& 0.5\end{array}\right]\phantom{\rule{thickmathspace}{0ex}}\mathbf{C}=\left[\begin{array}{cccc}5& 2& 1& 1\\ 2& 5& 1& 1\\ 1& 1& 5& 2\\ 1& 1& 2& 5\end{array}\right],$

then

$\begin{array}{c}{\mathbf{R}}^{\ast}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0.143& 0.143& 1& 0& 0& 0\\ 0.143& 0.143& 0& 1& 0& 0\\ 0.071& 0.071& 0& 0& 1& 0\\ 0.071& 0.071& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccccc}5& 2& 0& 0& 0& 0\\ 2& 5& 0& 0& 0& 0\\ 0& 0& 4.714& 1.714& 0.857& 0.857\\ 0& 0& 1.714& 4.714& 0.857& 0.857\\ 0& 0& 0.857& 0.857& 4.929& 1.928\\ 0& 0& 0.857& 0.857& 1.928& 4.929\end{array}\right]\\ \left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0.143& 0.143& 1& 0& 0& 0\\ 0.143& 0.143& 0& 1& 0& 0\\ 0.071& 0.071& 0& 0& 1& 0\\ 0.071& 0.071& 0& 0& 0& 1\end{array}\right]={\mathrm{\Lambda}}_{1}^{-1}{\mathbf{R}}_{1}{\mathrm{\Lambda}}_{2}^{-{1}^{\prime}}.\end{array}$

Now, we use a new Block LDL decomposition for ${\mathbf{R}}_{1}$ as:

$\begin{array}{c}{\mathbf{R}}_{1}=\left[\begin{array}{cc}{\mathbf{A}}_{1}& {{\mathbf{B}}^{\prime}}_{1}^{}\\ {\mathbf{B}}_{1}& {\mathbf{C}}_{1}\end{array}\right]=\left[\begin{array}{cc}\mathbf{I}& 0\\ {\mathbf{B}}_{1}{\mathbf{A}}_{1}^{-1}& \mathbf{I}\end{array}\right]\left[\begin{array}{cc}{\mathbf{A}}_{1}& 0\\ 0& {\mathbf{C}}_{1}-{\mathbf{B}}_{1}{\mathbf{A}}_{1}^{-1}{{\mathbf{B}}^{\prime}}_{1}^{}\end{array}\right]\left[\begin{array}{cc}\mathbf{I}& {\mathbf{A}}_{1}^{-1}{{\mathbf{B}}^{\prime}}_{1}^{}\\ 0& \mathbf{I}\end{array}\right]=\\ {\mathrm{\Lambda}}_{2}^{-1}{\mathbf{R}}_{2}{\mathrm{\Lambda}}_{2}^{-{1}^{\prime}},\end{array}$

where

${\mathbf{A}}_{1}=\left[\begin{array}{cccc}5& 2& 0& 0\\ 2& 5& 0& 0\\ 0& 0& 4.714& 1.714\\ 0& 0& 1.714& 4.714\end{array}\right]\phantom{\rule{thickmathspace}{0ex}}{\mathbf{B}}_{1}=\left[\begin{array}{cccc}0& 0& 0.857& 0.857\\ 0& 0& 0.857& 0.857\end{array}\right]\phantom{\rule{thickmathspace}{0ex}}{\mathbf{C}}_{1}=\left[\begin{array}{cc}4.929& 1.928\\ 1.928& 4.929\end{array}\right]$

$\begin{array}{c}{\mathbf{R}}_{1}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0.133& 0.133& 1& 0\\ 0& 0& 0.133& 0.133& 0& 1\end{array}\right]\left[\begin{array}{cccccc}5& 2& 0& 0& 0& 0\\ 2& 5& 0& 0& 0& 0\\ 0& 0& 4.714& 1.714& 0& 0\\ 0& 0& 1.714& 4.714& 0& 0\\ 0& 0& 0& 0& 4.7& 1.7\\ 0& 0& 0& 0& 1.7& 4.7\end{array}\right]\\ \left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0.133& 0.133\\ 0& 0& 0& 1& 0.133& 0.133\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]={\mathrm{\Lambda}}_{2}^{-1}{\mathbf{R}}_{2}{\mathrm{\Lambda}}_{2}^{-{1}^{\prime}}.\end{array}$

Then, the residual (co)variance matrix with the RMM (

**R**) is$\mathbf{R}={\mathbf{R}}_{2}=\left[\begin{array}{cccccc}5& 2& 0& 0& 0& 0\\ 2& 5& 0& 0& 0& 0\\ 0& 0& 4.714& 1.714& 0& 0\\ 0& 0& 1.714& 4.714& 0& 0\\ 0& 0& 0& 0& 4.7& 1.7\\ 0& 0& 0& 0& 1.7& 4.7\end{array}\right].$

The ${\mathrm{\Lambda}}^{-1}$ matrix is

$\begin{array}{c}{\mathrm{\Lambda}}^{-1}={\mathrm{\Lambda}}_{1}^{-1}{\mathrm{\Lambda}}_{2}^{-1}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0.143& 0.143& 1& 0& 0& 0\\ 0.143& 0.143& 0& 1& 0& 0\\ 0.071& 0.071& 0& 0& 1& 0\\ 0.071& 0.071& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0.133& 0.133& 1& 0\\ 0& 0& 0.133& 0.133& 0& 1\end{array}\right]\\ \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0.143& 0.143& 1& 0& 0& 0\\ 0.143& 0.143& 0& 1& 0& 0\\ 0.071& 0.071& 0.133& 0.133& 1& 0\\ 0.071& 0.071& 0.133& 0.133& 0& 1\end{array}\right]\end{array}$

and

$\mathrm{\Lambda}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ -0.143& -0.143& 1& 0& 0& 0\\ -0.143& -0.143& 0& 1& 0& 0\\ -0.033& -0.033& -0.133& -0.133& 1& 0\\ -0.033& -0.033& -0.133& -0.133& 0& 1\end{array}\right].$

Therefore, the recursive effects between the traits recorded at time 1 and the traits recorded at time 2 are 0.143, the recursive effects between traits recorded at time 1 and traits recorded at time 3 are 0.033, and between traits recorded at time 2 with the recorded at time 3 are 0.133.

Finally, the additive (co)variance matrix with the RMM (

**G**) is calculated as:$\mathbf{G}=\mathrm{\Lambda}{\mathbf{G}}^{\ast}{\mathrm{\Lambda}}^{\prime}=$

$\left[\begin{array}{c}100000\\ 010000\\ -0.143-0.1431000\\ -0.143-0.1430100\\ -0.033-0.033-0.133-0.13310\\ -0.033-0.033-0.133-0.13301\end{array}\right]\left[\begin{array}{c}211111\\ 121111\\ 112111\\ 111211\\ 111121\\ 111112\end{array}\right]\left[\begin{array}{c}10-0.143-0.143-0.033-0.033\\ 01-0.143-0.143-0.003-0.033\\ 0010-0.133-0.133\\ 0001-0.133-0.133\\ 000010\\ 000001\end{array}\right]$

$=\left[\begin{array}{cccccc}2& 1& 0.571& 0.571& 0.633& 0.633\\ 1& 2& 0.571& 0.571& 0.633& 0.633\\ 0.571& 0.571& 1.551& 0.551& 0.352& 0.352\\ 0.571& 0.571& 0.551& 1.551& 0.352& 0.352\\ 0.633& 0.633& 0.352& 0.352& 1.482& 0.482\\ 0.633& 0.633& 0.352& 0.352& 0.482& 1.482\end{array}\right].$

## REFERENCES

- Testing model nesting and equivalence.
*Psychol. Methods.*2010; 15 (20515234): 111-123 - Cell biology symposium: Genetics of feed efficiency in dairy and beef cattle.
*J. Anim. Sci.*2013; 91 (23345557): 1594-1613 - Exploring causal networks of bovine milk fatty acids in a multivariate mixed model context.
*Genet. Sel. Evol.*2014; 46 (24438068): 2 - Beyond genomic selection: The animal model strikes back (one generation)!.
*J. Anim. Breed. Genet.*2017; 134 (28508480): 224-231 - Estimation and comparison of multiple change-point models.
*J. Econom.*1998; 86: 221-241 - Genetic evaluation including intermediate omics features.
*Genetics.*2021; 219 (34849886)iyab130 - Image analysis and computer vision applications in animal sciences: An overview.
*Front. Vet. Sci.*2020; 7 (33195522)551269 - Nonlinear genetic relationships between traits and their implications on the estimation of genetic parameters.
*J. Anim. Sci.*1997; 75 (9419984): 3119-3125 - Non-linearity in the genetic relationship between milk yield and type traits in Holstein cattle.
*Livest. Prod. Sci.*1998; 57: 41-47 - Sampling-based approaches to calculating marginal densities.
*J. Am. Stat. Assoc.*1990; 85: 398-409 - Theory and analysis of threshold characters.
*J. Anim. Sci.*1982; 54: 1079-1096 - Quantitative genetic models for describing simultaneous and recursive relationships between phenotypes.
*Genetics.*2004; 167 (15280252): 1407-1424 - Non-linear recursive models for growth traits in the Pirenaica beef cattle breed.
*Animal.*2014; 8 (24673822): 904-911 - Monte Carlo sampling methods using Markov chains and their applications.
*Biometrika.*1970; 57: 97-109 - Applications of Linear Models in Animal Breeding.University of Guelph, 1984
- Inferring relationships between health and fertility in Norwegian Red cows using recursive models.
*J. Dairy Sci.*2009; 92 (19307661): 1778-1784 - The growth of structural equation modeling: 1994–2001.
*Struct. Equ. Modeling.*2003; 10: 35-46 - An application of change-point recursive models to the relationship between litter size and number of stillborns in pigs.
*J. Anim. Sci.*2010; 88 (20675604): 3493-3503 - Application of Bayesian causal inference and structural equation model to animal breeding.
*Anim. Sci. J.*2020; 91 (32219948)e13359 - Inferring causal structures and comparing the causal effects among calving difficulty, gestation length, and calf size in Japanese Black cattle.
*Animal.*2017; 11 (28478794): 2120-2128 - Relationships between milk yield and somatic cell score in Canadian Holsteins from simultaneous and recursive random regression models.
*J. Dairy Sci.*2010; 93 (20172242): 1216-1233 - Short communication: Recursive model approach to traits defined as ratios: Genetic parameters and breeding values.
*J. Dairy Sci.*2017; 100 (28284690): 3767-3772 - Genomic evaluation for feed efficiency in Canadian Holsteins.
*Interbull Bull.*2021; 56: 153-161 - Recursive relationships between milk yield and somatic cell score of Canadian Holsteins from finite mixture random regression models.
*J. Dairy Sci.*2010; 93 (20965363): 5474-5486 - Alternative parameterizations of the multiple-trait random regression model for milk yield and somatic cell score via recursive links between phenotypes.
*J. Anim. Breed. Genet.*2011; 128 (21749472): 258-266 - Restricted selection indices.
*Biometrics.*1959; 15: 10-19 - Efficiency of feed use in beef cattle.
*J. Anim. Sci.*1963; 22: 486-494 - Genetic and phenotypic relationships among milk urea nitrogen, fertility, and milk yield in Holstein cows.
*J. Dairy Sci.*2008; 91 (18946143): 4372-4382 - Investigation of influence of growing pigs' positive affective state on behavioral and physiological parameters using structural equation modeling.
*J. Anim. Sci.*2020; 98 (31999319)skaa028 - Structural equation modeling and whole-genome scans uncover chromosome regions and enriched pathways for carcass and meat quality in beef.
*Front. Genet.*2018; 9 (30555508): 532 - Symposium review: Omics in dairy and animal science—Promise, potential, and pitfalls.
*J. Dairy Sci.*2019; 102 (30268604): 4741-4754 - Analysis of fertility and dystocia in Holsteins using recursive models to handle censored and categorical data.
*J. Dairy Sci.*2007; 90 (17369243): 2012-2024 - Exploring biological relationships between calving traits in primiparous cattle with a Bayesian recursive model.
*Genetics.*2009; 181 (18984571): 277-287 - Modeling relationships between calving traits: A comparison between standard and recursive mixed models.
*Genet. Sel. Evol.*2010; 42 (20100345): 1 - Utilizing trait networks and structural equation models as tools to interpret multi-trait genome-wide association studies.
*Plant Methods.*2019; 15 (31548847): 107 - Heritable environmental variance causes nonlinear relationships between traits: Application to birth weight and stillbirth of pigs.
*Genetics.*2015; 199 (25631318): 1255-1269 - Development of a structural growth curve model that considers the causal effect of initial phenotypes.
*Genet. Sel. Evol.*2019; 51 (31046678): 19 - Causality: Models, Reasoning, and Inference.2nd ed. Cambridge University Press, 2009
- Searching for causal networks involving latent variables in complex traits: Application to growth, carcass, and meat quality in pigs.
*J. Anim. Sci.*2015; 93: 4617-4623 - The relative contribution of causal factors in the transition from infection to clinical chlamydial disease.
*Sci. Rep.*2018; 8 (29891934)8893 - Inferring relationships between clinical mastitis, productivity, and fertility: A recursive model application including genetics, farm associated herd management, and cow-specific antibiotic treatments.
*Prev. Vet. Med.*2013; 112 (23859301): 58-67 - Breeding and genetics symposium: Inferring causal effects from observational data in livestock.
*J. Anim. Sci.*2013; 91 (23230107): 553-564 - Estimating causal effects of treatments in randomized and nonrandomized studies.
*J. Educ. Psychol.*1974; 66: 688-701 - Structural equation models to disentangle the biological relationship between microbiota and complex traits: Methane production in dairy cattle as a case of study.
*J. Anim. Breed. Genet.*2020; 137 (31617268): 36-48 - A dimensional reduction approach to modulate the core ruminal microbiome associated with methane emissions via selective breeding.
*J. Dairy Sci.*2021; 104 (33896632): 8135-8151 - Alternative strategies for genetic analyses of milk flow in dairy cattle.
*J. Dairy Sci.*2015; 98 (26364101): 8209-8222 - Cause and Correlation in Biology.Cambridge University Press, 2000
- Marginal maximum likelihood estimation of variance components in Poisson mixed models using Laplacian integration.
*Genet. Sel. Evol.*1993; 25: 305-319 - Causal relationships between milk quality and coagulation properties in Italian Holstein-Friesian dairy cattle.
*Genet. Sel. Evol.*2015; 47 (25968045): 45 - Searching for recursive causal structures in multivariate quantitative genetics mixed models.
*Genetics.*2010; 185 (20351220): 633-644 - Is structural equation modeling advantageous for the genetic improvement of multiple traits?.
*Genetics.*2013; 194 (23608193): 561-572 - Searching for phenotypic causal networks involving complex traits: An application to European quail.
*Genet. Sel. Evol.*2011; 43 (22047591): 37 - On Bayesian structural inference in a simultaneous equation model.in: Stigum B.P. Econometrics and the Philosophy of Economics. Princeton University Press, 2003: 642-682
- APPENDIX_A.docx. figshare. Journal contribution.
- A cross-specific multiplicative binomial recursive model for the analysis of perinatal mortality in a diallel cross among three varieties of Iberian pig.
*Sci. Rep.*2020; 10 (33273670)21190 - A genetic analysis of mortality in pigs.
*Genetics.*2010; 184 (19901070): 277-284 - Joint analysis of binomial and continuous traits with a recursive model: A case study using mortality and litter size of pigs.
*Genetics.*2014; 196 (24414548): 643-651 - Analysis of litter size and average litter weight in pigs using a recursive model.
*Genetics.*2007; 177 (17720909): 1791-1799 - A multiple-trait Bayesian variable selection regression method for integrating phenotypic causal networks in genome-wide association studies.
*G3 (Bethesda).*2020; 10 (33020191): 4439-4448 - Correlation and causation.
*J. Agric. Res.*1921; 20: 557-585 - Inferring relationships between somatic cell score and milk yield using simultaneous and recursive models.
*J. Dairy Sci.*2007; 90 (17582135): 3508-3521 - Exploration of lagged relationships between mastitis and milk yield in dairy cows using a Bayesian structural equation Gaussian-threshold model.
*Genet. Sel. Evol.*2008; 40 (18558070): 333-357 - Bayesian structural equation models for inferring relationships between phenotypes: A review of methodology, identifiability, and applications.
*J. Anim. Breed. Genet.*2010; 127 (20074182): 3-15 - An alternative interpretation of residual feed intake by phenotypic recursive relationships in dairy cattle.
*JDS Commun.*2021; 2: 371-375

## Article info

### Publication history

Published online: March 02, 2023

Accepted:
October 30,
2022

Received:
July 26,
2022

### Publication stage

In Press Corrected Proof### Identification

### Copyright

© 2023 The Author(s).

### User license

Creative Commons Attribution (CC BY 4.0) | How you can reuse

Elsevier's open access license policy

Creative Commons Attribution (CC BY 4.0)

## Permitted

- Read, print & download
- Redistribute or republish the final article
- Text & data mine
- Translate the article
- Reuse portions or extracts from the article in other works
- Sell or re-use for commercial purposes

Elsevier's open access license policy