## 1. Introduction

Reproductive traits of animals determine the efficiency of food production, and hence are economically important traits for genetic improvement in the livestock industry (Rocha *et al.*, Reference Rocha, Eisen, Siewerdt, Van Vleck and Pomp2004*b*). However, since the reproductive traits usually have a low heritability, their improvement has met with limited success (Hansen *et al.*, Reference Hansen, Freeman and Berger1983; Avalos & Smith, Reference Avalos and Smith1987; Bennett & Leymaster, Reference Bennett and Leymaster1989; Lamberson, Reference Lamberson and Young1990). Fortunately, along with the recent advances of genomic technologies and the collective efforts of statistical genetics, methodologies have been introduced to identify quantitative trait loci (QTLs) regulating reproductive traits, and ultimately enhance genetic processes by using QTL-complemented breeding and selection strategies in animal food production (Rathje *et al.*, Reference Rathje, Rohrer and Johnson1997; Rohrer *et al.*, Reference Rohrer, Ford, Wise, Vallet and Christenson1999; Wilkie *et al.*, Reference Wilkie, Paszek, Beattie, Alexander, Wheeler and Schook1999; Cassady *et al.*, Reference Cassady, Johnson, Pomp, Rohrer, Van Vleck, Spiegel and Gilson2001; King *et al.*, Reference King, Jiang, Gibson, Haley and Archibald2003; Holl *et al.*, Reference Holl, Cassady, Pomp and Johnson2004; Rocha *et al.*, Reference Rocha, Eisen, Siewerdt, Van Vleck and Pomp2004*b*).

Most QTL mapping studies for counted reproductive traits have used interval mapping or composite interval mapping, which focus on detecting marginal (main) effects of QTL (Kirkpatrick *et al.*, Reference Kirkpatrick, Mengelt, Schulman and Martin1998; Spearow *et al.*, Reference Spearow, Nutson, Mailliard, Porter and Barkley1999; Rocha *et al*., 2004 *b*; Sato *et al.*, Reference Sato, Atsuji, Saito, Okitsu, Sato, Komatsuda, Mitsuhashi, Nirasawa, Hayashi, Sugimoto and Kobayashi2006; Ren *et al.*, Reference Ren, Ren, Xing, Guo, Wu, Yang, Mao and Huang2009). Epistatic interaction effects (epistasis) have not been thoroughly examined, even though the genetic architecture of most quantitative traits in organisms from plants to animals including humans is quite complex, including inputs from multiple QTLs and epistasis (Doebley & Stec, Reference Doebley and Stec1991; Cheverund & Routman, Reference Cheverund and Routman1995; Lynch & Walsh, Reference Lynch and Walsh1998; Price & Courtois, Reference Price and Courtois1999; Sugiyama *et al.*, Reference Sugiyama, Churchill, Higgins, Johns, Makaritsis, Gavras and Paigen2001; Wade, Reference Wade2001; Kopp *et al.*, Reference Kopp, Graze, Xu, Carroll and Nuzhdin2003; Carlborg & Haley, Reference Carlborg and Haley2004). Furthermore, some studies have shown that epistatically acting QTLs play a larger role than do the marginal effect QTLs in genetic modulation and evolution of quantitative traits (Yu *et al.*, Reference Yu, Li, Xu, Tan, Gao, Li, Zhang and Saghai Maroof1997; Yi *et al.*, Reference Yi, Diamont, Chiu, Kim, Allison, Fisler and Warden2004*a*, Reference Yi, Chiu, Allison, Fisler and Warden*b*, Reference Yi, Zinniel, Kim, Eisen, Bartolucci, Allison and Pomp2006; Moore, Reference Moore2005; Valdar *et al.*, Reference Valdar, Solberg, Gauguier, Cookson, Rawlins, Mott and Flint2006). Therefore, to better understand the genetic architecture of counted reproductive traits, statistical analyses are required to simultaneously identify main effects and epistatic interactions of multiple QTLs.

