A SNP variation in an expansin (EgExp4) gene affects height in oil palm

Plant material and phenotype details

The same GT oil palm population that was previously reported (Somyong et al., 2020) was used in this study. The population samples were kindly provided by the Golden Tenera Company Limited, Krabi, Thailand. The GT population contained 180 individuals that were planted in 2008. The population resulted from 30 crosses between six female parents and five male parents. The female parents (dura fruit type) consisted of A 43/9D, A 1/2D, R 15/14D, R 8/9D, R 10/1D and R 10/5D while the male parents (pisifera fruit type) consisted of R 9/8P, R 5/21P, R 3/8P, KA 17/2P and R 16/7P. Height (HT) was recorded four times in three successive years from 8 to 10 year oil palm: in March 2016, October 2016, March 2017 and November 2018.

DNA extraction

DNA samples were extracted from oil palm leaves by using a DNA extraction kit, DNeasy Plant Mini Kit (QIAGEN, Germantown, MD, USA) and using a CTAB/Chloroform-Isoamyl alcohol DNA extraction technique (Cullings, 1992). A mixture of 1 liter CTAB solution contained 1M Tris pH 8.0, 5M NaCl, 0.5M EDTA, and 20 g of CTAB (Cetyl Trimethyl Ammonium Bromide). PVP and β-mercaptoethanol were freshly added to the CTAB solution before it was used for DNA extraction. Quality and concentration of DNA were evaluated by agarose gel electrophoresis and a NanoDrop™ 1000 Spectrophotometer (Thermo Fisher Scientific, Fitchburg, WI, USA).

Oil palm samples, PCR amplification and barcoding preparation for PacBio SMRT sequencing

Forty oil palm individuals of the GT population, including 20 short oil palm individuals and 20 tall oil palm individuals, were selected for PacBio SMRT sequencing, based on their height phenotypes at the 4 recorded times. To obtain full-length genomic DNA sequences of five targeted genes, including EgDELLA1, EgGRF1, EgGA20ox1, EgAPG1 and EgExp4, specific primers were designed at 5′ UTR and 3′ UTR sites, using Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/) (Table S1). The PCR amplification was performed by Phusion High-Fidelity DNA Polymerase (Thermo Fisher Scientific, Baltics UAB, Lithuania) and ran on a Veriti® 96-Well Thermo cycler (Applied Biosystems, Waltham, MA, USA). Two round PCR amplification reported in previous work (Pootakham et al., 2017), with modification, was performed to obtain the barcoded amplicons, which were used for SMRT PacBio sequencing. For the 1^st round of PCR amplification by specific primers tagged with M13F and M13R (Table S2), 10–20 ng of genomic DNA was used in a 20 µl PCR volume. The 20 µl PCR volume consisted of 1 unit (U) of Phusion High-Fidelity (HF) DNA Polymerase, 1× Phusion HF Buffer, 0.2 mM dNTPs, 1.5 mM MgCl₂, 0.1 μM of each primer, and distilled water. The PCR conditions for amplification were the following; initial denaturation at 98 °C for 30 s and then 35 cycles of denaturation at 98 °C for 10 s, annealing at 55 °C for 30 s, elongation at 72 °C for 2 min, and an ending step at 72 °C for 10 min. The PCR products from the 1^st round of PCR were used as the DNA template for the 2^nd round of PCR. For the 2^nd round of PCR by M13F and M13R primers tagged with 16-base PacBio barcodes (Table S3), 2 µl of the 1^st round PCR product was used as the template in 50 µl PCR volume. A final volume of 50 µl consisted of 1 U of Phusion HF DNA polymerase, 1× Phusion HF Buffer, 0.2 mM dNTPs, 1.5 mM MgCl₂, 0.1 μM of each primer, and the remaining volume of distilled water. The same PCR conditions explained above were used. The DNA size of 5–10 µl of the PCR products was determined by agarose gel electrophoresis. The remaining 40–45 µl of PCR products was purified with Agentcourt AMPure magnetic beads (Beckman Coulter, Indianapolis, IN, USA), and finally diluted in 10 µl of 1× TE. The concentration of barcoded amplicons was measured using a Qubit 2.0 fluorometer and a Qubit dsDNA BR assay kit (Thermo Fisher Scientific, Waltham, MA, USA). The purified DNA was pooled in equimolar concentration in at least 500 ng of the pooled DNA (a volume of 30–40 µl). The pooled DNA was used to construct the SMRTbell libraries. The barcoded amplicons were ligated into SMRTbell adapters. All libraries were sequenced on a PacBio RSII system (Pacific Biosciences, Menlo Park, CA, USA), using the P6-C4 chemistry with 360-min movie lengths.

