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} = \mu + e_{j} + c_{cj} + x_{ip} * \beta_{pj} + \underline{ge}_{icj} + \underline{\epsilon}_{icj} \quad [1]$$

where $\underline{y}_{icj}$ is the genotype mean (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_{ip}$ is the number of allele copies coming from parent $p$, 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 one parent (the most frequently used) set as the reference in all environments (e.g. the central or recurrent parent in NAM).

The term $ge_{icj}$ is the residual genetic variation and $\epsilon_{icj}$ the plot error term. In model 1 because genotype values are averaged over replication $\epsilon_{icj}$ and $ge_{icj}$ are confounded. The variance of the $(ge_{icj} + \epsilon_{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_{\epsilon}^{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_{\epsilon}^{2} & 0 & \sigma_{g}^{2} & 0 & \sigma_{g}^{2} & 0 \\ 0 & \sigma_{g}^{2} + \sigma_{\epsilon}^{2} & 0 & \sigma_{g}^{2} & 0 & \sigma_{g}^{2} \\ \sigma_{g}^{2} & 0 & \sigma_{g}^{2} + \sigma_{\epsilon}^{2} & 0 & \sigma_{g}^{2} & 0 \\ 0 & \sigma_{g}^{2} & 0 & \sigma_{g}^{2} + \sigma_{\epsilon}^{2} & 0 & \sigma_{g}^{2} \\ \sigma_{g}^{2} & 0 & \sigma_{g}^{2} & 0 & \sigma_{g}^{2} + \sigma_{\epsilon}^{2} & 0 \\ 0 & \sigma_{g}^{2} & 0 & \sigma_{g}^{2} & 0 & \sigma_{g}^{2} + \sigma_{\epsilon}^{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_{\epsilon}^{2} & 0 & \sigma_{g_{12}} & 0 & \sigma_{g_{13}} & 0 \\ 0 & \sigma_{g_{1}}^{2} + \sigma_{\epsilon}^{2} & 0 & \sigma_{g_{12}} & 0 & \sigma_{g_{13}} \\ \sigma_{g_{12}} & 0 & \sigma_{g_{2}}^{2} + \sigma_{\epsilon}^{2} & 0 & \sigma_{g_{23}} & 0 \\ 0 & \sigma_{g_{12}} & 0 & \sigma_{g_{2}}^{2} + \sigma_{\epsilon}^{2} & 0 & \sigma_{g_{23}} \\ \sigma_{g_{13}} & 0 & \sigma_{g_{23}} & 0 & \sigma_{g_{3}}^{2} + \sigma_{\epsilon}^{2} & 0 \\ 0 & \sigma_{g_{13}} & 0 & \sigma_{g_{23}} & 0 & \sigma_{g_{3}}^{2} + \sigma_{\epsilon}^{2} \\ \end{bmatrix}$$

We considered the UN model as the default option.

The calculation of an ‘exact’ mixed model at each marker (QTL) position of the genome requires a lot of computer power, which can quickly make detection unfeasible, especially with large populations. Therefore, to reduce the computational power needed 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. In the second option, different $\hat{V}$ should be formed by removing the cofactor information that is too close of a tested QTL position.

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]$, $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 ($(N_{par}-1) * N_{env}$) for the main QTL term or one for each specific within environment allele).

In model 1, QTL terms are considered as fixed. Therefore, the user needs to provide best linear unbiased estimates (BLUEs) as trait value. Beyond estimating a general QTL significance, model 1 allows the estimation of QTL allelic effect for each parents in each environment. One parent is set as reference. By default in populations like the nested association mapping (NAM) populations the central parent is considered as the reference. For other types of MPPs the reference parent can be specified using a dedicated function argument (ref_par). 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 QTLs showing a significant QTL by environment (QTLxE) effects, it is possible to investigate which specific dimension of the environment (environmental covariate - EC) explains the QTLxE effect. Therefore, we propose a strategy to estimate the sensitivity of the parental QTL allele to a specific EC (rain, temperature, photoperiod, etc.). For that purpose, we extended the QTL term of model 1 like that $x_{ip} * \beta_{pj} = x_{ip} * (\beta_{p} + EC_{e}*S_{p} + l_{pe})$ where $\beta_{p}$ is the main parental allelic effect across environments, $EC_{e}$ is the EC value in each environment, $S_{p}$ is the parental allele sensitivity to $EC_{e}$ change and $l_{pe}$ the residual effect. The $S_{p}$ determine the rate of change of the parental QTL allelic additive effect given an extra unit of EC. The significance of $S_{p}$ allows to assess the degree of sensitivity of the QTL to the EC. To avoid model saturation, we implemented a procedure that first estimates the QTL parental alleles effect across (main effect term) and between environment (QTLxE term). Then, the function compare the significance of those two terms. The QTL by EC (QTLxEC) is only estimated for the parental allelic effect for which the QTLxE term is more significant than the main effect term.