EM Algorithm for Mapping Quantitative Trait Loci in Multivalent Tetraploids

Abstract

Multivalent tetraploids that include many plant species, such as potato, sugarcane, and rose, are of paramount importance to agricultural production and biological research. Quantitative trait locus (QTL) mapping in multivalent tetraploids is challenged by their unique cytogenetic properties, such as double reduction. We develop a statistical method for mapping multivalent tetraploid QTLs by considering these cytogenetic properties. This method is built in the mixture model-based framework and implemented with the EM algorithm. The method allows the simultaneous estimation of QTL positions, QTL effects, the chromosomal pairing factor, and the degree of double reduction as well as the assessment of the estimation precision of these parameters. We used simulated data to examine the statistical properties of the method and validate its utilization. The new method and its software will provide a useful tool for QTL mapping in multivalent tetraploids that undergo double reduction.

1. Introduction

Genetic analysis in polyploids has received considerable interest in recent years because of the biological and economic importance [1–3]. Genetic linkage maps constructed from molecular markers have been published for several major polyploids [4–10]. Statistical models for linkage analysis and map construction that consider unique biological properties of polyploids have been developed [11–14]. For bivalent polyploids, Wu et al. [15, 16] incorporated the so-called chromosomal pairing preference [17] into the linkage analysis framework, to increase the biological relevance of linkage mapping models. There have been several statistical models developed to map quantitative trait loci (QTLs) in bivalent polyploids [18, 19].

There is also a group of polyploids, called multivalent polyploids, in which chromosomes pair among more than two homologous copies at meiosis, rather than only two copies as like in bivalent polyploids. The origin of multivalent polyploids is mostly from the duplication of similar genomes and, for this reason, they are called autopolyploids [20, 21]. The consequence of multivalent pairing in autopolyploids is the occurrence of double reduction, that is, two sister chromatids of a chromosome sort into the same gamete [22]. Fisher [23] proposed a conceptual model for characterizing the individual probabilities of 11 different modes of gamete formation for a quadrivalent polyploid in terms of the recombination fraction between two different loci and their double reductions. Wu et al. [24] used Fisher's model to derive the EM algorithm for the estimation of the linkage between fully informative markers. Wu and Ma [25] extended this model to analyze any type of markers, regardless of their informativeness and dominant or codominant nature. The significant advantage of the models by Wu and colleagues directly lies in their generality, flexibility, and robustness.

In this paper, we develop a statistical method for QTL mapping in multivalent tetraploids by considering Fisher's [23] 11 classifications of gamete formation. The method allows the estimation and test of not only the QTL-marker linkage, but also the extent of double reduction of the QTL. Because of the inherent complexity of classification analyses of gamete formation, we will focus on the modeling and analysis of one-marker/one-QTL associations. A two-stage hierarchical model is derived to estimate the probabilities of gamete formation modes and therefore double reduction in the upper hierarchy and estimate the marker-QTL recombination fraction in the lower hierarchy within the maximum likelihood context implemented with the EM algorithm. The method is used to analyze a simulated data set, with results demonstrating statistical properties of the method and its analytical and biological merits.

2. Method

2.1. Genetic Design

Consider a heterozygous multivalent tetraploid line crossed with a homozygous line to generate a so-called pseudotest backcross population. For such a population, the genotypes of progeny are consistent with the genotypes of gametes produced by the heterozygous parent and, therefore, the derivation of mapping models can be based on the segregation of gametes. Suppose there are 𝑛 individuals in the pseudotest backcross population. A panel of codominant markers is typed for each individual, with which a linkage map is constructed. All the pseudotest backcross individuals are phenotyped for a quantitative trait that is assumed to be controlled by QTLs.

2.2. Tetrasomic Co-Inheritance

To simplify our analysis, we assume that the QTLs underlying the trait are mapped with single markers. Let 𝑀₁,…,𝑀₄ be the four alleles at a marker 𝐌, and let 𝑄₁,…,𝑄₄ be the four alleles at a QTL linked with the marker. The marker and QTL are linked with a recombination fraction of 𝑟. Because of the double reduction, the multivalent tetraploid generates 10 diploid gametes for locus, which are arrayed as (𝑀₁𝑀₁, 𝑀₂𝑀₂, 𝑀₃𝑀₃, 𝑀₄𝑀₄, 𝑀₁𝑀₂, 𝑀₁𝑀₃, 𝑀₁𝑀₄, 𝑀₂𝑀₃, 𝑀₂𝑀₄, 𝑀₃𝑀₄) for the markers, and (𝑄₁𝑄₁, 𝑄₂𝑄₂, 𝑄₃𝑄₃, 𝑄₄𝑄₄, 𝑄₁𝑄₂, 𝑄₁𝑄₃, 𝑄₁𝑄₄, 𝑄₂𝑄₃, 𝑄₂𝑄₄, 𝑄₃𝑄₄) for the QTL. In each case, the first four gametes are derived from the double reduction, whereas the second six gametes are derived from the chromosome paring. Let𝛼 and𝛽 be the frequencies of double reduction at the marker and QTL, respectively. The frequency of double reduction is a constant for any given locus, with the value depending on its distance from the centromere.