Yi & Banerjee (Reference Yi and Banerjee2009) proposed a unified methodology for mapping multiple interacting QTLs based upon the hierarchical generalized linear model framework. The key to the approach is the use of a continuous prior distribution on coefficients that favours sparseness in the fitted models and facilitates computation. The method of Yi & Banerjee (Reference Yi and Banerjee2009) can analyse various continuous and discrete phenotypes, and can fit a large number of effects, including covariates, main effects of numerous loci and epistatic and gene–environment interactions. Yi & Banerjee (Reference Yi and Banerjee2009) illustrated the hierarchical generalized linear model framework with normal and probit models for continuous and binary traits. In this study, we developed hierarchical Poisson and binomial models for mapping multiple interacting QTLs for complex count traits in experimental crosses based on the approach of Yi & Banerjee (Reference Yi and Banerjee2009).

We applied our methods to data collected from a QTL mapping study of female reproductive traits in mice (Rocha *et al.*, Reference Rocha, Eisen, Van Vleck and Pomp2004*a*, Reference Rocha, Eisen, Siewerdt, Van Vleck and Pomp*b*). The mouse sample was derived from parental lines having undergone long-term selection for high and low body weight gain. The previous analysis of these mouse data adopted the composite interval mapping method and detected 15 QTLs on five chromosomes for six reproductive traits. However, epistatic effects on the traits were not evaluated. In addition, two discrete counted reproductive traits, live fetuses (*LF*) and dead fetuses (*DF*), were treated as continuous traits. Improved models for these two count phenotypes should be Poisson models with observed number of *corpora lutea* (ovulation rate (OR)) treated as the baseline or *offset*, or binomial models treating *LF* and *DF* as the numbers of ‘success’ out of the corresponding number of *corpora lutea* (trials). The observed number of OR provides us the information about the upper bound of *LF* and *DF.* Therefore, to comprehensively and properly explore the genetic architecture of these reproductive traits in mice, we reanalysed the data by using the hierarchical Poisson and the binomial regression models. The aims of the present study were to estimate the epistatic effects of interacting QTLs on reproductive traits, and to detect possible additional QTLs for these traits.

## 2. Materials and methods

The mouse population, phenotypes and marker genotyping have been presented in detail in Rocha *et al.* (Reference Rocha, Eisen, Van Vleck and Pomp2004*a*, Reference Rocha, Eisen, Siewerdt, Van Vleck and Pompb), but are briefly described here.

### (i) Mouse lines and crosses

A total of 439 F_{2} female mice were bred from two long-term selected lines, the high-growth M16i line and the low-body-weight L6 line. The former was derived from an outbred population, whereas the latter was derived from a cross of four inbred lines. M16i females were intercrossed with L6 males, and the resulting F_{1} mice were *inter se* mated (no full-sib pairings) in two consecutive replicates encompassing a total of 64 full-sib F_{2} families. These mice were reared at 21°C in a 12:12 h light:dark cycle and 55% relative humidity. Food and water were supplied *ad libitum*. All the animals were handled according to the Institutional Animal Care and Use Committee Guidelines.

### (ii) Phenotypes and environmental effects

To measure the reproductive traits of female mice, 10-week-old F_{2} females were exposed to unrelated F_{1} males (B6C3F1/J) until a copulatory plug was detected. Pregnant females were subsequently euthanized at day 17 of gestation to obtain three counted reproductive phenotypes: *LF*, *DF* and OR (also called *corpora lutea*) (Fig. 1). Body weights at 10 weeks of age (*WK10*) were also measured, which are significantly correlated with the reproductive phenotypes (Rocha *et al.*, Reference Rocha, Eisen, Siewerdt, Van Vleck and Pomp2004*b*). Thus, our models adjusted for the covariates *WK10* as well as the family indicator as environmental factors.

### (iii) Genotyping

A total of 63 fully informative microsatellite markers spanning the 19 autosomes were genotyped (Rocha *et al.*, Reference Rocha, Eisen, Van Vleck and Pomp2004*a*, Reference Rocha, Eisen, Siewerdt, Van Vleck and Pomp*b*). Marker genotypes were determined by PCR and agarose gel electrophoresis protocols. Segregation distortion was evaluated by Chi-square test. Detection and correction of genotyping errors were conducted with MAPMAKER. A linkage map was generated with MAPMAKER/EXP and QTL analysis was carried out after marker distances were estimated. The marker linkage map covered 1257·8 cM (Kosambi) with an average spacing of 30 cM. A total of 21 (0·04%) marker genotypes were missing.

