De novo transcriptome analyses of host-fungal interactions in oil palm (Elaeis guineensis Jacq.)

Basal stem rot (BSR) is a fungal disease in oil palm (Elaeis guineensis Jacq.) which is caused by hemibiotrophic white rot fungi belonging to the Ganoderma genus. Molecular responses of oil palm to these pathogens are not well known although this information is crucial to strategize effective measures to eradicate BSR. In order to elucidate the molecular interactions between oil palm and G. boninense and its biocontrol fungus Trichoderma harzianum, we compared the root transcriptomes of untreated oil palm seedlings with those inoculated with G. boninense and T. harzianum, respectively. Differential gene expression analyses revealed that jasmonate (JA) and salicylate (SA) may act in an antagonistic manner in affecting the hormone biosynthesis, signaling, and downstream defense responses in G. boninense-treated oil palm roots. In addition, G. boninense may compete with the host to control disease symptom through the transcriptional regulation of ethylene (ET) biosynthesis, reactive oxygen species (ROS) production and scavenging. The strengthening of host cell walls and production of pathogenesis-related proteins as well as antifungal secondary metabolites in host plants, are among the important defense mechanisms deployed by oil palm against G. boninense. Meanwhile, endophytic T. harzianum was shown to improve the of nutrition status and nutrient transportation in host plants. The findings of this analysis have enhanced our understanding on the molecular interactions of G. boninense and oil palm, and also the biocontrol mechanisms involving T. harzianum, thus contributing to future formulations of better strategies for prevention and treatment of BSR.


Background
Basal stem rot (BSR) is a serious and prevalent fungal disease in oil palm that causes significant economic loss to oil palm plantations by reducing oil yield of infected palms [1,2], causing loss of stands and shortening replanting cycle of new palms [1,2]. BSR is caused by Ganoderma spp. which are able to degrade plant cell wall components including lignin, and cause white rot to oil palm [3]. In advanced infection stage, fruiting bodies may form on oil palm trunks before the rotted infected palms eventually collapse. The disease symptoms of BSR at the foliar part of oil palm include stunted and unopened spear leaves, yellowing and senescence of upper fronds, reduced and "one-sided mottling" canopy [2,4].
In recent years, fungal DNA-, biochemical-, sensorand geographic information system -based screening methods for Ganoderma spp. in oil palms have been developed [5][6][7][8] for early detection of BSR. The mitigation actions practiced by oil palm plantations include sanitation before replanting, isolation of infected or diseased palms by elimination and burning of infected palms [9], and application of fungicides [10]. The application of beneficial microbes as biocontrol agents such as Trichoderma spp., mycorrhizae and plant growth promoting bacteria have also been demonstrated to be able to alleviate BSR incidences [11][12][13].
Oil palm lines that are tolerant to Ganoderma spp. have been reported previously [14,15]. However, the genes that are responsible for their tolerance are unknown. Genes and gene products that contribute to the virulence of Ganoderma spp., and defense response of oil palm to these pathogens are not well known although this information is crucial to strategize effective measures against the devastating disease.
In recent years, transcript sequences encoding pathogenesis related (PR) proteins such as chitinases, glucanases, defensins, protease inhibitors [16][17][18][19], and proteins that are potentially involved in oil palm defense [17,20,21] have been reported. However, the information is partial and insufficient to delineate a complete picture of oil palm defense in response to Ganoderma spp. due to a limited number of transcripts and proteins being analyzed.
The main objective of this study was to elucidate the molecular responses of oil palm roots colonized by G. boninense by comparing the transcriptomes of inoculated and uninoculated oil palm seedlings with this pathogenic fungus. In addition, the fungal transcripts that were present in the transcriptome of treated oil palm roots were identified to lend insight into the infection mechanisms of G. boninense. We also compared the differentially expressed genes (DEGs) in oil palm roots upon treatments with G. boninense and T. harzianum, respectively, to shed light on different responses of oil palm roots to pathogenic and beneficial fungi, and to provide evidence supporting the biocontrol mechanisms proposed for T. harzianum against G. boninense. A comprehensive understanding of oil palm defense responses opens up opportunities for developing strategies to eradicate the disease and to enhance the disease tolerance/resistance of this important crop, noting that it is one of the major producers of edible oil and oil related products in the world.

