QTL detection in multiparental populations characterized in multiple environments

Multiparental populations (MPPs) like the nested association mapping (NAM; McMullen et al. (2009)) or multi-parent advanced generation inter-cross (MAGIC; Cavanagh et al. (2008)) populations have emerged as central genetic resources for research and breeding purposes (Scott et al. 2020). The diallels (Blanc et al. 2006), the factorial designs or the collections of crosses developed in breeding programs (Würschum 2012) can also be considered as MPPs. We can classify MPPs given the use of intrecrossing. In a first type of MPPs, there is no intercrossing. Each genotype inherits its alleles from two parents. Those MPPs can be seen as a collection of crosses (e.g. NAM or factorial design). A second type of MPPs like MAGIC involve genotypes that are a mosaic of more than two parents or founders. In this package, we only consider the first type of MPPs without intercrossing or collections of crosses with shared parents.

Even though several statistical methodologies and software are available to perform QTL detection in MPPs, there is potentially few or no open-source solution for the detection of QTL in MPPs characterized in multiple environments (MPP-ME). In this extension of mppR, we propose a method to detect QTLs in MPP-ME data.

We start from the QTL detection model 3 described in (Garin, Malosetti, and Eeuwijk 2020):

\[\underline{y}_{icj} = E_{j} + C_{cj} + x_{i_{q}p} * \beta_{pj} + \underline{GE}_{icj} + \underline{e}_{icj} \quad [1]\]

where \(\underline{y}_{icj}\) is the genotype adjusted mean (e.g. BLUE) of genotype \(i\) from cross \(c\) in environment \(j\). \(E_{j}\) is a fixed environment term, \(C_{cj}\) a fixed cross within environment interaction, \(x_{i_{q}p}\) is the number of allele copies coming from parent \(p\) at QTL position \(q\), and \(\beta_{pj}\) represent the QTL effect of parental allele \(p\) in environment \(j\). Different definitions of the QTL allelic effect are possible (Garin et al. 2017). In this version of the model we assumed that each parent carries a different allele and that each of those alleles can have an environmental specific effect. The model is therefore saturated with the estimation of \(N_{env} * N_{par}\) effects with the most frequent allele set as the reference in all environments (e.g. the central or recurrent parent in NAM). In model 1, QTL terms are considered as fixed. Therefore, the user needs to provide best linear unbiased estimates (BLUEs) as trait value.

The term \(GE_{icj}\) is the residual genetic variation and \(e_{icj}\) the plot error term. In model 1 because genotype values are averaged over replication \(e_{icj}\) and \(GE_{icj}\) are confounded. The variance of the \((GE_{icj} + e_{icj})\) term can be modeled with different variance covariance (VCOV) structures. A first possibility is to assume a constant genotypic variance across environments \(\sigma_{G}^{2}\) (compound symmetry - CS) and a unique variance for the plot error term \(\sigma_{e}^{2}\):