When the marker and QTL are co-segregating in a multivalent tetraploid, a total of 136 diploid gamete formation mechanisms are generated although there are only 100 gamete genotypes that are observable. Based on the presence/absence of double reduction and the number of recombinant events, Fisher [23] classified these 136 formation mechanisms into 11 gamete modes. Of these 11 gamete modes, however, only nine can be observed each with a frequency denoted by 𝑔_ℎ(ℎ=1,…,9). These 9 observable gamete modes were rearranged by Wu et al. [24] in matrix form expressed as𝑄₁𝑄₁𝑄₂𝑄₂𝑄₃𝑄₃𝑄₄𝑄₄𝑄₁𝑄₂𝑄₁𝑄₃𝑄₁𝑄₄𝑄₂𝑄₃𝑄₂𝑄₄𝑄₃𝑄₄𝑀₁𝑀₁𝑀₂𝑀₂𝑀₃𝑀₃𝑀₄𝑀₄𝑀𝐠=¹𝑀₂𝑀₂𝑀₁𝑀₁𝑀₃𝑀₃𝑀₁𝑀₁𝑀₄𝑀₄𝑀₁𝑀₂𝑀₃𝑀₃𝑀₂𝑀₂𝑀₄𝑀₄𝑀₂𝑀₃𝑀₄𝑀₄𝑀₃14𝑔₁1𝑔12²1𝑔12²1𝑔12²1𝑔12⁵1𝑔12⁵1𝑔12⁵1𝑔12⁶1𝑔12⁶1𝑔12⁶1𝑔12²14𝑔₁1𝑔12²1𝑔12²1𝑔12⁵1𝑔12⁶1𝑔12⁶1𝑔12⁵1𝑔12⁵1𝑔12⁶1𝑔12²1𝑔12²14𝑔₁1𝑔12²1𝑔12⁶1𝑔12⁵1𝑔12⁶1𝑔12⁵1𝑔12⁶1𝑔12⁵1𝑔12²1𝑔12²1𝑔12²14𝑔₁1𝑔12⁶1𝑔12⁶1𝑔12⁵1𝑔12⁶1𝑔12⁵1𝑔12⁵1𝑔12³1𝑔12³1𝑔12⁴1𝑔12⁴16𝑔₇1𝑔24⁸1𝑔24⁸1𝑔24⁸1𝑔24⁸16𝑔₉1𝑔12³1𝑔12⁴1𝑔12³1𝑔12⁴1𝑔24⁸16𝑔₇1𝑔24⁸1𝑔24⁸16𝑔₉1𝑔24⁸1𝑔12³1𝑔12⁴1𝑔12⁴1𝑔12³1𝑔24⁸1𝑔24⁸16𝑔₇16𝑔₉1𝑔24⁸1𝑔24⁸1𝑔12⁴1𝑔12³1𝑔12³1𝑔12⁴1𝑔24⁸1𝑔24⁸16𝑔₉16𝑔₇1𝑔24⁸1𝑔24⁸1𝑔12⁴1𝑔12³1𝑔12⁴1𝑔12³1𝑔24⁸16𝑔₉1𝑔24⁸1𝑔24⁸16𝑔₇1𝑔24⁸1𝑔12⁴1𝑔12⁴1𝑔12³1𝑔12³16𝑔₉1𝑔24⁸1𝑔24⁸1𝑔24⁸1𝑔24⁸16𝑔₇, (1)where 𝑔₁ and 𝑔₂ are associated with double reductions at both the marker and QTL, 𝑔₃ and 𝑔₄ with double reductions only at QTL 𝐐, 𝑔₅ and 𝑔₆ with double reductions only at marker 𝐌, and 𝑔₇−𝑔₉ with nondouble reductions. From matrix (1), we see that there are no, one and two recombinant events in the cells (𝑔₁), (𝑔₃,𝑔₅), and (𝑔₂,𝑔₄,𝑔₆,𝑔₉), respectively. The cells (𝑔₇) and (𝑔₈) are each a mixture of two different gamete formation mechanisms or configurations (𝐴 and 𝐵), that is, 𝑔₇=𝑔_7𝐴+𝑔_7𝐵 and 𝑔₈=𝑔_8𝐴+𝑔_8𝐵, with relative proportions determined by 𝑟. Because different configurations contain different numbers of recombination events, the expected number of recombination events in each cell, that is, an observable gamete genotype, should be the weighted average of the number of recombination events for each configuration. Wu et al. [24] used a matrix form (e) to count the expected number of recombination events for each observable gamete genotype expressed as𝑄𝐜=₁𝑄₁𝑄₂𝑄₂𝑄₃𝑄₃𝑄₄𝑄₄𝑄₁𝑄₂𝑄₁𝑄₃𝑄₁𝑄₄𝑄₂𝑄₃𝑄₂𝑄₄𝑄₃𝑄₄𝑀₁𝑀₁𝑀₂𝑀₂𝑀₃𝑀₃𝑀₄𝑀₄𝑀₁𝑀₂𝑀₂𝑀₁𝑀₁𝑀₃𝑀₃𝑀₁𝑀₁𝑀₄𝑀₄𝑀₁𝑀₂𝑀₃𝑀₃𝑀₂𝑀₂𝑀₄𝑀₄𝑀₂𝑀₃𝑀₄𝑀₄𝑀₃,022211122220221221122202212121222022121111222𝜙1+𝜓1+𝜓1+𝜓1+𝜓212121+𝜓2𝜙1+𝜓1+𝜓21+𝜓12211+𝜓1+𝜓2𝜙21+𝜓1+𝜓21121+𝜓1+𝜓22𝜙1+𝜓1+𝜓21211+𝜓21+𝜓1+𝜓2𝜙1+𝜓221121+𝜓1+𝜓1+𝜓1+𝜓2𝜙 (2)where𝑟𝜙=²10𝑟²,𝑟−18𝑟+9𝜓=.3−2𝑟 (3) Based on matrices (1) and (2), the expressions for the frequencies of double reduction (𝛼 and 𝛽) and the recombination fraction 𝑟 can be expressed in terms of 𝑔_𝑖 as𝛼=𝑔₁+𝑔₂+𝑔₃+𝑔₄,𝛽=𝑔₁+𝑔₂+𝑔₅+𝑔₆,1𝑟=2𝑔₃+𝑔₅+2𝑔₂+𝑔₄+𝑔₆+𝑔₉+2𝜙𝑔₇+1+𝜓𝑔₈. (4)