Sequence data analysis

The PacBio raw reads of targeted genes, including EgDELLA1, EgGRF1, EgGA20ox1, EgAPG1 and EgExp4, processed by PacBio RSII, were analyzed by SMRT analysis software package (version 2.3). The Circular Consensus Sequencing (CCS) reads were obtained using RS_ReadsOfInsert protocol with –minPasses 5 –minPredictedAccuracy 0.99 –maxLength 3,500 into a FASTQ file. These files were first separated by using barcodes, according to each oil palm sample, and then were grouped for each gene. This file separation or classification was performed by using Python script and the MUSCLE program for data alignment. Next, the CCS reads of each oil palm sample for each gene were used to select the full-length genomic DNA sequences. The CCS reads of the full-length genomic DNA sequences of each gene were pooled according to short or tall oil palm group. The grouped full-length genomic DNA sequences were mapped to the reference gene sequences using minimap2 version 2.17-r943-dirty. SNP and InDel polymorphic sites among short and tall oil palm groups were identified by GATK version 4.1.1.0. The SNPs and InDels with depth coverage of more than 20× were used for developing molecular markers for oil palm breeding programs. Details of CCS reads and the full-length genomic DNA sequences of each gene can be found at the NCBI database under accession number: PRJNA760254.

SNP Marker designing for high through-put PACE™ SNP genotyping

This study focused only on SNP marker development. SNP primers were designed by 3CR Bioscience Co. Ltd. SNP primers including specific nucleotides at each SNP variation at the targeted genes and nucleotides required for fluorescent emission (Table S4). PCR Allelic Competitive Extension (PACE) is an allele-specific technology for SNP genotyping that includes 2 allele-specific forward primers and 1 common reverse primer. Universal PACE™ Genotyping Master Mix (Standard ROX) (Essex CM20 2BU; 3CR Bioscience, UK) contained solution required for PCR amplification and fluorescent detection, including FAM (blue emission) and HEX (red emission) specified for each SNP change. Details of reaction preparation followed the PACE™ User Guide v1.6 (www.3crbio.com) with modification. The reaction mixes were placed in a 384-well plate that included 1–10 ng DNA template per reaction well and prepared by the wet DNA method. Five µl reaction volume comprised of 2.5 µl of DNA template, 2.5 µl of universal PACE™ Genotyping Master Mix and 0.07 µl of PACE assay mix. Universal PACE™ Genotyping Master Mix contained Tag DNA polymerase, universal fluorescent reporting cassette, dNTPs, performance enhancers, MgCl₂, and 5-carboxyl-x-rhodamine, succinimidyl ester (ROX). The assay mix contained allele-specific primer1-FAM, allele-specific allele2-HEX and common reverse primer in the ratio of final concentration as 0.168 µM: 0.168 µM: 0.42 µM, respectively, for each reaction. The PCR amplification was run on QuantStudio 6 Flex real-time PCR system (Thermo Fisher Scientific, Waltham, MA, USA). The PCR reaction for PACE genotyping was the following. For the 1^st step, enzyme activation was conducted at 94 °C for 15 min. For the 2^nd step, 10 cycles of touch down PCR were conducted by denaturation at 94 °C for 20 s, and touch down annealing and extension starting at 61 °C to 55 °C for 60 s by decreasing 0.6 °C per cycle. For the 3^rd step, 25 cycles were conducted by denaturation at 94 °C for 20 s, annealing and extension at 57 °C for 60 s and an ending step at 37 °C for 60 s. QuantStudio™ Real-Time PCR Software v1.3 was used to analyze SNP genotypes, which were two separated homozygous genotypes emitting either FAM (blue) or HEX (red) and heterozygous genotype emitting both fluorescences (green).

Population structure, statistical and association analyses

