Pmeg isolates Gh-ER1334 and Ca-ZTHO145, and Ppal isolates Gh-ER1349 and IC-SBR112.9 used in this study were isolated from BPR infected cacao in Ghana, Cameroon, and Côte d'Ivoire as described previously (Ali et al., 2016) and maintained on CV8 agar plate at 18°C. Zoospore inoculum was produced for various infection studies essentially as described by Lawrence (1978) and modified by Ali et al. (2016).

To identify the cacao responses to Pmeg and Ppal infection, various plant inoculation studies were carried out. For all the pod husk and leaf disc inoculation assays, Catongo cacao trees susceptible to Phytophthora; Crouzillat et al., 2000) were maintained in a USDA-ARS, Beltsville, MD greenhouse at 60% relative humidity and a photoperiod of 12 h light at 25 ± 3°C and 12 h dark at 21 ± 3°C. Ambient light was supplemented with 400 W high-pressure sodium lamps, as needed, to obtain 250 μmol.m⁻²s⁻¹. Automatic retractable shade cloths were used to limit light to a maximum of 1000 μmol.m⁻²s⁻¹. Plants were drip-irrigated two times daily and fertilized regularly with a controlled-release complete fertilizer (Nutricote Total 18-6-8 Type 80, Arysta Life Science North America Corp., NC) at a rate of 670 g.m⁻² every 12 weeks to maintain adequate levels of soil moisture and nutrition.

Unwounded pod husk inoculation assays were carried out to assess the virulence of Ppal (isolate IC-SBR112.9 and Gh-ER1349) and Pmeg (isolate Ca-ZTHO145 and Gh-ER1334) essentially as described by Ali et al. (2016). Briefly, 4–5 month old pods were cut into 2.5 × 2.5 cm pieces; the inner core materials removed and dipped in sterile 0.7 mM benzimidazole (senescence retardant) solution to prevent browning until used. Husk pieces were surface-sterilized with 6% (vv⁻¹) regular bleach (Clorox, USA) for 90 s followed by three rinses with sterile distilled water. Pieces were placed in sterile plastic containers (20 × 10 × 6 cm) lined with sterile tissue paper soaked in 0.7 mM benzimidazole solution at the bottom. For inoculation, 40 μl zoospore solutions were placed in the middle of the exterior part of each husk piece. Control husk pieces were treated with 40 μl sterile water. Containers were covered and incubated at 25°C with 50% RH and under 12 h light (200 lx) and dark cycles. Samples were photographed at 24 h intervals. The area under necrosis was quantified using Image J software (Abràmoff et al., 2004) and expressed as cm². Observations were taken up to 4 days and the area under disease progress curve (AUDPC) was calculated (Shaner and Finney, 1977). (1)AUDPC=∑​i = 1n[(Yi + 1+Yi)−2][ti + 1−ti)] Where, Y_i is an assessment of a disease progression (cm²) at the ith observation, t_i is time (day) at the ith observation, and n is the total number of observations. Each experiment was repeated independently twice with three replicating husk pieces per isolate per experiment.

For wounded pod husk inoculation assays, husk pieces were treated and placed in sterile plastic containers as mentioned above. Husk surface was punctured with a no. 1 cork borer (5 mm diameter) and 2–3 mm deep holes were created which were subsequently filled with 40 μl zoospore solutions of Ppal (isolate IC-SBR112.9 and Gh-ER1349) and Pmeg (isolate Ca-ZTHO145 and Gh-ER1334). Controls were filled with 40 μl sterile water. Containers were covered and incubated at 25°C with 50% RH and under 12 h light (200 lx) and dark cycles. Samples were photographed at 1 day intervals and the progression of necrotic area was quantified using Image J software (Abràmoff et al., 2004) and expressed as cm². Observations were taken up to 3 days and the AUDPC was calculated as mentioned above. Each experiment was repeated independently twice with three replicating husk pieces per isolate per experiment.

The secondary inoculums loads were estimated based on the surface sporangia count of the infected pod husk pieces. At 5 days post inoculation (dpi), each infected pod husk piece was washed with 1 ml Tween-20 solution (2% v/v). After serial dilution, sporangia solutions were plated on 20% V8-PARP agar (Ferguson and Jeffers, 1999). From each replicating husk pieces, 24 single colonies were picked and transferred to a 24-well plate containing 2 ml liquid CV8 and grown at room temperature for 4–5 days. Mycelia from each well were picked with sterile tooth picks and transferred to 2 ml safe lock tube (Eppendorf, USA), followed by flash freezing in liquid nitrogen and freeze drying. DNA from each sample was extracted as mentioned before by Ali et al. (2016) and subjected to a Pmeg and Ppal specific PCR assay as described by Ali et al. (2016) to determine the % of secondary inoculums.