2.3. Quantitative Genetic Model

For a given QTL, there are 10 different QTL gamete genotypes in the multivalent tetraploid, whose values can be partitioned into additive and dominance genetic effects of different types, expressed as𝜇₁₁=𝜇+𝑎₁,for𝑄₁𝑄₁,𝜇₂₂=𝜇+𝑎₂,for𝑄₂𝑄₂,𝜇₃₃=𝜇+𝑎₃,for𝑄₃𝑄₃,𝜇₄₄=𝜇−𝑎₁−𝑎₂−𝑎₃,for𝑄₄𝑄₄,𝜇₁₂=𝜇+𝑎₁+𝑎₂+𝑑₁₂,for𝑄₁𝑄₂,𝜇₁₃=𝜇+𝑎₁+𝑎₃+𝑑₁₃,for𝑄₁𝑄₃,𝜇₁₄=𝜇−𝑎₂−𝑎₃+𝑑₁₄,for𝑄₁𝑄₄,𝜇₂₃=𝜇+𝑎₂+𝑎₃+𝑑₂₃,for𝑄₂𝑄₃,𝜇₂₄=𝜇−𝑎₁−𝑎₃+𝑑₂₄,for𝑄₂𝑄₄,𝜇₃₄=𝜇−𝑎₁−𝑎₂+𝑑₃₄,for𝑄₃𝑄₄, (5) where𝜇 is the overall mean, 𝑎₁, 𝑎₂, and 𝑎₃ are the additive genetic effects of alleles 𝑄₁, 𝑄₂, and 𝑄₃ relative to allele 𝑄₄, and 𝑑₁₂, 𝑑₁₃, 𝑑₁₄, 𝑑₂₃, 𝑑₂₄, and 𝑑₃₄ are the dominant genetic effects due to interactions between different alleles 𝑄₁ and 𝑄₂, 𝑄₁ and 𝑄₃, 𝑄₁ and 𝑄₄, 𝑄₂ and 𝑄₃, 𝑄₂ and 𝑄₄, and 𝑄₃ and 𝑄₄, respectively.