Results
Transcriptomes of oil palm roots in response to G. boninense and T. harzianum In this experiment, we aimed to identify the genes that are involved in the molecular responses of oil palm roots to the artificial inoculation of G. boninense and T. harzianum. The G-treated oil palm root tissues (whole root sections) collected were not separated based on duration of inoculation due to the nonsynchronous nature of infection of G. boninense on the roots from the same plant at a particular time point. During the course of Ganoderma(G)-treatment, lesions and mycelia were observed on the root surface of G-treated seedlings at 6 and 12 weeks post inoculation (wpi), indicating that G. boninense has successfully colonized the roots. However, no other disease symptoms such as yellowing of leaves and formation of fruiting bodies were observed. The mRNA samples of oil palm roots from nine biological replicates in Gtreatment from 3-12 wpi, covering infection stages from recognition of pathogen to tissue necrosis, were pooled and sequenced to reveal the transcripts involved. The roots from an equal number of biological replicates of untreated oil palms were also sampled from 3-12 wpi, for mRNA isolation and sequencing. In total, 111,832 transcripts were assembled from 37,784,896 processed and filtered reads (approximately 2.79 × 10 9 bases) obtained from mRNA-seq analysis of untreated roots compared to 50,025 transcripts assembled from 14,495,512 processed and filtered reads (approximately 5.45 × 10 8 bases) which were obtained from mRNA-seq analysis of G-treated roots.
In addition, the transcriptome of oil palm inoculated with T. harzianum was also profiled. Throughout the 12 week-treatment, the spore counts of T. harzianum in the soil around the Trichoderma (T)-treated seedlings were maintained above 10 4 colonies per unit (cfu) per g of soil to ensure the treatment was effective while the spore counts in the untreated soil were estimated to be 10 2 cfu per g of soil. The mRNA samples from oil palm roots in T-treatment, from 3-12 wpi were pooled and sequenced. The processed and filtered reads (16,665,379 reads which total up to approximately 6.33 x 10 8 bases) were obtained from mRNA-seq analysis of T-treated roots and assembled into 45,968 transcripts.
The assembled sequences (transcripts) obtained from untreated, G-and T-treated oil palm roots were then clustered to reduce the number of overlapping sequences in each transcriptome for functional annotation. Approximately 50 % of the assembled sequences (unigenes) have significant matches (E values ≤10 −5 ) to the non-redundant protein database at National Center for Biotechnology Information (NCBI). About 84 % of these unigenes which have significant matches to protein encoding sequences from Viridiplantae were used for identification of differentially expressed genes (DEGs) in pairwise samples.
Biological averaging where pooled biological replicates instead of individual biological replicates was employed in DEG analysis in this study. Although the results in this study may not provide the same statistical resolution by replicated mathematical averaged experiments, biological averaging for DEG analysis is neither new nor unprecedented in both cDNA microarray and mRNA-seq [22][23][24][25][26][27][28]. Meanwhile, DEG analysis by DESeq [29] has been reported to be more conservative and predicted less number of DEGs in comparison with another method which has been developed for mRNA-seq analysis without replicates [22]. In this study, the variance was estimated by treating the samples being compared as if they were replicates of the same condition. We noted that without individual biological replicates, the statistical significance could only represent differences in G-or T-treated sample and untreated sample being compared. Nevertheless, the gene expression measurements represented averages across nine oil palm seedlings since roots from nine oil palm seedlings were pooled to generate the sequencing libraries.
The expression levels of DEGs in 21 pairwise samples were further verified by quantitative reverse-transcription (qRT)-PCR (Table 1) on the same RNA samples used for sequencing, and approximately 86 % of them support the results predicted in silico whereby most of the genes were shown to have consistent expression patterns (up-or down-regulation) by both DEG analysis and qRT-PCR. We noted differences in the magnitude of gene expression measured by these two approaches which are most probably caused by differences in sensitivity and normalization methods, however the gene expression patterns measured by both approaches are largely similar.

DEGs of oil palm roots during Ganoderma colonization
In total, 1,996 and 1,219 oil palm genes were found to be up-and down-regulated 4-fold and more in G-treated oil palm roots, respectively (Additional file 1: Table S1) compared to the untreated oil palm roots. These oil palm genes could either form part of the host defense against fungal invasion, or a consequence of the suppression of host defense system by the pathogen upon G-treatment. The up-regulated genes in oil palm treated with G. boninense were found to be enriched in 85 Gene Ontology (GO) Slim terms: 40 in biological process (BP), 33 in molecular function (MF), and 12 in cellular compartment (CC) ( Table 2). Among these GOSlim terms are those related to defense response, cellular response to salicylic acid (SA) stimulus, nitric oxide mediated signal transduction, 1-aminocyclopropane-1-carboxylate (ACC) oxidase activity, naringenin-chalcone synthase and others ( Table 2). On the other hand, the down-regulated DEGs were found to be enriched in 27 GOSlim terms (8 BP, 6 MF, 13 CC; Table 2). The enriched GOSlim terms among the downregulated DEGs in G-treated oil palm roots include those related to generation of precursor metabolites and energy, photosynthesis, electron transport chain, ribosome biogenesis, xyloglucan:xyloglucosyl transferase activity and others (Table 2).
Among the GOSlim terms that were enriched in the pool of up-regulated genes in G-treated oil palm roots,   In T-treatment only 17 of them (9 BP, 6 MF, 1 CC) were also found to be enriched among the up-regulated DEGs in T-treated oil palm roots, including those related to response to karrikin, generation of precursor metabolites and energy, translational elongation, amine transmembrane transporter activity, nitrate reductase (NADH) activity and others (Table 2). A total of 10 GOSlim terms that were enriched among the down-regulated genes in G-treated oil palm roots (4 BP, 4 MF, 2 CC) were also found to be enriched among the down-regulated DEGs in T-treated oil palm roots, including those related to ribosome biogenesis, electron transport chain, generation of precursor metabolites and energy ( Table 2). These GOSlim terms may reflect the general responses of oil palm roots to fungi irrespective of whether they are pathogenic or beneficial, and their modes of infection. True enough, 1,324 and 868 DEGs were up-and down-regulated by both G-and T-treatments while 672 and 351 DEGs were uniquely up-and down-regulated by G-treatment only ( Fig. 1). Figure 2 shows the DEGs in G-treated oil palm roots compared to the untreated oil palm roots that are involved in the defense response including signal perception and transduction, phytohormone biosynthesis and signaling, transcription factor, generation and scavenging of reactive oxygen species (ROS), production of secondary metabolites, cell wall biosynthesis and modification enzymes and others. It also compares the number of DEGs in T-treated oil palm seedlings in relative to that of the untreated oil palm seedlings. Snapshots of the DEGs in G-and T-treated oil palm roots over the main metabolic pathways and stress-related pathways are also displayed in Fig. 3.