Stage 3 cacao leaves (light green) (Bailey et al., 2005) were harvested from Catongo cacao trees. Leaf discs with 1.5 cm diameter were cut out with No. 24 cork borer and placed on 90 mm Petri plate lined with Whatman no. 2 filter paper soaked with 0.7 mM benzimidazole solution. For inoculation, a 20 μl zoospore suspension of Ppal (isolate Gh-ER1349) and Pmeg (Gh-ER1334) was placed at the center of the leaf disc. Controls were treated with water alone. Petri plates were covered and incubated at 25°C with 50% RH and under 12 h light (200 lx) and dark cycles.

Infected pod husk and leaf discs were harvested at 1, 2, and 3 dpi and flash frozen in liquid nitrogen. For RNA isolation, samples were ground finely in mortar and pestle and transferred to a 50 mL centrifuge tube containing 15 mL of 65°C extraction buffer (Bailey et al., 2005). Additional extraction methods were conducted as previously described (Bailey et al., 2013).

RNA was extracted from unwounded pod husk inoculated with Ppal isolate Gh-ER1349 and Pmeg isolate Gh-ER1334. RNA-Seq analysis was carried out by the National Center for Genome Resources (Santa Fe, NM, USA). cDNA was generated using the RNA library preparation TruSeq protocol developed by Illumina Technologies (San Diego, CA). Using the kit, mRNA was first isolated from total RNA by performing a polyA selection step, followed by construction of single end sequencing libraries with an insert size of 160 bp. Single-end sequencing was performed using the Illumina HiSeq2000 platform. Samples were multiplexed with unique six-mer barcodes generating filtered (for Illumina adapters/primers and PhiX contamination) 1 × 50 bp reads. The sequences acquired by RNA-Seq were verified by comparison to the Pmeg and Ppal genomes assembled in this study. RNA reads from RNA-Seq libraries ranging from 50 to 70 million reads in fastq format were aligned using memory-efficient short-read aligner Bowtie-2–2.1.0 (Langmead and Salzberg, 2012) to the coding sequences (CDS) of the cacao genome (Argout et al., 2010). Tabulated raw counts of reads to each CDS, respectively, were obtained from the bowtie alignment. Estimation of normalized RPKM (Reads Per Kilobase of transcript per Million mapped reads) and statistical analysis of expression level using count data of each gene with three replicates for each library was performed using the DEseq package (Anders, 2010) and R x64 2.15.2 program (http://www.r-project.org/). Genes were grouped into four categories, either up or down regulated (≥2-fold) for each species. KEGG (Kyoto Encyclopedia of Genes and Genomes) Automatic Annotation Server (KAAS) was used to obtain KEGG Orthology and KEGG pathways involving the differentially expressed genes by BLAST comparisons against the manually curated KEGG GENES database (Moriya et al., 2007). Gene ontology (GO) analysis was carried out using the program Blast2GO (Conesa et al., 2005).

Procedures for RT-qPCR and analysis were as described by Bailey et al. (2006). Pathogen load in each sample was calculated by comparing expression of Pmeg/Ppal reference genes (PmTP/PpTP) to expression of three cacao reference genes (TcRefs) (TcACT, TcGAPDH, and TcACP) (See Supporting Excel File S1, Sheet A). The geometric mean (C_T•TcRef) (Vandesompele et al., 2002; Ali et al., 2014) was calculated for cycle times of TcRefs (C_{T•TcRef(1−3)}] and the ΔC_T was calculated for PmTP/PpTP (C_{T•PmTP/PpTP}–C_T•TcRefs).

A total of 36 cacao genes were chosen for analysis by RT-qPCR. Nineteen cacao genes monitoring plant response during Pmeg/Ppal infection (See Supporting Excel File S1, Sheet A) were chosen based on their putative function and to represent a range of differential expression values between Pmeg, Ppal, and control treatments as observed in the RNA-Seq analysis to verify the RNA-Seq analysis. Another set of 17 genes was selected based on their putative participation in the plant hormone signal transduction pathway and plant-pathogen interaction pathway based on KEGG analysis. All 36 genes were analyzed across the three replicating sets of RNA from the Ppal isolate Gh-ER1349 and Pmeg isolate Gh-ER1334 infected and untreated control samples and analyzed for temporal and tissue specific transcriptions in response to Pmeg/Ppal infection in leaves and wounded and unwounded surfaces of pod husk pieces. The treatments include unwounded pod husks at 1 and 3 dpi, wounded pod husks at 1 dpi, and leaf disc at 1 and 2 dpi. Each treatment had three replicated samples for untreated controls, Pmeg, and Ppal infection. To obtain the relative transcript levels, the threshold cycle (C_T) values for all genes of interest (C_T•GOI) were normalized to the geometric mean (Vandesompele et al., 2002; Ali et al., 2014) C_T value of the three TcRefs as ΔC_T = (C_T•GOI) − (C_T•TcRefs). Expression levels for each cacao gene in each sample were calculated as % cacao Ref = 100^*[(E)^ΔCT] where E equals the primer efficiency.

Tabulated raw counts of reads to each CDS were obtained from the bowtie alignment. Read counts were normalized using the DEseq package (Anders, 2010) and R x64 2.15.2 program (http://www.r-project.org/). Normalized RNA-Seq count data were log₁₀ transformed and neural network self-organizing map (SOM) analysis (Kohonen, 1982, 1990) was employed (using MATLAB Neural Network Toolbox™) to objectively provide an unsupervised multi-variate and non-linear classification of the transcriptome into 24 different classes. The number of classes was chosen based on the number of treatments (control, Pmeg, and Ppal infection), number of samples, and allowing for differences in gene expression levels. The classification was done using a multivariate input vector describing the transcription response of each gene. GenAlEx 6.5 was then used to perform Principal Coordinates Analysis (PCoA) based on the covariance matrix with data standardization of log₁₀ transformed RNA-Seq count data. Results from 24 of the most highly expressed genes from each of the 24 SOM classes were used to generate a 3-D scatter plot with an Excel macro add-in.

To establish the similarity between the two publicly available cacao genomes, the Criollo cacao genome (Argout et al., 2010) was subjected to bidirectional BLAST using the BLASTn search against the Forastero cacao genome (Motamayor et al., 2013).

The data sets for the area under necrosis progression and pathogen load were normally distributed, as determined using the Shapiro–Wilk Test (Shapiro and Wilk, 1965) within the Statistical Package for the Social Sciences (IBM SPSS Statistics 22). The homogeneity of data sets across replicate experiments was confirmed within SPSS by two-tailed correlation analysis using mean data values (Pearson product moment; r ≥ 0.81; P = 0.01; Snedecor and Cochran, 1980). Therefore, data sets from the replicate experiments were pooled for the purposes of further statistical analysis. The significance of treatment effects was analyzed within SPSS by one-way ANOVA with Post-Hoc pair wise Least Significance Difference (LSD) comparisons (P = 0.05; Snedecor and Cochran, 1980).

Fisher's Exact Test between the four categories of genes, either up or down regulated (≥2-fold) for each species was carried out using the program Blast2GO (Conesa et al., 2005).

Relative expression data of cacao gene expression as % TcRefs was LOG-transformed to linearize the data (Rieu and Powers, 2009). Two-way RM ANOVA with Tukey's multiple comparisons test (P = 0.05) was carried out using GraphPad Prism version 7.0. A heat map of relative expression of the individual cacao genes was created using CIMminer (http://discover.nci.nih.gov/cimminer).

When zoospores of Pmeg and Ppal were applied to unwounded surfaces of cacao pod husks, the disease development, as measured by AUDPC, was dependent on the timing of the individual species isolate inoculated (Figure 1). When applied alone, the two Ppal isolates (IC-SBR112.9 and Gh-ER1349) had lower AUDPC values than the 2 Pmeg isolates (Ca-ZTHO145 and Gh-ER1334). When the Ppal zoospores were simultaneously inoculated along with Pmeg zoospores, the observed disease reactions were that of Pmeg (higher AUDPC values) and Pmeg was most often recovered from the sporangia produced at 5 dpi. When Pmeg zoospores were applied 24 h before Ppal zoospores, Pmeg was again recovered and the disease reaction was again that of Pmeg. Even when Ppal zoospores were applied to pod pieces 24 h before the Pmeg, Pmeg sporangia were recovered but the disease reactions were similar to that of the Ppal isolate included (a lower AUDPC resulted).

When zoospores of Ppal (IC-SBR112.9 and Gh-ER1349) and Pmeg (Ca-ZTHO145 and Gh-ER1334) were applied to wounded surfaces of cacao pod husk pieces, the disease development (AUDPC) and species specific recovery of sporangia was again dependent on the timing of the individual species isolates inoculated but results differed compared to those observed for pod pieces with unwounded surfaces (Figure 2). When applied alone, the Ppal isolates had AUDPC values similar to that of the Pmeg isolates and higher than observed on unwounded surfaces. When Ppal zoospores were co-inoculated with Pmeg zoospores onto wounded surfaces, Ppal was most often recovered from the sporangia produced. When Ppal zoospores were applied 24 h before the Pmeg zoospores were applied, again Ppal was most often recovered. When Pmeg zoospores were applied 24 h before Ppal zoospores, Pmeg was recovered.

Pathogen load was estimated as the percent of the expression levels of Pmeg/Ppal reference genes (PmTP/PpTP) relative to the expression of three cacao reference genes (detailed in experimental procedures). Pathogen load showed significant temporal variation for all three tissue types, leaf disks, wounded pod, and unwounded pod husk (Figure 3). There was a significant difference between Pmeg and Ppal load in infected leaf discs (p < 0.05) with Pmeg biomass being higher at both 1 and 2 dpi, but there was no significant difference in Pmeg and Ppal biomass (p > 0.05) in wounded or unwounded pod husks (Figure 3).

The RNA-Seq analysis approach used Pmeg (isolate Gh-ER1334) and Ppal (isolate Gh-ER1349) infected pod husk samples (3 dpi) from the susceptible clone Catongo. The pathogen loads were around 15% (% of PmTP/PpTP expression relative to TcRefs expression), a level of pathogen expected to give reasonable detection of Pmeg/Ppal transcripts based on previous experiments (Ali et al., in press). The RNA-Seq analysis identified 19,022, 19,392, and 19,393 possible cacao transcripts (with >10 RPKM) for the libraries of untreated control, Pmeg and Ppal infected samples respectively (see Supporting Excel File S2, Sheet A). After normalization, the base mean number of reads for each cacao reference gene varied among samples as follows: TcACP1 (4,459–4,847 reads), TcGAPDH (2,104–3,980 reads), and TcACT (63,879–80,794 reads). The geometric mean range for the three cacao reference genes as calculated for each clone varied from 8,573 to 11,227 reads, a difference of 23% or less between treatments after normalization.

There were 4,482 and 5,264 differentially expressed (>2-fold) transcripts between untreated control and Pmeg and Ppal infected samples, respectively, based on a p ≤ 0.05 (Benjamini and Hochberg, 1995). There were 3,181 and 3,521 transcripts induced and 1,301 and 1,743 transcripts repressed in Pmeg and Ppal compared to the untreated control, respectively (see Supporting Excel File S2, Sheet A). There were 463 transcripts induced and 580 transcripts repressed in response to Pmeg infection compared to Ppal infection although in most instances the trends were the same for both species with only the relative fold change varying. Among all the differential expression transcripts, 3,898 transcripts encode proteins with putative functions. For a better identification of the transcripts we established the similarity between two publicly available cacao genomes (Argout et al., 2010; Motamayor et al., 2013). There are 38,609 homologous genes (E < 1e-10) shared between the Criollo and Forastero cacao genomes (see Supporting Excel File S2, Sheet B).

To compare global cacao gene expression profiles in response to the two Phytophthora species, we used neural network software to produce SOM (Kohonen, 1982, 1990) from the RNA-Seq data. This enabled us to identify common patterns of gene expression in response to Pmeg/Ppal infection without matching genes one-to-one between treatments. Cacao transcripts (with >10 RPKM) for the libraries of untreated control, Pmeg, and Ppal infected samples were classified into 24 multi-variate and non-linear SOM classes (Figure 4A and see Supporting Excel File S1, Sheet B). Among the 24 SOM classes generated, classes 1, 2, 3, 7, 8, 13, 14, and 19 were highly expressed of which, classes 14 and 19 largely showed constitutive expression across all the RNA-Seq libraries with expression level being the primary divider. More than 80% of the genes belonging to classes 2, 3, 4, 5, 6, 8, and 9 were significantly induced in response to Pmeg/Ppal infection, whereas, between 20 and 30% of the genes belonging to classes 1, 6, 15, 16, and 17 were also significantly induced in response to Pmeg/Ppal infection (see Supporting Excel file S1, Sheet B). Genes repressed in response to Pmeg/Ppal infection were mostly grouped in class 12, 18, 22, 23, and 24. PCoA was carried out using log₁₀ transformed RNA-Seq count data of the 20 most highly expressed GeneMs from each of the 24 SOM classes. A 3D scatter plot PCoA comparison of coordinates 1, 2, and 3 is presented (Figure 4B). These coordinates account for 54.46, 15.37, and 11.26% of the data variation, respectively. As expected, genes from each of the SOM classes clustered separately, according to infection specificity and transcript level.

Gene ontology (GO) analysis using Fisher's Exact Test of the 4 categories of differentially expressed genes (>2-fold) responding to Pmeg/Ppal showed substantial differences (Figure 5 and Supporting Excel File S1, Sheet C). Genes belonging to GO classes corresponding to signaling, ion binding, hydrolase activity, cell wall, and signal transduction were more commonly observed in the induced categories than in repressed categories (FDR < 0.05 along with P < 0.05) in response to Pmeg/Ppal infection. Whereas, genes belonging to GO classes photosynthesis, cell division, and cell cycle were more commonly observed in the repressed categories than in the induced categories (FDR < 0.05 along with P < 0.05) in response to Pmeg/Ppal infection (Figure 5). Often genes belonging to defense associated GO classes (for example defense response, response to wounding, oxidative stress and biotic stimulus) were induced in response to Pmeg/Ppal infection but their GO classes did not meet the required level of significance (Figure 5). Most importantly, the Fisher's Exact Test in the blast2go suite indicated there were no differences for any GO class when directly comparing cacao gene categories induced by Pmeg/Ppal infection or repressed by Pmeg/Ppal infection (Figure 5 and Supporting Excel File S1, Sheet D).

Initially KEGG analysis was carried out on genes with the 8 highly expressed SOM classes (SOM class 1, 2, 3, 7, 8, 13, 14, and 19; see Figure 4) previously identified. Through this processes, candidate genes were identified for members of many pathways, often completing or nearly completing the pathways (see Supporting Excel File S1, Sheet E). Examples of particular relevance to this study include key pathways involved in the plant defense response including alpha-linolenic acid metabolism (jasmonate biosynthesis; ko00592), cysteine and methionine metabolism (ethylene biosynthesis; ko00270), plant-pathogen interaction (ko04626), endocytosis (ko04144), plant hormone signal transduction (ko04075), plant secondary metabolites including phenylpropanoid biosynthesis (ko00940), phenylalanine, tyrosine, and tryptophan biosynthesis (ko00400), and terpenoid backbone biosynthesis (ko00900). KEGG pathways analysis was then carried out incorporating 4 individual differentially expressed gene sets (>2-fold up or down regulated comparing Pmeg/Ppal infected to control), identifying 149 different metabolic pathways (see Supporting Excel File S1, Sheet F). The top 10 KEGG pathways containing the largest numbers of induced genes in response to Pmeg/Ppal infection were biosynthesis of amino acids, carbon metabolism, plant hormone signal transduction, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, endocytosis, plant-pathogen interaction, glycolysis/gluconeogenesis, protein processing in endoplasmic reticulum and carbon fixation in photosynthetic organisms (see Supporting Excel File S1, Sheet F). The top 10 KEGG pathways containing the largest numbers of repressed genes in response to Pmeg/Ppal infection were photosynthesis, cell cycle, starch and sucrose metabolism, carbon metabolism, biosynthesis of amino acids, plant hormone signal transduction, photosynthesis—antenna proteins, porphyrin, and chlorophyll metabolism, DNA replication and carbon fixation in photosynthetic organisms (see Supporting Excel File S1, Sheet F). Since both Pmeg and Ppal infection showed induction and repression of similar pathways, we carried out KEGG pathways analysis of 2 additional gene sets, induced by either or both Pmeg and Ppal infection and repressed by either or both Pmeg and Ppal infection (see Supporting Excel File S1, Sheet G). As expected, shifts in expression were observed in key pathways involved in the plant defense response mentioned above. There were 16 KEGG pathway ko04626 (plant-pathogen interaction) gene members showing more than two-fold induction in response to Pmeg/Ppal in infection (See Supplementary Table S1). The cacao genes induced include putative pathway recognition/signal components chitin elicitor receptor kinase 1 (CERK1), cyclic nucleotide gated channel (CNGC), and LRR receptor-like serine/threonine-protein kinase FLS2 (FLS2). In addition, at least one putative BAK1 (brassinosteroid insensitive 1-associated receptor kinase 1) candidate gene showed induction although less than the two-fold threshold used prior to KEGG analysis. Seventeen gene members associated with endocytosis (ko04144) showed up-regulation, with 7 of these being associated with late endosomes (See Supplementary Table S2). In addition, the pathway leading to jasmonic acid (JA) biosynthesis (ko00592 alpha-Linolenic acid metabolism) has 11 pathway steps with genes represented as induced in response to Pmeg/Ppal infection, again a near complete pathway (See Supplementary Table S1). Twenty two candidate members of KEGG pathway ko04075 (plant hormone signal transduction) showed induction in response to Pmeg/Ppal infection (See Supplementary Table S1). Multiple members of four primary genes in the ethylene biosynthesis pathway (within the cysteine and methionine metabolism pathway ko00270) were also induced (See Supplementary Table S1). Plant secondary metabolites including 12 candidate members of KEGG pathway ko00400 (phenylalanine, tyrosine, and tryptophan biosynthesis) showed induction in response to Pmeg/Ppal infection; 8 members of the terpenoid backbone biosynthesis pathway (ko00400) showed induction of two-fold or more with 2 additional members showing induction but less than two-fold and 8 genes in pathway terpenoid backbone biosynthesis (ko00900) showing induction (See Supplementary Table S3).

At least 240 genes putatively encoding Pr-related proteins were considered expressed (Table 1, mean expression ≥10 reads for at least one treatment) in the RNA-Seq analysis representing fourteen of the 17 characterized classes of Pr-proteins (Table 1). Infection repressed expression of 47 and 54 genes encoding Pr-related protein for Pmeg and Ppal, respectively while 81 and 87 genes encoding Pr-related proteins were induced by Pmeg and Ppal infection, respectively. Only 10 genes encoding Pr-related proteins were more highly expressed in response to Pmeg infection compared to Ppal infection while 24 genes encoding Pr-related proteins were more highly expressed in response to Ppal infection compared to Pmeg infection.

Using the 9 RNA samples used for RNA-Seq analysis, 17 out of the 19 cacao genes selected to show varied reactions to treatments had correlation coefficients >0.80 between the RNA-Seq and RT-qPCR expression results (13 with p ≤ 0.05) (see Supporting Excel File S1, Sheet H). In addition to these original 19 genes, another set of 17 cacao genes were selected based on their involvement in specific pathways including plant hormone signal transduction and plant-pathogen interactions and analyzed for temporal and tissue specific expression in response to Pmeg/Ppal infection (Figure 6 and see Supporting Excel File S1, Sheet I). Only 4 genes failed to show differential expression in response to Phytophthora infection in at least one time point/tissue combination (TcHyp2, TcHyp3, TcCop1, and TcCHD). Based on the euclidean distance of expression profile, the 36 genes can be categorized into 5 groups (Figure 6). Genes in group A were highly expressed in all the infected samples. These genes were repressed at the later time point in the control samples but showed steady or induced expression in infected samples. Genes in group B were highly induced in response to level of infection based on the pathogen load (Figure 6). Genes in group D were also induced in response to level of infection, but the level of induction was higher in infected pod husk (unwounded) compared to infected leaf disc. Genes in group C were either not differentially expressed or expression was altered in response to infection at a low level, and Group E genes were induced from very low constitutive expression levels (Figure 6 and see Supporting Excel File S1, Sheet I). Most of the genes studied showed coordinated regulation in response to infection by both Phytophthora species although in some cases the level and timing of induction was dependent on the Phytopthora species used (Figure 6 and see Supporting Excel File S1, Sheet I). The biggest difference was the more intense changes in expression for 14 of the 36 genes studied in Ppal inoculated wounded pod pieces.

PCoA was carried out using the RT-qPCR data from 36 genes for the 5 different treatment conditions (Figure 7). For each PCoA, the comparison of coordinates 1 and 2 is presented. In case of the leaf disc, a high amount of pathogen load was detected at 1 dpi and thus uninfected samples clustered separately from infected samples (Figure 7A). At the later time point (2 dpi), Ppal and Pmeg samples were further separated into upper and lower right quadrants (Figure 7B). In the case of unwounded pod husk, the level of infection was lower at 1 dpi and there was less separation of the samples between treatments but the Ppal treated unwounded pod husk samples were further separated from controls than the Pmeg samples (Figure 7C). At 3 dpi the Pmeg, Ppal and control samples clearly clustered separately in unwounded pod husk samples (Figure 7D). In wounded pod husks, only the Ppal treated samples clustered separately from control samples (Figure 7E).

Conventional RNA-Seq data analysis mainly focuses on differential induction or fold change but, it is also important to also consider mRNA expression level. Machine learning approaches, such as unsupervised SOM classification, allow us to consider not only fold change and expression level, but other unaccounted for variables. SOM is an unsupervised classification technique that reduces the dimension of data through the use of self-organizing neural networks, not imposing a priori class assumptions, but allowing the data to determine how genes should be classified. SOM-based gene expression analyses have been applied in diverse studies including studies on cell differentiation and development, organogenesis and tumor differentiation (Wirth et al., 2011). In this study, in addition to differential gene expression in response to treatment, a wide range in relative expression levels were demonstrated. In many cases, genes with “significant” differential expression associated with treatments were expressed at very low levels. This is especially true for many of the genes down-regulated in response to Phytophthora infection (SOM classes 12, 18, and 24). Often within a pathway, genes showing differential expression shared similar expression levels along with similar expression patterns. Clearly not all plant defense genes need to be induced to function. In some cases, individual members of a pathway failed to show induction but instead showed high constitutive expression. Furthermore, genes grouped together by SOM-based expression analysis can be used to set priorities for promoter analysis and other gene comparisons based on their similar regulation.

In a recent study using microarray technologies, Fister et al. (2016) identified 359 gene members within the cacao genome encompassing 15 of the 17 characterized classes of Pr-proteins. Approximately 50% of these genes are organized in tandem arrays within the cacao genome (Fister et al., 2016). Of the PR gene members, 67 were induced and 7 were repressed in cacao leaves in response to Ppal infection (Fister et al., 2016). In many cases leaves have been used as proxies due to the difficulty in acquiring pods in research settings. Since diseases like BPR primarily cause losses in yield by destroying the fruit directly, more studies focused on pods are needed. In the current study, all but 119 Pr-protein gene family members were found to be expressed (Table 1, RPKM ≥ 10 for at least one treatment) in pod which compares to the 126 gene members removed by the background filtration of the Fister et al. (2016) leaf microarray data set. Of the 67 Pr-protein gene family members induced in leaves infected by Ppal (Fister et al., 2016), 38 and 30 were also induced in cacao pod pieces inoculated with Ppal and Pmeg, respectively. In many cases, genes highly induced in leaves by Ppal were found to be expressed at high levels constitutively in control pod pieces and induced either at low levels or not at all in response to Phytophthora infection. Conversely, Pr-protein gene members induced at low levels in leaves often had low constitutive expression and were highly induced in pod pieces. On the surface this suggests basic differences in how genes encoding Pr-proteins respond to Phytophthora infection in pods and leaves but this would need to be verified using the same cacao germplasm. In total, 87 and 81 Pr-protein gene family members were induced in pod pieces inoculated with Ppal and Pmeg, respectively. In most cases, the genes observed uniquely induced in the current study also had low constitutive expression in pod pieces. Genes with high constitutive expression and limited induction tended to be placed in SOM class 1 whereas the genes with low constitutive expression and high induction tended to be placed in SOM classes 2 and 3. In our study, 54 and 47 Pr-protein genes were repressed in pod pieces responding to infection by Ppal and Pmeg, respectively, and these tended to fall into SOM classes 18, 22, 23, and 24. No genes encoding Pr-13 and 15 were found in cacao (Fister et al., 2016) and no genes encoding Pr-12 or Pr-14 were differentially expressed compared to control samples in response to Phytophthora infection.

Although, genes associated with many different KEGG pathways, some participating in basic aspects of pod physiology, were altered in response to Phytophthora infection, pathways associated with plant defense would seem to be most pertinent. For example, CERK1, CNGC, FLS2, and BAK1 modulate PAMP signal perception and transduction (Wu et al., 2014; Saand et al., 2015). Endocytosis is critical to the subsequent induction of plant innate immunity through FLS2 and BAK1 which participate together in this process (Robatzek et al., 2006; Chinchilla et al., 2007). FLS2 is targeted to the vacuole after PAMP recognition, phosphorylation, and ubiquination through the late endosome. Expression of several components associated with the FLS2 signaling pathway (or a parallel system in this case) was enhanced ending with the up-regulation of a constitutively highly expressed Pr-1 gene (Tc02_g002410). The most notable gene associated with endocytosis includes several Rab7 candidates, the encoded Rab7 protein being recruited to the extrahaustorial membrane (EHM) in interactions with Phytophthora (Bozkurt et al., 2015). BAK1 but not FLS2, is also recruited to the EHM (Bozkurt et al., 2015). Although, the inductions observed in the plant-pathogen interaction pathways were significant, the interactions studied herein were susceptible interactions. Moreover, the infected materials studied were in the necrotrophic phase. Therefore, the plant responses may have been limited compared to the reactions if immunity were induced.

Plant hormones are intimately intertwined with the induction of both susceptibility and resistance (Bari and Jones, 2009). Differential expression was noted for candidate genes associated with the actions of auxins, gibberellins, cytokinines, ABA, BRs, JA, and salicylic acid (SA). There were 14 candidate gene members associated with pathway ko04075 (plant hormone signal transduction) showing down-regulation and in many cases they represent orthologs of genes induced in the same pathway. Notably no members of the ethylene, JA or SA signal transduction pathways, pathways that are closely associated with plant defense, showed repression. The list of candidate genes showing induction in response to Phytophthora infection complete or nearly complete, respectively, the ethylene (Wang et al., 2002) and JA (Wasternack, 2007) biosynthesis pathways (See Supplementary Table S1). Though ethylene and JA play a major role in plant defense (Mersmann et al., 2010; Santino et al., 2013; Broekgaarden et al., 2015), induction of ethylene and JA related genes does not always mean involvement in resistance (Eshraghi et al., 2014).

From the literature, little support can be found for JA participating in penetration resistance to Phytophtora (a hemibiotroph) infection. Instead, SA is more often implicated (Halim et al., 2007; Loake and Grant, 2007). The SA signal transduction pathway was highlighted by the KEGG analysis. The non-expressor of PR1 (NPR1) plays a central role in regulating the transcriptional response to plant pathogens through the SA pathway (Fellbrich et al., 2002). NPR3, on the other hand, functions as a repressor of NPR1-mediated defense response (Shi et al., 2013a). Functional orthologs of NPR1 (Tc09_g007660) and NPR3 (Tc06_g011480) have been identified in cacao (Shi et al., 2010, 2013b) and the “NPR1” identified in our KEGG analysis as induced has been annotated and functionally proven to instead be NPR3 (Shi et al., 2013a). The functionally annotated NPR1 gene in cacao showed constitutive expression in the RNA-Seq analysis. NPR3's induction supports the contention that the Pr-1 gene induced in response to Phytophthora was more reasonably associated with the FLS2/BAK1signalling pathway (or rather a parallel system in the cacao/Phytophthora interaction). Along the same lines, although a candidate BAK1 gene showed up-regulation, a candidate BRI1 gene was down-regulated. This is of interest since BAK1 shares the ability to interact with FLS2 of the plant defense pathway and BRI1 of the brassinosteroid signaling pathway (Wang, 2012). In this case, the brassinosteroid signaling pathway appears more likely down regulated, something observed in a previous study of the interaction between cacao and Moniliophthora roreri, the causal agent of frosty pod rot (Ali et al., 2014).

Several component pathways associated with secondary metabolite production showed significant induction in response to Phytophthora infection. Although, genes identified were not always hugely induced or highly expressed, they make up almost the entire biosynthetic pathway from D-erythrose- 4-phosphate to phenylalanine. From phenylalanine, 13 members of the phenylpropanoid biosynthesis (PAL) pathway (ko00949) were induced, again largely completing the pathway. The PAL pathway is well studied and known to contribute to plant defense by its many off shoots including hydroxycinnamic acids, monolignols/lignin, flavonoids, isoflavonoids, and stilbenes (Dixon et al., 2002). A notable pathway gene member missing from this list is cinnamoyl-CoA reductase (EC:1.2.1.44), although a candidate cinnamoyl-CoA reductase gene (Tc06_g004070) found in the KEGG database showed high constitutive expression in all samples. Among secondary products, lignin and other compounds associated with the PAL pathway have been most commonly associated with resistance/tolerance to Phytophthora in cacao (Ndoumou et al., 1995; Omokolo et al., 2003; Boudjeko et al., 2007; Djocgoue et al., 2007; Louise et al., 2011).

The terpenoid backbone biosynthesis pathway (ko00900), like the PAL-pathway, showed significant induction. Among the secondary products dependent on this pathway are various terpenoids, alkaloids, carotinoids, and N-glycans. The link between JA and ethylene mediated signal transduction and terpene biosynthesis is well established (Dudareva et al., 2006). Terpenoid compound are commonly studied for their activity against herbivores (Mithöfer and Boland, 2012) but the antimicrobial activity of terpenoid compounds have not been studied in any detail in cacao.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.