From expression (5), we can solve the overall mean and additive and dominant effects as1𝜇=4𝜇₁₁+𝜇₂₂+𝜇₃₃+𝜇₄₄,𝑎₁=143𝜇₁₁−𝜇₂₂−𝜇₃₃−𝜇₄₄𝑎₂,=143𝜇₂₂−𝜇₁₁−𝜇₃₃−𝜇₄₄𝑎₃,=143𝜇₃₃−𝜇₁₁−𝜇₂₂−𝜇₄₄𝑑₁₂,=𝜇₁₂−143𝜇₃₃+3𝜇₄₄−𝜇₁₁−𝜇₂₂𝑑₁₃,=𝜇₁₃−143𝜇₂₂+3𝜇₄₄−𝜇₁₁−𝜇₃₃𝑑₁₄,=𝜇₁₄−143𝜇₂₂+3𝜇₃₃−𝜇₁₁−𝜇₄₄𝑑₂₃,=𝜇₂₃−143𝜇₁₁+3𝜇₄₄−𝜇₂₂−𝜇₃₃𝑑₂₄,=𝜇₂₄−143𝜇₁₁+3𝜇₃₃−𝜇₂₂−𝜇₄₄𝑑₃₄,=𝜇₃₄−143𝜇₁₁+3𝜇₂₂−𝜇₃₃−𝜇₄₄. (6)(7)(8)(9)(10)(11)(12)(13)(14)(15)

2.4. EM Algorithm

Ignoring the effects of other covariates, the phenotypic value, 𝑦_𝑖, for individual 𝑖 in the pseudotest backcross can be expressed in terms of the QTL effect and residual error as𝑦_𝑖=⁴_{𝑗1≤𝑗2=1}𝜉_{𝑗1𝑗2∣𝑖}𝜇_𝑗1𝑗2+𝑒_𝑖, (16) where 𝜉_{𝑗1𝑗2∣𝑖} is the indicator variable that is defined as 1 if individual 𝑖 has a QTL genotype 𝑗₁𝑗₂ (𝑗₁≤𝑗₂=𝑄₁,𝑄₂,𝑄₃,𝑄₄), and 0 otherwise, 𝜇_𝑗1𝑗2 is the genotypic value of QTL genotype 𝑗₁𝑗₂ as defined in (5), and 𝑒_𝑖 is the residual error assumed to be normally distributed with mean zero and variance 𝜎². We use Ω to denote the unknown vector (𝜇₁₁,𝜇₂₂,𝜇₃₃,𝜇₄₄,𝜇₁₂,𝜇₁₃,𝜇₁₄,𝜇₂₃,𝜇₂₄,𝜇₃₄,𝜎²).

For a QTL mapping experiment, marker genotypes are observable. Let 𝑛_𝑙1𝑙2 be the observation of marker genotype 𝑙₁𝑙₂ (𝑙₁≤𝑙₂=𝑀₁,𝑀₂,𝑀₃,𝑀₄). The likelihood of the phenotypic (𝐲) and marker data (𝐌) is constructed, within the mixture model framework, as𝐿(𝛀∣𝐲,𝐌)=^𝑛𝑀_𝑖=1⁴_𝑙1=𝑀1^𝑀4_𝑙2=𝑀1^𝑄4_𝑗1=𝑄1^𝑄4_𝑗2=𝑄1𝜋_{𝑗1𝑗2∣𝑙1𝑙2}𝑓_𝑗1𝑗2𝑦_𝑖,𝑙₁≤𝑙₂,𝑗₁≤𝑗₂, (17) where 𝜋_{𝑗1𝑗2∣𝑙1𝑙2} is the conditional probability of QTL genotype 𝑗₁𝑗₂ given marker genotype 𝑙₁𝑙₂, and 𝑓_𝑗1𝑗2(𝑦_𝑖) is assumed to follow a normal distribution with mean 𝜇_𝑗1𝑗2 and variance 𝜎². Prior conditional probability 𝜋_{𝑗1𝑗2∣𝑙1𝑙2} is calculated as the frequency of joint marker-QTL genotype 𝑙₁𝑙₂𝑗₁𝑗₂, expressed in terms of nine 𝑔 probabilities in matrix (1), divided by the frequency of marker genotype 𝑙₁𝑙₂. Marker genotype frequencies are 𝛼/4 for each of double reduction gametes 𝑀₁𝑀₁, 𝑀₂𝑀₂, 𝑀₃𝑀₃, and 𝑀₄𝑀₄, and (1−𝛼)/6 for each of nondouble reduction gametes 𝑀₁𝑀₂, 𝑀₁𝑀₃, 𝑀₁𝑀₄, 𝑀₂𝑀₃, 𝑀₂𝑀₄, and 𝑀₃𝑀₄.