DEGs of oil palm roots in response to T. harzianum
A total of 1,955 and 1,513 genes were found to be upand down-regulated 4-fold and more in T-treated oil palm roots, respectively, compared to the untreated oil palm roots (Additional file 1: Table S2). However, only 631 and 645 of these DEGs were found to be uniquely up-and down-regulated by T-treatment only (Fig. 1). Thirty five GOSlim terms (17 BP, 17 MF, 1 CC) were found to be enriched among the up-regulated DEGs in T-treated oil palm roots with 18 of them (8 BP, 10 MF) found to be enriched among the up-regulated DEGs in T-treated oil palm only but not among those in Gtreated oil palm roots ( Table 2). The majority of these GOSlim terms are related to transportation including amine transport, nitrogen compound transport and carboxylic acid transport (Table 2). On the other hand, the unique down-regulated DEGs in T-treated oil palm roots were found to be enriched in 17 GOSlim terms (6 BP, 9 MF, 2 CC) with only 7 of them being unique to Ttreated oil palm roots including response to chitin and anion binding ( Table 2).

Fungal transcripts in treated oil palm roots
The presence of fungal transcripts in G-treated oil palm roots allowed us to have an overview of the genes being expressed in G. boninense during its infection of host plant. The fungal transcripts in treated oil palm roots were identified among the unigenes as sequences with high identities to fungal sequences, or protein encoding sequences that showed high identities to those belonging to fungi. The number of reads that showed high identities to those of Ganoderma sp. 10597 SS1, Trichoderma harzianum CBS226.95, Trichoderma virens Gv29-8 and Trichoderma reesei were 3.1 %, 0.9 % and 2.4 % of the total mRNA-seq reads obtained from G-, T-treated and untreated oil palm roots, respectively. The number of protein encoding sequences that showed high identities to those belonging to fungi represents 13 % of the total unigenes that have significant matches to the nonredundant protein database at NCBI. The numbers of reads that contribute to these unigenes were 5.8 %, 0.4 % and 0.3 % of the mRNA-seq reads obtained from G-, Ttreated and untreated oil palm roots, respectively. Table 3 (and Additional file 1: Table S3) shows the fungal transcripts in treated and untreated oil palm roots. Among these transcripts are those encoding enzymes involved in lignin metabolism such as manganese peroxidase and laccases. Besides, there are a high number of fungal transcripts encoding glycoside hydrolases and glycosyltransferases from various families, carbohydrate binding modules, cellulases and polysaccharide lyase among the fungal transcripts in G-treated oil palm roots. On the other hand, transcripts encoding cell wall degrading/ modification enzymes were not present in the transcriptome of T-treated oil palm roots. The genes encoding phosphate transporters, sugar transporter and carbohydrate metabolism were found to be abundant among the fungal transcripts in T-treated oil palm roots (Table 3).

