Mixed model approach for IBD-based QTL mapping in a complex oil palm pedigree

Plant material

All palms analyzed belonged to families of the commercial oil palm (Elaeis guineensis Jacq.) breeding program of PalmElit, a CIRAD subsidiary and leading oil palm breeding company (www.palmelit.com). This breeding program is shared and conducted with partners, i.e. PT Socfin Indonesia (Indonesia) and INRAB (Benin). The 30852 palms with phenotypic evaluations belong to 478 crosses between 146 and 156 parents from heterotic groups A (GA) and B (GB), respectively. Most of the parents were planted in Pobè (INRAB, Benin), where the crosses were carried out. AxB crosses were planted in Aek Loba estate (PT Socfin Indonesia). Heterotic groups were based on the complementarity of yield components, bunch number and average weight, with GA palms having a low number of big bunches and GB palms having a high number of small bunches [26]. GA consisted of 132 palms from the “Deli” population derived from four oil palms planted in 1848 in Indonesia [27], while the remaining 11 palms originated from Angola. GB consisted of palms from La Mé (106, Côte d’Ivoire), Yangambi (23, Democratic Republic of the Congo), La Mé x Yangambi (5), La Mé x Sibiti (6, Democratic Republic of the Congo) and Nigeria (5).

Pedigree information was available for individuals in both heterotic groups (Additional file 1: Figure S1). However, for Deli individuals, the pedigree was known for a maximum of four generations. In this study, the Deli pedigree was reconstructed with MOLCOANC 3.0 software, as described in Cros et al. [28], to obtain kinship coefficients between all Deli palms.

Phenotypic data

Phenotypic evaluations were conducted in 26 trials on 30782 A × B oil palms that had been planted on 350 ha at Aek Loba (Indonesia, SOCFINDO estate) between 1995 and 2000 in a study lasting 11 years. The experimental designs of the trials involved randomized complete block designs (RCBD) with five or six blocks and balanced lattices of ranks four or five. The data used in this study was from oil palms which had reached a mature stage 6 to 10 years after planting. Mature bunch harvests were performed every 10 days, and the bunch numbers and weights were recorded. Based on these data, three production traits were analyzed, i.e. fresh fruit bunch yield (FFB, kg/year), bunch number (BN, bunch/year) and average bunch weight (ABW, kg), which is the ratio between FFB and BN estimated annually.

Phenotypic correlations (r _p ) were estimated in the A × B population via the Pearson correlation coefficient.

Molecular data and genetic map construction

The genotyped population consisted of palms from both heterotic groups which were parents of the crosses evaluated in the genetic trials. Palms were genotyped with 388 SSR markers developed in different studies [29–31].

Total genomic DNA was extracted from freeze-dried leaf samples of each progeny and parent using the commercial kit NucleoSpin 96 plants II (Macherey-Nagel, Germany) in accordance with the manufacturer’s protocol. Genomic DNA concentrations were estimated with a spectrometer (Infinite® 200 PRO NanoQuant®, Switzerland). Microsatellite fragment amplification was performed in a 384-well plate with 25 ng of DNA in a 10 μl final volume containing 1 U of Taq DNA polymerase, 1 μl of 10x buffer (10 mM Tris (pH 8.3), 50 mM KCl, 1.5 mM MgCl2, 0.001 % (w/v) glycerol), 2 mM dNTP, 0.6 mM MgCl₂, 0.08 μM of M13-tailed primer, 1 μM of the reverse primer, 1 μM of M13 primer-fluorescent dye VIC, PET, NED or FAM (Applied Biosystems, USA). The PCR conditions were as follows: initial denaturation at 94 °C for 5 min, followed by 35 cycles alternating 30 s of denaturation, 1 min 15 s of hybridation at annealing temperature, 1 min 30 s of extension, and ending with 30 min of final elongation at 72 °C. PCR products were pooled. We used 2 μl of the pooled PCR product with 10 μl of size standard GeneScan ™-600 LIZ ® at 0.0012 % in Hi-Di™ formamide for migration on a 3500xL Genetic Analyzer (Applied Biosystems, USA).

GeneMapper© V4.1 (Applied Biosystems, USA) software was used to determine allele sizes.

A consensus linkage map (Additional file 1: Figure S2) based on the pedigree of both groups A and B was constructed with CRIMAP version 2.4 [32], as described in [33]. The genetic map was drawn using MapChart software [34].

QTL mapping

