Genotyping by sequencing for the construction of oil palm (Elaeis guineensis Jacq.) genetic linkage map and mapping of yield related quantitative trait loci

Plant material and DNA extraction

Leaf tissue was collected from 112 palms and two parent palms. These 112 tenera full sibs were generated from a cross of Deli dura and a Serdang pisifera. The palms were managed and maintained by the Department of Agriculture (DOA) at a plantation in Quoin Hill, Tawau, Sabah, Malaysia (4°24′N, 118°1′E). The leaf samples were collected from individual trees and immediately submerged in liquid nitrogen (LN) before being stored at −80 °C freezer until DNA extraction was executed. One hundred milligrams of each individual leaf sample were used as a starting material for DNA extraction. The samples were frozen in the LN before being disrupted with beads using the TissueLyser II. The genomic DNA from the tissue samples was extracted using the DNeasy^® Plant Mini Kit (Qiagen, Hilden, Germany), following the manufacturer’s protocol, but with minor modifications to prolong the incubation time.

The genomic DNA concentration and purity were quantified and qualified using the Nanodrop-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The samples were normalized in 100 ng/µl using Buffer AE (Qiagen, Hilden, Germany). Then, the integrity of the extracted DNA was further examined on a 1% TAE agarose gel (Choice Care, Malaysia) for an hour at 100 V. The DNA ladder (Thermo Scientific, Waltham, MA, USA) was added to serve as an indicator for reference bands. After an hour, the agarose gel was stained in an ethidium bromide solution for half an hour before being viewed under an AlphaImager HP DE-500 (Alpha Innotech, San Leandro, CA, USA). DNA samples of all 112 F₁ and the two parents that passed the quality control checks with acceptable ratios of absorbance reading at A260/A280 and A260/A230 1.8 and 2.0–2.2 (Thermo Scientific, Waltham, MA, USA), respectively, were considered to have acceptable DNA purity. These DNA samples were stored in a 4 °C chiller (Thermo Scientific, Waltham, MA, USA) prior to the sequencing works.

GBS library construction

The sequencing works as well as the GBS library construction on the 112 progenies and their parents was outsourced to Genting Laboratory Service Sdn. Bhd. (GLS). In this study, the quality of genomic DNA was quantified using the Qubit^® dsDNA HS Assay Kit (Life Technologies, Carlsbad, CA, USA). GBS libraries were constructed using the double digestion methods with HindIII and TaqI enzymes (New England Biolabs, Ipswich, MA, USA) according to the standard GBS protocol (Elshire et al., 2011). HindIII is a type II endonuclease that recognizes a degenerate 6 bp sequence (5′…AAGCTT…3′) and cleavages between the AA of the sequence resulting in 5′ overhang. TaqI. on the other hand, recognizes a degenerate 4 bp sequence (5′…TCGA…3′) and cuts between the TC of the sequence to create a 5′ overhang. A total of 200 ng of genomic DNA from each sample was digested with the enzymes for 2 h at 75 °C. The digested samples were examined through gel electrophoresis (Scie-Plas, Gloucestershire, England) on a 1% TAE agarose gel before the ligation adapters were added. The ligation products were pooled, and those with sizes of 400–600 bp were selected using Pipin Prep (Sage Science, Beverly, MA, USA). Two different types of adapters, barcode and common adapters, were used to enable paired-end and multiplex sequencing on the Illumina HiSeq2000 platform (Illumina Inc., San Diego, CA, USA).

In this experiment, 12 unique adapters with specific barcodes were designed to carry out the multiplexing of 12 samples per Illumina flow cell lane. The library construction was performed using a manufacturer’s protocol. Finally, libraries without adapter dimers were retained for DNA sequencing using the Illumina HiSeq2000 sequencing platform following the standard protocol (Illumina Inc., San Diego, CA, USA).

SNP and InDel calling