### (iv) Statistical analyses

#### (a) Hierarchical Poisson and binomial models

Poisson and binomial models were used to detect QTL for the numbers of *LF* and *DF*. The Poisson model treated the number of *LF* and *DF* as counted responses relative to the number of *OR*, and the binomial model treated the number of *LF* and *DF* as the events of ‘success’ out of the corresponding ‘trials’ (the number of *OR*). Under each of these models, two approaches were applied: the first considering only main effects of QTL, and the second considering not only main effects but also their epistatic effects. All the models included *WK10* and the family indicator as continuous and categorical covariates, respectively.

Denoting the number of *LF* or *DF* by *y _{i}* and the number of

*OR*by

*T*, for mouse

_{i}*i*=1, …,

*n*, the Poisson epistatic QTL model can be expressed as

and the binomial epistatic QTL model can take the form

where Poisson(*T _{i}*θ

*) presents Poisson distribution with mean*

_{i}*T*θ

_{i}*, and Bin(*

_{i}*T*) represents binomial distribution with size

_{i},p_{i}*T*and probability

_{i}*p*;

_{i}

**X**_{E},

**X**_{G}and

**X**

_{GG}are design matrices for the environmental effects, main effects (additive and dominance) and epistatic effects (additive–additive, additive–dominance, dominance–additive and dominance–dominance), respectively; β

_{0}is the intercept, and

**β**

_{E},

**β**

_{G}and

**β**

_{GG}are vectors of the environmental, main and epistatic effects, respectively. In the Poisson model, θ

*is the underlying rate per unit of*

_{i}*OR*for the

*i*th mouse, and the logarithm of θ

*is expressed as a linear predictor of the environmental effects, main and epistatic effects of potential QTL. In the binomial model,*

_{i}*p*is the probability of each

_{i}*OR*being live or dead for the

*i*th mouse, and the logit transformation is used to relate

*p*to the environmental effects, main and epistatic effects of potential QTLs. The Poisson and binomial models are different but both reasonable for the present data analyses. It would be expected that these two models may each enable us to detect some different QTLs.

_{i} The genetic variables were coded by using a modified Cockerham genetic model (Yi & Banerjee, Reference Yi and Banerjee2009). The dominance variable was coded as −0·5 and 0·5 for homozygotes and heterozygote of each QTL, respectively. To achieve the same scale as that of the dominance contrast, the additive variable was coded as −2^{−0·5}, 0 and 2^{−0·5} for the three genotypes of F_{2} population, rather than the original Cockerham codes of −1, 0 and 1. The epistasis **X**_{GG} was constructed by multiplying two corresponding main-effect variables. For each environmental variable, the raw values were transformed to have a mean of 0 and a standard deviation of 0·5 (Gelman *et al.*, Reference Gelman, Jakulin, Pittau and Su2008; Yi & Banerjee, Reference Yi and Banerjee2009). This transformation standardized all the environmental effects to the scale of all the genetic main effects described above. For loci with missing genotypic values, the codes of contrasts were replaced by their conditional expectations given the observed marker data (Haley & Knott, Reference Haley and Knott1992). Although the proposed models can include loci between the observed markers, we describe our methods by considering observed markers as potential QTLs (e.g. Xu, Reference Xu2007; Yi & Banerjee, Reference Yi and Banerjee2009).

A total of 63 observed markers lead to a total of 7938 effects, including 63 additive effects, 63 dominance effects and 7812 epistatic effects. Most of these effects are expected to be zero or at least negligible. To incorporate this notion into the model and to make the large-scale model identifiable, an independent Student-*t* prior was assumed on coefficients β* _{j}* for

*j*=1, …,

*J*(Yi & Banerjee, Reference Yi and Banerjee2009). For the intercept β

_{0}, a weakly informative prior with ν

_{0}=1 and

*s*

_{0}=10 was used, and for main and epistatic effects, ν

