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.
A linear mixed model was designed to account for non-genetic (trial, block) and genetic (general combining ability (GCA) and specific combining ability (SCA)) effects and adjusted to the data:
$$ Y=X\beta +{Z}_1u+{Z}_2s+{Z}_A{g}_A+{Z}_B{g}_B+e $$
(1)
where Y is a n × 1 observation vector of production traits (n = 30852) averaged over 4 years in the A × B population, X is a n × m design matrix relating observations to trial fixed effects β with β being a m × 1 vector (m = 26), Z
1 is a n × b design matrix relating observations to block random effects u ~ N(0, Iσ
2
u
) with u being a b × 1 vector (b = 677), Z
2 is a n × c design matrix relating observations to SCA random effects s ~ N(0, Iσ
2
s
) with s being a c × 1 vector (c = 478), Z
A
and Z
B
are n × q
A
and n × q
B
design matrices relating observations to GCA random effects for heterotic groups A and B, g
A
~ N(0, A
A
σ
2
a
A
) and g
B
~ N(0, A
B
σ
2
a
B
) respectively, with g
A
and g
B
being q
A
× 1 and q
B
× 1 vectors, respectively (q
A
= 146 and q
B
= 156), and eis the n × 1 vector of residual effects ~ N(0, Iσ
2
e
). I is an identity matrix and A
A
and A
B
are the pedigree-based kinship matrices of heterotic groups A and B, respectively.
Calculation of heritabilities and genetic correlations were based on variances estimated by the linear mixed model (1). Narrow-sense heritabilities (h2) of each trait were obtained as the ratio of the variance of general combining abilities, for both groups A and B, to the total phenotypic variance of crosses according to the formula:
$$ {h}_{A/B}^2=\frac{\sigma^2{a}_{A/B}}{\sigma_b^2+{\sigma}_s^2+{\sigma}^2{a}_A+{\sigma}^2{a}_B+{\sigma}_e^2} $$
Genetic correlations (r
g
) between traits T1 and T2 in each heterotic group were calculated using covariances and variances estimated in a multivariate linear mixed model, based on model (1). The model had the form:
$$ \left[{}_{T_2}^{T_1}\right]=\left[{{}_0^X}_X^0\right]\kern0.5em \left[{}_{\beta_{T2}}^{\beta_{T1}}\right]+\left[{{}_0^{Z_1}}_{Z_1}^0\right]\kern0.5em \left[{}_{u_{T2}}^{u_{T1}}\right]+\left[{{}_0^{Z_2}}_{Z_2}^0\right]\left[{}_{S_{T2}}^{S_{T1}}\right]+\left[{{}_0^{Z_A}}_{Z_A}^0\right]\left[{}_{{g_A}_{T_2}}^{g_{A{T}_1}}\right]+\left[{{}_0^{Z_B}}_{Z_B}^0\right]\left[{}_{g_{B_{T2}}}^{g_{{B_T}_1}}\right]+\left[{}_{e_{T_2}}^{e_{T_1}}\right] $$
(2)
GCA were structured as:
$$ \left[\begin{array}{c}\hfill {g_A}_{T_1}\hfill \\ {}\hfill {g_A}_{T_2}\hfill \end{array}\right]\sim N\left(0,\left[\begin{array}{cc}\hfill {\sigma}^2{a}_{A_{T_1}}\hfill & \hfill Co{v}_A\left({T}_1,{T}_2\right)\hfill \\ {}\hfill Co{v}_A\left({T}_1,{T}_2\right)\hfill & \hfill {\sigma}^2{a}_{A_{T_2}}\hfill \end{array}\right]\otimes {A}_A\right) $$
$$ \left[\begin{array}{c}\hfill {g_B}_{T_1}\hfill \\ {}\hfill {g_B}_{T_2}\hfill \end{array}\right]\sim N\left(0,\left[\begin{array}{cc}\hfill {\sigma}^2{a}_{B_{T_1}}\hfill & \hfill Co{v}_B\left({T}_1,{T}_2\right)\hfill \\ {}\hfill Co{v}_B\left({T}_1,{T}_2\right)\hfill & \hfill {\sigma}^2{a}_{B_{T_2}}\hfill \end{array}\right]\otimes {A}_B\right) $$
where Cov
A
(T
1, T
2) and Cov
B
(T
1, T
2) are additive genetic covariances.
Residual effects were structured as:
$$ \left[\begin{array}{c}\hfill {e}_{T_1}\hfill \\ {}\hfill {e}_{T_2}\hfill \end{array}\right]\sim N\left(0,\left[\begin{array}{cc}\hfill {\sigma}^2{{}_e}_{T_1}\hfill & \hfill Co{v}_e\left({T}_1,{T}_2\right)\hfill \\ {}\hfill Co{v}_e\left({T}_1,{T}_2\right)\hfill & \hfill {\sigma}^2{{}_e}_{T_2}\hfill \end{array}\right]\otimes I\right) $$
where Cov
e
(T
1, T
2) is residual covariance.
Genetic correlations (r
g
) were calculated according to the formula:
$$ {r}_{gA/gB}=\frac{Co{v}_{A/B}\left({T}_1,{T}_2\right)}{\sigma_{A/B\ }{T}_1\times {\sigma}_{A/B\ }{T}_2} $$
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 MgCl2, 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.
IBD matrix computation
QTL mapping using a variance component approach requires identity-by-descent (IBD) relationship matrices as variance-covariance matrices for testing QTL effects at each genetic position. An IBD analysis using SimWalk2 software [35] was conducted to estimate the IBD matrices formed by empirical kinship coefficients for each pair of palm trees present in group A and B pedigrees. Simwalk2 implements a Bayesian framework using a Markov chain Monte Carlo (MCMC) algorithm to select the most likely genetic descent graph depicting inheritance patterns within pedigrees [36]. Empirical kinship coefficients were calculated at every marker position and over a 3 cM interval grid, leading to 922 evaluation points for QTL presence testing using IBD analysis in Simwalk2.
Statistical modeling
The model (1) (Phenotypic data section) was used as a null model for QTL mapping. Two other models were designed to test for the QTL presence in each heterotic group by adding a random QTL effect to model (1). A genome scan was performed by fitting models for groups A and B at k genetic positions used to compute IBD kinship matrices testing for the QTL presence. The models had the form:
$$ Y=X\beta +{Z}_1u+{Z}_2s+{Z}_A{g}_A+{Z}_B{g}_B+{Z}_A{v}_{A,k}+e $$
(3)
$$ Y=X\beta +{Z}_1u+{Z}_2s+{Z}_A{g}_A+{Z}_B{g}_B+{Z}_B{v}_{B,k}+e $$
(4)
where v
A,k
and v
B,k
are q
A
× 1 and q
B
× 1 vectors of QTL random effects at the kth genetic position in heterotic group A and B respectively (q
A
= 146 and q
B
= 156), with v
A,k
~ N(0, A
A,k
σ
2
a
A,k
) and v
B,k
~ N(0, A
B,k
σ
2
a
B,k
). A
A,k
and A
B,k
are the IBD kinship matrices at the kth genetic position, calculated in groups A and B, respectively.
The log-likelihood ratio test (LRT) was performed to compare models (3)/(4) and (1) at each k position. As the distribution of LRT is unknown under the null hypothesis [19], permutations were performed to estimate significance thresholds (see next paragraph). A QTL was declared significant if the LRT exceeded the threshold, and sets of putative QTLs in both heterotic groups were constituted.
According to a multiple QTL mapping strategy [37] to better account for the QTL variance of QTL elsewhere in the genome when testing a genetic position, a second genome scan was then performed with updated null (5) and test (6)/(7) models that incorporated putative QTL random effects for both heterotic groups:
$$ Y=X\beta +{Z}_1u+{Z}_2s+{Z}_A{g}_A+{Z}_B{g}_B + {\displaystyle \sum_{c\ in\ {C}_A\ }}{Z}_A{v}_{A,c}+{\displaystyle \sum_{c\ in\ {C}_B\ }}{Z}_B{v}_{B,c} + e $$
(5)
$$ Y=X\beta +{Z}_1u+{Z}_2s+{Z}_A{g}_A+{Z}_B{g}_B+{\displaystyle \sum_{c\ in\ {C}_A\ }}{Z}_A{v}_{A,c}+{\displaystyle \sum_{c\ in\ {C}_B\ }}{Z}_B{v}_{B,c}+{Z}_A{v}_{A,k}+e $$
(6)
$$ Y=X\beta +{Z}_1u+{Z}_2s+{Z}_A{g}_A+{Z}_B{g}_B+{\displaystyle \sum_{c\ in\ {C}_A\ }}{Z}_A{v}_{A,c}+{\displaystyle \sum_{c\ in\ {C}_B\ }}{Z}_B{v}_{B,c}+{Z}_B{v}_{B,k}+e $$
(7)
Where C
A
and C
B
are subsets of genetic positions where putative QTLs were identified in the first genome scan for groups A and B, respectively. To test the QTL presence at genetic positions within a 20 cM window centered on C
A
and C
B
positions, the random putative QTL effect at the considered positions was removed in the null and test models (5) and (6)/(7), respectively.
As for the first genome scan model, (6)/(7) and (5) were compared using LRT and QTLs were declared significant if the LRT exceeded the threshold determined by permutation. Approximate confidence intervals of QTL location were estimated by 1-LRT support interval calculation, i.e. the interval in which LRT is within 1 units of its maximum [38]. Final multi-QTL models were fitted for each production trait incorporating random effects for all QTLs identified in the second genome scan. Best linear unbiased predictors (BLUP) were obtained as predictors of additive value of each individual at a given QTL. BLUPs were plotted on the pedigree graph grouping all tested parents using Pedimap software [39].
All linear mixed effect models were estimated by REML using the R-ASReml package (Butler et al., 2009) for R [40].
LRT threshold calculation
LRT thresholds used to determine the presence or absence of QTLs at the tested genetic positions were set independently for each heterotic group, each production trait and for the first and second genome scans. The threshold calculation was based on the distribution of the LRT under the hypothesis of the absence of QTLs in the genome. Null distributions were obtained by permutations [41], conducting genome scans with de-correlated QTL effects from phenotypes without affecting other effects correlations. To do that, vectors of IBD coefficients constituting the IBD kinship matrices associated to QTL effects, A
A,k
and A
B,k
in the scan model (see statistical modeling section), were randomly shuffled before each replicate of genome scan, whereas all other components were unchanged comparing to the null model. 1000 and 500 replicates were obtained to determine each threshold for the first and second genome scans of QTL mapping process respectively, for each trait and heterotic group (Additional file 1: Figure S3). For each replicate, the maximum LRT value obtained was recorded to construct the distribution of maximum LRT value under the null hypothesis. The 95 % percentile was calculated based on this distribution to determine a 5 % LRT threshold. To save computing time, a reduced density of evaluation points was tested for each replicate genome scan, with 1 out of 25 and 50 positions for the first and second genome scan of QTL mapping process respectively, which led to 40,000 and 10,500 tests of null hypothesis, respectively, to obtain the null distribution. The testing frame was shifted for each permutation sample to ensure that all genetic positions were assessed.
Assessment of QTL pleiotropic effects
18 QTL regions were identified based on LRT thresholds and the approximate confidence intervals of QTL location. Linear mixed models of the form (5) including effects for the 18 QTL regions were adjusted for each production variable. Additive QTL effects were predicted for parents of corresponding heterotic groups and Pearson correlation coefficients between them were calculated for each pair of production variables.