QTL mapping was performed following the two-step variance component approach described in George et al. [19]. Van Eeuwijk et al. [21] presented an adapted approach to identify QTLs in a maize hybrid program that tested parents from heterotic groups like in oil palm. The approach presented here was extended to an outbreeding population for which genotyping data were available only on the parents of both heterotic groups A and B and phenotypic data only on A × B individuals. The first step involved computation of IBD kinship matrices over a grid of genetic positions. The second step consisted of QTL presence tests over the grid of genetic positions by comparison of various linear mixed models.

Variation for production traits in AxB population

Fresh fruit bunch yield (FFB) and underlying components bunch number (BN) and average bunch weight (ABW) were normally distributed within the hybrid population (30852 palms) derived from crosses between palms of heterotic groups A and B, and ranged from 6 to 411 kg/palm/year, 2 to 34 bunches/year and 1.2 to 31.4 kg, respectively (Fig. 1). Coefficients of variation were around 16 % for FFB and ABW and 22 % for BN. In the A × B population, BN and ABW were highly and negatively correlated and FFB was highly and positively correlated with BN, while it was negatively correlated with ABW (Fig. 1c). Genetic correlations showed the same pattern, with closer correlations found in the heterotic group A than B (Fig. 1c). Consequently, the isoproduction curves on Fig 1b show that the higher FFB values were mainly reached with high BN rather than high ABW.

Heritability based on estimation of variance components with only pedigree information ranged from 0.15 for FFB in heterotic group A to 0.36 for ABW in heterotic group B (Fig. 1d). Heritability was lower in heterotic group A than in heterotic group B for the three production traits, with a lower additive genetic variance found in group A (Fig. 1d).

QTL mapping for production traits

Figure 2 shows a heat map of the log-likelihood ratio test (LRT) for the second genome scan of QTL mapping process (see Material and Methods) in each heterotic group for the three production traits. Eighteen significant QTL regions were identified, i.e. seven in group A and eleven in group B. The QTL regions were distributed in 14 linkage groups, with linkage group 1 and 8 harboring two and four QTL regions, respectively (Fig. 2). QTL confidence intervals spanned 14.5 cM on average, ranging from 2 to 42 cM (Fig. 2).

A major QTL was found for group B on the bottom of linkage group 4, accounting for 21 and 44 % of the phenotypic variance of BN and ABW, respectively (Fig. 2). However, this region was not significant for FFB due to the opposite effects on BN and ABW. The strongest QTL found in group A was located on the top of linkage group 12 and it accounted for 7, 12 and 11 % of the phenotypic variance of FFB, BN and ABW, respectively (Fig. 2).

Pleiotropic effects of QTLs

Correlations obtained between production traits at the hybrid population level indicated that common genetic regions should control several traits. To assess the prevalence of these pleiotropic QTLs, correlations between additive values of each parent at a given QTL were calculated for the 18 QTL regions identified for the three production traits (see Material and methods). Figure 3a shows the correlation patterns among pairs of traits, heterotic groups and QTLs. In heterotic group A, all correlations were high and positive between FFB and BN, and negative between ABW and the other two traits (Fig. 3a). The same pattern was observed in eight of the eleven QTLs identified in heterotic group B (Fig. 3a). Two of the three remaining QTLs, i.e. B_8@24 and B_16@57, had opposite effects on BN and ABW, but the effects on ABW were positively correlated with the effects on FFB (Fig. 3a). B_5@69 was the only QTL that had positive effects on the three production traits (Fig. 3a), but no significant effect was found on BN and ABW (Table 1).

Figure 3b illustrates the different types of QTL found in heterotic group B. QTL B_1@159 had opposite and significant effects on BN and ABW, like the majority of QTLs found (Fig. 3b), which could explain the correlations obtained in the hybrid population. However, QTLs that had trait-specific effects were found, on BN exclusively (B_15@117, Fig. 3b) or with superior effects on ABW, while the effect on BN was non-significant (B_16@57, Fig. 3b).

Overlapping of QTLs between heterotic groups

To assess the complementarity of QTL patterns as a possible cause of the heterosis observed in the hybrid population, the findings of a log-likelihood ratio test (LRT) for the presence of additive QTLs in each heterotic group at all genomic evaluation points were plotted for group A versus group B (Fig. 4). For FFB, point locations mainly along the X and Y axes indicated that the QTL patterns differed between the two heterotic groups (Fig. 4a). On the contrary, for BN and ABW, points near the diagonal indicated overlapping of QTLs for these traits between the two heterotic groups (Fig. 4b and c). Among the 18 QTLs in the final model, these overlapping QTLs were found only on the bottom on linkage group 1 (1@165 for group A and 1@159 for group B, Fig. 2), but some others reached the significance threshold in only one of the two groups (linkage group 12, Fig. 2).

Variance component approach for QTL mapping