Discussion
What are the genes and biological processes involved in oil palm roots infected by hemibiotrophic G. boninense?
The DEGs in G-treated oil palm roots in comparison with the untreated oil palm roots enabled us to deduce the proteins and pathways involved in the defense response of host plant in response to G. boninense. In the initial stage of infection by G. boninense, recognition of pathogen by oil palm could be initiated by the binding of pathogen-associated molecular patterns (PAMPs) from G. boninense to pattern recognition receptors (PRRs) residing extracellular or at the plasma membrane. This is evident by the presence of DEGs encoding numerous receptors in the transcriptome of G-treated oil palm roots such as receptor-like kinases (RLKs) and receptor-like proteins (RLPs).
The components of fungal cell wall or cell membrane such as chitin, glucans and ergosterol have been reported as PAMPs in many pathosystems [30]. This is supported by our findings which revealed that various oil palm chitinases and glucanases were indeed up-regulated in oil palm roots in response to G-treatment. These cell wall degrading enzymes (CWDEs) may degrade fungal cell wall, releasing chitins and glucans. Alternatively, ergosterol which was reported to be important for the virulence of pathogen [30] may also elicit PAMP-triggered immunity (PTI) in oil palm. In addition, degraded plant cell wall or cell wall components (also known as damage-associated molecular patterns, DAMPs) such as oligogalacturonates may serve as endogenous elicitors that bind to PRRs [31]. Fungal transcripts encoding enzymes involved in the biosynthesis of ergosterol, and CWDEs such as polygalacturonases, laccases and many other glycosyl hydrolases were found to be in high abundance in G-treated oil palm roots (Table 3).
In addition to PTI, a rapid and robust plant response i.e. effector-triggered immunity (ETI) which is based on highly polymorphic resistance (R) proteins [32] can also be activated in plants upon the recognition of avirulence (Avr) effectors released by biotrophic pathogens. ETI often blocks further invasion of pathogens by inducing a hypersensitive response and other locally induced defence responses which confine pathogens at the infection site and limit their growth and supply of nutrients. So far, only a Toll/interleukin 1 receptor domain R protein has been reported to be associated with the resistance of Arabidopsis against the necrotrophic Leptospharia maculans [33]. Although it is premature to speculate the presence of R proteins that are associated with the hemibiotrophic Ganoderma spp. in oil palm, several oil palm genes encoding R proteins with yet to be known functions were found to be differentially regulated in this study (Fig. 2). Among them, the gene encoding RPM1-interacting protein 4 (RIN4) which is a mediator of recognition of an R protein (RPM1) and also a negative regulator of PTI [34], was found to be down-regulated in G-treated oil palm roots. However, the gene encoding RPM1 was not found to be differentially regulated in G-treated oil palm roots. Its association with the susceptibility of oil palm to Ganoderma spp. was not known.
Signaling, ion fluxes, respiratory burst and generation of reactive oxygen species (ROS) may occur in infected roots as evident by the presence of DEGs encoding putative ADP-ribosylation factor (ARF), ARF GTPase activating protein, Rac small GTPases, GTP binding protein, mitogen associated protein kinase (MAPK), calcium dependent protein kinase (CDPK), ion channels, calcium transporters, respiratory burst oxidase (Rboh) and superoxide dismutase (SOD) in G-treated roots. Plasma membrane dynamics play an important role in plant defense whereby the residing calcium and anion channels and associated ion fluxes are important signaling events that regulate host defence. Besides, Rboh or NADPH oxidase which resides in the plant plasma membrane produces superoxide and other ROS including hydrogen peroxide that can cause lipid peroxidation, cross-linking of cell wall proteins and lignification. Hydrogen peroxide can also travel into plant cytoplasm serving as signaling molecule of other defense response including ethylene (ET) biosynthesis and hypersensitive reaction (HR)-associated cell death [35].
Lipid peroxidation of plasma membrane followed by lipoxygenation of fatty acids may lead to the biosynthesis of oxylipins, including 12-oxophytodienoate (OPDA) which is the precursor of jasmonic acid (JA). A few genes encoding lipoxygenases were up-regulated in G-treated oil palm roots which may possibly lead to an increase in oxylipins, however the expression of gene encoding 13-lipoxygenase which is involved in the lipoxygenation of α-linolenic acid leading to the production of 12-OPDA was not differentially expressed in G-treated oil palm roots. Lipoxygenase-catalysed reactions can alternatively produce toxic volatile and non-volatile fatty acid-derived secondary metabolites that can directly attach to pathogens [36]. In addition, the genes encoding 12-OPTA reductase (OPR) and JA O-methyltransferase were down-and up-regulated in G-treated oil palm roots, respectively; indicating that the production of JA could be reduced while the pre-existing JA could be methylated to methyl-JA (an inactivated form of JA and a signal molecule). Taken together, it is plausible that the JA-induced defence against G. boninense in G-treated oil palm was downregulated by the pathogen.
The genes encoding coronatine-insensitive 1 (COI) was up-regulated in G-treated oil palm roots. COI is a component of SCF (COI1) E3 ubiquitin ligase complexes which play crucial roles in regulating response to JA and regulating plant gene expression during plant-pathogen interactions [37]. Mutants coi1-1 to coi1-14 were found to have enhanced resistance to Pseudomonas syringae  atropurpurea [38]. Thus, the up-regulation of this gene could possibly reduce the resistance of oil palm roots to G. boninense due to suppression of JA-induced defence response. The decrease of JA in G-treated oil palm roots is further corroborated by the down-regulation of a JAinducible gene encoding phenylalanine ammonia lyase (PAL). PAL is a key enzyme for the biosynthesis of cinnamic acid (CA) precursor which can be channeled to the biosynthetic pathways of monolignols, flavone/ flavonoid, isoflavone phytoalexin, and SA. The downregulation of an oil palm PAL gene (EgPAL) in G-treated oil palms was also supported by a previous study [21]. The accumulation of transcripts related to the methylation of SA to its activated form (methyl-SA) in Gtreated oil palm suggests the involvement of methyl-SA in hormone signaling and possibly in systemic response in neighboring uninfected cells. True enough, the genes encoding NPR1 and NPR interactor which are involved in the signal transduction pathway that leads to SAmediated systemic acquired resistance (SAR) [39], were also up-regulated in G-treated oil palm roots. The level of ET in G-treated oil palm tissues could be modulated by the expression of genes encoding ACC oxidase-like proteins. Different isoforms of ACC oxidases appear to be the principal targets whereby different gene sequences encoding this enzyme were up-and down-regulated in Gtreated oil palm roots, complicated the deduction of the role of ET in G-treated oil palm roots. There could be "tug-of-war" between the pathogen and host to regulate cell death through ET in G-treated oil palm roots. ET was demonstrated to inhibit symptom development in necrotrophic pathogen infection but enhance the cell death caused by other type of pathogen infection [40]. Since G. boninense is a hemibiotroph, the regulation of ET in oil palm could be of central importance to determine the host resistance to this pathogen.
Our results suggested that the interplay of JA, SA and ET signaling pathways may partly determine the downstream molecular response in G-treated oil palm. SA accumulation was reported to increase the resistance of host plant to hemibiotrophic fungi but promote the susceptibility to necrotrophic pathogens [41,42]. The accumulation of SA and SA signaling may negatively regulate the gene expression of ascorbate oxidase and ascorbate peroxidase which are involved in the scavenging of hydrogen peroxide [43] in G-treated oil palm roots. Hydrogen peroxide has been reported to be involved in HR which limits disease progression caused by biotroph but promotes the invasion of necrotroph. Since pathogenic G. boninense is a hemibiotroph, the modulation of cellular hydrogen peroxide in infected host cells could be critical in determining the disease tolerance/resistance of the host. In this study, a gene encoding hypersensitiveinduced reaction (HIR) protein and two genes encoding autophagy related proteins were up-regulated; while a gene encoding defender against apoptotic death was down-regulated in G-treated oil palm roots. The oil palm HIR may have similar function as the rice HIR1 (OsHIR1) which can trigger hypersensitive cell death [44]. This may partly contribute to the susceptibility of oil palm roots to the necrotrophic G. boninense.
The downstream oil palm defense responses may involve the activation of various transcription factors (TFs) mainly AP2/ERE, bZIP, Zn finger, NAC, MYB and WRKY (Figs. 2 and 3) by concerted signals. As a result, various pathogenesis-related (PR) proteins were accumulated as evident by the up-regulation of their gene expression. The putative transcripts encoding PR proteins (including glucanases, chitinases, type 2 ribosome inactivating proteins, Bowman Birk serine protease inhibitor, type 2 proteinase inhibitor proteinases, cysteine and subtilisin proteinases), cell wall modifying enzymes such as pectinesterase and its inhibitor, pectate lyase, callose synthase, cellulose synthase and laccases, were upregulated in G-treated oil palms. The production of PR proteins and cell wall modification (including callose deposition) are some of the defense mechanisms reported to be operating in host plants in response to pathogenic invasion [45]. Callose deposition causes blockage of plasmodesmata impeding cell-to cell movement of pathogens while methyl esterified cell wall pectin is less susceptible to hydrolysis of fungal polygalacturonases [46]. It is noted that many genes encoding xyloglucan endotransglucosylase/hydrolase and xyloglucan glycosyltransferases were down-and up-regulated, respectively, in the G-treated oil palm roots. Xylose was found to be able to enhance the production of fungal laccase which is a lignolytic enzyme [47]. Infected oil palm could possibly strengthen the cell wall by reducing the accumulation of xyloglucan degrading enzymes (thus reducing fungal laccase), by modifying xyloglucans in the host cell wall. The amount of lignin may also increase in G-treated oil palms due to the accumulation of transcripts related to monolignol biosynthesis especially trans-cinnamate 4monooxygenase and caffeic acid 3-O-methyltransferase. As the first line of barrier to pathogen, cell wall and its biosynthesis and modification may play an important role in oil palm in its defense against G. boninense.
The up-regulation of genes encoding isoflavone 2'-hydroxylase and isoflavone 3'-hydroxylase in G-treated oil palms suggested an accumulation of unknown isoflavones. However, the down-regulation of a gene encoding isoflavone reductase in G-treated oil palms may dampen the accumulation of phytoalexin in G-treated oil palm roots. A rice isoflavone reductase-like gene, OsIRL, was demonstrated to be positively regulated by JA and negatively regulated by SA [48]. It is possible that the fungal elicitor from G. boninense could suppress JA and the JAinduced defence in G-treated oil palms including the down-regulation of isoflavone reductase. Furthermore, a high number of transcripts encoding sesquiterpene synthases, cytochrome P450s and a few other enzymes (in less number of transcripts) that are involved in terpene, sterol and steroid biosynthesis were also differentially expressed suggesting the accumulation of biochemical compounds in G-treated oil palms. A complex regulation of terpenoid pathway may lead to the production of sesquiterpenoid phytoalexins with fungitoxicity and also phytosterols that may restrict cell membrane permeability [49]. However, the types and roles of these secondary metabolites in oil palm defense await further elucidation.
What are the expressed genes from G. boninense that are involved in colonization of oil palm roots?
As a white rot fungus, G. boninense is able to degrade and mineralize all components of plant cell walls, including the recalcitrant lignin. Transcriptome sequencing of Gtreated oil palm roots also allowed us to have an overview of the fungal genes being expressed in G-treated oil palm roots, including genes encoding enzymes involved in lignin degradation (manganese peroxidase and laccases). Besides, there are also fungal transcripts encoding many glycoside hydrolases and glycosyltransferases from various families, carbohydrate binding modules, cellulases, pectinesterases, polysaccharide lyase etc., suggesting their active roles in degrading other plant cell wall components. Fungal transcripts encoding glyoxal oxidase which is required for filamentous growth and pathogenicity in Ustilago maydis [50] were also found among the fungal transcripts in G-treated oil palm roots. Glyoxal oxidase was also found to have a role in peroxide production for diverse oxidative reactions during wood decay [51].
Fungal gene sequences that are involved in ergosterol biosynthesis were found among transcripts in G-treated oil palms exclusively. As mentioned earlier, ergosterol which is prevalent in the membrane of G. boninense may be perceived by oil palm PRP as a PAMP which triggers PTI in the host plant. Ergosterol was shown to be able to trigger differential changes in the metabolome and variation in the biosynthesis of secondary metabolites in tobacco cells [49]. It also induces expression of genes encoding PR proteins and PAL [52]. In addition, immunomodulatory protein and allergens from G. boninense may also function in triggering host defense. The presence of a fungal gene encoding necrosis-and ethyleneinducing (NEP) 1 protein in G-treated oil palm roots is of particular interest as this protein has been found to be able to induce the formation of necrotic lesions and expression of genes encoding PR proteins, as well as to trigger the production of ROS in host plants [53]. NEP1 may cause phytotoxicity to oil palm roots through cytolysis and increased membrane permeability [54].
On the other hand, the presence of fungal transcripts encoding catalases, toxin detoxification proteins, and multidrug transporters, may assist the pathogen surviving the defense mechanism of host plants. Furthermore, two different fungal transcripts encoding SA hydroxylases were found in G-treated oil palm root and untreated roots respectively. The questions as to whether these enzymes were secreted into host cells and able to degrade host SA are unknown. A few transcripts encoding fungal elicitors that have been reported in other pathosystems were discovered among the fungal transcripts, however, their roles await further investigation. All in all, the gene products mentioned above may contribute to the virulence of G. boninense.
How did the avirulent symbiont T. harzianum evade host defense?
Approximately 20 % and 37 % of the GOSlim terms that are enriched among the up-and down-regulated genes in G-treated oil palm roots, was found to be also enriched among those in oil palm roots treated with T. harzianum, respectively. This shows that a common set of oil palm genes were up-and down-regulated by both fungi irrespective of whether they are pathogenic or avirulent fungi. However, different individual genes categorized under the same GOSlim may play different roles, shaping the final response adapted to the challenging pathogens (Table 2).
In plants colonized by Trichoderma sp., the endophytic fungus was found to be limited to a few root cortical cell layers in plant roots [55]. T. harzianum was able to colonize oil palm roots without causing any lesions. Since the T-treated oil palm seedlings were asymptomatic, the host may be able to limit the colonization of Trichoderma by reinforcing cell wall and producing antimicrobial compounds and ROS [56][57][58][59][60]. This is further substantiated by the differential expression of genes encoding enzymes involved in the production of secondary metabolites (mainly alkaloids, glycoside and ketone with unknown functions), cell wall biosynthesis, cell wall modification and respiratory burst in T-treated oil palm roots. T. harzianum may remain accommodated by the host as an avirulent symbiont by suppressing plant defense through the secretion of still unknown fungal effectors [61].
The enzymes involved in the biosynthesis of ergosterol was absent in T-treated oil palm roots indicating that this fungal cell membrane sterol may not serve as a PAMP in T-treated oil palm roots that triggers the downstream defense response. Indeed, the ergosterol-inducible genes encoding sesquiterpene synthases [62] were found to be not differentially expressed in T-treated oil palm roots.
The interaction between the oil palm and T. harzianum must have been finely tuned with promising benefits to both partners [61]. The expression of three genes encoding ACC oxidases were found to be down-regulated, while another one was found to be up-regulated in T-treated oil palm roots compared to those in untreated oil palms. Since many of these genes were different from those encoding the same enzyme in G-treated oil palms, it is obvious that these genes were elegantly regulated by different signals. The genes encoding OPR and COI were down-and up-regulated in T-treated oil palm roots, respectively, indicating the suppression of JA biosynthesis and JA-mediated defense response in T-treated oil palm roots. It is not surprising since T. harzianum may not be causing necrosis or wounding to oil palm roots. The up-regulation of three genes encoding NPR1-like proteins suggests that an active SA signaling pathway is present in T-treated oil palms. The cross talks of hormones may decide the final defense response in the host, allowing the survival of T. harzianum but effectively containing its spread in plant host.
Although the genes encoding Rbohs were up-regulated and SA signaling could negatively regulate the gene expression of enzymes which are involved in the scavenging of hydrogen peroxide in T-treated oil palm roots (Fig. 2), the gene expression of two genes encoding defender against apoptotic death were up-regulated, suggesting that their gene products may serve as negative regulators of cell death [63] in T-treated oil palm roots. It is logical that T. harzianum could avoid cell death in host cells as one of its survival strategies in host plants.
Unlike G-treated oil palm roots, the gene encoding polygalacturonase was down-regulated in T-treated oil palm roots, implicating differences in cell wall modification mechanisms and possibly also the types of DAMP induced by the two fungi. On the other hand, fewer fungal genes encoding cell wall degrading/modification enzymes were expressed in the transcriptome of T-treated oil palm roots. Transcripts encoding laccases and peroxidases were not found among the fungal genes in T-treated oil palm roots. It is possibly important for the fungus to avoid from generating DAMPs that may trigger host defense system in its attempt to establish a symbiotic (or at least a neutral) relationship with its host plant. In addition, T. harzianum may avoid over accumulation of SA in the host plants by producing SA hydroxylase coinciding with the presence of a transcript encoding this enzyme among the fungal transcripts in the T-treated oil palm seedlings.
What are the fungal genes from T. harzianum that are involved in biological control of G. boninense?
Several fungal elicitors from T. virens were reported to be able to activate plant basal immunity such as the ET-inducing xylanase [64], the proteinaceous nonenzymatic elicitor Sm1 [65,66] and the 18 mer peptaibols [67]. Although these genes were not found among the genes in T-treated oil palm roots, we do not exclude the possibility that they are regulated by other mechanisms at translation or post-translational levels. The biological control mechanisms by Trichoderma spp. could also be achieved by direct effect of Trichoderma on pathogen through competition for space and nutrients, parasitisation of other fungi by producing antimicrobial compounds and antifungal enzymes (reviewed in [68]). Indeed, the genes encoding for phosphate transporters, hexose transporter, monosaccharide transporters and carbohydrate metabolism were found to be abundant among the genes in T-treated oil palm roots, in favor of better competition for nutrients. Fungal genes from T. harzianum that are potentially involved in parasitisation of G. boninense were unreported while those involved in the production of antimicrobial compounds and antifungal enzymes were not identified. Our results demonstrated that T. harzianum is more likely to strengthen the host plants by improving the nutrient status of host plants.