The estimates of unknown parameters that maximize the likelihood (17) can be obtained by implementing the EM algorithm. In step E, we calculate the posterior probability of a QTL genotype given a specific marker genotype of individual 𝑖 byΘ_{𝑗1𝑗2∣𝑙1𝑙2𝑖}=𝜋_{𝑗1𝑗2∣𝑙1𝑙2}𝑓_𝑗1𝑗2𝑦_𝑖∑^𝑄4_{𝑗′1=𝑄1}∑^𝑄4_{𝑗′2=𝑄1}𝜋_{𝑗12𝑗′2∣𝑙1𝑙2}𝑓_{𝑗′1𝑗′2}𝑦_𝑖. (18)

In step M, we calculate the frequencies of nine observable gamete modes based on the calculated posterior probabilities using the following: 𝑔₁=12𝑛₁₁+𝑛₂₂+𝑛₃₃+𝑛₄₄×^𝑛11_𝑖=1Θ_11∣11𝑖+^𝑛22_𝑖=1Θ_22∣22𝑖+^𝑛33_𝑖=1Θ_33∣33𝑖+^𝑛44_𝑖=1Θ_44∣44𝑖,𝑔₂=12𝑛₁₁+𝑛₂₂+𝑛₃₃+𝑛₄₄×^𝑛11_𝑖=1Θ_22∣11𝑖+Θ_33∣11𝑖+Θ_44∣11𝑖+^𝑛22_𝑖=1Θ_11∣22𝑖+Θ_33∣22𝑖+Θ_44∣22𝑖+^𝑛33_𝑖=1Θ_11∣33𝑖+Θ_22∣33𝑖+Θ_44∣33𝑖+^𝑛44_𝑖=1Θ_11∣44𝑖+Θ_22∣44𝑖+Θ_33∣44𝑖,𝑔₃=12𝑛₁₂+𝑛₁₃+𝑛₁₄+𝑛₂₃+𝑛₂₄+𝑛₃₄×^𝑛12_𝑖=1Θ_11∣12𝑖+Θ_22∣12𝑖+^𝑛13_𝑖=1Θ_11∣13𝑖+Θ_33∣13𝑖+^𝑛14_𝑖=1Θ_11∣14𝑖+Θ_44∣14𝑖+^𝑛23_𝑖=1Θ_22∣23𝑖+Θ_33∣23𝑖+^𝑛24_𝑖=1Θ_22∣24𝑖+Θ_44∣24𝑖+^𝑛34_𝑖=1Θ_33∣34𝑖+Θ_44∣34𝑖,𝑔₄=12𝑛₁₂+𝑛₁₃+𝑛₁₄+𝑛₂₃+𝑛₂₄+𝑛₃₄×^𝑛12_𝑖=1Θ_33∣12𝑖+Θ_44∣12𝑖+^𝑛13_𝑖=1Θ_22∣13𝑖+Θ_44∣13𝑖+^𝑛14_𝑖=1Θ_22∣14𝑖+Θ_33∣14𝑖+^𝑛23_𝑖=1Θ_11∣23𝑖+Θ_44∣23𝑖+^𝑛24_𝑖=1Θ_11∣24𝑖+Θ_33∣24𝑖+^𝑛34_𝑖=1Θ_11∣34𝑖+Θ_22∣34𝑖,𝑔₅=12𝑛₁₁+𝑛₂₂+𝑛₃₃+𝑛₄₄×^𝑛11_𝑖=1Θ_12∣11𝑖+Θ_13∣11𝑖+Θ_14∣11𝑖+^𝑛22_𝑖=1Θ_12∣22𝑖+Θ_23∣22𝑖+Θ_24∣22𝑖+^𝑛33_𝑖=1Θ_13∣33𝑖+Θ_23∣33𝑖+Θ_34∣33𝑖+^𝑛44_𝑖=1Θ_14∣44𝑖+Θ_24∣44𝑖+Θ_34∣44𝑖,𝑔₆=12𝑛₁₁+𝑛₂₂+𝑛₃₃+𝑛₄₄×^𝑛11_i=1Θ_23∣11𝑖+Θ_24∣11𝑖+Θ_34∣11𝑖+^𝑛22_𝑖=1Θ_13∣22𝑖+Θ_14∣22𝑖+Θ_34∣22𝑖+^𝑛33_𝑖=1Θ_12∣33𝑖+Θ_14∣33𝑖+Θ_24∣33𝑖+^𝑛44_𝑖=1Θ_12∣44𝑖+Θ_13∣44𝑖+Θ_23∣44𝑖,𝑔₇=12𝑛₁₂+𝑛₁₃+𝑛₁₄+𝑛₂₃+𝑛₂₄+𝑛₃₄×^𝑛12_𝑖=1Θ_12∣12𝑖+^𝑛13_𝑖=1Θ_13∣13𝑖+^𝑛14_𝑖=1Θ_14∣14𝑖+^𝑛23_𝑖=1Θ_23∣23𝑖+^𝑛24_𝑖=1Θ_24∣24𝑖+^𝑛34_𝑖=1Θ_34∣34𝑖,𝑔₈=12𝑛₁₂+𝑛₁₃+𝑛₁₄+𝑛₂₃+𝑛₂₄+𝑛₃₄×^𝑛12_𝑖=1Θ_13∣12𝑖+Θ_14∣12𝑖+Θ_23∣12𝑖+Θ_24∣12𝑖+^𝑛13_𝑖=1Θ_12∣13𝑖+Θ_14∣13𝑖+Θ_23∣12𝑖+Θ_34∣12𝑖+^𝑛14_𝑖=1Θ_12∣14𝑖+Θ_13∣14𝑖+Θ_24∣14𝑖+Θ_34∣14𝑖+^𝑛23_𝑖=1Θ_12∣23𝑖+Θ_13∣23𝑖+Θ_24∣23𝑖+Θ_34∣23𝑖+^𝑛24_𝑖=1Θ_12∣24𝑖+Θ_14∣24𝑖+Θ_23∣24𝑖+Θ_34∣24𝑖+^𝑛34_𝑖=1Θ_13∣34𝑖+Θ_14∣34𝑖+Θ_23∣34𝑖+Θ_24∣34𝑖,𝑔₉=12𝑛₁₂+𝑛₁₃+𝑛₁₄+𝑛₂₃+𝑛₂₄+𝑛₃₄×^𝑛12_𝑖=1Θ_34∣12𝑖+^𝑛13_𝑖=1Θ_24∣13𝑖+^𝑛14_𝑖=1Θ_23∣14𝑖+^𝑛23_𝑖=1Θ_14∣23𝑖+^𝑛24_𝑖=1Θ_13∣24𝑖+^𝑛34_𝑖=1Θ_12∣34𝑖 (19) which lead to the estimates of the frequencies of double reduction as𝛼=𝑔₁+𝑔₂+𝑔₃+𝑔₄=1𝑛𝑛₁₁+𝑛₂₂+𝑛₃₃+𝑛₄₄,𝛽=𝑔₁+𝑔₂+𝑔₅+𝑔₆,̂1𝑟=2𝑔₃+𝑔₅+2𝑔₂+𝑔₄+𝑔₆+𝑔₉+2𝜙𝑔₇+1+𝜓𝑔₈. (20)