*=0·01 and*

_{j}*s*=10

_{j}^{−4}were assigned. These shrinkage priors give each coefficient a high probability of being near zero while still allowing for occasionally large effects (Yi & Xu, Reference Yi and Xu2008; Yi & Banerjee, Reference Yi and Banerjee2009). To facilitate the estimation of parameters, the

*t*distribution was expressed as a mixture of normal distributions with mean 0 and variance distributed as scaled inverse-χ

^{2}:

#### (b) Model fit and search algorithm

The computational algorithm of Yi & Banerjee (Reference Yi and Banerjee2009) as implemented in the freely available R package qtlbim (Yandell *et al.*, Reference Yandell, Mehta, Banerjee, Shriner, Venkataraman, Moon, Neely, Wu, von Smith and Yi2007) was used to fit the Poisson and binomial models by estimating the posterior modes of the coefficients. Yi & Banerjee (Reference Yi and Banerjee2009) developed a procedure to fit generalized linear models with the Student-*t* prior by incorporating an expectation–maximum (EM) algorithm into the usual iteratively weighted least squares (IWLS). The IWLS algorithm approximates the generalized linear model by a weighted normal linear regression (Gelman *et al.*, Reference Gelman, Carlin, Stern and Rubin2003). At each iteration, one calculates pseudo-datum *z _{i}* and pseudo-variance σ

_{i}^{2}for each observation

*i*based on the current estimates of parameters, approximates the generalized linear model likelihood

*p*(

*y*

_{i}|**X**

_{i}**β**,φ) by the normal likelihood

*N*(

*z*

_{i}|**X**

_{i}**β**,σ

_{i}^{2}), and then updates the parameters β

*by a weighted linear regression. The EM algorithm of Yi & Banerjee (Reference Yi and Banerjee2009) treats variances τ*

_{j}

_{j}^{2}as missing data and replaces them by their posterior means at each E-step. Given the variances τ

_{j}^{2}, the parameters β

_{j}are estimated by including the priors of β

*τ*

_{j}|

_{j}^{2}~

*N*(0,τ

_{j}^{2}) as additional ‘data points’ in the normal likelihood

*N*(

*z*

_{i}|**X**

_{i}**β**,σ

_{i}^{2}). Briefly, the algorithm starts with initial values for each τ

_{j}^{2}and β

*, and then proceeds as follows: (1) based on the current values of β*

_{j}*, calculate pseudo-datum*

_{j}*z*and pseudo-variance σ

_{i}

_{i}^{2}; (2) E-step: replace each variance τ

_{j}^{2}by its conditional posterior expectation; (3) M-step: perform the weighted least square regression based on the normal likelihood approximation to obtain estimates and (4) repeat steps 1–3 until convergence.

The hierarchical model can simultaneously handle many correlated variables (Xu, Reference Xu2007; Yi & Banerjee, Reference Yi and Banerjee2009). However, fitting a model with thousands of variables requires large memory and intensive computation. Thus, the model search strategy proposed by Yi & Banerjee (Reference Yi and Banerjee2009) was used to build a parsimonious model. The approach sets two threshold values *t* _{1} (say 10^{−10}) and *t* _{2} (say 0·005) to control the effect size and the *P*-values, respectively, and then proceeds to identify significant effects as follows: (1) searching for main effects: for each chromosome *c* (*c*=1, 2, …, 19), simultaneously add all possible main effects for markers on chromosome *c* into the current model, fit the model, and then delete main effects satisfying ; (2) searching for epistatic effects among the included main effects: simultaneously add all possible interactions among the included main effects into the current model, fit the model, and then delete epistatic effects satisfying ; (3) searching for epistatic effects between the excluded and the included main effects: for chromosome *c* (*c*=1, 2, …, 19), simultaneously add all possible interactions between the remaining main effects (not included in the current model) on chromosome *c* and the included main effects into the current model, fit the model, and then delete epistatic effects ; (4) removing main and epistatic effects with *P*-values larger than *t* _{2} from the current model.