\[V\begin{bmatrix} y_{i11} \\ y_{i'21} \\ y_{i12} \\ y_{i'22} \\ y_{i13} \\ y_{i'23} \end{bmatrix} = \begin{bmatrix} \sigma_{G}^{2} + \sigma_{e}^{2} & 0 & \sigma_{G}^{2} & 0 & \sigma_{G}^{2} & 0 \\ 0 & \sigma_{G}^{2} + \sigma_{e}^{2} & 0 & \sigma_{G}^{2} & 0 & \sigma_{G}^{2} \\ \sigma_{G}^{2} & 0 & \sigma_{G}^{2} + \sigma_{e}^{2} & 0 & \sigma_{G}^{2} & 0 \\ 0 & \sigma_{G}^{2} & 0 & \sigma_{G}^{2} + \sigma_{e}^{2} & 0 & \sigma_{G}^{2} \\ \sigma_{G}^{2} & 0 & \sigma_{G}^{2} & 0 & \sigma_{G}^{2} + \sigma_{e}^{2} & 0 \\ 0 & \sigma_{G}^{2} & 0 & \sigma_{G}^{2} & 0 & \sigma_{G}^{2} + \sigma_{e}^{2} \\ \end{bmatrix}\]

The CS model requires the estimation of a single term (\(\sigma_{G}^{2}\)) to describe genetic the covariance between all pairs of environments. This model assumes that the background polygenic effect (effect of all genome positions except the tested QTL and cofactors) is the same in all environments (Van Eeuwijk et al. 2001).

A more realistic option called unstructured (UN) model allows the genetic covariance to be different in every pairs of environments. From a genetic point of view, it means that a set of genes have a different effect in each environments (Van Eeuwijk et al. 2001). In the unstructured model, each pair of environment get a specific genetic covariance term \(cov(y..j, y..j') = \sigma_{G_{jj'}}^{2}\)

\[V\begin{bmatrix} y_{i11} \\ y_{i'21} \\ y_{i12} \\ y_{i'22} \\ y_{i13} \\ y_{i'23} \end{bmatrix} = \begin{bmatrix} \sigma_{G_{1}}^{2} + \sigma_{e}^{2} & 0 & \sigma_{G_{12}} & 0 & \sigma_{G_{13}} & 0 \\ 0 & \sigma_{G_{1}}^{2} + \sigma_{e}^{2} & 0 & \sigma_{G_{12}} & 0 & \sigma_{G_{13}} \\ \sigma_{G_{12}} & 0 & \sigma_{G_{2}}^{2} + \sigma_{e}^{2} & 0 & \sigma_{G_{23}} & 0 \\ 0 & \sigma_{G_{12}} & 0 & \sigma_{G_{2}}^{2} + \sigma_{e}^{2} & 0 & \sigma_{G_{23}} \\ \sigma_{G_{13}} & 0 & \sigma_{G_{23}} & 0 & \sigma_{G_{3}}^{2} + \sigma_{e}^{2} & 0 \\ 0 & \sigma_{G_{13}} & 0 & \sigma_{G_{23}} & 0 & \sigma_{G_{3}}^{2} + \sigma_{e}^{2} \\ \end{bmatrix}\]

We considered the UN model as the default option.

The QTL detection procedure is composed of the following steps:

1. SIM scan
2. Selection of cofactors
3. CIM scan
4. Selection of QTLs
5. Estimation of the QTL allele effects
6. Estimation of the QTL R2 contribution
7. Plot of the QTL profile
8. Visualisation of the allelic significance along the genome

The function mppGE_proc() perform the whole procedure detailed above. For further details, see the examples below.

The calculation of an ‘exact’ mixed model at each marker (QTL) position of the genome is computationally demanding, which can quickly make detection unfeasible, especially with large populations. Therefore, to reduce the computation demand for a genome scan, we implemented an approximate mixed model computation similar to the generalized least square strategy implemented in Kruijer et al. (2015) or Garin et al. (2017). The procedure consists of estimating a general VCOV (\(\hat{V}\)) using model 1 without the tested QTL position for the simple interval mapping (SIM) scan. For the composite interval mapping (CIM) model, \(\hat{V}\) is estimated including the cofactors as fixed effect without the QTL position. For the CIM scan, there will be as many \(\hat{V}\) as the number of cofactors combinations given cofactor exclusion around the tested position.

There are two option for the estimation of \(\hat{V}\). In the first option we estimate a unique \(\hat{V}\) taking all cofactors into consideration (VCOV_data = "unique"). In the second option, different \(\hat{V}\) should be formed by removing the cofactor information that is too close of a tested QTL position (VCOV_data = "minus_cof").

The statistical significance of the tested QTL position and the allelic effects is obtained by substituting \(\hat{V}\) to get the following Wald statistic

\[W_Q = \beta^{T} [V(\beta)]^{-1} \beta \quad [2]\]

where,

\[\beta = (X^{T} \hat{V}^{-1} X)^{-1} X^T \hat{V}^{-1} y \quad [3]\]

and

\[V(\beta) = (X^{T} \hat{V}^{-1} X)^{-1}\]

\(X\) represents the fixed effect matrix including the cross effect, the cofactors, and the QTL position, and \(y\) is the vector of genotypic values. \(W_Q \sim \chi_{d}^{2}\) with degree of freedom (\(d\)) equal to the number of tested QTL allelic effects (e.g. \((N_{par}-1) * N_{env}\) for a QTL term with environment-specific allelic effects).

It is possible to estimate the QTL parental allelic effect using model 1. Model 1 allows the estimation of QTL allelic effect for each parents in each environment. One parent is set as reference. By default, it is the most frequent parental allele. Therefore, in populations like the nested association mapping (NAM) populations, the central parent is considered as the reference. The reference parent can be specified by the user through the ref_par argument. The parental QTL allelic effect must be interpreted as the extra effect of one allele copy with respect to the reference parent in the considered environment.

For the parental allelic effects showing a significant interaction with the environment in model 4, it is possible to extend the model to test the sensitivity of the QTL allele to a specific environmental covariate (\(EC_j\), e.g. temperature or humidity). For that purpose, we can rewrite the QTL effect term \(\beta_{pj}\) of the parents with a significant environmental interaction (pxE). We can replace \(\beta_{pj}\) with \(EC_j*S_p+l_{p\epsilon}\). \(EC_j\) represents the EC value in environment j associated with the sensitivity term \(S_p\). The \(S_{p}\) determines the rate of change of the parental QTL allelic additive effect given an extra unit of EC. Finally, \(l_{p\epsilon}\) is a residual effect. The fitted model becomes:

\[\underline{y}_{icj} = E_{j} + C_{cj} + \sum_{q=1}^{n_{QTL}} x_{i_{q}p} * (\alpha_p + \beta_{pj}) + x_{i_{q}pxE} * (\alpha_p + EC_j*S_p+l_{p\epsilon}) + \underline{GE}_{icj} + \underline{e}_{icj} \quad [5]\] The significance of the sensitivity term \(S_p\) can be evaluated with the Wald test. Contrary to the model fitted for the QTL detection scan, model 1 for QTL effect estimation, models 4 and 5 are fitted using an ‘exact’ REML solution using the R pakage nlme.