The genotypic value of QTL genotype 𝑗₁𝑗₂ and residual variance are estimated by𝜇_𝑗1𝑗2=∑^𝑛_i=1∑^𝑀4_𝑙1=𝑀1∑^𝑀4_𝑙2=𝑀1Θ_{𝑗1𝑗1∣𝑙1𝑙2𝑖}𝑦_𝑖∑^𝑀4_𝑙1=𝑀1∑^𝑀4_𝑙2=𝑀1Θ_{𝑗1𝑗1∣𝑙1𝑙2𝑖},𝑙₁≤𝑙₂;𝑗₁≤𝑗₂,𝜎²=1𝑛^𝑛𝑀_𝑖=1⁴_𝑙1=𝑀1^𝑀4_𝑙2=𝑀1^𝑄4_𝑗1=𝑄1^𝑄4_𝑗2=𝑄1Θ_{𝑗1𝑗1∣𝑙1𝑙2𝑖}𝑦_𝑖−𝜇_𝑗1𝑗2²,𝑙₁≤𝑙₂;𝑗₁≤𝑗₂. (21)(22) The iteration is repeated between the E step, (3) and (18), and M step, (19)–(22), until stable estimates are obtained. The stable estimates are the maximum likelihood estimates (MLEs) of parameters.

2.5. Hypothesis Testing