The contribution of targeted SNP markers to the phenotype traits was analyzed by ANOVA and association analysis. Descriptive statistics of height phenotype data (HT) were analyzed by SPSS 11.5. This statistic package was also used to analyze preliminary relationships between the polymorphic loci of targeted genes with HT by comparing mean height using One-Way ANOVA. The significant association of the polymorphic loci with the traits was then analyzed by TASSEL 2.1 (http://www.maizegenetics.net/#!tassel/c17q9) by GLM with Q model with 10,000 permutations that included Q-matrix. The required information for the association analysis included genotype, height phenotype and Q-matrix information. Inferred ancestry of individuals of optimal K value from STRUCTURE output was used as Q-matrix information. STRUCTURE 2.3.4 (http://pritchardlab.stanford.edu/structure.html) and STRUCTURE harvester (http://taylor0.biology.ucla.edu/structureHarvester/) were used to analyze population structure and determine optimal K value, respectively. The inferred ancestry of individuals was used as Q-matrix information by setting its value as covariance in the association analysis of the targeted markers with height. Details of STRUCTURE analysis of the GT population were already explained in our previous work (Somyong et al., 2019). Briefly, STRUCTURE was run three times at this setting: Length of Burnin Period = 50,000, number of MCMC Reps after Burnin = 50,000 and K setting = 1–10. The optimal K was determined by using STRUCTURE Harvester (Earl & vonHoldt, 2012) and the optimal K calculated using the Evanno method (Evanno, Regnaut & Goudet, 2005). The optimal K of the GT population (180 individuals) was 3 because the highest Delta K was found for this value (Somyong et al., 2019). For marker information, previous markers (Somyong et al., 2019; Somyong et al., 2020), including mEgDELLA1-1, mEgDELLA1-11, mEgACCO-pr2, mEgSSRffb10-8, mEgGRF1-3 and mEgGA20ox1-4, and an additional marker, mEgExp4_SNP118 from this study, were used in TASSEL analysis. The p value applied to determine significance of marker-trait association was less than 0.05.

Details for height phenotype of the GT population

Height information of HT-1, HT-2, HT-3 and HT-4 for the GT population (180 individuals) was explained in our previous work (Somyong et al., 2020). HT-1, HT-2 and HT-3 were recorded in 6-month intervals in 8–9 year oil palm plants while HT-4 was recorded in 10-year oil palm plants. The mean height of the GT population was 192 cm for HT-1, 231 cm for HT-2, 263 cm for HT-3, and 380 cm for HT-4, with a 75 cm average increase within a year. To identify nucleotide variations of height-related genes by PacBio SMRT sequencing, 40 individuals were selected based on height phenotype. Height phenotype was classified as short and tall groups, including 20 individuals of the shortest oil palm plants and 20 individuals of the tallest oil palm plants from the GT population. Height phenotype of the same short and tall oil palm individuals was recorded for all 4 height recordings. The height distribution of the short and tall oil palm groups is illustrated in Fig. S1. For HT-1, average height was 132 cm for the short group and 254 cm for the tall group. For HT-2, average height was 164 cm for the short group and 300 cm for the tall group. For HT-3, average height was 190 cm for the short group and 342 cm for the tall group. For HT-4, average height was 289 cm for the short group and 471 cm for the tall group. The height difference between short and tall oil palm groups was 122 cm for HT-1, 136 cm for HT-2, 152 cm for HT-3 and 182 cm for HT-4. Height between short and tall groups was significantly different for the 4 height recordings with a p value of 0.000 by using t-test.

Primer testing for amplification of full-length genomic DNA sequences of the height-related genes

The first step of this study was primer design and full-length amplification testing in short and tall oil palm groups from the GT population. Five candidate genes, EgDELLA1, EgGRF1, EgGA20ox1, EgAPG1 and EgExp4, were targeted in this study. The targeted sites of primer design were 5′ UTR and 3′ UTR sites of the genes. The full-length genomic DNA sequences of these genes were obtained from oil palm draft sequences of the Malaysian Palm Oil Board (MPOB) (http://genomsawit.mpob.gov.my/index.php?track=30). For full-length gene amplification, two primer pairs were designed for each gene from 5′ UTR and 3′ UTR sites (Table S1). The selected primers that amplified the full-length gene products of expected sizes, from 2,516 to 3,015 bp, included EgDELLA1-P1, EgGRF1-P2, EgGA20ox1-P2, EgAPG1-P1 and EgExp4-P1 (Fig. S2).

Sample preparation for PacBio SMRT sequencing of the height-related genes

The 40 short and tall oil palm samples were prepared for PacBio SMRT sequencing of the EgDELLA1, EgGRF1, EgGA20ox1, EgAPG1 and EgExp4 genes. The full-lengths of these genes were successfully amplified by the M13-tagged EgDELLA1-P1, EgGRF1-P2, EgGA20ox1-P2, EgAPG1-P1 and EgExp4-P1 primers (Table S2) in the 1^st round of PCR with bands close to the expected sizes of 3,015 bp, 2,516 bp, 2,759 bp, 2,917 bp and 2,586 bp, respectively. The PCR products from the 1^st round of PCR were used as a DNA template for the 2^nd round of PCR. These M13-tagged PCR products were then carried to the second amplification with several combinations of barcode-tagged M13 primer sets (Table S3). Each barcode-tagged M13 primer set consisted of forward and reverse primers that were used in the separation of each oil palm sample for each gene amplification. The results of amplification by M13-tagged primers and barcode-tagged M13 primers also had bands close to the expected sizes of 3,015 bp, 2,516 bp, 2,759 bp, 2,917 bp and 2,586 bp respectively, as shown in Fig. S3.

Analysis of sequence variations among short and tall oil palm groups

The full-length genomic DNA sequences of EgDELLA1, EgGRF1, EgGA20ox1, EgAPG1 and EgExp4 genes of short and tall oil palm samples were revealed by PacBio SMRT sequencing technology (Table 1). For EgDELLA1, 37 of the 40 total oil palm plants, including 18 short and 19 tall oil palms, have been successfully sequenced, representing 1,166 full-length genomic DNA sequences with a size range of 2,959–3,246 bp. For EgGRF1, 37 of the 40 total oil palm plants, including 18 short and 19 tall oil palms, have been successfully sequenced, representing 909 full-length genomic DNA sequences with a size range of 1,933–2,607 bp. For EgGA20ox1, 36 of the 40 total oil palm plants, including 18 short and 18 tall oil palms, have been successfully sequenced, representing 1,494 full-length genomic DNA sequences with a size range of 2,138–2,913 bp. For EgAPG1, all 40 total oil palm plants have been successfully sequenced, representing 387 full-length genomic DNA sequences with a size range of 2,060–2,952 bp. For EgExp4, all 40 total oil palm plants have also been successfully sequenced, representing 5,384 full-length genomic DNA sequences with a size range of 2,493–2,705 bp. The results showed that the highest number of full-length genomic DNA sequences was found from EgExp4 while the lowest number of sequences was found from EgAPG1. The variations of the number of full-length genomic DNA sequences results from amplification ability during the two rounds of PCR amplification.

The results of sequence variation, position and variation type of all these genes, among short and tall oil palm groups, are listed in Table 2. Variations were found for 4 of the 5 genes. No variation was found for the EgAPG1 gene. For the EgDELLA1 gene, there was a total of 1,166 sequences, representing short and tall oil palm groups with 578 and 677 sequences, respectively. We found 3 variations, containing 1 insertion and 2 SNPs. The position of the insertion was at 312 (T/TA), and the SNP variations were at 2,100 (T/A) and 2,248 (G/A), all of which are illustrated on the EgDELLA1 reference gene sequence, which has a full-length of 3,015 bp (Fig. S4). For the EgGRF1 gene, there was a total of 909 sequences, representing short and tall oil palm groups with 468 and 441 sequences, respectively. We found 4 variations, including 2 deletions and 2 insertions. The positions of deletions were at 1,044 (TATA/T) and 1,553 (ATCTCTC/A), and insertions were at 1,100 (G/GT) and 1,553 (A/ATCTC), all of which are illustrated on the EgGRF1 reference gene sequence, which has a full-length of 2,516 bp (Fig. S5). For the EgGA20ox1 gene, there was a total of 1,494 sequences, representing short and tall oil palm groups with 767 and 726 sequences, respectively. We found 3 variations, including 2 insertions and 1 SNP. The positions of insertions were at 1,428 (G/GAA) and 1,428 (G/GA, GAAA), and the SNP variation was at 1,468 (T/G), all of which are illustrated on the EgGA20ox1 reference gene sequence, which has a full-length of 2,759 bp (Fig. S6). For the EgExp4 gene, there was a total of 5,384 sequences, representing short and tall oil palm groups with 3,249 and 2,135 sequences, respectively. We found 2 variations containing 1 deletion and 1 SNP. The position of the deletion was at 989 (TC/T), and the SNP variation was at 118 (T/C), both of which are illustrated on the EgExp4 reference gene sequence, which has a full-length of 2,586 bp (Fig. S7). No variation was found for the EgAPG1 gene, which is also illustrated on the EgAPG1 reference gene sequence, which has a full-length of 2,917bp (Fig. S8). This may be due to either no-true variation or an inadequate amount of the full-length genomic DNA sequences that is found between short and tall groups for this gene.

High through-put genotyping by PACE™ SNP genotyping

This work has targeted only SNP marker development. PACE™ SNP genotyping was used in this study. Four SNP primer sets were designed from three target genes, including EgDELLA1, EgGA20ox1 and EgExp4, and are listed in Table S4. These SNP primer sets included mEgDELLA1_SNP2100, mEgDELLA1_SNP2248, mEgGA20ox1_SNP1468 and mEgExp4_SNP118. Three of the SNP primer sets, including mEgDELLA1_SNP2100, mEgDELLA1_SNP2248 and mEgExp4_SNP118, were polymorphic in the GT population while mEgGA20ox1_SNP1468 was monomorphic in the same population. The allelic discrimination plots of these polymorphic markers are illustrated in Fig. 1. Genotypes of mEgDELLA1_SNP2100 included A/A (10 individuals), T/A (65 individuals) and T/T (103 individuals) (Fig. 1A). Genotypes of mEgDELLA1_SNP2248 included G/A (51 individuals) and G/G (124 individuals) (Fig. 1B). Genotypes of mEgExp4_SNP118 included C/C (25 individuals), T/C (78 individuals) and T/T (66 individuals) (Fig. 1C). The genotypes of mEgExp4_SNP118 were found to have the most distribution of allele frequencies in the two homozygous groups and heterozygous group when compared with the other two markers.

ANOVA and association analysis of mEgExp4_SNP118 with height

ANOVA analysis shows that SNP changes of mEgExp4_SNP118 in the GT population, from nucleotide T to C, affected height significantly while that of mEgDELLA1_SNP2100 and mEgDELLA1_SNP2248 did not. Height details and the height distribution of mEgExp4_SNP118 genotypes are shown in Table S5. We suggest that the EgExp4 gene contributes to the height trait in the GT population. ANOVA analysis confirmed that mEgExp4_SNP118 genotypes, C/C, T/C and T/T, have significant height differences between them, for all 4 height-recordings, with p values of 0.002–0.047 (Table 3A). Oil palm individuals with genotype C/C were significantly taller than individuals with genotype T/C and genotype T/T in all 4 height-recordings, by 17–39 cm. In addition, genotypes T/C and T/T did not have significant height differences between them.

For HT-1, individuals with genotypes C/C were significantly taller than individuals with genotypes T/C and T/T by 17 cm and 25 cm, respectively. For HT-2, individuals with genotype C/C were significantly taller than individuals with genotypes T/C and T/T by 18 cm and 26 cm, respectively. For HT-3, individuals with genotype C/C were significantly taller than individuals with genotypes T/C and T/T by 15 cm and 21 cm, respectively. For HT-4, individuals with genotype C/C were significantly taller than individuals with genotypes T/C and T/T by 29 cm and 39 cm, respectively. Even though, the individuals with genotype T/C had intermediate height, they were not significantly different from the individuals with genotype T/T. We suggest that a SNP change from T to C, in position 118 of EgExp4, is involved in increasing height, while T acts in the opposite way by decreasing height. The contribution of EgExp4 to height was further confirmed by association analysis using TASSEL.

To perform association analysis by TASSEL, besides the genotype and phenotype data, a Q-matrix from the GT population was needed. The Q-matrix was determined by STRUCTURE analysis and was shown in our previous study (Somyong et al., 2019, 2020). The trait information consisted of HT-1, HT-2, HT-3 and HT-4. After association analysis, with the permutation value set as 10,000, we found a significant association between mEgExp4_SNP118 and height, as shown in Table 3B. The mEgExp4_SNP118 marker was found to be significantly associated with height in 2 from 4 height-recordings with p values of 0.0383 for HT-1 and 0.0263 for HT-4. This result supports data from the ANOVA analysis. Based on our previous work and this work, the three markers that have been shown to associate with height include mEgDELLA1-1 (Somyong et al., 2019), mEgGRF1 (Somyong et al., 2020) and mEgExp4_SNP118 that target DELLA, growth regulating factor and expansin genes, respectively.