After these steps, a final model would be obtained with both the preset environmental effects and the genetic effects that are associated with the phenotype. The estimates of these effects and the corresponding *P*-values indicate how strongly they influence the phenotype. We present estimate, standard error and *P*-value for each effect included in the final model, and calculated two measures of model fit and comparison, the deviance and the Akaike information criterion (AIC).

## 3. Results

### (i) LF

#### (a) Poisson models

Using the main-effect model, five main effects were detected (the left panel of Fig. 2). Among them, two additive effects were located on chromosomes 2 and 10, and three dominance effects on chromosomes 1, 2 and 10, respectively. The strongest main effect on *LF* was on chromosome 2, at position 65·3 cM, and it accounted for the highest proportion of phenotypic variance (about 3%). The *P*-values for the main effects were all less than 0·05.

The Poisson epistatic model detected five main effects and seven epistatic effects (the right panel of Fig. 2). The five main effects were identical to those of the main effect model. Among the seven epistatic effects, there were two additive–additive, three additive–dominance, one dominance–additive and one dominance–dominance effects on *LF* located on some pairs of the seven chromosomes (2, 4, 10, 11, 14, 17 and 19). The epistatic effects involved three QTLs identified in the main effect model and seven novel QTLs, suggesting that the epistasis model can identify QTLs with weak main effects but strong epistasis. The strongest epistatic effect on *LF* was observed between chromosome 10 and chromosome 14. The absolute values of estimated epistatic effects ranged from 0·13 to 0·16, which were about two times greater than those of the main effects, and the total proportion of phenotypic variances explained by the epistatic effects was approximately 9%, which was higher than that of the main effects. This revealed evidence that QTL affect *LF* mainly through their epistatic effects and the epistasis plays a more important role than does the main effect in controlling the genetic variation of *LF*. The *P*-values for the main and epistatic effects were all less than 0·05.

For comparison purposes, the Poisson main-effect model for *LF* was refitted without model search (i.e. always including all of 126 additive and dominance effects in the model). Figure 3 shows that the main effects detected in this saturated model coincided with those in the parsimonious model (the left panel of Fig. 2). The majority of genetic effects were shrunk to zero, indicating that the proposed model with shrinkage prior can capture the notion that most genetic effects influencing *LF* are very weak. We also tried to simultaneously fit all the main effects and epistatic interactions, but the attempt failed because of requirement of large memory. This demonstrated that a model search strategy is necessary to build a parsimonious model by seeking significant genetic effects when the number of effects is huge.

#### (b) Binomial models

The binomial main-effect model identified seven main effects (the left panel of Fig. 4). Among these main effects, four additive effects were located on chromosomes 2, 6 and 9, and three dominance effects on chromosomes 1, 2 and 10, respectively. The strongest main effect on *LF* was on chromosome 9, at position 0 cM, and it contributed the highest proportion (about 6%) to the phenotypic variance. The *P*-values for the main effects were all less than 0·0005.

Under the epistatic model, five main effects and 11 epistatic effects were identified (the left panel of Fig. 4). Among the five main effects, four effects overlapped those in the binomial main-effect model. The strongest main effect on *LF* was on chromosome 2, at position 65·3 cM, and it explained the highest proportion (about 5%) of phenotypic variance. Among the 11 epistatic effects, there were one additive–additive, two additive–dominance, five dominance–additive and three dominance–dominance effects. The epistatic effects involved three main effects that were identified and 16 effects that were not identified in the main-effect model. This shows that the model can identify QTLs which have a weak main effect but a strong interaction with other QTLs. The strongest epistatic effect was observed between chromosome 13 and chromosome 16. The absolute values of estimated epistatic effects ranged from 0·4 to 0·9, two times greater than those of main effects, and the total proportion of phenotypic variance explained by the epistatic effects was approximately 23%, about two times greater than that of main effects. All the *P*-values for the main and epistatic effects in the model were less than 0·001.

### (ii) DF

#### (a) Poisson models

The Poisson main-effect model detected nine main effects on *DF* (the left panel of Fig. 5). Among these main effects, six additive effects were on chromosomes 1, 3, 6, 8, 10 and 11, and three dominance effects on chromosomes 1, 5 and 14, respectively. The strongest main effect on *DF* was on chromosome 6, at position 0 cM, and it contributed the highest proportion (about 2%) to the overall phenotypic variance. The *P*-values for the main effects were all less than 0·05.