QTL mapping is interesting when assessing perennial species [4], but biological and economic issues can hamper its use even for major crop species such as oil palm for which genetic studies are scarce. As proposed by Parisseaux and Bernardo [17], in silico mapping, i.e. based on data from current breeding programs, has many combined advantages. No extra costs are required to develop and phenotype specific mapping populations and this approach enables screening a high number of individuals and alleles, increasing relevance for the breeding process. Linear mixed model approaches used to map QTLs in such genetic material are flexible and powerful [42]. In the present study, this allowed us to test for the presence of QTLs in two heterotic groups, as proposed in maize by van Eeuwijk et al. [21], and to account for confounding nongenetic effects and polygenic effects, which can be important in such large trials with different genetic backgrounds. Moreover, one of the interesting benefits specific to this study was the possibility of linking phenotypes that were recorded exclusively in A × B individuals to genotypes available only for parents and grand-parents in some cases. Thus it enabled testing the presence of QTLs in both heterotic groups in a single model without first estimating breeding values of genotyped individuals. By genotyping only parents and grandparents, recombination events at the AxB individuals were not taken into account. This could limit the accuracy of the position of QTL and their confidence intervals which are mainly determined by the number of recombination events analyzed. However, the number of individuals genotyped, more than 300, is greater than the population sizes commonly found in oil palm and allowed, with moderate genotyping cost, to conduct QTL mapping based on a very high amount of phenotypic data, thus enhancing the detection power. In future studies, the flexibility of a mixed model framework should allow us to easily extend this approach to test environmental responses and dominance or epistasis effects.

In a two-step variance component approach, QTL mapping is influenced by the precision of IBD kinship matrices that highlight allele identities between individuals and recombination events between genetic positions. Marker and pedigree information were combined in multipoint analyses using Simwalk2 Bayesian software to estimate these matrices [35], which allowed us to handle such large and complex pedigrees. The pedigree information quality is a key factor in this step. Unlike the group B pedigree, which is well known for the different genetic origins that compose it, the group A pedigree is not fully known in higher generations. MOLCOANC software [28] was used to reconstruct the pedigree, and QTL mapping using the reconstructed pedigree identified two additional QTLs compared to the raw pedigree (Additional file 1: Figure S4). One of them on linkage group 12 is highly supported by a QTL found in the same location for the same traits in a multiparental population involving related palm trees as parents [6].

Comparison to previous results in oil palm

Few other studies have aimed at mapping QTLs for the same production traits, despite frequent limitations for collecting data for a number of individuals that would be optimal for QTL mapping. Rance et al. [5] first attempted to map QTLs based on 84 individuals from an F2 population derived from a single palm selfing progeny. A saturated genetic map could not be obtained, and three putative QTLs were identified, one for each of three production traits, i.e. FFB, BN and ABW. More recently two studies investigated production traits in F1 populations based on 208 individuals [10] and 52 individuals [7]. The latest involved very few individuals and identified two QTLs for FFB and one for ABW, but none for BN. Ukoskit et al. [10], using a greater number of individuals, identified three QTLs for FFB. Three, two and five QTLs were found for FFB, BN and ABW, respectively, in a multiparental population pooling 299 individuals from an incomplete factorial design involving two parents of each heterotic group [6]. The two-step variance component approach used in this study identified eight, ten and nine QTLs for FFB, BN and ABW, respectively, and some other suggestive QTLs that were just below LRT threshold. The higher detection power achieved in our study compared to previous works was likely due to the higher genetic variation screened through pedigree based analysis. Analysis in F1 from single crosses, especially in oil palm for which our commercial elite genotypes are partially inbreed, screen a very restricted genetic basis and raises the risk that palms used as parents had homozygous positive alleles at many QTLs. Although the multi-parent design overcame this issue, it was less efficient than pedigree-based analysis. Moreover, the difference in detection power could be linked to the higher number of individuals in the mapping population, at least compared to studies that have fewer. The number of genotypes upon which the QTL mapping was based was actually not very high (around 150 in each heterotic group) as only the parents were genotyped, but the general combining abilities of these individuals were estimated very precisely via numerous progeny tests, with an accuracy of 0.9 [43].

CIRAD microsatellite marker resources have been widely used in oil palm genetics since the publication of a reference genetic map [29] that facilitated comparison of QTL mapping results between different studies. Among the 18 QTL regions identified in this study, eight were found by other authors for the same traits. Billotte et al. [6] identified the set of QTLs closest to our, with five shared regions as well as three with suggestive QTLs. This was not surprising as the parents derived from their factorial design were founders in the pedigrees of the group A and B individuals used in our study. The co-localization of pedigree-based QTLs with those found in previous studies in related or unrelated populations confirms the ability of the two-step variance component approach to identify regions that control production traits.