Following parameter estimation, several hypotheses should be tested. The hypothesis about the presence of a QTL segregating in the pseudotest backcrosses is formulated as𝐻₀:𝑎₁=𝑎₂=𝑎₃=𝑑₁₂=𝑑₁₃=𝑑₁₄=𝑑₂₃=𝑑₂₄=𝑑₃₄𝐻=0,₁:atleastoneofthemisnotequaltozero. (23) The difference between the log-likelihood functions under the null and alternative hypotheses are calculated. But the distribution of this log-likelihood ratio (LR) is not known because of the violation of regularity conditions for the mixture model (1). For this reason, a commonly used empirical approach based on permutation tests by reshuffling the relationships between the marker genotypes and phenotypes [26] is used to determine the critical threshold, in order to judge whether there is a QTL for the trait.

After a significant QTL is detected, the next hypothesis is about the additive genetic effect of the QTL. This can be tested by formulating the null hypothesis, 𝐻₀:𝑎₁=𝑎₂=𝑎₃=0, (24) under which the estimates of genotypic values of QTL genotypes can be obtained with the EM algorithm as described above, but posing three constraints derived from (7), (8), and (9). Similarly, the dominant genetic effects can be tested with the null hypothesis, 𝐻₀:𝑑₁₂=𝑑₁₃=𝑑₁₄=𝑑₂₃=𝑑₂₄=𝑑₃₄=0, (25) with estimates of genotypic values under the constraints derived from (10)–(15). All these genetic effects can be tested individually.

2.6. Application to Simulated Data

A pseudotest backcross for a multivalent tetraploid was hypothesized, in which a marker is assumed to be linked with a QTL that affects a quantitative trait. Marker and QTL genotypes were simulated for the pseudotest backcross of different sample sizes (𝑛=100,200,400) based on a range of double reduction (0.05, 0.15, 0.30) and recombination fraction (0.05, 0.25). We assume the same frequency of double reduction between the marker and QTL. The phenotypic value of an individual is expressed as the summation of genotypic values of a QTL genotype carried by this individual and a normally distributed error. The genotypic values of a QTL genotype are calculated by (5), where the overall mean𝜇 is assigned as 1, and the additive and dominant effects assigned as 𝑎₁=𝑎₂=𝑎₃=0.6 and 𝑑₁₂=𝑑₁₃=𝑑₁₄=𝑑₂₃=𝑑₂₄=𝑑₃₄=0.5. The error variance is determined according to the heritability of 𝐻²=0.1 and0.4, respectively.

In this simulation study, fully informative markers and QTL are assumed and, thus, the double reduction at the marker can be estimated analytically. The estimates of the parameters converge to stable values at a rapid rate given that there are closed forms for parameter estimators in the EM framework. We evaluate the estimation of the other parameters related to QTL segregation, effects, and position. The means of the MLEs of the QTL-related parameters and their standard errors based on 1000 simulation replicates are illustrated in Tables 1, 2, and 3. With a small sample size (100), the double reduction of the QTL was accurately estimated, with the precision of estimation relatively independent of the magnitude of heritability and the degree of QTL-marker linkage (Table 1). The most significant factor that affected the estimate of QTL position (in terms of its recombination with the marker) was the heritability, followed by sample size and the degree of QTL-marker linkage. In general, at least a sample size of 200 was required to reliably estimate the QTL position for a major gene that explains about 20%–30% of the phenotypic variance.

Table 1 

Maximum likelihood estimates of the double reduction, recombination fraction, overall mean, additive effects, and dominant effects under different simulation scenarios in a mapping population of size 100. Numbers in parentheses are the standard errors of the estimates.

Table 2 

Maximum likelihood estimates of the double reduction, recombination fraction, overall mean, additive effects, and dominant effects under different simulation scenarios in a mapping population of size 200. Numbers in parentheses are the standard errors of the estimates.

Table 3 

Maximum likelihood estimates of the double reduction, recombination fraction, overall mean, additive effects, and dominant effects under different simulation scenarios in a mapping population of size 400. Numbers in parentheses are the standard errors of the estimates.