Using the Poisson epistatic model, seven main effects and nine epistatic effects were identified (the right panel of Fig. 5). Among the seven main effects, five were identical to those identified in the main-effect model. The strongest main effect on *DF* was the same as that in the main-effect model. Among the nine epistatic effects, there were three additive–additive, two additive–dominance, three dominance–additive and one dominance–dominance effects. The epistatic effects involved only one main effect that was identified and 15 effects that were not identified in the main-effect model. The strongest epistatic effect on *DF* was observed between chromosome 2 and chromosome 13, accounting for the highest proportion of phenotypic variance (about 3%). The absolute values of estimated epistatic effects ranged from 0·5 to 0·9, which were about three times greater than those of main effects, and the total proportion of phenotypic variances explained by the epistatic effects was approximately 15%, about two times greater than that of main effects. The *P*-values for the main and epistatic effects were less than 0·001.

#### (b) Binomial models

The binomial main-effect model identified nine main effects (the left panel of Fig. 6). Seven additive effects were located on chromosomes 1, 3, 6, 8, 9, 10 and 11, and two dominance effects on chromosomes 14 and 18, respectively. The strongest main effect on *DF* was on chromosome 6, at position 0 cM, and it contributed the highest proportion (about 2%) to the phenotypic variance. The *P*-values for the main effects were all less than 0·01.

The epistatic model detected seven main effects and ten epistatic effects (the right panel of Fig. 6). Among the seven main effects, five coincided with those identified in the binomial main-effect model. The strongest main effect on *DF* was located on chromosome 1, at position 26·1 cM, and it accounted for the highest proportion of phenotypic variance (about 2%). There were three additive–additive, three additive–dominance, three dominance–additive and one dominance–dominance effects on *DF* Among these epistatic effects, six were also identified by the preceding Poisson epistatic model for *DF*. The epistatic effects involved three main effects that were identified and 15 effects that were not identified in the main-effect model. The strongest epistatic effect on *DF* was between chromosome 2 and chromosome 6. The absolute values of estimated epistatic effects ranged from 0·6 to 1·0, which were about two times greater than those of main effects, and the total proportion of phenotypic variances explained by the epistatic effects was approximately 20%, two times greater than that of main effects. All the *P*-values for the main and epistatic effects in the model were less than 0·001.

### (iii) Model comparison

We used two summary measures, the deviance and the AIC, to compare different models. The deviance, defined as −2 times the log-likelihood, is a statistical summary of model fit; lower deviance means better fit to data. The AIC, defined as deviance +2 (number of predictors), measures the predictive power; a model is estimated to reduce out-of-sample prediction error if the AIC decreases. In all the analyses, the epistatic model had lower deviance and AIC than the corresponding main-effect model (Figs 2–6). This indicated that inclusion of the significant epistatic interactions improved the fit of the model to data and reduced out-of-sample prediction error.

### (iv) Treating discrete *LF* and *DF* as normally continuous traits

To investigate whether or not count phenotypes can be analysed by methods for continuous traits, we performed Bayesian methods by treating the count phenotypes *LF* and *DF* or the ratios *LF*/*OR* and *DF*/*OR* as normally continuous traits. We also tried to use the Box–Cox transformation to these traits and obtained similar results displayed here. Our analyses used the same priors and the model search method as in the analyses of our Poisson and Binomial models. Figure 7 displays the estimates, standard errors and *P*-values for the effects detected in the final models. Most of the effects detected by the normal models were also detected in the Poisson and the binomial models. However, the normal models missed several strong effects that were detected by both the Poisson and the Binomial models. The results indicated that the proposed methods were more appropriate for analysing the count phenotypes. This is expected from the general framework of generalized linear models (e.g. Gelman *et al.*, Reference Gelman, Carlin, Stern and Rubin2003).

## 4. Discussion