Determinism of heterotic effects on fresh fruit bunch yield

The implementation of the large genetic trial used in this study follows the RRS with progeny testing of parents, motivated by the observations of increased FFB in A × B palms compared to the maximal value of intra-group crosses [3]. This heterosis effect might be explained by the better combination of underlying production traits, i.e. bunch number and average weight, that exhibited marked differences between heterotic groups, with group A palms having a low number of big bunches and group B palms having a high number of small bunches. Heterosis is an important phenomenon in plant breeding, especially in crops that exploit it through commercial hybrid production, but the details of its molecular determinism are still not well known. Dominance through complementation of inferior alleles and over dominance are the main mechanisms proposed without being mutually exclusive [44]. In maize, a QTL mapping approach on large pedigree of heterotic groups showed that complementation of positive alleles in the hybrids could be an important factor [21]. The low occurrence of collocations of QTLs between heterotic groups noted in the present study, only one among 18 significant QTL regions, confirms that both heterotic groups are genetically distant and suggests that complementation of favorable alleles would also be part of the phenomenon in oil palm. However, the experimental design used in this study was not adequate to identify QTL fixed for alternative alleles in heterotic groups and could have a major role on heterosis, because QTL effects were tested within heterotic groups and were identified only if the QTL was polymorphic within the group.

Inclusion of QTL information in oil palm breeding programs

Our QTL mapping approach based on data from current breeding programs overcame usual ineffective analysis problems by pooling a sufficient amount of data, and the identified QTLs are highly relevant for integrating information for breeding. Moreover, the incomplete factorial design allowed us to obtain estimates of QTL allele effects from a wide range of genetic backgrounds, thus enhancing the transferability to breeding programs. However, as mentioned by Würschum [18], a validation step is necessary in marker-assisted selection to avoid the use of false-positive or overestimated QTLs, a problem that can arise due to biases in the statistical design. QTLs can be validated in independent oil palm pedigrees in other genetic trials to assess the stability of QTLs in different genetic or environmental settings, as has been done in other perennial species [45, 46]. An alternative approach would be to use the A × B population, described in this study, to validate QTLs at the progeny level by genotyping a high number of A × B individuals with markers identified at QTL locations. In addition, genotyping directly A × B individuals in several related progenies would offer two advantages: to give the statistical power for testing QTL effects conditionally of genetic background, and so detect dominance or epistasis effects between heterotic group A and B alleles, and to increase the number of recombination events in the pedigree analyzed in regions of interest, which combined with an increased genotyping density of these regions would reduce the confidence interval around QTLs. The latter would be interesting for QTLs identified in less well covered regions (e.g. top of linkage group 14) caused by difficulties in genetic map construction [33].

For the QTL validation or marker-assisted selection process, target QTL genotypes of pedigree members must be known in order to be able to design efficient crosses that test and combine segregating QTLs. The variance component approach modeled QTL effects as a random term, and estimates of the additive value of each individual at a given QTL were obtained with BLUPs. Despite the normal distribution of BLUPs for QTL effects, the segregation in full-sib families and alleles carried by pedigree members could be inferred qualitatively, e.g. by looking at the pedigree with projected BLUP values (Additional file 1: Figure S5). The Bayesian approach implemented in FlexQTL software offers an interesting way to predict the QTL genotype [23]. The presence of biallelic QTLs in the genome are tested simultaneously with QTL alleles carried by founders by inference on QTL genotype probabilities for any members of the pedigree based on pedigree, marker and phenotype data.

Such robust information about QTLs and their segregation in genetic backgrounds used in a breeding scheme could be integrated in the RRS to favor complementation of QTLs for production traits and thus to enhance the heterotic effects. As shown in Additional file 1: Figure S5, QTL segregation in heterotic groups differed depending on the genetic origins, e.g. La Mé and Yangambi, and recombination between these genetic origins could lead to improvements within heterotic groups by combining favorable alleles. Conventional pyramiding of interesting QTL approaches could be a way to obtain in few generations improved parents in each heterotic group. Although pyramiding could be efficient in the case of traits with simple genetic determinism, genomic selection is becoming an attractive approach for traits with genetic architecture following an infinitesimal model, especially in species with long breeding cycles and biological constraints. Various authors have proposed the inclusion of genetic architecture information to prioritize genomic regions in the estimation of breeding value [47–49]. The benefits of a genomic selection approach was evaluated using data from the same oil palm breeding program [43] and our results provide an opportunity to test, based on empirical data, a combined QTL and genomic selection approach to achieve efficient marker-assisted selection.