Conclusions
Investigation on the transcriptomes of oil palm roots in response to G. boninense has enabled us to dissect the defense mechanisms involved in the establishment of BSR: 1. Ergosterol, fungal cell wall components, host cell wall components (degraded by fungal CWDEs) may elicit the defense responses of oil palm roots; 2. JA biosynthesis, signaling and JA-mediated defense response may be suppressed while the SA biosynthesis, signaling and SAmediated defense could be activated in G-treated oil palm roots by the pathogen; 3. G. boninense and the host plant may compete to control the establishment of disease symptom through the transcriptional regulation of ET biosynthesis, hydrogen peroxide production and scavenging, apoptosis and autophagy; 4. Colonized oil palm roots may produce PR proteins and secondary metabolites that have antifungal activities against G. boninense; and 5. Colonized oil palm may strengthen its cell wall through cell wall modifying enzymes as its defense response.
The endophytic T. harzianum could possibly evade the defense response of oil palm roots by suppressing both JA and ET biosynthesis, signaling and defense responses and avoid establishment of cell-death. However, oil palm could limit the spread of this avirulent symbiotic fungus by producing tailored made secondary metabolites from different groups. Analyses of the transcriptome of oil palm roots colonized with T. harzianum shed light on the possible biocontrol mechanisms employed by T. harzianum against G. boninense which include improved nutrition status of host plants through mobilization of nutrients. The findings of this analysis enhance our understanding on the molecular interactions of G. boninense and oil palm, and also the biocontrol mechanisms involving T. harzianum, thus contributing to future formulations of better strategies to eradicate BSR.

