Transcriptome analysis of oil palm inflorescences revealed candidate genes for an auxin signaling pathway involved in parthenocarpy

Plant materials

An oil palm variety, commercially named as Surat Thani 2, was selected for our study because of its long bunch stalks, allowing for clear female flower visibility. A total of 15 6-year-old oil palm plants with visible female inflorescences during an anthesis stage were selected. Oil palm, Surat Thani 2, is a cross between Deli dura (mother plant) and Lame (father plant), developed by Department of Agriculture, Ministry of Agriculture and Cooperatives, Thailand (http://www.doa.go.th/main/). The oil palm trees were grown in Saraburi, Thailand by the Oil Palm Technology Development for commercial bio diesel industry in newly planted area project (OPTD), Department of Agronomy, Faculty of Agriculture, Kasetsart University (Bangken Campus).

Hormonal treatment and collection of oil palm pistils

Synthetic auxin, 2, 4, 5-tri-chlorophenoxy propionic acid (2, 4, 5-TP), which has been reported to induce parthenocarpy in oil palm (Thomas et al., 1973), was used for spraying the female inflorescences. One gram of 2, 4, 5-TP was first dissolved in 100 ml methanol. The solution was poured into a five l pressure sprayer and then filled with distilled water to make a total volume of five l. About 250 ml (about 50 mg of 2, 4, 5-TP) of the solution was sprayed on each inflorescence, using the pressure sprayer. The treatments were performed on three consecutive days during different pollination stages and again after 7 and 14 days. The oil palm plants with visible female inflorescences selected for the treatment had female inflorescences in different stages of anthesis. These anthesis stages included the stage with unopened stigma, assigned as day after pollination (DAP) = 0, stigma opened for 1 day (DAP = 1), stigma opened for 2 days (DAP = 2) and 3 days (DAP = 3). So, the stages with the first day of 2, 4, 5-TP spraying included DAP = 0, DAP = 1, DAP = 2 and DAP = 3. In summary, 13 female inflorescences (named as Inflo. 1T-13T) of 13 oil palm trees of four different pollination stages were subjected to auxin treatment and four female inflorescences (Inflo. 1C-4C) of four oil palm trees were used as the controls, without auxin treatment. Details of the oil palm trees used for both treated and untreated inflorescences and the first day of treated stages are listed in Table S1. The parthenocarpy phenotype was observed 2 weeks, 6 weeks and at the ripening stage (18 weeks) after pollination. For the auxin treated oil palm plants, five female flowers were collected for each inflorescence before the first auxin spraying and again 24 and 48 h after the spraying, for RNA extraction. The collected flowers were immediately frozen in liquid nitrogen during collection and stored at −80 °C until the extraction was performed. For the control oil palm plants, without the auxin treatment, the pistil samples were collected at the same stages as were the treatments.

Summary of RNA samples and pooling strategy

For expression studies, oil palm plants that were treated with auxin on DAP = 0 and DAP = 1 were selected because spraying during these stages appeared to result in a higher degree of parthenocarpy. Details of the three control and three auxin-treated inflorescences, used for pistil collection at the given stages, are represented on Table S2. The RNA samples used as the auxin treated group included Inflo.1T/WA, Inflo.6T/WA and Inflo.4T/WA. The RNA samples used as the control group included Inflo.1T/NA, Inflo.1C/NA, Inflo.2C/NA1, Inflo.2C/NA2 and Inflo.4C/NA. Since the total RNA for each pollination stage of some samples was not enough for mRNA isolation, some RNA samples were pooled. The pooled samples included RNA samples, Inflo.1T/WA (24 and 48 h auxin-treated samples were pooled), Inflo.6T/WA (24 and 48 h auxin-treated samples were pooled), Inflo.4T/WA (24 and 48 h auxin-treated samples were pooled), Inflo.1C/NA (samples that were not auxin-treated at DAP = 0 of day 2 and 3 collection were pooled), Inflo.2C/NA1 (samples that were not auxin-treated at DAP = 1 and 2 were pooled) and Inflo.4C/NA (samples that were not auxin-treated at DAP = 2 and 3 were pooled) (Table S2).

RNA extraction, mRNA isolation and mRNA purification

RNA extraction was performed on a large scale by the CTAB method. A total RNA of about 75 μg was needed for mRNA isolation. Five oil palm pistils for each inflorescence were ground in liquid nitrogen and transferred to five two ml tubes containing 1,000 μl extraction buffer (2% CTAB, 1.4 M NaCl, 2% PVP, 20 mM EDTA pH 8, 100 mM Tris–HCl pH 8) and 20 μl of 20% SDS. After adding β-mercaptoethanol, the extraction was cleaned with chloroform: Isoamyl (24:1) for five times. The final supernatant was precipitated with 1/3 volume of 8M LiCl resulting in a final concentration of 2M LiCl. The pellet was washed with 70% ethanol, then air dried and resuspended in 50–100 μl of DEPC water. The quality and concentration was determined by 1% agarose gel electrophoresis and NanoDrop^®ND-1000.

Isolation of mRNA was performed using the Dynabeads mRNA purification kit (Thermo Fisher Scientific Baltics UAB, Vilnius, Lithuania). About 75 μg of total RNA was needed and dissolved in DEPC water to make final volume of 100 μl. If the RNA sample was less than 75 μg diluted in more than 100 μl of solution, whole volume was used. After purification by the kit, the mRNA was eluted by adding 25 μl of 10 mM Tris–HCl pH 7.5. The quality and concentration were determined on a Fragment Analyzer (Advanced Analytical Technologies, Inc., Ames, IA, USA), using the DNF-472M33 kit. The mRNA was diluted to about 1,000–2,500 pg/μl and two μl of mRNA was used for the determination. Purification of mRNA was performed by using the Magnetic Bead Cleanup Module (Thermo Fisher Scientific, Waltham, MA, USA). Five μl of 37 °C pre-warmed nuclease-free water was used to elute purified mRNA and the solution was kept for preparing the barcoded RNA library.

Barcoded RNA library preparation

The barcoded RNA libraries were made following the Ion Total RNA-Seq Kit v2 (Life Technologies, Carlsbad, CA, USA) protocol. Briefly, all five μl of mRNA from the previous step was hybridized and ligated with the Ion Adaptor Mix v2. The ligation reaction was performed on a Veriti 96 well Thermal Cycler (Applied Biosystems, Foster City, CA, USA), at 30 °C for 1 h. Then, reverse transcription was performed using SuperScript^® III Enzyme Mix and incubated at 42 °C for 30 min. The cDNA was purified by a magnetic bead cleanup module. A total of 12 μl of 37 °C pre-warmed nuclease-free water was used to elute the purified cDNA. Next, the barcodes were added in the amplifying cDNA step. Six μl of purified cDNA was used. The barcoded libraries were made using Platinum^® PCR SuperMix High Fidelity, Ion Xpress™ RNA 3′ Barcode primer and Ion Xpress™ RNA-Seq Barcode BC primer. In this study, the BC primers included BC1-BC8 and BC12. The reaction was run on a thermal cycler with the following setting. The reaction was first incubated at 94 °C for 2 min, then performed for two cycles at 94 °C for 30 s, 50 °C for 30 s, 68 °C for 30 s, then performed for 16 cycles at 94 °C for 30 s, 62 °C for 30 s, 68 °C for 30 s and ended at 68 °C for 5 min. The amplified cDNA was also purified using a magnetic bead cleanup. A total of 15 μl of 37 °C pre-warmed nuclease-free water was used to elute the purified amplified cDNA.

The quality and concentration were determined on a Fragment Analyzer (Advanced Analytical Technologies, Inc., Ames, IA, USA), using the DNF-474-33 kit (HS NGS fragment 1-6000 bp). The cDNA was diluted to about 2,500–5,000 pg/μl and two μl of cDNA was also used for the determination.

Ion proton sequencing and transcriptional analysis

Details of RNA sequencing are explained in the Ion 540™ Kit-OT2 User guide. Briefly, each barcoded cDNA library from the previous step was diluted to 100 pM. Then, the libraries were pooled and six μl of the pooled cDNA was used for template preparation. The template-positive Ion 540™ ion sphere particles (ISPs), containing amplified cDNA (about 200 bp), were prepared using the Ion 540™ kit-OT2 kit and the reaction was performed in the One Tough™ 2 Instrument. The template-positive ISPs were enriched using Ion One Tough™ ES. Next, the enriched template-positive ISPs were loaded on an Ion 540™ Chip and sequenced on an Ion S5™ XL sequencer.

For the expression analysis, raw reads were first mapped to the oil palm genome (https://www.ncbi.nlm.nih.gov/genome/2669), which is E. guineensis assembly EG5, by using bowtie2-2.2.3. Reads in gene regions were counted by HT-seq count version 0.9.1 .The raw read count was normalized using the R package DESeq2. The expression fold-change of each gene was calculated between treatment and control conditions with the threshold for DEGs set to p-value < 0.05.

Reverse transcription polymerase chain reaction analysis and quantitative real-time RT-PCR of candidate genes controlling parthenocarpy

To confirm the RNA seq expression of candidate genes, the same amplified cDNA samples of mRNA used for RNA seq were subjected to RT-PCR. The RT-PCR reaction was performed in three replicates, for both target and reference genes. Cyclophilin 2 (renamed as EgCyp2) was chosen for the reference gene (Yeap et al., 2014). The primer sequences were 5′ CTCGTCTGATGTCGTCTA 3′ as the forward primer and 5′ CTGCTGGTACTCTGGTAA 3′ as the reverse primer. The RT-PCR reaction was conducted in a 20 μl solution, containing three μl of cDNA (one ng/μl of the amplified cDNA), 1× High Fidelity buffer, two mM MgSO₄, 400 μM dNTP, 0.5 μM of each primer and 0.5 unit of Platinum^®Tag DNA polymerase High Fidelity (Invitrogen by Thermo Fisher Scientific, Waltham, MA, USA). The amplification was performed on the Veriti 96 Well Thermocycler (Applied Biosystems by Thermo Fisher Scientific, Waltham, MA, USA). The PCR reaction was set as the following; 5 min at 94 °C of initiation PCR activation, 35 cycles of 1 min at 94 °C, 0.30 min at 55 °C and 1 min at 68 °C. The reaction ended in a 7 min-step at 68 °C. The products were analyzed on 1% agarose gel, containing Invitrogen SYBR safe DNA gel stain.

Quantitative real-time RT-PCR was performed using Fast Start SYBR Green Master (Roche Applied Science, Penzberg, Germany). The same cDNA and reference gene, EgCyp2, as RT-PCR was used in this qRT-PCR. The reaction was conducted in a 20 μl solution, containing three μl of cDNA (one ng/μl of the amplified cDNA), 1× FastStart SYBR Green Master, 0.5 μM of each primer. The qRT-PCR reaction was performed in three replicates, for both target and reference genes. The amplification was performed on the CFX 96 Real-Time PCR detection system (Bio-Rad Laboratories, Inc., Hercules, CA, USA) at this setting; one cycle of 3 min at 95 °C, 50 cycles of 0.10 min at 95 °C, 0.15 min at 55 °C and 0.30 min at 72 °C, then the reaction was ended with 95 °C for 0.10 min. This was followed with the melting steps from 65 °C to 95 °C by increments of 0.5 °C for 0.05 min in each step. The fold expression was calculated using the normalized expression (ΔΔCq) method, normalized to EgCyp2 expression by CFX manager software. Settings of the graph data were X-axis as target, Y-axis as linear, scaling as unscaled and error type as standard error of the mean.

Characteristics of stigmas during different pollination stages

The characteristics of stigma during pollination indicate the occurrence of pollen-pistil interaction. If a stigma is still white and closed, it means that pollen-pistil interaction has not occurred. We assigned stage DAP = 0 as the stage in which the stigma is still white and closed and pollination has not taken place. The pollination usually occurs during DAP = 1 to DAP = 2 stages because the white stigma is still opened for day 1 and day 2, respectively. The stigma has a light pink color during stage DAP = 3, a dark pink color during DAP = 4, and then a dry black color by DAP = 10. Auxin treatment may be most effective for inducing parthenocarpy if pollination has not occurred. This work started spraying auxin in four different stages, DAP = 0, DAP = 1, DAP = 2 and DAP = 3, as shown in Fig. S1.

The effect of auxin treatment on parthenocarpy before the ripening stage

Before the ripening stage, the effect of auxin on parthenocarpy was observed 1, 2 and 6 weeks after pollination. Several characteristics were recorded, including stigma characteristics, such as color, opening or closing, dryness, the fertility of fruits, fruit shape and the presence of kernels. In control inflorescences, normal fruits developed 6 weeks after pollination. The fruits contained soft shells and liquid endosperm, which is part of the kernel. Black dry stigmas were still attached at the top of the fruits. However, there were several changes found in the treated fruits, when compared with the control fruits. Stigmas of the treated fruits (first treated at DAP =1, DAP = 2 and DAP = 3) were still fresh and red for 2 weeks while stigmas of the control fruits were dry and black since DAP = 5, as in Fig. S2A. Moreover, stigmas of Inflo.1T, first treated with auxin at DAP = 0, were still fresh and red, and the freshness and redness lasted for 6 weeks after pollination, as shown in Fig. S2B. The development of treated oil palm fruits of Inflo.4T was delayed, causing fruits to be smaller than those of controls. Some fruits, including some from Inflo.5T and Inflo.8T, failed to develop. The shape of fruits produced by treated plants was more oval than that of the controls, which produced more rounded fruits. However, the majority of treated oil palm fruits developed to the ripening stage successfully. In addition, auxin had an effect on oil palm plants that were treated on DAP = 3 because the fruits developed kernels. The development, however, was slower than that of the controls.

The effect of auxin treatment on parthenocarpy at the ripening stage

The effect of auxin on parthenocarpy was observed about 4.5 months after pollination (ripening stage) that seedless fruits was the most visible by cutting fruits in half. All fruits of the control oil palm plants (Inflo.1C-4C) contained kernels, hard shells and oily mesocarp. Auxin treatment had the greatest parthenocarpic effect when auxin was sprayed on DAP = 0, when stigmas had not yet opened. This experiment was performed on two oil palm trees (Inflo.1T and Inflo.2T). The treated bunches were smaller than the controls (Fig. 1A) and the treated fruits showed 100% parthenocarpy (n = 30), without seeds (no white kernel), as shown in Fig. 1B. When auxin was first sprayed on DAP = 1, the effect of auxin on parthenocarpy was reduced because some pistils were pollinated. The experiment was conducted on four oil palm plants (Inflo.3T-Inflo.6T). The fruit set of Inflo.5T failed, while the other three oil palms showed varying degrees of parthenocarpy, with Inflo.3T, Inflo.6T and Inflo.4T having 27%, 50% and 90% parthenocarpic fruits from 20 randomly selected fruits for each oil palm plant, as shown in Fig. 2. Moreover, the effect of auxin on parthenocarpy was reduced further on DAP = 2 because pistils were pollinated more 2 days after stigma opening. This auxin treatment was performed on DAP = 2 in four oil palm trees (Inflo.7T-Inflo.10T). The levels of parthenocarpy were variable with Inflo.10T, Inflo.7T, Inflo.8T and Inflo.9T having 5%, 6%, and 10% parthenocarpic fruits from 20 randomly selected fruits for each oil palm, respectively. Even though the pollen-pistil interaction was expected to stop at DAP = 3, we still found that auxin treatment induced some parthenocarpic fruits but to a smaller degree. These oil palms including Inflo.12T, Inflo.11T and Inflo.13T, which contained 0%, 5% and 10% parthenocarpic fruits from 20 randomly selected fruits for each oil palm, respectively. The above information indicates that hormone treatment can induce parthenocarpic development throughout fruit progression to the ripening stage. The hormonal treatment caused 100%, 55%, 8% and 5% parthenocarpy when the inflorescences were auxin-treated on DAP = 0, DAP = 1, DAP = 2 and DAP = 3, respectively.

RNA sequencing and data analysis

Eight barcoded RNA libraries from the above eight RNA samples were used for Ion Proton sequencing. Raw reads ranged from 8,425,859 to 11,811,166 reads, with average sizes ranging from 99 to 137 bp, as shown in Table 1. The raw data was deposited in the NCBI database (https://www.ncbi.nlm.nih.gov/) (SRA accession: PRJNA490783). After comparing with the reference oil palm genome (assembly EG5), the mapped reads ranged from 8,179,948 to 11,320,799 reads, representing 95.85–98.01% of the map. The DEGs of the auxin-treated oil palms were compared with those of the control oil palms. There are two groups that were separated for comparisons. For oil palms that were first treated on DAP = 0, the data for comparison includes RNA expression data from RNA libraries, Inflo.1T/NA, Inflo.1T/WA and Inflo.1C/NA. For oil palms that were first treated on DAP = 1, the data for comparison includes RNA expression data from RNA libraries, Inflo.6T/WA, Inflo.4T/WA, Inflo.2C/NA1, Inflo.2C/NA2 and Inflo.4C/NA. The expression-fold-change of each gene was calculated between treatment and control oil palm plants with the threshold for DEGs set to a p-value <0.05. To reveal DEGs, eight pair comparisons were performed between treated and control samples, listed in Table S3. The number of DEGs found was from 120 to 246 genes, for each library. About 9% of the DEGs are categorized by KEGG (Kyoto Encyclopedia of Genes and Genomes; http://www.genome.jp/kegg/) (Fig. S3). A total of 83 from a total of 136 KEGG categorized DEGs are nonredundant genes. The KEGG pathway that contains the highest number of down-regulated DEGs is amino sugar and nucleotide sugar metabolism. These genes included chitinase 1-like, endochitinase EP3-like, GDP-mannose dehydratase 1-like, and UDP-glucose 6-dehydrogenase 5-like. The KEGG pathway that contains the most up-regulated DEGs is drug metabolisms, including cytochrome P450, metabolism of xenobiotics by cytochrome P450 and glutathione metabolism. These genes included glutathione S-transferase U8-like, probable carboxylesterase 18 and probable glutathione S-transferase GSTF1. Some DEGs that are not categorized by KEGG are in the auxin signaling pathway, which is the pathway that we are focusing on in this experiment. These genes are the potential candidate genes involved in parthenocarpy in oil palm.

Candidate genes in the auxin signaling pathway controlling parthenocarpy

We proposed candidate genes involved in the auxin signaling pathway by the eight pair comparisons explained earlier. Five out of eight comparisons between the treatments and controls revealed six candidate genes in the auxin signaling pathway, showing differential expression during auxin treatment on DAP = 0 and DAP = 1, as shown in Table 2. In the treatment at DAP = 0 group, two candidate genes were found, including probable flavin-containing monooxygenase 1 (cds26665) and auxin-responsive protein SAUR71-like (cds32638). In the treatment at DAP = 1 group, four candidate genes were found, including probable indole-3-acetic acid (IAA)-amido synthetase GH3.8 (cds27914), probable IAA-amido synthetase GH3.1 (cds30335), IAA-induced protein ARG7-like (cds31450) and tryptophan aminotransferase-related protein 3-like (cds7351).

From plants treated with auxin on DAP = 0, we found flavin-containing monooxygenase 1 (EgFMO1) to be the strongest candidate gene for parthenocarpy because it showed differential expression, when compared with the two controls, Inflo.1T/NA and Inflo.1C/NA. Another gene, auxin-responsive protein SAUR71-like (EgSAUR71) was found to be differentially expressed, compared to Inflo.1C/NA. In plants treated on DAP = 1, probable IAA-amido synthetase GH3.8 (EgGH3.8) is the strongest candidate for parthenocarpy because it was differentially expressed, when compared with the three controls, Inflo.2C/NA1, Inflo.2C/NA2 and Inflo.4C/NA. The next potential candidate for parthenocarpy is probable IAA-amido synthetase GH3.1 (EgGH3.1) because it was differentially expressed, compared to the two controls, Inflo.2C/NA1 and Inflo.2C/NA2. IAA-induced protein ARG7-like (EgARG7) and tryptophan aminotransferase-related protein 3-like (EgTAA3) were differentially expressed compared with one control. Inflo.2C/NA1 and Inflo.4C/NA, respectively.

Expression confirmation of the candidate genes in the auxin signaling pathway by reverse transcription polymerase chain reaction and quantitative real-time RT-PCR

We performed RT-PCR and qRT-PCR experiments confirming RNA seq results by focusing on the auxin signaling genes. The sequences of the six candidate genes that were used to design primers for RT-PCR and qRT-PCR are available on the NCBI database. The accession numbers and details of primers that were used to amplify these six genes are listed in Table S4. Confirming by RT-PCR and qRT-PCR, five candidate genes were found to be clearly differentially expressed between treated and control samples, as shown in Figs. 3–6. EgCyp2 was used as the reference gene with an expected size of 163 bp. These candidate genes included probable IAA-amido synthetase GH3.8 (EgGH3.8), IAA-amido synthetase GH3.1 (EgGH3.1), IAA induced ARG7 like (EgARG7), tryptophan amino transferase-related protein 3-like (EgTAA3) and flavin-containing monooxygenase 1 (EgFMO1).

The first RT-PCR comparison was between Inflo.1T/WA and Inflo.1C/NA, by searching for differences in EgFMO1 and EgSAUR71 amplifications. The EgFMO1 expression was not clearly different between the treatment and controls (Fig. 3A). The amplification of EgSAUR71 by the two given primers (Table S4) failed. However, after confirming by qRT-PCR, the EgFMO1 expression was found more in the control than it was in the treatment, as shown in Fig. 4. This result correlates with RNA seq that shows expression of EgFMO1 in the control to be about 2.5 times greater than it is in the treatment. The second RT-PCR and qRT-PCR comparisons were between Inflo.6T/WA and Inflo.2C/NA1, by searching for differences in EgGH3.8, EgGH3.1 and EgARG7 amplifications (Figs. 3B and 5). The EgGH3.8 expression was only found in the treatment. This result correlates with RNA seq that shows expression of EgGH3.8 in the treatment to be about 6.2 times greater than it is in the control. EgGH3.1 expression was also only found in the treatment. This result also correlates with RNA seq that shows expression of EgGH3.1 in the treatment to be about 5.8 times greater than it is in the control. EgARG7 expression in the control was more than in the treatment. This result correlates with RNA seq that shows expression of EgARG7 in the control to be about 6.9 times greater than in the treatment. The third RT-PCR and qRT-PCR comparisons were between Inflo.6T/WA and inflo.4C/NA, by searching for differences in EgTAA3 amplification (Figs. 3C and 6). The EgTAA3 expression was found more in the control than it was in the treatment. This result correlates with RNA seq that shows expression of EgTAA3 in the control to be about 4.6 times greater than in the treatment. EgCyp2 was used as the reference gene in all three comparisons.