The estimation precision of the QTL effects depended on the heritability, sample size, and degree of QTL-marker linkage. As heritability, sample size, and linkage degree increased, the estimates of various QTL effects were more precise. As compared with the dominant genetic effects, the estimates of the additive genetic effects required a larger sample size, more precise phenotypic measurements (leading to a higher heritability), and a denser linkage map (with a stronger degree of QTL-marker linkage). We found that the estimates of QTL effects were influenced by the frequency of QTL double reduction. At low frequencies of double reduction, the effects of QTL were more accurately estimated than at higher frequencies. For a QTL undergoing a strong double reduction (say 𝛼=0.3), a sample size of at least 400 is required even if the QTL explains a large proportion of the phenotypic variance (0.4). For a modest-sized QTL, a much larger sample size was required.

We performed a simulation study to test how the misspecification of double reduction affects the estimate of QTL-related parameters. This was done by using traditional mapping models (without considering double reduction) to analyze the simulated data of QTL genotypes with different degrees of double reduction. When a QTL undergoes double reduction, traditional models that do not consider double reduction provided misleading results about the estimates of QTL effects and position (data not shown). Furthermore, increasing heritabilities and sample sizes did not improve the estimates. In this case, the power of QTL detection was reduced.

3. Discussion

A statistical method for genetic mapping of quantitative trait loci (QTLs) in a multivalent tetraploid undergoing a double reduction process is described. As an important cytological characteristic of polyploids, double reduction may play a significant role in plant evolution and maintenance of genetic polymorphism in natural populations. Also, because double reduction affects the result of linkage analysis through the crossing-over events between different chromosomes [24, 25], it is important to incorporate double reduction into a QTL mapping framework. This method provides a powerful tool for QTL mapping and understanding the genetic control of a quantitative trait in a multivalent tetraploid.

The method capitalizes on 11 different classifications of two-locus gamete formations, derived by Fisher [23], during multivalent tetraploid meiosis and has proven to be powerful for simultaneous estimation of the frequencies of double reduction and the recombination fraction between different loci. Although a couple of statistical approaches have been proposed to map multivalent tetraploid QTLs [26, 27], this method has for the first time incorporated Fisher's tetrasomic inheritance into the mapping framework, thus enhancing the cytological relevance of QTL detection. Results from simulation studies showed that the method can be used to map QTLs in a controlled cross of multivalent tetraploids when the mapping population is adequately large (say 400). When a QTL undergoes double reduction, traditional mapping approaches will incorrectly estimate the position and effects of the QTL, proportional to the degree of double reduction. The new method can estimate the double reduction of a QTL, an important parameter related to the genetic diversity and evolution of polyploids [28, 29].

Because of the high complexity of the mixture model implemented with tetrasomic inheritance, we only considered a one-marker model for QTL mapping. Interval mapping, which localizes a QTL with two flanking markers, has proven to be more advantageous in parameter estimation over the one-marker model [30]. It will be worthwhile to integrate components of our model into the interval mapping framework to fully explore the statistical merits of interval mapping for QTL mapping in multivalent tetraploids. Furthermore, the model proposed in this article assumes the segregation of fully informative codominant loci, each with 10 distinct genotypes, in a controlled cross of multivalent tetraploids. For partially informative codominant markers, a two-stage hierarchical mixture model will be needed to model the different allelic configurations for a phenotypically identical genotype. Although molecular marker technologies have improved in recent years, dominant markers may still be used in genetic mapping projects of some underrepresentative species including polyploids. Thus, it is also important to extend our model to map QTLs with dominant markers. For partially informative loci, the number of QTL genotypes may be unknown and, thus, a model selection procedure should be incorporated to determine the optimal number of genotypes at a QTL.

The genetic mapping of polyploids is complex because of their complex inheritance modes. Sophisticated statistical models are required to tackle genetic problems hidden in the polysomic inheritance of polyploids. Currently, there are some debates on the optimal modeling of tetrasomic inheritance in linkage analysis [13, 25] and QTL mapping [18, 31] partly because of our limited knowledge about these fascinating species. Before a detailed understanding of the cytological mechanisms for meioses in multivalent polyploids is obtained, this type of debate will continue. In any case, the development of powerful statistical models for polyploid mapping continues to be a pressing need. The application of these models to real-world data will not only test their usefulness, but also provide an unprecedented opportunity to understand the genetic differentiation among polyploid genomes and characterize the genetic architecture of quantitatively inherited traits for this unique group of species. Software for the method described is available at http://statgen.psu.edu/.

Acknowledgments

This work is supported by Joint grant DMS/NIGMS-0540745. Additional support is through the Office of Sciences (BER), U.S. Department of Energy, Interagency Agreement no. DE-A102-07ER64453.