After the sequencing run, the sequencing reads were filtered using a FASTQ file with a quality confidence threshold of ≥25. The sequencing reads with the minimum read coverage of five for both parents and progenies were mapped and aligned to the EG5 genome build (Singh et al., 2013), which was obtained from the GenomeSawit webpage (http://genomsawit.mpob.gov.my), using the Burrows-Wheeler Aligner (Li & Durbin, 2009). Then, SNP and InDel callings were performed using the genome analysis toolkit, GATK ver. 3.3.0, with default parameters (McKenna et al., 2010).

SNP and InDel filtering

Using sequence data with a minimum depth reading of 5X and a Q-value of 25 for each sample, markers across the progeny palms with less than 10% missing calls and major alleles, as well as a genotype frequency of more than 95% were excluded. The filtering task was performed using VCFtools ver. 0.1.10 (Danecek et al., 2011) before the final versions of the SNP and InDel genotyping calls were tabulated into PLINK format (Purcell et al., 2007). In the preparation of the input file for the construction of a genetic linkage map, the allele calls for the parental palms were also investigated. Markers with genotypes that were different from the expected Mendelian segregation patterns between the progenies and the parents for F₁ segregation were filtered out.

Genotype data preparation

Markers with major allele and genotype frequencies of more than 95% were eliminated using Microsoft Excel 2010. Subsequently, the heterozygous markers were coded based on JoinMap ver. 4.1 for a cross pollinated (CP) population. Markers coded with <lmxll> and <nnxnp> represent markers with one heterozygous parents, and <hkxhk>, <efxeg> and <abxcd> represent markers with both parents being heterozygous with the presence of two alleles, three alleles and four alleles, respectively.

The segregation ratios for each marker were examined using the chi-square goodness-of-fit test, and the threshold p-value was set at 0.001. Segregation distortion markers with a p-value of less than 0.001 were excluded. Subsequently, the marker similarities were checked for reducing marker redundancy, as two identical markers would be located at the same position in the linkage map. If a pair of markers has a similarity value that exactly equals 1, one of the markers in the pair would be excluded. After removing markers with the filtering parameters as described above, the data set was ready for the construction of genetic linkage maps.

Genetic linkage map construction

Genetic linkage maps of the female parent, D200, and the male parent, QP447, were constructed using the significant SNP and InDel markers of the 112 F₁ population. The analysis was performed using the JoinMap 4.1 software (Van Ooijen, 2006) and the maps were drawn using the MapChart V2.2 software (Voorrips, 2002), as described below. Further data clean-up was performed prior to the construction of the genetic linkage map, as described in the following details. In order to capture the population structure in this study, 2,291 informative SNP markers were used to plot a principal component analysis (PCA) using PLINK ver. 1.90 (Weeks, 2010) and visualization using Microsoft Excel 2010 to access the F₁ mating design.

Maternal and paternal maps were constructed using the filtered markers sets. To determine the significant LOD threshold, the start value was set at 2.0 and the end value at 10.0 with a step size of 1.0 (Van Ooijen, 2006). The LOD score was calculated for the recombination frequency and was performed using the software based on the G2 statistic, G2 = 2 ∑ O log (O/E) where, O is the observed value and E is the expected number of individuals in a cell, log is the natural algorithm and ∑ is the sum over all cells. Then, the values were multiplied by 0.217, which is derived from 0.5*log10(e) to obtain the normal LOD score (Van Ooijen, 2006). The grouping parameter was selected at a LOD-value of more than 4.0 which indicates the likelihood of two loci being linked is greater than 10,000:1. The maximum likelihood (ML) mapping algorithm was selected to compute the mutual distance of loci in the respective groups.

The map order was further improved by eliminating markers exhibiting a Nearest Neighbor Stress value (N.N. stress) of more than 3 cM. The N.N. stress measurement assesses the extent to which markers are located outside of their expected region (Van Ooijen, 2006). After removing markers with a N.N. stress of more than 3 cM in each of the linkage groups, the maps were recalculated until the optimal map order was achieved. Finally, the genetic map distances were calculated using the Haldane mapping function. Once the linkage group maps were effectively constructed, both parental maps were drawn and visualized using MapChart 2.2 (Voorrips, 2002), and the genetic map information was used in the QTL analysis.

Phenotypic traits and analyses

In this study, the FFB yield of the individual palm from five consecutive years (1984–1988) was recorded, and the average FFB yield was calculated. With a minimum of at least three bunches per palm, O/B, O/DM, O/WM, M/F, K/F, S/F and F/B were analyzed as per standard industry practice (Blaak, Sparnaaij & Menedez, 1963). The OY was calculated from the multiplication of the 5-year average FFB yield and the average O/B. A total of nine traits (Supplemental Data) was compiled for descriptive statistics. Out of a total of 112 genotyped palms, only 111 palms were collected for FFB yield and 97 palms with yield related components. The phenotype data collection was performed by the Department of Agriculture (DOA), Sabah. During the field data recording, a few palms died, but their leaves were already sampled earlier for DNA sequencing and molecular analysis. As a result, fourteen samples were excluded from the bunch analysis due to the insufficient number of bunches recorded. A Pearson correlation analysis and the significant level between traits were also computed using SPSS statistical software. The R-package heritability was used to estimate narrow-sense heritabilities for the nine traits (Kruijer, 2019).

QTL analysis

To declare the presence of a significant QTL, the threshold LOD values were estimated at the genome-wide (GW) and chromosome-wide (CW) levels. For both cases, the acceptable error level of the permutation test with 1,000 iterations was 5%. The interval mapping (Lander & Botstein, 1989) analysis results were used to identify the nearest marker position to a detectable QTL region that surpasses the LOD threshold, which was determined in the permutation test at the GW level. Besides these two analyses, the Kruskal-Wallis (KW) test was also performed to detect significant marker-trait associations at p < 0.05. The proportion of phenotypic variance explained (PVE) by a single QTL was calculated by the square of the partial correlation coefficient (r²) using MapQTL ver. 6 (Van Ooijen, 2009). Based on the results of the QTL analysis, the genotype effects of significant QTLs for the traits were investigated using the analysis of variance (ANOVA) and t-test. Based on the physical position of the shortlisted markers on the EG5 genome, gene prediction was mined from MPOB’s predicted gene list (http://genomsawit.mpob.gov.my).

SNP and InDel markers identification

A total of 2,684,179 raw reads were obtained. These showed variances from the reference genome with a minimum base quality confidence threshold of 25 and a minimum read depth coverage of five reads in each sample. These raw reads are comprised of 2,530,632 SNPs (94.3%) and 153,547 InDels (5.7%), with a file size of 1.6 GB and 110 MB sequence data, respectively.

For SNP marker discovery, 2,472,925 SNPs that had missing genotype calls of more than 10% across the progeny samples were excluded. The remaining 57,707 SNPs, which represented 2.33% of the total discovered SNP variances, were retained. For InDel markers, 3,854 markers were retained after the elimination process. Of these, 2,154 InDels (55.9%) fit into the genome build, and 1,700 InDels (44.1%) resided on the scaffold.

Among the 57,707 SNPs with a minimum genotype calls of 90%, there were 49,756 SNPs (86.2%) that were filtered out, and 7,951 informative SNPs were left. Ultimately, only 665 InDels (17.3%) were retained for the subsequent analysis. In total, 3,113 SNPs and 225 InDels were eliminated after evaluating by x² test for goodness-of-fit to the Mendelian segregation. As a result, 4,838 SNPs and 440 InDels were informative (Table S1) and were retained for the construction of the genetic linkage maps as shown in Table 1. Based on the results of the PCA analysis (Fig. 1), using a subset of 2,291 informative SNPs, it seems that the progenies belonged to the same cross.

The dura (maternal) and pisifera (paternal) linkage maps

Among the 5,278 shortlisted markers, 3,224 SNPs (66.6%) and 300 InDels (68.2%) were heterozygous in the male parent, whereas 1,605 SNPs (33.2%) and 125 InDels (28.4%) were heterozygous in the female parent. Nine SNPs and 13 InDels showed the segregation patterns of being heterozygous in both parents with two alleles present. There were only two heterozygous InDels with three alleles. Based on the results, the heterozygosity markers were found to be more abundant in the paternal oil palm compared to the maternal oil palm (Table S2). Subsequently, the shortlisted markers, comprising 4,838 SNPs (91.7%) and 440 InDels (8.3%) were used to construct the genetic linkage map for Deli dura x Serdang pisifera. Of these, only 1,754 markers, consists of the segregation patterns lmxll, hkxhk and efxeg were used for the maternal map construction. Whereas 3,548 markers, comprising markers with the segregation patterns of nnxnp, hkxhk and efxeg were used for the paternal map construction. Of these, the redundant markers having identical genotype calls were excluded; 280 for the maternal and paternal map. In total, 236 markers for the maternal map and 109 markers for the paternal map were excluded due to of unsuccessful assignment to any grouping nodes.

Table 2 shows the markers that were assigned to the groups at a LOD score of 5.0 for the genetic linkage map. The major grouping nodes were confidently scored and assigned to 16 chromosomes corresponding with the physical map of the oil palm genome (Singh et al., 2013). Finally, 16 linkage groups for both parental maps were constructed using 1,239 markers for the maternal map (D200) and 2,918 markers for the paternal map (QP447), as shown in Table 3. The linkage groups were named according to the physical map (Singh et al., 2013) for easy traceability. The LOD threshold was set at 5.0 to minimize erroneous grouping assigned for each of the markers. The LOD value of 5.0 was the best synteny to form the grouping with the chromosome from the reference genome based on the physical map information.

Furthermore, markers with a N.N. stress value of more than 4.0 were excluded in order to retain only markers that were correctly assigned to the respective groups (Van Ooijen, 2006). The remaining 1,071 and 2,437 markers for the maternal and paternal maps, respectively, were fully informative in constructing the genetic linkage maps in this CP population. The marker distributions are presented in Fig. 2A for the maternal map (D200) and Fig. 2B for the paternal map (QP447). The maternal (D200, Deli dura) and paternal (QP447, Serdang pisifera) maps are 2,737.6 cM and 4,571.6 cM long, respectively. The average marker density for the maternal map is 2.9 cM; whereas for the paternal map, it is 2.0 cM. Likewise, the recombination rate for the maternal palm is 1.52 cM/Mb which is lower than the paternal palm’s 2.54 cM/Mb.

The SNP and InDel markers used in the construction of the genetic linkage map were further investigated with respect to the physical map’s position. The results in both maternal and paternal maps showed that the genetic linkage map would be able to position the markers on the scaffold to the respective linkage map as shown in Fig. 3.

Verification of map orders

The collinearity of the markers in the genetic linkage and physical map was investigated as shown in Fig. 4. Only the markers that fell on the pseudo-chromosome were used in the investigation as these markers are located in the oil palm reference genome with the physical marker location (Singh et al., 2013). Figure 4 shows that the genetic map constructed in this study was highly correlated to the physical map as indicated by the R-squared (r²) value, especially for linkage groups (LG) LG01, LG02, LG07, LG09, LG10, LG13 and LG14 for both the maternal and paternal maps with r² greater than 0.9. For maternal maps, all the groups showed high collinearity with an r² value of more than 0.7 between the linkage map and the physical map, except for LG06. Similarly, all the groups from the paternal map showed moderate to high collinearity, except for LG16 which showed very low collinearity with an r² value of 0.18.

Distribution, correlation and heritability of the phenotypic data

In the Deli dura and Serdang pisifera cross, the coefficient of variations (CVs) of the phenotypic data were less than 30% for all nine analyzed traits, ranging from 2.2% to 23.3%. As shown in Table 4, the CVs were 17.7% for FFB, 20% for OY, 11.4% for O/B, 2.2% for O/DM, 7.2% for O/WM, 5.2% for M/F, 20.0% for K/F, 23.2% for S/F and 4.9% for F/B. Table 5 shows the results of the pairwise Pearson correlations and the significance levels. The OY and other traits, except for F/B, showed significant correlations above 0.80. On the other hand, K/F and S/F showing a strong negative correlation with M/F is not surprising. A similar correlation between M/F, and K/F and S/F has been reported in previous reports (Ithnin et al., 2017; Jeennor & Volkaert, 2014). Negative correlations were also observed between K/F and S/F and the oil factors, including OY, O/B, O/DM, and O/WM. Contrary to K/F and S/F, M/F showed a significant correlation with low to moderate positive correlations with the oil factors.

Heritability estimation values for the nine phenotypic traits are tabulated in Table 6. The heritability estimation (h²) for FFB, OY, O/DM and O/WM was very low, ranging from 0.47 to 9.29%. The heritability value for F/B was slightly higher at 18.31%. In contrast, a moderate heritability estimation for O/B was recorded at 39.07%. In this population study, the estimation value for M/F was the highest at 97.03%, followed by K/F and S/F which were obtained at 88.92% and 83.96%, respectively.

Identification of significant QTLs on the Deli dura (D200) and Serdang pisifera (QP447)

A total of eight significant markers associated with M/F, K/F and S/F were identified on LG13, which originated from the female parent (D200). Of these, two markers for M/F, five markers for K/F and one marker for S/F were discovered. All associated markers were located on chromosome 13 of the physical map (Singh et al., 2013).

One of the QTLs associated with M/F was located near SNP_sc00239_385519 with a phenotypic variation (PVE) of 20.6%. The SNP marker was located at 35.9 cM on the constructed genetic map and was located on the scaffold of the physical map (Fig. 5). Another QTL detected for M/F was detected on SNP_sc06366_4201 at 60.6 cM. This marker explained 23.9% of PVE and is located on the scaffold of the physical map (Fig. 5). Another QTL peak was detected on the scatter plot, as shown in Fig. 5, at 19.5 cM, but this marker was excluded due to its insignificant value in the KW test.

Meanwhile, five other QTLs were identified for K/F at a GW significant threshold level of 4.1 as shown in Fig. 6. These markers gave a PVE value ranging from 19.4% to 25.6%. Out of the five markers, two of them were InDels namely IND_Chr13_2331610 and IND_sc00391_495837. The former was located on the chromosome while, the latter was positioned on the scaffold of the physical map. Other than InDels, three other SNP markers associated with K/F namely, SNP_Chr13_2672021, SNP_sc00239_385519 and SNP_sc00391_495539 were identified. Their PVE ranged between 19.5% and 25.6%. Of these, only the first SNP was located on the chromosome and the rest were on the scaffold. The result showed that the QTLs for K/F, M/F and S/F were found on the maternal map instead of the paternal map. These traits could possibly be genetically controlled by the maternal line.

Besides QTLs for M/F and K/F, a QTL was also detected on SNP_sc06366_4201 for S/F at the GW’s significant threshold level of 4.2 and PVE which was at 18.1% (Fig. 7). The SNP marker was located on the scaffold and at 60.6 cM in the genetic map. The QTL detected on S/F overlapped with K/F (Fig. 8). Unfortunately, there was no QTL region with a LOD value greater than the GW threshold identified for FFB, OY, O/B, O/DM, O/WM and F/B.

In this study, the shortlisted markers influencing M/F, K/F and S/F were compared with the allele segregation pattern contributing to the traits, as shown in Fig. 9. In general, the shortlisted markers showed a significant contribution to the variation of the phenotypes as depicted in the statistical analysis (ANOVA and t-test) for all the traits. The mean M/F for oil palms with genotype AG was higher compared to those with genotype AA for SNP_sc00239_385519 (Fig. 9A). Nevertheless, palms with genotype AG of SNP_sc06366_4201 had a lower mean M/F than those with genotype GG (Fig. 9B). Palms carrying heterozygous genotypes for these two markers possessed a higher M/F than palms with homozygous genotypes.

For K/F traits, palms with homozygous genotypes AA of SNP_Chr13_2672021, SNP_sc00239_385519 and TT of SNP_sc00391_495539 exhibited higher means of K/F than those with heterozygous genotypes (Fig. 9C–9E). For InDel markers, oil palms with genotypes C/C showed a higher mean K/F than those with genotypes CAT/C of IND_sc00391_495837 (Fig. 9F). Meanwhile, palms with genotype CG/CG showed a higher mean compared to those with genotypes C/CG (Fig. 9G). For S/F, a higher mean was observed when oil palms exhibited genotype AG compared to GG for SNP_sc06366_4201 (Fig. 9H). These effects on the mean phenotypes for genotype segregations are summarized in Table 7.