Preparation of Ganoderma boninense-inoculated rubber wood blocks
Rubber wood blocks (6 cm x 6 cm x 12 cm each) were sterilized at 121°C for 30 min before they were added with 100 ml potato sucrose agar (PSA) per rubber wood block and autoclaved. The autoclaved wood blocks were inoculated with Ganoderma boninense PER71 (Malaysian Palm Oil Board, Bangi, Malaysia) cultures that were grown on PSA at 28°C for 7 d, and incubated in the dark at room temperature for 8 weeks. The fungal strain selected for this study is a reference strain which has been reported to be pathogenic to oil palm [12,13,15].

Preparation of surface mulch and conidial suspension from Trichoderma harzianum
For the preparation of conidial suspension, Trichoderma harzianum Rifai strain T32 (Faculty of Science, Universiti Putra Malaysia, Serdang, Malaysia) was grown on potato dextrose agar (PDA; BD, France) at 25°C for 7 d. The spores were then harvested and the concentration of conidia was determined using a haemocytometer (Neubauer, Marienfeld, Germany). To prepare surface mulch from T. harzianum, approximately 150 g palm pressed fibers (PPF) of empty fruit bunch (EFB) were rinsed with water and autoclaved at 121°C for 45 min in a plastic bag. The autoclaved PPF and EFB were inoculated with 20 ml of T. harzianum conidial suspension containing 1-9 x 10 8 spores/ml, sealed and incubated at 25°C for two weeks in the dark.