We have developed hierarchical Poisson and binomial models for mapping multiple interacting QTLs for count phenotypes based upon the unified generalized linear model framework of Yi & Banerjee (Reference Yi and Banerjee2009). Our method can fit a large number of effects, including covariates, main effects of numerous loci, and epistatic and gene–environment interactions, and can accommodate the correlation among the variables. Many complex traits, including reproductive phenotypes in mice, were measured as Poisson or binomial data. However, statistical methods for mapping multiple interacting QTLs for such phenotypes have not been fully developed previously.

To better characterize the genetic architecture of reproductive traits in mice, we applied our methods to two counted reproductive traits, *LF* and *DF*, in an F_{2} female mouse population. Since the observed number of OR provides important information about the variation of *LF* and *DF*, the numbers of *LF* and *DF* were modelled using the Poisson model with *OR* as an offset, and as the events of ‘success’ out of total ‘trials’ (*OR*) using the binomial model. As described earlier, the Poisson and binomial models capture different properties of the trait, and allow us to detect some different QTLs. In contrast, since the Poisson distribution can be derived as a limiting case to the binomial distribution as the number of trials goes to infinity and the expected number of successes remains fixed, the two models can lead to similar results (Casella & Berger, Reference Casella and Berger2001). In the present study, however, since the sample size (*n*=439) is not large enough and the proportions of *LF* and *DF* (0·75 and 0·10, respectively) are not small, it is not surprising that the QTLs detected in the Poisson models are not in full agreement with those detected in the binomial models. Furthermore, since the proportion of *DF* is lower than that of *LF*, there is more coincidence between two models for *DF* than that for *LF*.

More than 10 QTLs involved in the main and epistatic effects were identified for *LF* and *DF*, respectively, which exhibit a complex pattern of genetic influence on *LF* and *DF*. Most QTLs show a very weak main effect but a strong epistatic effect on *LF* and *DF*. Compared to the initial study (Rocha *et al.*, Reference Rocha, Eisen, Siewerdt, Van Vleck and Pomp2004*b*), in which only 3 and 1 QTLs were identified for *LF* and *DF*, respectively, and no epistasis was evaluated, the current study not only identified additional QTLs but also provided new information about the genetic architecture of *LF* and *DF* through the epistasis. These results also demonstrated that the models incorporating analysis of epistasis can identify QTLs, which might have a weak main effect but a strong epistatic effect with other QTLs. Moreover, based on the number of QTLs involved in the epistatic effects, the absolute values of estimated epistatic effects, and the total proportion of phenotypic variances accounted for by the epistatic effects, one conclusion can be drawn that the epistasis plays a more crucial role than does the main effect in regulating the genetic variation of *LF* and *DF*. Given the importance of epistasis in the genetic architecture of complex traits, appropriate statistical analyses should accommodate epistatic effects (Manolio *et al.*, Reference Manolio, Collins, Cox, Goldstein, Hindorff, Hunter, McCarthy, Ramos, Cardon, Chakravarti, Cho, Guttmacher, Kong, Kruglyak, Mardis, Rotimi, Slatkin, Valle, Whittemore, Boehnke, Clark, Eichler, Gibson, Haines, Mackay, McCarroll and Visscher2009).

Several main and epistatic effects are shared by *LF* and *DF*, which suggests that the pleiotropy plays an important role in the reproductive traits in this particular context of F_{2} mice. Among all of the chromosomes on which some QTLs were detected for *LF* and *DF* in the present study, the most active one is chromosome 2, on which about seven QTLs involved in the main and epistatic effects were detected. Furthermore, the QTL on chromosome 2 has the strongest main or epistatic effects and contribute the highest proportion to the overall phenotypic variance in some fitted models, which suggests that chromosome 2 has potentially biological relevance to *LF* and *DF*. The result is consistent with that in the initial study (Rocha *et al.*, Reference Rocha, Eisen, Siewerdt, Van Vleck and Pomp2004*b*). Other frequently involved chromosomes include chromosomes 1, 6, 9 and 10.

This work was supported by National Institutes of Health (NIH) Grants R01 GM069430 to N. Y., and P30DK056336. The opinions expressed herein are those of the authors and do not necessarily represent those of the NIH or any other organization with which the authors are affiliated.