Treatments of oil palm seedlings
A total of 27 five-month-old oil palm Elaeis guineensis GH 500 Series (Dura x Pisifera) seedlings were purchased from Sime Darby Seeds and Agricultural Services Sdn. Bhd. (Banting, Malaysia) for the following treatments: 1. artificial inoculation with G. boninense (Gtreatment); 2. artificial inoculation with T. harzianum (T-treatment); and as controls (untreated oil palm seedlings). The genotypes of these oil palm seedlings have not been reported to be tolerant against Ganoderma spp. In fact, previous studies have found them to be susceptible to G. boninense [21].
For oil palm seedlings in the G-treatment, the roots of each oil palm seedling were placed in contact with a rubber wood block inoculated with G. boninense and covered with a mix of top soil and sand in a ratio of 2:1. For oil palm seedlings in the T-treatment, 300 g of T. harzianum surface mulch was placed at the basal part of each seedling, and 500 ml conidial suspension of T. harzianum was applied directly to each seedling every fortnight. The seedlings were arranged in a complete randomized design in a green house. Three seedlings were harvested from respective treatments at 3, 6 and 12 week post inoculation (wpi) with three untreated oil palm seedlings as controls. The whole roots were collected for analysis.

RNA preparation
Total RNA from the roots of oil palm was isolated using a modified CTAB method as described by Wang et al. [69] and treated with DNase I (New England Biolabs, Hitchin, UK). The quality and integrity of RNA were evaluated using a spectrophotometer (Eppendorf, Hamburg, Germany) and 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The RNA samples of nine biological replicates from each treatment or control were pooled for RNA sequencing and qRT-PCR.
Massive RNA sequencing and de novo assembly of RNA-seq RNA sequencing was conducted on RNA samples from untreated, G-, and T-treated oil palm roots, with the Illumina technology (Illumina, San Diego, CA, USA). The mRNA-seq data have been deposited at European Nucleotide Archive (ENA) under the accession number PRJEB7252. The short reads obtained from the roots were converted to FASTQ format and analyzed with Fastxtool (http://hannonlab.cshl.edu/fastx_toolkit/) for quality check with a minimum Phred value of 20. To remove the sequences that have high identities to fungal sequences, the sequence reads were mapped to the genome sequences of Ganoderma sp. 10597 SS1, T. harzianum CBS226.95, T. virens Gv29-8 and T. reesei (downloaded from Department of Energy Joint Genome Institutes; JGI; http://genome.jgi-psf.org/), with Bowtie 0. 12.7 (Langmead et al. 2009). After the removal of probable fungal sequences, de novo assembly of RNAseq was performed with Velvet (version 1.1.05, [70]) and Oases assembler (version 0.1.22, [71]) at an optimized k-mer size for each data set. The assembled sequences from the root samples were then clustered with TGICL [72] and CAP3 [73] with a minimum overlap length of 100 bp and an identity at 97 % to form unigenes.

Annotation of trancripts
The unigenes were translated into six reading frames and compared with non-redundant (nr) database at the National Center for Biotechnology Information (NCBI) with BLASTX [74]. Matches with E-values equal or less than 10 −5 were treated as 'significant matches'. Sequences with no hits or matches with E-values more than 10 −5 were classified as 'non-significant matches'. Unigenes with their most significant matches originating from bacteria, fungi, animal and other non-plant organisms, were separated from plant sequences. The unigenes of plant and fungal origins were mapped to Gene Ontology (GO) with Blast2GO [75] (www.blast2go.com/). GOSlim terms were assigned to each unigene.

Identification of differentially expressed genes (DEGs)
The sequence reads from each sample (without the fungal sequences) were mapped to the unigenes originating from plants with Tophat 2.02 [76]. The number of reads that mapped to each assembled sequence was calculated using HTseq-count (http://www-huber.embl.de/ users/anders/HTSeq/). DESeq [29] was used to identify differentially expressed genes (DEGs) between two samples. The variance was estimated by treating the samples being compared as if they were replicates of the same condition [29]. Genes that have a log 2 [Number of reads in treated oil palm roots/Number of reads in untreated oil palm roots] ≥ |2.0| (corresponding to 4fold or more for up-down-regulation, respectively) with a P-value ≤ 0.05, were defined as significant DEGs in the two samples being compared. These oil palm DEGs were matched to Arabidopsis coding sequences (CDS) TAIR10 build 29 (ftp://ftp.ensemblgenomes.org/pub/plants/release-29/fasta/arabidopsis_thaliana/cds/) by tblastx with a maximum E-values equal or less than 10 −5 , and displayed by MapMan software (http://mapman.gabipd.org/) [77].

Enrichment of GO Slim terms among DEGs
The enrichment of GO Slim terms among the DEGs was tested with Fisher's Exact Test with Multiple Testing Correction of false discovery rate (FDR) [78]. GOSlim term with a FDR which is equal or less than 0.05 was considered to be significantly "enriched". The enriched GO terms were summarized by REVIGO [79] by removing redundant GO terms and presented by Circos software [80].

Verification of DEGs
Verification of DEGs was conducted using the same plant materials and RNA samples used for mRNA-seq. RNA treated with DNaseI (1.2 μg) was reversetranscribed with Affinity Script QPCR cDNA Synthesis Kit (Stratagene, La Jolla, CA, USA) according to the instructions of the manufacturer. The genes and primers used for real-time qRT-PCR are listed in Additional file1: Table S4. Real-time quantitative reverse transcription (qRT)-PCR was performed in BIO-RAD iQ5 iCycler Thermal Cycler (Bio-Rad, Berkeley, CA, USA) in quadruplicates. Each PCR reaction (20 μl) contained 200 ng of the first strand cDNA, 200 nM forward and reverse primers, 30nM reference dye, and 1× master mix from Brilliant III Ultra-Fast SYBR® Green QPCR Master Mix (Stratagene, La Jolla, CA, USA). The PCR condition was 95°C for 3 min; 40 cycles of 95°C for 5 s and 60°C for 10 s. Two endogenous controls encoding actin and cyclophilin from oil palm were amplified in parallel for normalization and quantification. The amplification efficiency was estimated according to the equation: E=[(10 -1/y )-1]×100 [81]. The relative transcript levels were calculated based on a comparative C T method using multiple reference genes for normalization [82].

Availability of supporting data
The raw reads for this project have been deposited in the European Nucleotide Archive (ENA) under the accession number PRJEB7252.