- 1Hans Merensky Chair in Avocado Research, University of Pretoria, Pretoria, South Africa
- 2Department of Biochemistry, Genetics and Microbiology, Faculty of Natural and Agricultural Sciences, University of Pretoria, Pretoria, South Africa
- 3Forestry and Agricultural Biotechnology Institute, Faculty of Natural and Agricultural Sciences, University of Pretoria, Pretoria, South Africa
The hemibiotrophic plant pathogen Phytophthora cinnamomi Rands is the most devastating pathogen of avocado (Persea americana Mill.) and, as such, causes significant annual losses in the industry. Although the molecular basis of P. cinnamomi resistance in avocado and P. cinnamomi virulence determinants have been the subject of recent research, none have yet attempted to compare the transcriptomic responses of both pathogen and host during their interaction. In the current study, the transcriptomes of both avocado and P. cinnamomi were explored by dual RNA sequencing. The basis for partial resistance was sought by the inclusion of both susceptible (R0.12) and partially resistant (Dusa®) rootstocks sampled at early (6, 12 and 24 hours post-inoculation, hpi) and late time-points (120 hpi). Substantial differences were noted in the number of differentially expressed genes found in Dusa® and R0.12, specifically at 12 and 24 hpi. Here, the partially resistant rootstock perpetuated defense responses initiated at 6 hpi, while the susceptible rootstock abruptly reversed course. Instead, gene ontology enrichment confirmed that R0.12 activated pathways related to growth and development, essentially rendering its response at 12 and 24 hpi no different from that of the mock-inoculated controls. As expected, several classes of P. cinnamomi effector genes were differentially expressed in both Dusa® and R0.12. However, their expression differed between rootstocks, indicating that P. cinnamomi might alter the expression of its effector arsenal based on the rootstock. Based on some of the observed differences, several P. cinnamomi effectors were highlighted as potential candidates for further research. Similarly, the receptor-like kinase (RLK) and apoplastic protease coding genes in avocado were investigated, focusing on their potential role in differing rootstock responses. This study suggests that the basis of partial resistance in Dusa® is predicated on its ability to respond appropriately during the early stages following P. cinnamomi inoculation, and that important components of the first line of inducible defense, apoplastic proteases and RLKs, are likely to be important to the observed outcome.
Introduction
Avocado (Persea americana Mill.) is an economically important fruit crop belonging to Lauraceae. This family is arguably the most economically important family in Laurales, consisting of over 2,500 species with a wide tropical and subtropical distribution (Kostermans, 1957; Gentry, 1988; Li, 1995; van der Werff and Richter, 1996; Chanderbali et al., 2001; Song et al., 2016). Of these, avocado contributes significantly to the global economy with production surpassing 7 million tons in 2019 and an estimated gross production value of $6.56 billion (constant 2014–2016, int. $; http://www.fao.org/faostat/en/#home). In recent years, avocados have gained popularity because of their nutritional profile coupled with numerous health benefits associated with their consumption (Hurtado-Fernández et al., 2018). However, the production of avocado is severely threatened by the persistence of several pathogens in major avocado-growing regions of the world (van den Berg et al., 2021).
One of the most significant threats to avocado production is Phytophthora root rot (PhRR), caused by the oomycete Phytophthora cinnamomi Rands, a plant pathogen of global significance due to its effects on wild and cultivated plants, which presents a serious threat to the diversity and structure of natural ecosystems(Hardham and Blackman, 2018). P. cinnamomi is a hemibiotrophic oomycete, implying that its interaction with the host is initially biotrophic, following host colonization, after which a necrotrophic phase is invoked. In avocado, P. cinnamomi invades the sub-surface feeder roots where, in susceptible rootstocks, the root structure is rapidly colonized and functionally damaged, limiting water and nutrient uptake. Control and management of PhRR spread remain a challenging endeavor because of P. cinnamomi's broad host range coupled with the ability of tolerant or partially resistant plants to remain asymptomatic. Chemical control of P.cinnamomi has predominantly been conducted using phosphites in the form of sprays or trunk injections (Pegg et al., 1985; Hardy et al., 2001). However, decreased sensitivity to treatment with phosphites and the questionable sustainability of clearing infested orchards point to integrated approaches as the only viable method for PhRR control and management (Dobrowolski et al., 2008; Wolstenholme and Sheard, 2010). These integrated approaches include methods such as use of tolerant or partially resistant rootstocks such as Dusa®, together with chemical control and efficient orchard management (Hardham, 2005; Wolstenholme and Sheard, 2010). Nonetheless, the impact of PhRR remains because of reduced fruit size, quality and shelf life, decreasing yield and, in the case of severe infection, killing of trees (Dann et al., 2013).
To augment the use of partially resistant rootstocks as a control measure for future management of PhRR, a comprehensive understanding of the P. americana-P. cinnamomi interaction is essential (van den Berg et al., 2021). One way of understanding host-pathogen interactions is quantifying gene expression in both the pathogen and the host during infection. This would invariably lead to a better understanding of P. cinnamomi molecular responses during infection, paving the way for the development of novel targets for precise and viable control measures. Meanwhile, understanding host responses could aid in the development or selection of novel PhRR-resistant rootstocks.
The availability of high-quality omics data to study the P. americana-P. cinnamomi pathosystem has improved significantly over the past few years. The newest P. cinnamomi genome represents a significant improvement over previous genome assemblies (Engelbrecht et al., 2021). Similarly, the P. americana West-Indian pure accession (WI) rootstock genome from the Avocado Genome Consortium boasts a chromosome-level reference assembly (Avocado Genome Consortium, article in preparation). These resources are invaluable for future studies on the P. americana-P. cinnamomi interaction at the genetic level. There are, however, only a few resources available for understanding the expression of genes for this pathosystem on a global scale (Mahomed and van den Berg, 2011; Reeksting et al., 2014, 2016; Reitmann et al., 2017; van den Berg et al., 2018b; Engelbrecht et al., 2021). Notably, insights into this pathosystem could be significantly improved by high-throughput transcriptome sequencing such as RNA-seq. This technology profiles transcriptomes using deep-sequencing technologies, enabling quantification of gene expression and providing a novel understanding of changing RNA expression levels that occur in response to external stimuli (Martin et al., 2013). In particular, dual RNA sequencing data, which allow for investigations on the interaction of both host and pathogen during the process of infection (Westermann et al., 2012; Wolf et al., 2018), could further aid in understanding this important pathosystem.
To understand any pathosystem, one needs to only consider the evolutionary arms race that exists between a pathogen and a host. Generally, plants recognize pathogens through one of two methods, the first is by recognition of pathogen-associated molecular patterns (PAMPs; Davis and Hahlbrock, 1987; Zipfel and Felix, 2005). These molecular patterns are recognized by membrane-bound pattern recognition receptors (PRRs), leading to activation of the innate immune-response, pattern-triggered immunity (PTI). Notably, most PRRs are either transmembrane receptor-like kinases (RLKs) or receptor-like proteins (RLPs) and, as such, constitute an essential part of the pathogen-sensing machinery in plants (Böhm et al., 2014; Couto and Zipfel, 2016; Yu et al., 2017). Here, plant proteases have important roles in the signaling process, both prior to and following the recognition of PAMPs or host-derived damage-associated molecular patterns (DAMPs; Hou et al., 2018; Godson and van der Hoorn, 2021). Intriguingly, pathogens have been known to evade recognition by PRRs through the actions of effector proteins, some of which seem to trick the host plant into activating RLKs, which are involved in growth and development (Dou and Zhou, 2012). In turn, pathogen effectors are often recognized either directly or indirectly by various resistance (R) proteins, culminating in the second stronger inducible immune-response, effector-triggered immunity (ETI; Jones and Dangl, 2006; Monteiro and Nishimura, 2018). Typically, ETI leads to a specialized form of localized programmed cell death (PCD), known as the hypersensitive response (HR), to arrest pathogen spread (Cui et al., 2015). As such, ETI is known to be effective against biotrophic or hemibiotrophic pathogens but not against necrotrophs (Glazebrook, 2005). In plants, the largest group of R proteins are that of cytoplasmic nucleotide-binding leucine-rich repeat (NLR) proteins (McDowell and Woffenden, 2003). Notably, the full complement of avocado NLRs was recently described using some of the RNA-seq data presented in this study (Fick et al., 2022).
To gain a better understanding of the P. americana- P. cinnamomi interaction at a molecular level, a comprehensive root transcriptomic analysis was performed using dual RNA sequencing technology. The experimental approach included the use of both a susceptible (R0.12) and a partially resistant avocado rootstock (Dusa®), inoculated with P. cinnamomi and harvested at 6, 12, 24, and 120 h post-inoculation (hpi). This approach allows for the simultaneous investigation of molecular components that govern susceptibility and/or resistance in avocado roots over time, as well as pathogen virulence strategies, at both the biotrophic and necrotrophic life stages and the transition between the two. We utilized the reference genomes for P. americana WI pure accession and P. cinnamomi GKB4 for RNA-seq read alignment. Differentially expressed genes (DEGs) were identified, and the investigation implicated several P. americana genes encoding for pathogenesis-related (PR) genes, receptor-like kinases (RLKs) and proteases, and others in the defense response against P. cinnamomi. Correspondingly, several P. cinnamomi effector coding genes were identified and linked to the infection process including but not limited to RxLRs, elicitins, extra-cellular protease inhibitors (EPIs), and xyloglucan-specific endo-beta-1,4-glucanases (XEGs). Most importantly, the analyses conducted in this study showed that Dusa® exhibited drastic responsive gene expression differences during the biotrophic-necrotrophic transition phase of infection when compared to the susceptible host. This observation suggests that a sweeping responsive change in gene expression is the molecular basis of partial resistance.
This is the first study that conducted dual RNA-seq to examine the global gene expression of both P. cinnamomi and two inoculated avocado rootstocks of varying susceptibility during infection. The study provides the most comprehensive overview of avocado defense response pathways to date and a glimpse at the arsenal of virulence factors utilized by P. cinnamomi to promote disease. Besides candidates for functional validation, the results generated by this study will provide resources for future identification of markers for improved rootstock breeding and screening.
Materials and Methods
Plant and Pathogen Material
Both in vitro and in planta materials were used for this research. For the in vitro material, P. cinnamomi (isolate GBK4) was grown on 5% V8 agar for 5 days at 25°C in the dark. Mycelia were harvested and snap-frozen in liquid nitrogen before being stored at −80°C until RNA extraction. For the in planta material, clonal avocado plantlets approximately 2 years old were obtained from Westfalia Technological Services (Tzaneen, Limpopo, South Africa). A zoospore suspension of P. cinnamomi (isolate GKB4) was produced as described in Engelbrecht and van den Berg (2013). Both partially resistant (Dusa®) and susceptible (R0.12) avocado rootstocks were inoculated using the zoospore suspension (1.4 x 105 zoospores/ml) or dH20 (mock-inoculation), submerging the roots and lower stem for 2 h. Following inoculation, all the plantlets were re-planted in 1.5-L plastic bags containing a 1:1 perlite:vermiculite mixture. Root samples were harvested for three replicate samples, with three plantlets per replicate at each time-point: 6, 12, 24, and 120 hpi. Due to resource limitations, control plants were only harvested at 24 hpi. The harvested materials were snap-frozen in liquid nitrogen before being stored at −80°C until RNA extraction.
RNA Extraction and Sequencing
All the sample materials were processed to a fine powder using the IKA®Tube Mill control(IKA®, Staufen, Germany).A modified CTAB RNA extraction protocol was used to extract total RNA from all the samples as described in Engelbrecht and van den Berg (2013). RNA concentration was quantified using a Nanodrop ND−1000 spectrophotometer (Nanodrop Technologies Inc., Montchanin, DE, United States). RNA integrity was assessed on a 2% TAE agarose gel under non−denaturing conditions. Potential contaminating DNA was removed from total RNA using RNase−free DNase (Fermentas Life Sciences, Hanover, MD, United States) according to the manufacturer′s guidelines. DNAse−treated RNA was then purified using an RNeasy®Mini−Elute® cleanup kit (Qiagen, Valencia, CA, United States) according to the manufacturer′s guidelines. Purified total RNA quality was assessed on Agilent 2100 Bioanalyzer (Agilent Technologies, CA, United States). Finally, 1 μg purified RNA from each biological replicate was combined to constitute the samples used for each of the respective libraries (Supplementary Table 1). Purified RNA was shipped to Macrogen (Macrogen Inc., Seoul, GG, South Korea) for sequencing on an Illumina HiSeq platform using the PE150 mode with a minimum raw read count of 80 million per sample.
RNA Sequencing Quality Control and Genome Alignment
Raw RNA-seq data were filtered using Trimmomatic v0.39 [trimming parameters: ILLUMINACLIP: TruSeq3-PE.fa:2:30:10:2:true LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36] to remove low-quality bases and Illumina sequencing adapters (Bolger et al., 2014). The read quality before and after trimming was assessed using FastQC v0.11.9 (Andrews, 2010) and summarized using MultiQC v1.9 (Ewels et al., 2016). The P. americana WI pure accession genome was obtained from the Avocado Genome Consortium (Avocado Genome Consortium, article in preparation). The P. cinnamomi GKB4 genome was obtained from Engelbrecht et al. (2021). The genomes were concatenated and used for RNA-seq data alignment to eliminate alignment bias from conserved sequences. Trimmed RNA-seq reads were aligned to the combined P. americana-P.cinnamomi genome using HISAT2 v2.0.6 (Kim et al., 2019). The resultant alignment files were processed and combined using samtools (Li et al., 2009) and used to obtain high-confidence splice junctions from Portcullis v1.2.0 (Mapleson et al., 2018). A second HISAT2 alignment was performed using the high confidence splice junctions obtained from Portcullis. Alignment statistics were calculated using samtools idxstats and flagstat. Transcripts per million (TPMs) were calculated using StringTie v2.1.4 (Pertea et al., 2015).
Differential Gene Expression Analyses
Read counts were quantified at the gene level using featureCounts v2.0.1 (Liao et al., 2014), informed by the genome annotation file of either P. americana or P. cinnamomi in order to separate the reads from each of the respective genomes. Read data were imported into R 4.0.4 (R Core Team, 2021) and analyzed using DeSeq2 v1.32.0 (Love et al., 2014). Data for transcripts in each library that fell below 10 reads were removed, and transcripts without any read data were removed from the data set. Mock-inoculated P. americana libraries for R0.12 or Dusa® were set as reference libraries for the P. americana gene sets against the inoculated P. americana R0.12 or Dusa® rootstock libraries, respectively. The mycelial libraries were set as the reference for the P. cinnamomi gene sets against the inoculated P. americana rootstock libraries. DEGs were identified by Wald test, and multiple hypotheses testing correction was conducted using the Benjamini-Hochberg (BH) false discovery rate (FDR) method. Significantly DEGs were defined using an FDR cutoff (adjusted p-value) of <0.05 and a log2(fold change) of more than 1 (upregulated) and < −1 (downregulated). The varianceStabilizingTransformation (VST) function in the DESeq2 package was used to transform the DESeq2 normalized data. The VST data were used for principal component analyses (PCAs) of the top 500 variable genes using plotPCA. Hierarchical clustering of 50 genes with the highest variation across the samples was conducted using hclust and dendextend (Galili, 2015), followed by visualization using heatmap.2 in gplots v3.1.1 (Warnes et al., 2009). Additionally, a list of the top 100 variable genes from the P. americana dataset was exported for further analyses.
Additional time course analyses were also conducted for both P. americana and P. cinnamomi using TCseq v1.16.0 (Wu and Gu, 2021) for downstream temporal pattern analyses. Raw count data (featureCounts) for both P. americana and P. cinnamomi were imported into R separately and divided into two datasets, one for each rootstock. Differential analyses were conducted using the DBanalysis function in TCseq, filtering out sample data points with <10 reads and any genes in which less than two samples contain read counts. A time course table was generated using Reads Per Kilobase Million (RPKM) normalized read counts. Significantly DEGs were defined using an FDR cutoff of <0.05 and a log2(fold change) of more than 1 and < −1. The Z-scores of DEGs were calculated and used for clustering analyses using fuzzy c-means (Futschik and Carlisle, 2005) with k = 5, as implemented by the timeclust function in TCseq. Clusters were visualized using the timeclustplot function.
Gene Ontology Enrichment Analyses
P. americana and P. cinnamomi significant DEG lists from DESeq2 were imported into R 4.0.4 and split into up- and downregulated sets. Gene ontology enrichment was performed using GoSeq v1.44.0 (Young et al., 2010). Over- and under-representation of GO terms among DEGs in each gene set were calculated using the Wallenius approximation, and overrepresented GO terms with p-values <0.05 (BH FDR) were selected to constitute enriched GO term lists. Gene lists from each cluster obtained following TCseq clustering analyses were also subjected to GO enrichment using the same parameters. Lists of GO terms were reduced using REVIGO (Supek et al., 2011) and visualized in Tableau v2021.2.1.
Significant Pathway, Gene Family, and Virulence Factor Identification
Predicted proteins from the P. americana WI pure accession genome were assigned to MapMan bins using both Mercator v3.6 and Mercator v4.0 applying default parameters and all available datasets (Lohse et al., 2014; Schwacke et al., 2019). Significant DEGs consequent to DeSeq2 analyses, from each time point and rootstock, were input into MapMan v3.6.0RC1 (Schwacke et al., 2019) to visualize and compare pathways and gene families for which expression differs in the rootstocks over time and between the rootstocks at defined time points following inoculation with P. cinnamomi. Alternatively, DEGs were visualized using heatmaps generated using pheatmap v1.0.12 (Kolde, 2012) in RStudio v1.4.1106 (RStudio Team, 2021). Localization prediction for predicted proteins of interest was performed using a lightweight deep neural network web service (Stärk et al., 2021). In order to identify potential P. cinnamomi virulence factors, predicted protein sequences from the corresponding P. cinnamomi TCseq cluster gene lists were compared to the Pathogen-Host Interactions (PHI) database v4.12 (Urban et al., 2019). This comparison was conducted using blastp in the BLAST+ suite (Camacho et al., 2009), with an Expect value (E) cutoff of 1.0e−5 and returning only the top match for each query.
Results
Transcript Expression Analyses
Following quality control, a total of 4.99 billion RNA sequencing reads were generated, with approximately 2.54 billion from the Dusa® libraries, 1.98 billion from the R0.12 libraries, and 475 million from the P. cinnamomi mycelial libraries (Supplementary Table 2). Thus, an average of 150 million filtered reads were obtained per library. Overall, 81.11% of the reads in the Dusa® libraries and 74.29% of the reads in the R0.12 libraries aligned to the P. americana genome. The read percentage from both Dusa® and R0.12 inoculated libraries that aligned to the P. americana genome was substantially lower at 120 h hpi in comparison to the earlier time points, 62.32% vs. 81.54%, coinciding with a sharp rise in the average read percentage aligning to the P. cinnamomi genome at 120 h hpi, 3.07% vs. 0.16% at earlier time points. An average of 0.89% of the reads in the inoculated sample libraries and 94.36% in the mycelial libraries aligned to the P. cinnamomi genome. On average, 0.004% of the reads from the mock-inoculated (dH2O) R0.12 libraries and 0.044% of the reads from the mock-inoculated Dusa® libraries mapped back to the P. cinnamomi genome.
Gene expression data were obtained for up to 32,383 Dusa®, 32,588 R0.12, and 16,985 P. cinnamomi genes, depending on the sample. When comparing the 6, 12, 24, and 120 hpi Dusa®samples to the mock-inoculated samples, 7,916, 7,342, 6,159 and 4,626 P. americana DEGs were observed at each time point, respectively (Figure 1A; Supplementary Table 3). In the R0.12 6, 12, 24, and 120 hpi samples, 8,148, 721, 97, and 4,003 P. americana DEGs were observed, respectively, when compared to the mock-inoculated samples (Figure 1A; Supplementary Table 4). When considering P. cinnamomi, 4,821, 3,288, 5,011, and 6,779 genes were differentially expressed (DE) in the 6, 12, 24, and 120 hpi Dusa® inoculated samples compared to the mycelial controls (Figure 1B; Supplementary Table 5). Similarly, when compared to the mycelial controls, the 6, 12, 24, and 120 hpi R0.12 samples yielded 5,733, 1,755, 3,533, and 6,902 P. cinnamomi DEGs, respectively (Figure 1B; Supplementary Table 6). The differences in log2(fold change) between rootstocks, compared at each time point, were also obtained for both P. americana (Supplementary Table 7) and P. cinnamomi genes (Supplementary Table 8).
Figure 1. Graphical representation of differentially expressed genes (DEGs) over time in both a susceptible (R0.12) and a partially resistant (Dusa®) rootstock inoculated with Phytophthora cinnamomi. (A) Total number of significantly up- [log2(fold change) > 1, p = 0.05] and downregulated avocado genes [log2(fold change) < −1, p = 0.05] in P. cinnamomi-inoculated rootstocks as compared to mock-inoculated controls, shown at 6, 12, 24, and 120 h post inoculation (hpi) in both partially resistant and susceptible avocado rootstocks. (B) Total number of significantly up- [log2(fold change) > 1, p = 0.05] and downregulated P. cinnamomi genes [log2(fold change) < −1, p = 0.05] in inoculated avocado rootstocks as compared to in vitro mycelial controls, over time in both Dusa® and R0.12.
Principal component analyses (PCAs) were performed using the VST-normalized expression data of the top 500 variable genes from either P. americana or P. cinnamomi. The PCAs were performed to visualize and group the sample libraries in a way conducive to gaining a broad perspective on similarities and differences. For P. americana top variable genes, the PCAs identified three distinct sample groups (Figure 2A). The first group contained representative points for mock-inoculated samples from both rootstocks and inoculated R0.12 samples at 12 and 24 hpi. The second group was composed exclusively of early time-point samples from both Dusa® (6, 12, and 24 hpi) and R0.12 (6 hpi). The third group contained only late time-point samples (120 hpi) from both R0.12 and Dusa®. Variance across the first principal component axis, which separates groups 1 and 3 from group 2, represented the largest variation between samples at 64%, while variance across the second principal component axis, separating groups 1 and 2, was only 21%. For P. cinnamomi top variable genes, The PCAs resulted in four distinct sample groups (Figure 2B). Group 1 represented in vitro mycelia samples only, while group 2 contained all but one sample of R0.12 at 12 and 24 hpi. Group 3 contained representatives from only the late time-point samples (120 hpi) in both rootstocks. Group 4 contained the early time-point samples for Dusa® (6, 12, and 24 hpi) and R0.12 (6 hpi) along with one 12 hpi R0.12 sample. Variance across the first principal component axis accounted for the most variation between groups at 74%, exemplified by groups 1 and 4 occupying opposite ends of the axis. Variance across the second principal component axis accounted for only 8% of the difference between samples, largely distinguishing group 3 from the rest.
Figure 2. Principal component analysis (PCA) representing variation in gene expression. Expression values were obtained following differential gene expression analyses and variance stabilized transformation. Principal component 1 (PC1) and PC2 were arbitrarily assigned to either axis and represent the total representative variation between samples along that axis. Samples that cluster together were grouped and labeled accordingly. (A) Persea americana genes over time, following inoculation with Phytophthora cinnamomi. The top 500 variably expressed P. americana genes across all samples were used for PCA. Samples include both partially resistant (Dusa®) and susceptible (R0.12) P. americana rootstocks that were either mock-inoculated (dH2O) or inoculated with P. cinnamomi and harvested at 6, 12, 24, or 120 h post inoculation (hpi). (B) Phytophthora cinnamomi genes over time, following inoculation of Persea americana. The top 500 variably expressed P. cinnamomi genes across all samples were used for PCA. Samples include in vitro P. cinnamomi mycelial controls and in planta samples, partially resistant (Dusa®) and susceptible (R0.12) P. americana rootstocks that were inoculated with P. cinnamomi and harvested at 6, 12, 24, or 120 h post inoculation (hpi).
Hierarchical clustering of the top 50 variable genes in both P. americana and P. cinnamomi was also performed to broadly visualize the differences in gene expression between samples. The dendrogram of the top 50 variable genes in P. americana showed two distinct groups of genes with approximately opposite patterns of expression (Figure 3A). The first smaller group contained genes that were expressed to a greater degree in the mock-inoculated samples, 12, 24, and 120 hpi R0.12 samples, and 120 hpi Dusa® samples. Genes in the second larger group displayed greater levels of expression in the 6 hpi R0.12 and 6, 12 and 24 hpi Dusa® samples. In P. cinnamomi, four groups of genes with varying expressions were observed (Figure 3B). The first group contained genes that were expressed highest at 120 hpi in both rootstocks, while the second group was expressed highest at 24 and 120 hpi, although to a greater degree in Dusa®. The third group's genes were expressed to the greatest degree in mycelial samples and at 12 and 24 hpi in R0.12 but not Dusa®. Lastly, genes in group 4 were expressed largely at 6 and 120 hpi in R0.12 and 6 and 12 hpi in Dusa®.
Figure 3. Hierarchical clustering of the top 50 variably expressed Persea americana and Phytophthora cinnamomi genes. Expression values were obtained following differential gene expression analyses and variance stabilized transformation. The green color represents upregulation, while red indicates downregulation, the scale of which is indicated by the intensity of the color. Gene identifiers are located on the right of the heatmap. Dual RNA-seq library references are indicated along the bottom of the heatmap. Sample clustering is represented by a dendrogram on the left of the heatmap that is color-coded and labeled accordingly. (A) Top 50 variably expressed P. americana genes across all samples were visualized by hierarchical clustering. Samples include both partially resistant (Dusa®) and susceptible (R0.12) P. americana rootstocks that were either mock-inoculated (dH2O) or inoculated with P. cinnamomi and harvested at 6, 12, 24, or 120 h post inoculation (hpi), indicated along the top of the heatmap. (B) Top 50 variably expressed P. cinnamomi genes across all samples were visualized by hierarchical clustering. The green color represents upregulation, while red indicates downregulation, the scale of which is indicated by the intensity of the color. Samples include in vitro P. cinnamomi mycelial controls and in planta samples, partially resistant (Dusa®) and susceptible (R0.12) P. americana rootstocks that were inoculated with P. cinnamomi and harvested at 6, 12, 24, or 120 h post inoculation (hpi), indicated along the top of the heatmap.
Temporal Cluster Analyses
P. americana
Temporal cluster analyses were performed to separate the P. americana genes based on their expression in mock-inoculated samples (0 hpi), and over time following P. cinnamomi inoculation. Five clusters with distinct temporal patterns were obtained for both Dusa® and R0.12 (Figures 4A, 5A; Supplementary Table 9). In Dusa®, cluster 1 represented genes with reduced expression at 6, 12, and 24 hpi followed by further reduction at 120 hpi. Overrepresented biological process (BP) GO terms resulting from this cluster included lignin, pectin, and chitin catabolic processes, several terms related to growth and development such as response to auxin, root development, nitrate import, and several terms related to oxidative stress (Figure 4B; Supplementary Table 10). The expression of genes in cluster 2 increased significantly at 6 hpi, reducing gradually at each subsequent time point until the expression at 120 hpi closely resembled that of the mock-inoculated samples (0 hpi). Overrepresented BP terms from genes in this cluster included transcription, translation, RNA modification, intracellular protein transport, and exocytosis. Cluster 3 contained genes that were expressed at their lowest levels at 6, 12, and 24 hpi, while the expression at 0 and 120 hpi seemed similar. Gene ontology enrichment of genes in this cluster only yielded three BP terms, DNA-templated regulation of transcription, protein folding, and translation. Cluster 4 contained genes that increased in expression during the early time points, peaking at 24 hpi, followed by a drop slightly below baseline expression at 120 hpi. Enriched BP terms from this cluster included phosphorylation, defense response, DNA repair, RNA-templated transcription, and several terms related to metabolism and photosynthesis. Expression of genes in cluster 5 remained largely stable across the mock-inoculated and early time-point samples, followed by a definitive increase at 120 hpi. Overrepresented BP terms from genes in this cluster included phosphorylation, negative regulation of endopeptidase activity, immune response, plant type hypersensitive response, cellular aromatic compound metabolic process, response to auxin, response to chitin, wounding and fungi, and several terms related to oxidative stress.
Figure 4. Temporal clustering and Gene Ontology enrichment of Persea americana differentially expressed genes (DEGs) in Dusa® following inoculation with Phytophthora cinnamomi. (A) Sample libraries from the partially resistant P. americana rootstock Dusa® following inoculation with P. cinnamomi or mock inoculation (dH2O) were used in clustering analyses. DEG Z-scores, based on reads per kilobase million (RPKM) normalized read counts, were calculated and used for clustering analyses using the fuzzy c-means (Futschik and Carlisle, 2005) in the R package TCseqv1.16.0 (Wu and Gu, 2021). Inoculated samples included 6, 12, 24, and 120 h post inoculation (hpi) and are indicated on the x-axis, with the mock-inoculated samples indicated as 0 hpi. Membership scores produced by the fuzzy c-means algorithm are color-coded and indicated to the right of each cluster. (B) Member genes from each cluster were subjected to Gene Ontology (GO) enrichment, and significantly overrepresented GO terms, with adjusted p-value < 0.05 Benjamini-Hochberg false discovery rate, were reduced in REVIGO (Supek et al., 2011). The cluster from which the terms originate is indicated on the leftmost side of the Figure, followed by the biological process GO term. The overrepresented p-value of each GO term is indicated along the x-axis as –log10 (p-adjusted).
Figure 5. Temporal clustering and Gene Ontology enrichment of Persea americana differentially expressed genes (DEGs) in R0.12 following inoculation with Phytophthora cinnamomi. (A) Sample libraries from the susceptible P. americana rootstock R0.12 following inoculation with P. cinnamomi or mock inoculation (dH2O) were used in clustering analyses. DEG Z-scores, based on reads per kilobase million (RPKM) normalized read counts, were calculated and used for clustering analyses using the fuzzy c-means (Futschik and Carlisle, 2005) in the R package TCseqv1.16.0 (Wu and Gu, 2021). Inoculated samples included 6, 12, 24, and 120 ho post inoculation (hpi) and are indicated on the x-axis, with the mock-inoculated samples indicated as 0 hpi. Membership scores produced by the fuzzy c-means algorithm are color-coded and indicated to the right of each cluster. (B) Member genes from each cluster were subjected to Gene Ontology (GO) enrichment and significantly overrepresented GO terms, with adjusted p-value < 0.05 Benjamini-Hochberg false discovery rate, were reduced in REVIGO (Supek et al., 2011). The cluster from which the terms originate is indicated on the leftmost side of the Figure, followed by the biological process GO term. The overrepresented p-value of each GO term is indicated along the x-axis as –log10 (p-adjusted).
In the susceptible rootstock R0.12, temporal cluster analyses yielded five clusters that were objectively different to those obtained for Dusa® (Figure 5A; Supplementary Table 9). Genes in cluster 1 showed substantially decreased expression at 6 hpi, however, at 12, 24, and 120 hpi, their expression could be matched to the expression in the mock-inoculated samples. GO enrichment of genes in this cluster resulted in a limited list of BP terms, DNA-templated regulation of transcription, positive regulation of catalytic activity, and response to osmotic stress (Figure 5B; Supplementary Table 11). Cluster 2 showed a trend opposite to that of cluster 1, with genes from this cluster displaying increased expression at 6 hpi, followed by a return to baseline at 12 hpi. Overrepresented BP terms obtained from genes in this cluster included phosphorylation, defense response, several terms related to transcription, translation, growth and development, metabolism, and photosynthesis. The expression of genes in cluster 3 remained relatively unchanged in the mock-inoculated and early time-point samples, with a distinctive increase at 120 hpi. Several enriched BP terms were obtained from genes in this cluster, including phosphorylation, response to chitin, wounding and fungi, response to auxin, cellular response to nitrogen starvation, negative regulation of endopeptidase activity, and several terms related to metabolism and oxidative stress. Cluster 4 contained genes that experienced a decrease in expression at 6 hpi, an increase back to baseline levels at 12 hpi, and a subsequent progressive decrease at 24 and 120 hpi. Biological process terms significantly enriched in this group included translation, auxin activated signaling pathway, lignin catabolic process, plant cell wall organization, and phloem/xylem histogenesis. Lastly, the expression of genes in cluster 5 was equivalent at 0 and 6 hpi, noticeably increased at 12 hpi, and decreased to well below baseline levels at 120 hpi. Enriched BP terms in this cluster included gibberellic acid-mediated signaling pathway, response to auxin, several terms related to oxidative stress and, interestingly, many terms related to growth and metabolism. The results of GO enrichment analyses for P. americana genes at individual time points are also available for both Dusa® and R0.12 (Supplementary Figures 1, 2; Supplementary Tables 12–15).
The member genes, expression profiles, and GO profile of each cluster in both rootstocks were also compared to one another. Two clusters from each rootstock were highly similar, cluster 2 from both Dusa® and R0.12 with 2,015 shared genes and clusters 5 from Dusa® and 3 from R0.12 that shared 1,652 genes (Supplementary Table 16). The expression of genes in cluster 2 from both rootstocks peaked at 6 hpi, although in R0.12, the expression returned to mock-inoculated levels at 12 hpi, but in Dusa® the expression only returned to mock-inoculated levels at 120 hpi. Interestingly, the number of shared genes between clusters 2 from R0.12 and 4 from Dusa® was also substantial, but in Dusa®, the average gene expression rose steadily from 6- to 24 hpi. Likely the most similar were clusters 5 from Dusa® and 3 from R0.12, in which the gene expression remained mostly unchanged until 120 hpi. Additionally, in both clusters, enriched BP terms were near identical. Cluster 1 from Dusa®seemed most like cluster 4 from R0.12 because of the similarity of some BP terms, expression at 0, 6, and 120 hpi, and shared member gene identity. Still, some similarities were seen when comparing cluster 1 from Dusa® to cluster 5 from R0.12. However, the expression of genes from cluster 1 in Dusa® only decreased over time, whereas in R0.12, the expression peaked at 12 hpi and decreased thereafter. Lastly, cluster 3 in Dusa® seemed to resemble both clusters 1 and 4 in R0.12 to some degree. In all three of these clusters, the gene expression decreased at 6 hpi while remaining low at 12 and 24 hpi in Dusa® but not in R0.12.
P. cinnamomi
Temporal cluster analyses were performed to separate P. cinnamomi genes based on their expression in in vitro cultured mycelia, and over time following inoculation of a partially resistant (Dusa®) and a susceptible (R0.12) P. americana rootstock (Figures 6A, 7A; Supplementary Table 17). Following inoculation of Dusa®, the expression of P. cinnamomi genes belonging to cluster 1 remained similar to that of mycelia (0 hpi) at early time-points, while at 120 hpi, the expression increased considerably. Overrepresented BP terms resulting from genes in this cluster included transmembrane transport, modulation by symbiont of host defense-related PCD, and several terms related to metabolism (Figure 6B; Supplementary Table 18). The expression of genes in cluster 2 showed a single peak at 24 hpi, and only one enriched BP term was obtained for this cluster, glutamine metabolic process. Genes in cluster 3 showed lower expression at 6, 12, and 24 hpi when compared to mycelia, and at 120 hpi, and biological process terms associated with genes from this cluster were mostly related to translation. The expression of genes in cluster 4 was higher at 6 hpi than at other time points or that of mycelia; however, no enriched BP terms were recovered for this cluster. Cluster 5 contained genes that were expressed at much lower levels in planta than in mycelia, and BP terms associated with this cluster included autophagy, transcription, phospholipid biosynthesis, and positive regulation of ATPase activity.
Figure 6. Temporal clustering and Gene Ontology enrichment of Phytophthora cinnamomi differentially expressed genes (DEGs) following inoculation of the Persea americana rootstock Dusa®. (A) Sample libraries from in vitro mycelial controls and the partially resistant P. americana rootstock Dusa®inoculated with P. cinnamomi were used in clustering analyses. DEG Z-scores, based on reads per kilobase million (RPKM) normalized read counts, were calculated and used for clustering analyses using the fuzzy c-means (Futschik and Carlisle, 2005) in the R package TCseqv1.16.0 (Wu and Gu, 2021). Inoculated samples included 6, 12, 24, and 120 h post inoculation (hpi) and are indicated on the x-axis, with the mock-inoculated samples indicated as 0 hpi. Membership scores produced by the fuzzy c-means algorithm are color-coded and indicated to the right of each cluster. (B) Member genes from each cluster were subjected to Gene Ontology (GO) enrichment and significantly overrepresented GO terms, with adjusted p-value < 0.05 Benjamini-Hochberg false discovery rate, were reduced in REVIGO (Supek et al., 2011). The cluster from which the terms originate is indicated on the left most side of the Figure, followed by the biological process GO term. The overrepresented p-value of each GO term is indicated along the x-axis as –log10 (p-adjusted).
Figure 7. Temporal clustering and Gene Ontology enrichment of Phytophthora cinnamomi differentially expressed genes (DEGs) following inoculation of the Persea americana rootstock R0.12. (A) Sample libraries from in vitro mycelial controls and the susceptible P. americana rootstock R0.12 inoculated with P. cinnamomi were used in clustering analyses. DEG Z-scores, based on reads per kilobase million (RPKM) normalized read counts, were calculated and used for clustering analyses using the fuzzy c-means (Futschik and Carlisle, 2005) in the R package TCseqv1.16.0 (Wu and Gu, 2021). Inoculated samples included 6, 12, 24, and 120 h post inoculation (hpi) and are indicated on the x-axis, with the mock-inoculated samples indicated as 0 hpi. Membership scores produced by the fuzzy c-means algorithm are color-coded and indicated to the right of each cluster. (B) Member genes from each cluster were subjected to gene ontology (GO) enrichment, and significantly overrepresented GO terms, with adjusted p-value < 0.05 Benjamini-Hochberg false discovery rate, were reduced in REVIGO (Supek et al., 2011). The cluster from which the terms originate is indicated on the left most side of the Figure, followed by the biological process GO term. The overrepresented p-value of each GO term is indicated along the x-axis as –log10 (p-adjusted).
Temporal cluster analyses of R0.12 yielded five clusters with profiles similar to those seen in Dusa®. In R0.12, cluster 1 included genes that are expressed highest in the in vitro mycelia samples rather than in planta (Figure 7A; Supplementary Table 17). The only enriched BP term obtained from genes in this cluster was pyruvate metabolic process (Figure 7B; Supplementary Table 19). Cluster 2 included genes that displayed high expression at 24 hpi, while at all the other time points, levels were equivalent. Overrepresented BP terms associated with genes from this cluster included rRNA processing, protein folding, and several terms related to translation and metabolism. Genes that formed part of cluster 4 were expressed at much higher levels at 120 hpi than at any other time-point or in mycelia. Enriched BP terms resulting from genes in cluster 4 included transmembrane transport and several terms related to metabolism. The expression of genes in cluster 5 was substantially lower at 6 hpi when compared to mycelia; at 12 hpi, levels were equivalent to mycelia followed by slightly lower levels at 24 and 120 hpi. The only BP term enriched using genes in this cluster was translation. The results of GO enrichment analyses for P. cinnamomi genes at individual time points following inoculation of both Dusa® and R0.12 are also available (Supplementary Figures 3, 4; Supplementary Tables 20–23)
The member genes, expression profiles, and GO terms for each P. cinnamomi temporal cluster was compared between the rootstocks (Supplementary Table 24). Although near-identical expression profiles could be identified for each cluster from either rootstock, member gene lists were only highly similar for two clusters from each rootstock, cluster 1 from Dusa® and cluster 4 from R0.12, and cluster 5 from Dusa® and cluster 1 from R0.12, which shared 1,120 and 1595 genes, respectively. Cluster 3 in Dusa® shared 471 genes with cluster 5 in R0.12; they also shared enriched BP terms related to translation. Meanwhile, cluster 4 in Dusa® and cluster 3 in R0.12 shared 552 genes, neither of which contained any enriched BP terms. Cluster 2 from both rootstocks contained genes that were expressed highest at 24 hpi, although they only shared 231 genes. For perspective, cluster 2 from R0.12 shared a similar number of genes with clusters 1, 3, and 5 from Dusa®.
P. cinnamomi Virulence Factors
Following cluster analyses, the corresponding predicted proteins from P. cinnamomi genes in each cluster were compared to the PHI database to identify important virulence determinants at defined time points (Supplementary Tables 25, 26). On average, 33% of P. cinnamomi proteins from described clusters, following the inoculation of Dusa® and R0.12, matched to the PHI database sequences with an E-value < 1.0e−5 (Table 1). Overall, 19 different types of effectors were identified, 10 of which were well-documented in Phytophthora species; RxLR effector proteins were predominant, followed by elicitins, necrosis, and ethylene-inducing peptide 1 (Nep1)-like proteins (NLPs),) glucanase inhibitor proteins (GIPs), and EPIs, among others (Table 2). Of the 40 unique RxLRs identified, only 16 were identified in both rootstocks. A similar observation was made for all identified Phytophthora effectors except for XEG1 and Phytophthora cactorum-Fragaria (PcF), indicating a distinctive compliment of effectors was used during the infection of each rootstock.
Table 1. Summary statistics following comparison of temporal cluster member protein sequences to the Pathogen Host Interactions (PHI) database.
Table 2. Number and types of known Phytophthora effectors identified following the comparison of temporal cluster member protein sequences to the Pathogen-Host Interactions (PHI) database.
When only considering matches to the effector proteins from the PHI database, cluster 4 from Dusa® and cluster 3 from R0.12, which both peak at 6 hpi, contain several notable matches. This comparison identified several P. cinnamomi genes, following the inoculation of either rootstock, which encoded for proteins with similarity to CBEL and RxLR effector proteins from other Phytophthora species (Supplementary Tables 25, 26). In both rootstocks, more unique predicted CBELs, four in Dusa® and three in R0.12, were identified in the 6 hpi peak clusters than in any of the other clusters. Contrastingly, the lowest number of predicted proteins with similarity to RxLRs was identified in these clusters as compared to the 24 and 120 hpi peak clusters. Here, only one predicted RxLR was identified in R0.12 (jg4048), whereas seven were identified in Dusa®. Interestingly, in inoculated Dusa®, there were also matches to a crinkling and necrosis (CRN) protein from Phytophthora sojae. In R0.12, similarities to an elicitin from Phytophthora infestans and an NLP from Phytophthora capsici were found.
The expression of P. cinnamomi genes in cluster 2 from both Dusa® and R0.12 peaked at 24 hpi. In these clusters, as might be expected, far fewer predicted proteins with similarity to CBEL were identified. However, following the inoculation of both rootstocks, several genes that code for proteins with similarity to GIP2 in P. sojae and EPIs from P. infestans were identified. It was noted, however, that the compliment of EPIs and GIPs present in either rootstock was distinctly different, with less than half of the proteins being shared between clusters. Further differences between the rootstocks included the presence of three proteins with similarity to elicitins that were present in R0.12 but not in Dusa®. Likewise, 19 proteins with similarity to RxLRs from P. sojae and P. infestans were identified in R0.12, while only nine were identified in Dusa®, contrasting with the observation made at 6 hpi. Three predicted NLP-like proteins were also identified following the inoculation of R0.12, while only one was identified in Dusa®. Surprisingly, the predicted protein sequence jg108, which showed similarity to the phytotoxic protein PcF, was found in cluster 2 of R0.12 but not in Dusa®.
Otherwise, the repertoire of effector proteins identified in cluster 4 from R0.12, in which the gene expression peaked at 120 hpi, was significantly less diverse when compared to R0.12 cluster 2; cluster 4 contained only six predicted RxLRs, three EPIs, one GIP, and four XEGs. In the analogous cluster from Dusa®, cluster 1, there were 14 predicted RxLRs, six EPIs, seven GIPs, and six XEGs. Elicitins were also identified in Dusa® for the first time, totaling 12, whereas seven were identified in R0.12. In Dusa® and R0.12, the total number of predicted NLP effector proteins also increased to nine and six, respectively. Notably, jg1926, a predicted NLP protein with high similarity to necrosis-inducing Phytophthora protein 1 (NPP1) from Phytophthora parasitica, was identified in the cluster from R0.12 but not in Dusa®. Intriguingly, the same predicted CBEL, jg11698, was identified in both Dusa® and R0.12 at this time point. In contrast to the observations made on clusters 2 from both rootstocks, PcF was identified in Dusa® cluster 1 and not in R0.12 cluster 4.
P. americana Defense Responses
Mercator was used to assign the predicted P. americana proteins to MapMan bins. Mercator3 successfully annotated 54.2% of the 47,687 predicted P. americana protein sequences, assigning 32.39% to the bins and occupying 94.25% of the available bins (Supplementary Table 27). Mercator4 successfully annotated 47.52% of the predicted P. americana predicted sequences, assigning 33.59% to the bins and occupying 94.34% of the available bins (Supplementary Table 28). Mercator 3 annotation conducted for further investigations because of the higher number of annotated P. americana proteins, 25,845 vs. 22,661. DEGs were then visualized in MapMan to assist in pathway identification and comparison.
However, because of the significant contrast in the number of DEGs between Dusa® and R0.12 at 12 and 24 hpi (13,501 and 818, respectively), the list of top 100 variably expressed genes was first used to filter potential candidate genes and pathways for further investigation (Supplementary Table 29). To further reduce this list, only genes known to play a direct role in plant defense or stress responses were kept, while genes involved in biological processes, e.g., photosynthesis and electron transport, were omitted. The filtered list included 35 genes involved in the biosynthesis and regulation of some phytohormones, PR protein coding genes, protease inhibitors, proteases, several extensins, and a wall-associated kinase (WAK), among others (Table 3). Given the limited knowledge of receptor-like kinases and apoplastic proteases in P. americana, we decided to focus the further analysis on these two groups of proteins.
Table 3. Filtered list of the most variably expressed defense- and stress-related genes between R0.12 and Dusa® following P. cinnamomi inoculation.
Receptor-Like Kinases
When visualized using MapMan, clear differences in the number of DE receptor-like kinases (RLKs) were apparent, even at 6 and 120 hpi, where R0.12 and Dusa® displayed less differences overall (Table 4). In Dusa®, far more L-type lectin (L-lectin), thaumatin-like protein kinase (TLPK), receptor-like kinase in flowers 3 (RFK3-like), crinkly4-like (CR4-like), proline-rich extensin-like receptor protein kinase (PERK-like), and leaf rust 10 disease-resistance locus receptor-like protein kinase (LRK10-like) coding genes were significantly downregulated at 6, 12, and 24 hpi than those that were upregulated. While at 120 hpi the opposite was observed - far more were significantly upregulated for each than were downregulated. A related distinction could be made between 6 and 120 hpi for L-lectin, TLPK, RFK3, CR4-like, PERK-like, and LRK10-like coding genes in R0.12. However, fewer of these genes, except for CR4-like, were significantly upregulated in R0.12 at 120 hpi when compared to Dusa®.
Table 4. Number of predicted receptor-like kinase (RLK) family proteins encoded for by genes that are significantly up- and downregulated in Dusa® or R0.12 following inoculation with Phytophthora cinnamomi.
Domain of unknown function 26 (DUF26) and S-domain RLKs were among the largest families of identified RLKs in P. americana (Table 4). A substantial portion of these were significantly upregulated throughout all the time-points in Dusa®, following P. cinnamomi inoculation, and at 6 and 120 hpi in R0.12. Interestingly, more DUF26 and S-domain RLKs were significantly upregulated in R0.12 at 6 hpi when compared to Dusa®, while the opposite was observed at 120 hpi. Another large family of identified RLKs were wall-associated kinases (WAKs); of these, many were significantly upregulated at 6 hpi in both Dusa® and R0.12, while a smaller but notable proportion was downregulated. In Dusa®, the downregulated portion of WAKs decreased steadily over time, while the upregulated portion decreased at 12 and again at 24 hpi followed by an increase at 120 hpi. Interestingly, the vast majority of DE extensins were downregulated in both rootstocks.
Leucine-rich repeat RLKs were, by far, the largest family of identified RLKs in P. americana (Table 4). In contrast to the RLKs reported earlier, a large percentage LRR-RLKs, specifically of types I, VIII, and XII were significantly upregulated in Dusa® at the earliest time point but decreased thereafter. Similarly, a large proportion of these LRR-RLKs were upregulated at 6 hpi in R0.12, while less was observed at 120 hpi. The share of significantly upregulated type X and XIII LRR-RLKs was also high at 6 hpi but low at 120 hpi in both rootstocks. However, in Dusa®, the largest percentage of upregulated type X and XIII LRR-RLKs were seen at 24 hpi. Opposingly, the majority of DE types III, VII, and XI LRR-RLKs were downregulated at all the time-points in Dusa® and at 6 and 120 hpi in R0.12. The downregulated portion of DE types II and V LRR-RLKs remained consistent at all the time points in Dusa®, while the upregulated portion of these RLKs peaked at 12 hpi. In R0.12, a larger percentage of type II was downregulated at 6 and 120 hpi when compared to Dusa®. Lastly, the majority of DE types IV and VI LRR-RLKs were significantly downregulated at 120 hpi in both rootstocks. A small percentage of type VI was upregulated at 24 and 120 hpi in Dusa® and at 120 hpi in R0.12. It is worth noting, however, that at 12 and 24 hpi, there was almost no up- or downregulation of these or any other putative RLK coding genes in R0.12. If not for the 12 and 24 hpi time points, very few differences would be observed between R0.12 and Dusa®.
Apoplastic Proteases
Predicted proteases were extracted from the Mercator 3 annotation file. A total of 62 aspartic proteases, 77 cysteine proteases, 50 metalloproteases, and 229 serine proteases were recovered (Supplementary Table 30). Of the 229 serine proteases identified, 108 were subtilases (SBTs) and separated during subsequent analyses. This list was filtered to further exclude any proteases that did not show significant differential expression at any time or in any rootstock when compared to the appropriate mock-inoculated control. The subcellular localization of these DE proteases was predicted using an online web service, available at (http://embed.protein.properties). Twenty-five aspartic proteases were DE, of which 20 were predicted to localize to the extracellular space and two to the cell membrane. Interestingly, only 18 of the predicted cysteine proteases were DE, of which only three were expected to localize extra-cellularly. Altogether, 23 metalloproteases were DE; however, only five and two were predicted to localize to the extracellular space and the cell membrane, respectively. Forty-one non-SBT serine proteases were DE, of which 18 were expected in the extracellular space, and only one was predicted to localize to the cell membrane. In total, 37 SBTs were DE; surprisingly all the 37 were predicted to exist in the extracellular space. Only proteases predicted to localize n the extracellular space or in the cell membrane, henceforth collectively referred to as apoplastic proteases, were retained for further investigation.
Only a small number of genes encoding for apoplastic aspartic proteases were significantly upregulated in either rootstock at any given time (Figure 8). Three groups were assigned based on the clustering of genes on the dendrogram obtained from Dusa®. In group 1, genes were generally upregulated following P. cinnamomi challenge. One, Peame105C05g001950, was significantly upregulated at the three earliest time points in Dusa® when compared to the mock-inoculated control. However, in R0.12, this gene was not upregulated at any time-point but instead tended toward downregulation, although not significant, at 12 hpi. The remaining two were upregulated at 6 and 120 hpi in both Dusa® and R0.12. Most of the genes in group 2 were downregulated at 6-, 12, and 24 hpi in Dusa® and 6 hpi in R0.12. Three of these were significantly upregulated at 120 hpi in Dusa®, while only one reached significance in R0.12. Genes that fall in group 3 were generally downregulated at most time points in Dusa® and at 6 and 120 hpi in R0.12. Interestingly, two of these, Peame105C07g010840 and Peame105C03g002900, were significantly upregulated in R0.12 at 6 and 12 hpi, respectively.
Figure 8. Heatmap and dendrogram illustrating the expression of genes encoding for apoplastic aspartic proteases in Persea americana following inoculation with Phytophthora cinnamomi. Differentially expressed genes (DEGs) encoding for apoplastic aspartic proteases from both partially resistant (Dusa®) and susceptible (R0.12) P. americana rootstocks are indicated. Differential expression at 6, 12, 24, and 120 h post inoculation (hpi), following inoculation with P. cinnamomi, was calculated using mock-inoculated (dH2O) sample libraries as the control. Significantly upregulated [log2(fold change) > 1, p = 0.05] genes are indicated by a red color, while significantly downregulated [log2(fold change) < −1, p = 0.05] genes are indicated in blue, the intensity of which represents scale. Significant DEGs are indicated with a black dot. The expression of aspartic protease coding genes was broadly clustered into three groups based on the dendrogram from Dusa® and were indicated numerically on the dendrogram. Only aspartic proteases that were differentially expressed in one or more samples and were predicted to localize to the cell membrane or extracellular space are shown.
Cysteine proteases and metalloproteases represented the two smallest groups of DE apoplastic proteases identified in this study. Of the cysteine proteases described, two were significantly upregulated at 120 hpi in Dusa®, while only one reached significance in R0.12 (Figure 9A). Intriguingly, Peame105C08g003910 was significantly downregulated at 6 and 12 hpi in R0.12 but not in Dusa®. The heatmap describing the metalloproteases was divided into two groups based on the dendrogram obtained from Dusa® (Figure 9B). Metalloproteases in the first group contained four genes that were downregulated at 6 and 12 hpi in Dusa®, two of which were also downregulated at 24 hpi. The same two were downregulated in R0.12 but only at 6 hpi, while the remaining two were not DE at 6, 12, or 24 hpi. Conversely, all four were significantly upregulated in both Dusa® and R0.12 at 120 hpi. The second group of metalloproteases was mostly downregulated at various early time points in Dusa® and R0.12; one of these was also downregulated at 120 hpi in R0.12. Based on the avocado genome annotation, all except one DE apoplastic metalloproteinase, Peame105C02g015270, were predicted to encode for matrix metalloproteases (MMPs).
Figure 9. Heatmap and dendrogram illustrating the expression of genes encoding for apoplastic cysteine proteases and metalloproteases in Persea americana following inoculation with Phytophthora cinnamomi. Differentially expressed genes (DEGs) encoding for (A) apoplastic cysteine proteases and (B) apoplastic metalloproteases from both partially resistant (Dusa®) and susceptible (R0.12) P. americana rootstocks are indicated. Differential expression at 6, 12, 24, and 120 h post inoculation (hpi), following inoculation with P. cinnamomi, was calculated using mock-inoculated (dH2O) sample libraries as the control. Significantly upregulated [log2(fold change) > 1, p = 0.05] genes are indicated by a red color, while significantly downregulated [log2(fold change) < −1, p = 0.05] genes are indicated in blue, the intensity of which represents scale. Significant DEGs are indicated with a black dot. The expression of metalloprotease coding genes was broadly clustered into two groups based on the dendrogram from Dusa® and were indicated numerically on the dendrogram. Only cysteine proteases and metalloproteases that were differentially expressed in one or more samples and were predicted to localize to the cell membrane or extracellular space are shown.
The heatmaps depicting the expression of genes encoding for non-SBT serine proteases could be divided into three groups based on the dendrogram from Dusa® (Figure 10). The avocado genome annotation classified all but one, Peame105C06g005480, as encoding for serine carboxypeptidase-like proteins (SCLPs). Genes in the first group were not DE at 6, 12, or 24 hpi in either Dusa® or R0.12. However, all except one from each rootstock were significantly downregulated at 120 hpi. In comparison, genes from group 2 were significantly downregulated during various early time points in Dusa® and at 6 hpi in R0.12. However, neither rootstock displayed up- or downregulation for any group 2 apoplastic serine proteases at 120 hpi. For group 3, there were two genes that were significantly upregulated at 120 hpi in both rootstocks but none at any other time point. Three of the remaining group 3 serine protease coding genes were significantly upregulated in Dusa® at 12 hpi, one of which was significantly downregulated at 120 hpi. Only one of these, Peame105C03g003010, was significantly upregulated at an early time-point in R0.12. The last group 3 serine protease, Peame105C04g026050, was upregulated at 6 and 24 hpi in Dusa® while only being downregulated in R0.12 at 120 hpi. Although the same gene was inclined toward downregulation in Dusa® at 120 hpi, it did not reach significance.
Figure 10. Heatmap and dendrogram illustrating the expression of genes encoding for apoplastic serine proteases in Persea americana following inoculation with Phytophthora cinnamomi. Differentially expressed genes (DEGs) encoding for apoplastic serine proteases from both partially resistant (Dusa®) and susceptible (R0.12) P. americana rootstocks are indicated. Differential expression at 6, 12, 24, and 120 h post inoculation (hpi), following inoculation with P. cinnamomi, was calculated using mock-inoculated (dH2O) sample libraries as the control. Significantly upregulated [log2(fold change) > 1, p = 0.05] genes are indicated by a red color, while significantly downregulated [log2(fold change) < −1, p = 0.05] genes are indicated in blue, the intensity of which represents scale. Significant DEGs are indicated with a black dot. The expression of serine protease coding genes was broadly clustered into three groups based on the dendrogram from Dusa® and were indicated numerically on the dendrogram. Only serine proteases that were differentially expressed in one or more samples and were predicted to localize to the cell membrane or extracellular space are shown.
Surprisingly, all significant DE SBTs were predicted to exist in the extracellular space; as such, they constitute the largest group of apoplastic proteases identified in P. americana. The expression of SBT coding genes could be broadly clustered into four groups based on the dendrogram from Dusa® (Figure 11). The genes found in group 1 from Dusa® were not DE at 6, 12, or 24 hpi; however, the majority was significantly downregulated at 120 hpi. In R0.12, the expression of genes in group 1 was largely comparable to Dusa®, except for three. The first, Peame105C009g008400, was significantly downregulated at 6 and 12 hpi; the second, Peame105C09g024080, was significantly upregulated at 12 hpi; the third, Peame105C08g006250, tended toward upregulation at 6, 12, and 24 hpi, achieving significance at 120 hpi. The second group of SBTs was generally downregulated at most time-points in Dusa® and at 6 and 120 hpi in R0.12. In contrast, group 3 contained genes that were significantly upregulated in Dusa® at several early time-points, especially 6 hpi. In R0.12, only two of these genes were significantly upregulated and only at 6 hpi. Peame105C01g045590 was particularly interesting, as it was significantly upregulated in Dusa® at all the time-points; however, in R0.12 it was not DE at any. Based on Blastp, the predicted protein coding sequence for Peame105C01g045590 and Peame105C01g045560 closely resembled that of SBT3.3 in Cinnamomum micranthum. Lastly, group 4 consisted primarily of genes that were downregulated at 6 hpi in either Dusa® or R0.12 but not at any of the subsequent time points. Group 4 also contained one gene that was upregulated at 6 and 120 hpi in both rootstocks. Interestingly, the last member of group 4, Peam105C08g006080, was upregulated at 24 and 120 hpi in Dusa® but tended toward downregulation in R0.12. Here, the Blastp analysis revealed that the predicted protein sequence of Peam105C08g006080 showed high similarity to SBT3 like proteins in C. micranthum.
Figure 11. Heatmap and dendrogram illustrating the expression of genes encoding for apoplastic subtilases (SBTs) in Persea americana following inoculation with Phytophthora cinnamomi. Differentially expressed genes (DEGs) encoding for apoplastic SBTs from both partially resistant (Dusa®) and susceptible (R0.12) P. americana rootstocks are indicated. Differential expression at 6, 12, 24, and 120 h post inoculation (hpi), following inoculation with P. cinnamomi, was calculated using mock-inoculated (dH2O) sample libraries as the control. Significantly upregulated [log2(fold change) > 1, p = 0.05] genes are indicated by a red color, while significantly downregulated [log2(fold change) < −1, p = 0.05] genes are indicated in blue, the intensity of which represents scale. Significant DEGs are indicated with a black dot. The expression of SBT coding genes were broadly clustered into four groups based on the dendrogram from Dusa® and were indicated numerically on the dendrogram. Only SBTs that were differentially expressed in one or more samples and were predicted to localize to the cell membrane or extracellular space are shown.
Discussion
The molecular basis of P. cinnamomi resistance in avocado has been the topic of a recent discussion (van den Berg et al., 2021). This review provided an up-to-date justification and classification of avocado rootstocks as resistant, partially resistant, tolerant, or susceptible to P. cinnamomi. The current study attempted to contextualize the basis of partial resistance by temporal transcriptomic analyses of representatively susceptible (R0.12) and partially resistant (Dusa®) avocado rootstocks inoculated with P. cinnamomi (Engelbrecht and van den Berg, 2013; van den Berg et al., 2018b). Using dual RNA sequencing technology, this study provides a holistic overview of both avocado and P. cinnamomi gene expression over time, during their interaction. Notably, sample libraries were collected for the characterization of gene expression at both early, 6, 12, and 24 hpi, and late, 120 hpi, responses. As expected, the relative abundance of transcripts that mapped back to the P. cinnamomi genome increased substantially between the early and late sample libraries in both rootstocks, indicating host colonization and spread of P. cinnamomi. Contrastingly, the number of avocado transcripts that mapped back to the P. americana genome decreased substantially at the late time points. This last observation was somewhat expected, given that root necrosis is often observed following P. cinnamomi inoculation of avocado rootstocks at later time points (Engelbrecht and van den Berg, 2013; Engelbrecht et al., 2013). Supportively, a previous study provided evidence that P. cinnamomi switches to a necrotrophic life stage at around 24 hpi in the rootstock Dusa®(van den Berg et al., 2018b).
In this study, the initial analyses revealed overwhelmingly significant differences in the number of avocado DEGs between rootstocks, specifically at 12 and 24 hpi. At these time points, the number of DEGs in the susceptible rootstock paled in comparison to DEGs in the partially resistant rootstock. In fact, in the susceptible rootstock, slightly more than 800 DEGs were identified when combining both time points, while over 13,000 were found in Dusa®. In contrast, at 6 and 120 hpi, both rootstocks displayed a similar number of up- or downregulated genes. The PCAs and hierarchical clustering of avocado DEGs further supported the variation in the transcriptomes of R0.12 and Dusa® at 12 and 24 hpi. In fact, the R0.12 12 and 24 hpi sample libraries were indistinguishable from the mock-inoculated sample libraries following PCA, while a clear difference was seen among mock, early, and late responses in Dusa®. Together, these observations provide the first indication of major transcriptomic differences between Dusa® and R0.12, specifically at 12 and 24 hpi, following inoculation with P. cinnamomi.
Following temporal cluster analyses of avocado DEGs, GO enrichment confirmed the activation of pathways related to defense responses in the gene sets from clusters that were upregulated at 6 hpi in both Dusa® and R0.12. Surprisingly, however, in R0.12, the expression of these genes seemed to revert to near mock-inoculated levels at 12 and 24 hpi. Instead, the cluster analyses and subsequent GO enrichment revealed that R0.12 was activating pathways related to growth and development at 24 hpi. In contrast, similar pathways were identified following the GO enrichment of the cluster representing genes that were downregulated at all post-inoculation time points in Dusa®. Instead, in Dusa®, the activation of pathways related to defense continued at 12 and 24 hpi. Clusters representing gene sets that were upregulated at 120 hpi in both Dusa® and R0.12 gave rise to similarly enriched GO terms. Here, the pathways related to phosphorylation were the most significantly enriched in both rootstocks. Notably, phosphorylation is known to play an essential role in plant immune responses, constituting a requisite part of PRR-mediated defense response signaling (Park et al., 2012). Otherwise, common defense-related terms such as response to chitin, oxidative stress, and wounding, were also enriched in these clusters. Together, this evidence points toward a high likelihood that most of the variation in resistance between R0.12 and Dusa® lies in the early response to P. cinnamomi challenge.
P. cinnamomi
Interestingly, the increased expression of P. cinnamomi elicitin homologs, a class of effectors known for inducing HR (Derevnina et al., 2016), occurred far earlier in R0.12 than in Dusa® (Table 2). Similarly, the P. cinnamomi homologs of the phytotoxic P.cactorum-Fragaria (PcF) protein and xyloglucan-specific endo-beta-1,4-glucanase 1 (XEG1) were expressed at 24 hpi in R0.12 but not in Dusa®. Notably, XEGs are an important group of cell wall-degrading enzymes (CWDEs); Phytophthora sojae XEG1 (PsXEG1) was shown to be an important determinant of virulence in infected soybean (Ma et al., 2015, 2017; Xia et al., 2020). Meanwhile, PcF application was revealed to cause necrosis in both tomato and strawberry (Orsomando et al., 2001). Together, these observations suggest that the necrotrophic life stage of P. cinnamomi may initiate earlier in R0.12 following inoculation when compared to Dusa®. This hypothesis should constitute the focus of a separate study given the possible implications for the industry.
The effector proteins used by P. cinnamomi showed further variance when comparing Dusa® and R0.12, exemplified by the upregulation of putative RxLR effector coding genes. RxLRs are an important class of cytoplasmic effectors in Phytophthora species and the subject of extensive study (Win et al., 2007; Jiang et al., 2008; McGowan and Fitzpatrick, 2017; Reitmann et al., 2017; Hardham and Blackman, 2018; Engelbrecht et al., 2021; Joubert et al., 2021). Although a large proportion of the upregulated RxLR coding genes in either Dusa® or R0.12 were shared between the rootstocks, each was host to a unique compliment. Interestingly, in R0.12, the largest compliment of upregulated RxLR coding genes was found at 24 hpi. Meanwhile, the number and diversity of upregulated RxLR coding genes in Dusa® increased steadily as the infection progressed. Importantly, Phytophthora RxLRs are known to play a fundamental role in the evasion of host defense signaling and immune responses (Kelley et al., 2010; Vetukuri et al., 2017; Dalio et al., 2018; Naveed et al., 2019). Thus, it is reasonable to assume that the unique compliment of upregulated RxLR coding genes in R0.12 is at least partly responsible for the inappropriate response seen in this rootstock at 12 and 24 hpi. However, the upregulated RxLR coding genes in Dusa® should not be overlooked given that this rootstock is not fully resistant to P. cinnamomi. For more details on the expression of the putative P. cinnamomi RxLR effector coding gene repertoire, using some of the data from this study, readers are referred to the study performed by Joubert et al. (2021).
Although fewer in number, several other important groups of effectors were upregulated in both Dusa® and R0.12, and showed a varied expression between the rootstocks over time. The first of these, the NLP effectors, represent a large family of microbial effectors that contribute significantly to virulence and host defense response elicitation (Pemberton and Salmond, 2004; Seidl and van den Ackerveken, 2019). Intriguingly, only half of the upregulated putative NLPs were shared between Dusa® and R0.12. Additionally, more upregulated NLP coding genes were identified in R0.12 at 6 and 24 hpi than in Dusa®, while more were identified in Dusa® at 120 hpi. Together, these observations would suggest that the NLPs likely contribute to differences in the observed virulence of P. cinnamomi on Dusa® and R0.12.
Entirely unique compliments of upregulated CRN effector coding genes were seen when comparing P. cinnamomi gene expression in Dusa® and R0.12 (Table 2). CRNs are a large, functionally diverse family of modular cytoplasmic effectors that have been described in a wide array of plant-pathogenic oomycetes (Tyler et al., 2006; Kamoun, 2007; Haas et al., 2009; Schornack et al., 2010; Stam et al., 2013; McGowan and Fitzpatrick, 2020). These effectors were first identified in Phytophthora infestans for their ability to induce host leaf-crinkling and necrosis as well as induce defense responses (Torto et al., 2003). Since then, CRNs have also been implicated in the suppression of host defense responses (Chen et al., 2013; Mafurah et al., 2015; Xiang et al., 2021). It is worth noting that 45 CRNs were previously identified from the P. cinnamomi genome (Engelbrecht et al., 2021). However, only five upregulated CRNs were identified in the current study using the PHI database. Given the apparent importance of this group of effectors, it seems unlikely that the current study accurately described the full extent of CRN expression in P. cinnamomi. Nonetheless, differences were still seen in the number and temporal association of expressed CRNs, warranting further fine-grained analyses and characterization.
Interestingly, substantially more upreguated GIPs were observed in Dusa® at 120 hpi than in R0.12. In studied Phytophthora spp., GIPs bind to and inhibit plant secreted endo-β-1,3-glucanases (EGases; Damasceno et al., 2008). Notably, EGases act on an abundant component of Phytophthora cell walls, β-1,3 glucan, directly damaging the pathogen while releasing PAMPs known to elicit host defense responses (Erwin et al., 1983; van Loon et al., 2006). Markedly, β-1,3-glucanase has been linked to increased P. cinnamomi resistance in the partially resistant rootstock R0.06 (van den Berg et al., 2018a). Thus, it is possible that the EGases produced by Dusa® and not R0.12 might have led to the consequent upregulation of P. cinnamomi GIPs at later time points, to directly counteract cell wall damage and the release of related PAMPs.
Similarly, more upregulated EPI effector coding genes were seen at both 24 and 120 hpi in P. cinnamomi-inoculated Dusa® than in R0.12. Notably, EPI effectors are known to inhibit plant-derived extracellular proteases. For example, the P. infestans effectors EPI1 and EPI10 interact with and inhibit the apoplastic serine protease P69B (Tian et al., 2004, 2005). Similarly, the homolog from Phytophthora palmivora, EPI10, inhibits the extracellular subtilase HbSPA from Hevea brasiliensis (Ekchaweng et al., 2017). Inhibition of plant proteases seems to be an advantageous strategy that has developed in diverse plant pathogenic species such as Cladosporium fulvum, which secretes the effector Avr2 (Rooney et al., 2005). This effector has been shown to inhibit several apoplastic cysteine proteases in tomato, including Required for Cladosporium resistance 3 (Rcr3) and Phytophthora-inhibited protease 1 (Rooney et al., 2005; Shabab et al., 2008; van Esse et al., 2008; Pip1). Similarly, the P. infestans cystatin-like EPIs EPIC1 and EPIC2B are also known to target apoplastic proteases, including Rcr3 and Pip1 (Tian et al., 2007; Song et al., 2009). It would be reasonable to assume, given the larger compliment of upregulated EPIs observed following the inoculation of Dusa®, that P. cinnamomi might be responding to a larger or more devastating set of secreted proteases from the partially resistant rootstock.
Finally, several CBELs were upregulated in both Dusa® and R0.12, most prominently at 6 hpi. This is in keeping with the primary function of CBELs, which allow Phytophthora mycelium to adhere to cellulosic structures such as the cell wall of plants (Mateos et al., 1997; Gaulin et al., 2002). However, CBELs are also known to elicit a strong host defense response, acting like PAMPs (Gaulin et al., 2006). As such, there has been some interest in production of recombinant CBEL proteins for agricultural applications as activators of plant immune responses (Larroque et al., 2011). Thus, the CBELs identified here might be deemed novel targets for future research aimed at understanding the early determinants of P. cinnamomi virulence in avocado and increasing resistance to this devastating pathogen.
P. americana
Given that PRRs constitute the starting point for inducible defense responses in plants, it is conceivable that the overwhelming lack of defense responses at 12 and 24 hpi in R0.12 might stem from inappropriate PRR signaling or activation. To this end and given the prevalence of GO terms related to phosphorylation in both Dusa® and R0.12, we sought to characterize differences in the expression of major classes of RLKs in avocado. When considering the representative percentages of up- or downregulated RLK families at 6 and 120 hpi in Dusa® and R0.12, very few substantiative differences were noted between the rootstocks. This fact limited the usefulness of a direct comparison between the rootstocks, given the broad approach used in this study. However, when considering only the DE of RLKs in Dusa®, some interesting patterns emerged. For instance, some RLKs, such as those encoding for extensins and LRR-RLKs belonging to types III, IV, VII, and XII protein families, were substantially downregulated across several or all time points. Although RLKs are often considered PRRs, it is important to note that many have diverse roles in a myriad of non-pathogen-related functions (Breiden and Simon, 2016). Thus, it is likely that a vast number of RLKs in these families might either be linked to growth and development or are otherwise disadvantageous during P. cinnamomi challenge. There were also RLKs that were mostly upregulated during the earlier time points, such as types I, VIII, X, XII, and XIII LRR-RLK gene families, while contrastingly, genes encoding for L-lectin, TLPK, RFK3-like, PERK-like, LRK10-like, and CR4-like RLK proteins were, by and large, upregulated at later time-points. Of note, RLKs have diversly described roles against a multitude of biotrophic, hemibiotrophic, and necrotrophic pathogens (Llorente et al., 2005; Veronese et al., 2005; Roux et al., 2011; Mosher et al., 2013; Zhang et al., 2013a,b). Thus, it is interesting to broadly posit roles for many of these RLK gene families in either early (biotrophic) or late (necrotrophic) responses to P. cinnamomi based on their temporal affinity. It might also be worthwhile to determine the roles of RLK coding genes, such as S-domain, DUF26, WAK-like, and types II, V, and VII LRR-RLKs, of which a large percentage were upregulated at all time-points following P. cinnamomi inoculation.
The advent of PTI is also often associated with increased callose deposition at the site of infection (Voigt, 2014). Interestingly, in the partially resistant avocado rootstock R0.06, challenged by P. cinnamomi, confocal microscopy noted rapid deposition of callose as early as 6 hpi near the site of infection (van den Berg et al., 2018a). In contrast, some scant callose deposition was only noted at 96 hpi in the susceptible rootstock R0.12. A similar observation has been made for Dusa®, in which P. cinnamomi inoculation leads to rapid accumulation of callose near the site of infection (unpublished study, personal communication). Surprisingly, DE callose synthase genes were not identified from transcriptomic data (van den Berg et al., 2018b). Nevertheless, it has been suggested that the regulation of callose biosynthesis might not occur primarily at the transcriptional level but rather at the translational level or by the involvement of proteases (Nakashima et al., 2003; Saheed et al., 2009). In support of this hypothesis, Arabidopsis thaliana plants overexpressing Vitis quinquangularis aspartic protease 13 (VqAP13) showed enhanced callose deposition at the site of Pseudomonas syringae pv. tomato (Pst) DC3000 inoculation when compared to wild-type plants (Guo et al., 2016). Similarly, overexpression of the Oryza sativa matrix metalloprotease OsMMP1 in tobacco led to increased callose and cellulose deposition (Das et al., 2018). Thus, it is entirely possible that enhanced callose deposition in some partially resistant avocado rootstocks may, to some extent, be regulated by proteases. In the present study, several proteases were significantly upregulated at 6 hpi in Dusa® but not in R0.12, including the aspartic protease Peame105C05g001950, the serine protease Peame105C04g026050, and the subtilase Peame105C01g045590. Given the early onset (6 hpi) of enhanced callose deposition in partially resistant rootstocks, these proteases may deserve further study in this regard along with the non-apoplastic proteases that were not looked at in this study. Nonetheless, elucidating the role of these and other avocado proteases in callose deposition as well as other aspects of PTI should constitute the topic of future investigations.
Otherwise, several subtilases in either rootstock were significantly upregulated at 6 hpi when compared to the mock-inoculated control plants, and while the susceptible rootstock R0.12 increased the expression of three subtilases at 6 hpi, four were upregulated in Dusa®. Several subtilases were also significantly upregulated at 12 and 24 hpi in Dusa® but not in R0.12. These increases occurred during the theoretical biotrophic phase of P. cinnamomi's life cycle; as such it is tempting to posit that they might play a role in inducing ETI and subsequent HR. HR, characterized by rapid localized cell death at the initial site of pathogen infiltration, constitutes an effective strategy for defense against biotrophic pathogens and early interactions with hemibiotrophic pathogens (Balint-Kurti, 2019). In animals, caspases (cysteine-containing aspartate-specific proteases) play an important part in the control of PCD (Boatright and Salvesen, 2003). Interestingly, several proteases, specifically apoplastic subtilases, have been implicated in caspase-like activity in plants, promoting cell death during pathogenic challenge (Coffeen and Wolpert, 2004; Chichkova et al., 2010; Vartapetian et al., 2011; Fernández et al., 2012; Zimmermann et al., 2016). In tomato, the subtilase P69B was induced by both biotrophic and hemibiotrophic pathogens, as well as salicylic acid (SA) application (Fischer et al., 1989; Tornero et al., 1996; Jordá et al., 1999; Jord and Vera, 2000). In an elegant study, Zimmermann et al. (2016) showed that P69B promoted cell death in tomato, although the overexpression of P69B alone could not. Thus, control over PCD is an undoubtedly important aspect of defense responses, and apoplastic subtilases form an integral part of its control, although additional proteases also participate (Godson and van der Hoorn, 2021). Investigations aimed at determining whether these subtilases display caspase-like activity would be highly informative.
It was worth noting that two tomato MMPs, Sl2-MMP and Sl3-MMP, act antagonistically upstream of P69B, cleaving it to prevent cell death (Zimmermann et al., 2016). Interestingly, in the current study, all the identified apoplastic metalloproteases were downregulated in Dusa® at 6 hpi, while two in R0.12, Peame105C10g014010 and Peame105C10g013220, were not. The same metalloproteases were then significantly upregulated in both rootstocks at 120 hpi. This pattern of expression would come as no surprise if one were to first presume, as was the case for Sl2-MMP and Sl3-MMP, that these metalloproteases suppress cell death. It, therefore likely that further investigations might reveal similar functions for the avocado metalloproteases described in this study, and that the initiation of HR in Dusa® is more strongly induced as a result.
Some extracellular proteases have been implicated in the regulation of systemic acquired resistance (SAR) and priming of host immune responses (Xia et al., 2004; Prasad et al., 2009; Ramírez et al., 2013; Breitenbach et al., 2014). Pertinently, extracellular aspartic proteases in both tobacco and tomato have been associated with degradation of PR proteins (Rodrigo et al., 1989, 1991). Thus, it was previously suggested that an extracellular aspartic protease may be responsible for the cleavage of PR-1b (Hou et al., 2018). The resulting C-terminal PR-1b fragment, CAP-Derived Peptide 1 (CAPE1), has been defined as a DAMP and linked to the activation of jasmonic acid (JA) and SA defense response pathways, as well as SAR activation (Chen et al., 2014). Additionally, the aspartic protease Constitutive Disease Resistance 1 (CDR1) has been implicated in inducing SAR both locally and systemically (Xia et al., 2004; Prasad et al., 2009), while others such as apoplastic, EDS1-dependent 1 (AED1) serve to limit SAR signaling (Breitenbach et al., 2014). In A. thaliana, the subtilase SBT3.3 was found to participate in chromatin remodeling of SA defense-related genes, activating them following pathogen challenge and priming future immune response (Ramírez et al., 2013). Meanwhile, the V. quinquangularis aspartic protease coding gene VqAP13 has been positively associated with the SA defense response pathway (Guo et al., 2013, 2016). Thus, the differences seen between R0.12 and Dusa® regarding differential expression of aspartic proteases and subtilases cannot be overlooked. For instance, it is conceivable that inappropriate expression of an AED1 homolog in R0.12 could lead to suppression of defense responses in favor of growth. In contrast, upregulation of avocado AP13, CDR1, or SBT3.3 homolog in Dusa® might lead to increased SAR or priming of SA defense response genes systemically. Congruently, it was previously shown that the SA defense response is an important part of the early immune response in Dusa® challenged by P. cinnamomi (van den Berg et al., 2018b). Interestingly, the predicted protein coding sequences for Peame105C01g045560 and Peame105C01g045590 closely resemble SBT3.3 in C. micranthum. Similarly, the product of Peame105C08g006080 showed high similarity to several SBT3-like proteins in C. micranthum. Based on the near opposite regulatory patterns seen for two of these genes, Peame105C01g045590 and Peame105C08g006080, in Dusa® and R0.12, further investigation is warranted. It would be interesting to determine whether the overexpression of these genes has the potential to enhance disease resistance and prime immune responses in a model system such as A. thaliana challenged by the oomycete Hyaloperonospora arabidopsidis or bacterium Pst DC3000.
Besides being one of the most well-studied plant proteases, SBTs are most often secreted (Schaller et al., 2018). Quite notably, the predicted proteins of all DE avocado apoplastic SBTs were expected to localize to the extracellular space. Additionally, all but one DE avocado apoplastic serine protease coding gene, excluding SBTs, was predicted to code for serine carboxypeptidase-like proteins (SCLPs). Markedly, SCLPs are known to contain signal peptides and are quite often found in the apoplast (Fraser et al., 2005; Sueldo et al., 2014; Grosse-Holz et al., 2018). Similarly, MMPs from the M10 family in plants have all been described as containing signal peptides and are the only metalloproteases in plants that are classed as apoplastic (Godson and van der Hoorn, 2021). Therefore, it is interesting to note that all but one DE avocado apoplastic metalloprotease coding gene was predicted to belong to the M10 MMP family. Together these observations provide support for the accuracy of the web service used to predict the subcellular localization of the proteases in this study (Stärk et al., 2021) and the relevance of the identified DE avocado apoplastic proteases for future characterization studies.
Conclusion
In the current study, the basis of P. cinnamomi partial resistance in avocado was explored by dual RNA sequencing of a susceptible (R0.12) and a partially resistant rootstock (Dusa®). The initial analyses revealed significant differences in the number of DEGs between the rootstocks, particularly at 12 and 24 hpi. The GO enrichment indicated that while Dusa® activated various defense responses during these times, R0.12 activated pathways related to growth and development. P. cinnamomi DEGs reflected this difference, indicating that the use of effector proteins may either be causing or responding to the lack of defense responses in R0.12 at 12 and 24 hpi. In either case, the lack of appropriate defense responses from R0.12 at these intermediate time points coincided with the expected window in which P. cinnamomi switches between its biotrophic and necrotrophic life stages. In an attempt to uncover the root cause of such strikingly different responses, we also described broadly the expression of several avocado genes encoding RLKs and apoplastic proteases. To the best of our knowledge, no literature currently describes two plants of the same species, differing in their susceptibility to a pathogen, with such overwhelmingly different responses. As such, understanding the root cause of such a sweeping change in gene expression could provide novel insights into the avocado-P. cinnamomi interaction and the basis of partial resistance, and should preface investigations going forward. It should be noted that this study only used a single mock-inoculated avocado control time point because of resource restraints. Additionally, the use of a single susceptible and partially resistant rootstock limited the potential to highlight additional defense response pathways. Future investigations should endeavor to include controls at each time point to strengthen statistical analyses and use several rootstocks varying in their susceptibility to P. cinnamomi to provide a more holistic view of avocado defense responses following inoculation. Finally, it should be noted that low read counts were obtained for P. cinnamomi at the earliest time points. This limitation and potential cause for high inter-sample variation should be noted when considering the conclusions made in this study on P. cinnamomi.
Data Availability Statement
The raw data for all RNA-sequencing sample libraries used in this study are available from the Sequence Read Archive of NCBI Genbank under accession number PRJNA675400.
Author Contributions
NB, JE, and RB: conceptualized and reviewed the manuscript. JE: performed all the wet work. RB: performed all the data analyses and drafted the manuscript. All authors contributed to and approved the final version of the manuscript.
Funding
This study was supported by the University of Pretoria and the Forestry and Agricultural Biotechnology Institute (FABI). Funding was provided by the Hans Merensky Foundation.
Conflict of Interest
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.
Publisher's Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We are grateful to the Avocado Genome Consortium for providing early access to their genome prior to publication.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.928176/full#supplementary-material
References
Andrews, S. (2010). FastQC: a quality control tool for high throughput sequence data [Online]. Babraham Bioinformatics: Babraham Institute. Available online at: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed February 11, 2021).
Balint-Kurti, P. (2019). The plant hypersensitive response: concepts, control and consequences. Mol. Plant Pathol. 20, 1163–1178. doi: 10.1111/mpp.12821
Boatright, K. M., and Salvesen, G. S. (2003). Mechanisms of caspase activation. Curr. Opinion Cell Biol. 15, 725–731. doi: 10.1016/j.ceb.2003.10.009
Böhm, H., Albert, I., Fan, L., Reinhard, A., and Nürnberger, T. (2014). Immune receptor complexes at the plant cell surface. Curr. Opinion Cell Biol. 20, 47–54. doi: 10.1016/j.pbi.2014.04.007
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Breiden, M., and Simon, R. (2016). QandA: how does peptide signaling direct plant development? BMC Biol. 14, 58. doi: 10.1186/s12915-016-0280-3
Breitenbach, H. H., Wenig, M., Wittek, F., Jordá, L., Maldonado-Alconada, A. M., Sarioglu, H., et al. (2014). Contrasting roles of the apoplastic aspartyl protease apoplastic, enhanced disease susceptibility1-dependent1 and legume lectin-like protein1 in arabidopsis systemic acquired resistance. Plant Physiol. 165, 791–809. doi: 10.1104/pp.114.239665
Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., et al. (2009). BLAST+: architecture and applications. BMC Bioinform. 10, 421. doi: 10.1186/1471-2105-10-421
Chanderbali, A. S., van der Werff, H., and Renner, S. S. (2001). Phylogeny and historical biogeography of Lauraceae: evidence from the chloroplast and nuclear genomes. Annal. Missouri Botanic. Garden 1, 104–134. doi: 10.2307/2666133
Chen, X. R., Xing, Y. P., Li, Y. P., Tong, Y. H., and Xu, J. Y. (2013). RNA-Seq reveals infection-related gene expression changes in Phytophthora capsici. PLoS ONE 8, e74588. doi: 10.1371/journal.pone.0074588
Chen, Y. L., Lee, C. Y., Cheng, K. T., Chang, W. H., Huang, R. N., Nam, H. G., et al. (2014). Quantitative peptidomics study reveals that a wound-induced peptide from PR-1 regulates immune signaling in tomato. Plant Cell 26, 4135–4148. doi: 10.1105/tpc.114.131185
Chichkova, N. V., Shaw, J., Galiullina, R. A., Drury, G. E., Tuzhikov, A. I., Kim, S. H., et al. (2010). Phytaspase, a relocalisable cell death promoting plant protease with caspase specificity. EMBO J. 29, 1149–1161. doi: 10.1038/emboj.2010.1
Coffeen, W. C., and Wolpert, T. J. (2004). Purification and characterization of serine proteases that exhibit caspase-like activity and are associated with programmed cell death in Avena sativa. Plant Cell 16, 857–873. doi: 10.1105/tpc.017947
Couto, D., and Zipfel, C. (2016). Regulation of pattern recognition receptor signalling in plants. Nat. Rev. Immunol. 16, 537–552. doi: 10.1038/nri.2016.77
Cui, H., Tsuda, K., and Parker, J. E. (2015). Effector-triggered immunity: from pathogen perception to robust defense. Ann. Rev. Plant Biol. 66, 487–511. doi: 10.1146/annurev-arplant-050213-040012
Dalio, R., Maximo, H., Oliveira, T., Dias, R., Breton, M., Felizatti, H., et al. (2018). Phytophthora parasitica effector PpRxLR2 suppresses Nicotiana benthamiana immunity. Mol. Plant-Microbe Interact. 31, 481–493. doi: 10.1094/MPMI-07-17-0158-FI
Damasceno, C. M., Bishop, J. G., Ripoll, D. R., Win, J., Kamoun, S., and Rose, J. K. (2008). Structure of the glucanase inhibitor protein (GIP) family from Phytophthora species suggests coevolution with plant endo-beta-1,3-glucanases. Mol. Plant-Microbe Interact. 21, 820–830. doi: 10.1094/MPMI-21-6-0820
Dann, E., Whiley, T., Pegg, K., Dean, J., Chandra, K., and Smith, L. (2013). Evaluation of fruit yields and tree health of Hass trees grafted to 8 rootstocks in high Phytophthora root rot conditions at Childers, QLD. Talking Avo. 24, 34–37. Available online at: https://avocado.org.au/wp-content/uploads/2016/10/Spring_v24_n3_2013.pdf
Das, P. K., Biswas, R., Anjum, N., Das, A. K., and Maiti, M. K. (2018). Rice matrix metalloproteinase OsMMP1 plays pleiotropic roles in plant development and symplastic-apoplastic transport by modulating cellulose and callose depositions. Scientific Rep. 8, 2783. doi: 10.1038/s41598-018-20070-4
Davis, K. R., and Hahlbrock, K. (1987). Induction of defense responses in cultured parsley cells by plant cell wall fragments. Plant Physiol. 84, 1286–1290. doi: 10.1104/pp.84.4.1286
Derevnina, L., Dagdas, Y. F., De la Concepcion, J. C., Bialas, A., Kellner, R., Petre, B., et al. (2016). Nine things to know about elicitins. New Phytol. 212, 888–895. doi: 10.1111/nph.14137
Dobrowolski, M., Shearer, B., Colquhoun, I., O'brien, P., and Hardy, G. S. (2008). Selection for decreased sensitivity to phosphite in Phytophthora cinnamomi with prolonged use of fungicide. Plant Pathol. 57, 928–936. doi: 10.1111/j.1365-3059.2008.01883.x
Dou, D., and Zhou, J. M. (2012). Phytopathogen effectors subverting host immunity: different foes, similar battleground. Cell Host Microbe 12, 484–495. doi: 10.1016/j.chom.2012.09.003
Ekchaweng, K., Evangelisti, E., Schornack, S., Tian, M., and Churngchow, N. (2017). The plant defense and pathogen counterdefense mediated by Hevea brasiliensis serine protease HbSPA and Phytophthora palmivora extracellular protease inhibitor PpEPI10. PLoS ONE 12, e0175795. doi: 10.1371/journal.pone.0175795
Engelbrecht, J., Duong, T. A., Prabhu, S. A., Seedat, M., and van den Berg, N. (2021). Genome of the destructive oomycete Phytophthora cinnamomi provides insights into its pathogenicity and adaptive potential. BMC Genom. 22, 1–15. doi: 10.1186/s12864-021-07552-y
Engelbrecht, J., Duong, T. A., and van den Berg, N. (2013). Development of a nested quantitative real-time PCR for detecting Phytophthora cinnamomi in Persea americana rootstocks. Plant Dis. 97, 1012–1017. doi: 10.1094/PDIS-11-12-1007-RE
Engelbrecht, J., and van den Berg, N. (2013). Expression of defence-related genes against Phytophthora cinnamomi in five avocado rootstocks. South Afric. J. Sci. 109, 1–8. doi: 10.1590/sajs.2013/20120058
Erwin, D. C., Bartnicki-Garcia, S., and Tsao, P. H. (1983). Phytophthora : Its Biology, Taxonomy, Ecology, and Pathology. St. Paul, Minn.: American Phytopathological Society.
Ewels, P., Magnusson, M., Lundin, S., and Käller, M. (2016). MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32, 3047–3048. doi: 10.1093/bioinformatics/btw354
Fernández, M. B., Daleo, G. R., and Guevara, M. G. (2012). DEVDase activity is induced in potato leaves during Phytophthora infestans infection. Plant Physiology and Biochemistry 61, 197–203. doi: 10.1016/j.plaphy.2012.10.007
Fick, A., Swart, V., Backer, R., Bombarely, A., Engelbrecht, J., and van den Berg, N. (2022). Partially resistant avocado rootstock dusa® shows prolonged upregulation of nucleotide binding-leucine rich repeat genes in response to Phytophthora cinnamomi infection. Front. Plant Sci. 13, 644. doi: 10.3389/fpls.2022.793644
Fischer, W., Christ, U., Baumgartner, M., Erismann, K. H., and Mösinger, E. (1989). Pathogenesis-related proteins of tomato: II. biochemical and immunological characterization. Physiologic. Mol. Plant Pathol. 35, 67–83. doi: 10.1016/0885-5765(89)90008-8
Fraser, C. M., Rider, L. W., and Chapple, C. (2005). An Expression and bioinformatics analysis of the arabidopsis serine carboxypeptidase-like gene family. Plant Physiol. 138, 1136–1148. doi: 10.1104/pp.104.057950
Futschik, M. E., and Carlisle, B. (2005). Noise-robust soft clustering of gene expression time-course data. J. Bioinform. Computat. Biol. 3, 965–988. doi: 10.1142/S0219720005001375
Galili, T. (2015). dendextend: an R package for visualizing, adjusting and comparing trees of hierarchical clustering. Bioinformatics 31, 3718–3720. doi: 10.1093/bioinformatics/btv428
Gaulin, E., Drame, N., Lafitte, C., Torto-Alalibo, T., Martinez, Y., Ameline-Torregrosa, C., et al. (2006). Cellulose binding domains of a Phytophthora cell wall protein are novel pathogen-associated molecular patterns. Plant Cell 18, 1766–1777. doi: 10.1105/tpc.105.038687
Gaulin, E., Jauneau, A., Villalba, F., Rickauer, M., Esquerre-Tugaye, M. T., and Bottin, A. (2002). The CBEL glycoprotein of Phytophthora parasitica var. nicotianae is involved in cell wall deposition and adhesion to cellulosic substrates. J. Cell Sci. 115, 4565–4575. doi: 10.1242/jcs.00138
Gentry, A. H. (1988). Changes in plant community diversity and floristic composition on environmental and geographical gradients. Annal. Missouri Botanic. Garden 11, 1–34. doi: 10.2307/2399464
Glazebrook, J. (2005). Contrasting mechanisms of defense against biotrophic and necrotrophic pathogens. Ann. Rev. Phytopathol. 43, 205–227. doi: 10.1146/annurev.phyto.43.040204.135923
Godson, A., and van der Hoorn, R. A. L. (2021). The front line of defence: a meta-analysis of apoplastic proteases in plant immunity. J. Experiment. Bot. 72, 3381–3394. doi: 10.1093/jxb/eraa602
Grosse-Holz, F., Kelly, S., Blaskowski, S., Kaschani, F., Kaiser, M., and van der Hoorn, R. A. L. (2018). The transcriptome, extracellular proteome and active secretome of agroinfiltrated Nicotiana benthamiana uncover a large, diverse protease repertoire. Plant Biotechnol. J. 16, 1068–1084. doi: 10.1111/pbi.12852
Guo, R., Tu, M., Wang, X., Zhao, J., Wan, R., Li, Z., et al. (2016). Ectopic expression of a grape aspartic protease gene, AP13, in Arabidopsis thaliana improves resistance to powdery mildew but increases susceptibility to Botrytis cinerea. Plant Sci. 248, 17–27. doi: 10.1016/j.plantsci.2016.04.006
Guo, R., Xu, X., Carole, B., Li, X., Gao, M., Zheng, Y., et al. (2013). Genome-wide identification, evolutionary and expression analysis of the aspartic protease gene superfamily in grape. BMC Genom. 14, 1–18. doi: 10.1186/1471-2164-14-554
Haas, B. J., Kamoun, S., Zody, M. C., Jiang, R. H., Handsaker, R. E., Cano, L. M., et al. (2009). Genome sequence and analysis of the Irish potato famine pathogen Phytophthora infestans. Nature 461, 393–398. doi: 10.1038/nature08358
Hardham, A. R. (2005). Phytophthora cinnamomi. Mol. Plant Pathol. 6, 589–604. doi: 10.1111/j.1364-3703.2005.00308.x
Hardham, A. R., and Blackman, L. M. (2018). Phytophthora cinnamomi. Mol. Plant Pathol. 19, 260–285. doi: 10.1111/mpp.12568
Hardy, G. S. J., Barrett, S., and Shearer, B. (2001). The future of phosphite as a fungicide to control the soilborne plant pathogen Phytophthora cinnamomi in natural ecosystems. Austral. Plant Pathol. 30, 133–139. doi: 10.1071/AP01012
Hou, S., Jamieson, P., and He, P. (2018). The cloak, dagger, and shield: proteases in plant-pathogen interactions. Biochem. J. 475, 2491–2509. doi: 10.1042/BCJ20170781
Hurtado-Fernández, E., Fernández-Gutiérrez, A., and Carrasco-Pancorbo, A. (2018). “Avocado fruit—Persea americana,” in Exotic Fruits. (New York, NY: Elsevier), p. 37–48.
Jiang, R. H., Tripathy, S., Govers, F., and Tyler, B. M. (2008). RXLR effector reservoir in two Phytophthora species is dominated by a single rapidly evolving superfamily with more than 700 members. Proceed. Nat. Acad. Sci. 105, 4874–4879. doi: 10.1073/pnas.0709303105
Jones, J. D., and Dangl, J. L. (2006). The plant immune system. Nature 444, 323–329. doi: 10.1038/nature05286
Jordá, L., Coego, A., Conejero, V., and Vera, P. (1999). A genomic cluster containing four differentially regulated subtilisin-like processing protease genes is in tomato plants. J. Biologic. Chemistr. 274, 2360–2365. doi: 10.1074/jbc.274.4.2360
Jordá, L., and Vera, P. (2000). Local and systemic induction of two defense-related subtilisin-like protease promoters in transgenic Arabidopsis plants. Luciferin induction of PR gene expression. Plant Physiol. 124, 1049–1058. doi: 10.1104/pp.124.3.1049
Joubert, M., Backer, R., Engelbrecht, J., and van den Berg, N. (2021). Expression of several Phytophthora cinnamomi putative RxLRs provides evidence for virulence roles in avocado. PLoS One 16, e0254645. doi: 10.1371/journal.pone.0254645
Kamoun, S. (2007). Groovy times: filamentous pathogen effectors revealed. Curr. Opinion Cell Biol. 10, 358–365. doi: 10.1016/j.pbi.2007.04.017
Kelley, B. S., Lee, S. J., Damasceno, C. M., Chakravarthy, S., Kim, B. D., Martin, G. B., et al. (2010). A secreted effector protein (SNE1) from Phytophthora infestans is a broadly acting suppressor of programmed cell death. Plant J. 62, 357–366. doi: 10.1111/j.1365-313X.2010.04160.x
Kim, D., Paggi, J. M., Park, C., Bennett, C., and Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907–915. doi: 10.1038/s41587-019-0201-4
Larroque, M., Ramirez, D., Lafitte, C., Borderies, G., Dumas, B., and Gaulin, E. (2011). Expression and purification of a biologically active Phytophthora parasitica cellulose binding elicitor lectin in Pichia pastoris. Prot. Express. Purificat. 80, 217–223. doi: 10.1016/j.pep.2011.07.007
Li, H. (1995). The origin and evolution of Litsea genera group (Laureae) in Lauraceae. Acta Botanica Yunnanica 17, 251–254.
Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352
Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930. doi: 10.1093/bioinformatics/btt656
Llorente, F., Alonso-Blanco, C., Sánchez-Rodriguez, C., Jorda, L., and Molina, A. (2005). ERECTA receptor-like kinase and heterotrimeric G protein from Arabidopsis are required for resistance to the necrotrophic fungus Plectosphaerella cucumerina. Plant J. 43, 165–180. doi: 10.1111/j.1365-313X.2005.02440.x
Lohse, M., Nagel, A., Herter, T., May, P., Schroda, M., Zrenner, R., et al. (2014). Mercator: a fast and simple web server for genome scale functional annotation of plant sequence data. Plant Cell Environ. 37, 1250–1258. doi: 10.1111/pce.12231
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 1–21. doi: 10.1186/s13059-014-0550-8
Ma, Z., Song, T., Zhu, L., Ye, W., Wang, Y., Shao, Y., et al. (2015). A Phytophthora sojae glycoside hydrolase 12 protein is a major virulence factor during soybean infection and is recognized as a PAMP. Plant Cell 27, 2057–2072. doi: 10.1105/tpc.15.00390
Ma, Z., Zhu, L., Song, T., Wang, Y., Zhang, Q., Xia, Y., et al. (2017). A paralogous decoy protects Phytophthora sojae apoplastic effector PsXEG1 from a host inhibitor. Science 355, 710–714. doi: 10.1126/science.aai7919
Mafurah, J. J., Ma, H., Zhang, M., Xu, J., He, F., Ye, T., et al. (2015). A virulence essential CRN effector of Phytophthora capsici suppresses host defense and induces cell death in plant nucleus. PLoS One 10, e0127965. doi: 10.1371/journal.pone.0127965
Mahomed, W., and van den Berg, N. (2011). EST sequencing and gene expression profiling of defence-related genes from Persea americana infected with Phytophthora cinnamomi. BMC Plant Biol. 11, 167. doi: 10.1186/1471-2229-11-167
Mapleson, D., Venturini, L., Kaithakottil, G., and Swarbreck, D. (2018). Efficient and accurate detection of splice junctions from RNA-seq with Portcullis. GigaScience 7, giy131. doi: 10.1093/gigascience/giy131
Martin, L., Fei, Z., Giovannoni, J., and Rose, J. K. C. (2013). Catalyzing plant science research with RNA-seq. Front. Plant Sci. 4, 66. doi: 10.3389/fpls.2013.00066
Mateos, F. V., Rickauer, M., and Esquerre-Tugaye, M. T. (1997). Cloning and characterization of a cDNA encoding an elicitor of Phytophthora parasitica var. nicotianae that shows cellulose-binding and lectin-like activities. Mol. Plant-Microbe Interact. 10, 1045–1053. doi: 10.1094/MPMI.1997.10.9.1045
McDowell, J. M., and Woffenden, B. J. (2003). Plant disease resistance genes: recent insights and potential applications. Trends Biotechnol. 21, 178–183. doi: 10.1016/S0167-7799(03)00053-2
McGowan, J., and Fitzpatrick, D. A. (2017). Genomic, network, and phylogenetic analysis of the oomycete effector arsenal. MSphere 2, e00408–00417. doi: 10.1128/mSphere.00408-17
McGowan, J., and Fitzpatrick, D. A. (2020). Recent advances in oomycete genomics. Adv. Genetics 105, 175–228. doi: 10.1016/bs.adgen.2020.03.001
Monteiro, F., and Nishimura, M. T. (2018). Structural, functional, and genomic diversity of plant nlr proteins: an evolved resource for rational engineering of plant immunity. Ann. Rev. Phytopathol. 56, 243–267. doi: 10.1146/annurev-phyto-080417-045817
Mosher, S., Seybold, H., Rodriguez, P., Stahl, M., Davies, K. A., Dayaratne, S., et al. (2013). The tyrosine-sulfated peptide receptors PSKR1 and PSY1R modify the immunity of Arabidopsis to biotrophic and necrotrophic pathogens in an antagonistic manner. Plant J. 73, 469–482. doi: 10.1111/tpj.12050
Nakashima, J., Laosinchai, W., Cui, X., and Malcolm Brown, R. (2003). New insight into the mechanism of cellulose and callose biosynthesis: proteases may regulate callose biosynthesis upon wounding. Cellulose 10, 369–389. doi: 10.1023/A:1027336605479
Naveed, Z. A., Bibi, S., and Ali, G. S. (2019). The Phytophthora RXLR Effector Avrblb2 modulates plant immunity by interfering with ca2+ signaling pathway. Front. Plant Sci. 10, 374. doi: 10.3389/fpls.2019.00374
Orsomando, G., Lorenzi, M., Raffaelli, N., Dalla Rizza, M., Mezzetti, B., and Ruggieri, S. (2001). Phytotoxic protein PcF, purification, characterization, and cDNA sequencing of a novel hydroxyproline-containing factor secreted by the strawberry pathogen Phytophthora cactorum. J. Biologic. Chemistr. 276, 21578–21584. doi: 10.1074/jbc.M101377200
Park, C. J., Caddell, D. F., and Ronald, P. C. (2012). Protein phosphorylation in plant immunity: insights into the regulation of pattern recognition receptor-mediated signaling. Front. Plant Sci. 3, 177–177. doi: 10.3389/fpls.2012.00177
Pegg, K., Whiley, A., Saranah, J., and Glass, R. (1985). Control of Phytophthora root rot of avocado with phosphorus acid. Austral. Plant Pathol. 14, 25–29. doi: 10.1071/APP9850025
Pemberton, C. L., and Salmond, G. P. C. (2004). The Nep1-like proteins-a growing family of microbial elicitors of plant necrosis. Mol. Plant Pathol. 5, 353–359. doi: 10.1111/j.1364-3703.2004.00235.x
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122
Prasad, B. D., Creissen, G., Lamb, C., and Chattoo, B. B. (2009). Overexpression of rice (Oryza sativa L.) OsCDR1 leads to constitutive activation of defense responses in rice and Arabidopsis. Mol. Plant Microbe Interact. 22, 1635–1644. doi: 10.1094/MPMI-22-12-1635
R Core Team (2021). R: A language and environment for statistical computing [Online]. Vienna, Austria: R Foundation for Statistical Computing. Available: https://www.r-project.org/ (accessed February 17, 2021).
Ramírez, V., López, A., Mauch-Mani, B., Gil, M. J., and Vera, P. (2013). An extracellular subtilase switch for immune priming in arabidopsis. PLoS Pathogens 9, e1003445. doi: 10.1371/journal.ppat.1003445
Reeksting, B. J., Coetzer, N., Mahomed, W., Engelbrecht, J., and van den Berg, N. (2014). De novo sequencing, assembly, and analysis of the root transcriptome of Persea americana (Mill.) in response to Phytophthora cinnamomi and flooding. PLoS One 9, e86399. doi: 10.1371/journal.pone.0086399
Reeksting, B. J., Olivier, N. A., and van den Berg, N. (2016). Transcriptome responses of an ungrafted Phytophthora root rot tolerant avocado (Persea americana) rootstock to flooding and Phytophthora cinnamomi. BMC Plant Biol. 16, 205. doi: 10.1186/s12870-016-0893-2
Reitmann, A., Berger, D. K., and van den Berg, N. (2017). Putative pathogenicity genes of Phytophthora cinnamomi identified via RNA-Seq analysis of pre-infection structures. Euro. J. Plant Pathol. 147, 211–228. doi: 10.1007/s10658-016-0993-8
Rodrigo, I., Vera, P., and Conejero, V. (1989). Degradation of tomato pathogenesis-related proteins by an endogenous 37-kDa aspartyl endoproteinase. Euro. J. Biochemistr. 184, 663–669. doi: 10.1111/j.1432-1033.1989.tb15064.x
Rodrigo, I., Vera, P., Van Loon, L. C., and Conejero, V. (1991). Degradation of tobacco pathogenesis-related proteins: evidence for conserved mechanisms of degradation of pathogenesis-related proteins in plants. Plant Physiol. 95, 616–622. doi: 10.1104/pp.95.2.616
Rooney, H. C., Van't Klooster, J. W., van der Hoorn, R. A., Joosten, M. H., Jones, J. D., and de Wit, P. J. (2005). Cladosporium Avr2 inhibits tomato Rcr3 protease required for Cf-2-dependent disease resistance. Science 308, 1783–1786. doi: 10.1126/science.1111404
Roux, M., Schwessinger, B., Albrecht, C., Chinchilla, D., Jones, A., Holton, N., et al. (2011). The Arabidopsis Leucine-Rich Repeat Receptor–Like Kinases BAK1/SERK3 and BKK1/SERK4 are required for innate immunity to hemibiotrophic and biotrophic pathogens. Plant Cell 23, 2440–2455. doi: 10.1105/tpc.111.084301
Saheed, S. A., Cierlik, I., Larsson, K. A. E., Delp, G., Bradley, G., Jonsson, L. M. V., et al. (2009). Stronger induction of callose deposition in barley by Russian wheat aphid than bird cherry-oat aphid is not associated with differences in callose synthase or β-1,3-glucanase transcript abundance. Physiologia Plantar. 135, 150–161. doi: 10.1111/j.1399-3054.2008.01180.x
Schaller, A., Stintzi, A., Rivas, S., Serrano, I., Chichkova, N. V., Vartapetian, A. B., et al. (2018). From structure to function - a family portrait of plant subtilases. New Phytol. 218, 901–915. doi: 10.1111/nph.14582
Schornack, S., van Damme, M., Bozkurt, T. O., Cano, L. M., Smoker, M., Thines, M., et al. (2010). Ancient class of translocated oomycete effectors targets the host nucleus. Proceed. Nat. Acad. Sci. 107, 17421–17426. doi: 10.1073/pnas.1008491107
Schwacke, R., Ponce-Soto, G. Y., Krause, K., Bolger, A. M., Arsova, B., Hallab, A., et al. (2019). MapMan4: a refined protein classification and annotation framework applicable to multi-omics data analysis. Mol. Plant 12, 879–892. doi: 10.1016/j.molp.2019.01.003
Seidl, M. F., and van den Ackerveken, G. (2019). Activity and phylogenetics of the broadly occurring family of microbial nep1-like proteins. Ann. Rev. Phytopathol. 57, 367–386. doi: 10.1146/annurev-phyto-082718-100054
Shabab, M., Shindo, T., Gu, C., Kaschani, F., Pansuriya, T., Chintha, R., et al. (2008). Fungal effector protein AVR2 targets diversifying defense-related cys proteases of tomato. Plant Cell 20, 1169–1183. doi: 10.1105/tpc.107.056325
Song, J., Win, J., Tian, M., Schornack, S., Kaschani, F., Ilyas, M., et al. (2009). Apoplastic effectors secreted by two unrelated eukaryotic plant pathogens target the tomato defense protease Rcr3. Proceed. Nat. Acad. Sci. 106, 1654–1659. doi: 10.1073/pnas.0809201106
Song, Y., Yao, X., Tan, Y., Gan, Y., and Corlett, R. T. (2016). Complete chloroplast genome sequence of the avocado: gene organization, comparative analysis, and phylogenetic relationships with other Lauraceae. Canad. J. For. Res. 46, 1293–1301. doi: 10.1139/cjfr-2016-0199
Stam, R., Jupe, J., Howden, A. J., Morris, J. A., Boevink, P. C., Hedley, P. E., et al. (2013). Identification and Characterisation CRN Effectors in Phytophthora capsici shows modularity and functional diversity. PLoS One 8, e59517. doi: 10.1371/annotation/90bd45cb-33a7-426f-a928-9ddc351b08cc
Stärk, H., Dallago, C., Heinzinger, M., and Rost, B. (2021). Light attention predicts protein location from the language of life. Bioinform. Adv. 1, 1. doi: 10.1093/bioadv/vbab035
Sueldo, D., Ahmed, A., Misas-Villamil, J., Colby, T., Tameling, W., Joosten, M. H., et al. (2014). Dynamic hydrolase activities precede hypersensitive tissue collapse in tomato seedlings. New Phytol. 203, 913–925. doi: 10.1111/nph.12870
Supek, F., Bošnjak, M., Škunca, N., and Šmuc, T. (2011). REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS One 6, e21800. doi: 10.1371/journal.pone.0021800
Tian, M., Benedetti, B., and Kamoun, S. (2005). A Second Kazal-like protease inhibitor from Phytophthora infestans inhibits and interacts with the apoplastic pathogenesis-related protease P69B of tomato. Plant Physiol. 138, 1785–1793. doi: 10.1104/pp.105.061226
Tian, M., Huitema, E., Da Cunha, L., Torto-Alalibo, T., and Kamoun, S. (2004). A Kazal-like extracellular serine protease inhibitor from Phytophthora infestans targets the tomato pathogenesis-related protease P69B. J. Biologic. Chemistr. 279, 26370–26377. doi: 10.1074/jbc.M400941200
Tian, M., Win, J., Song, J., van der Hoorn, R., van der Knaap, E., and Kamoun, S. (2007). A Phytophthora infestans cystatin-like protein targets a novel tomato papain-like apoplastic protease. Plant Physiol. 143, 364–377. doi: 10.1104/pp.106.090050
Tornero, P., Conejero, V., and Vera, P. (1996). Primary structure and expression of a pathogen-induced protease (PR-P69) in tomato plants: similarity of functional domains to subtilisin-like endoproteases. Proceed. Nat. Acad. Sci. 93, 6332–6337. doi: 10.1073/pnas.93.13.6332
Torto, T. A., Li, S., Styer, A., Huitema, E., Testa, A., Gow, N. A., et al. (2003). EST mining and functional expression assays identify extracellular effector proteins from the plant pathogen Phytophthora. Genome Res. 13, 1675–1685. doi: 10.1101/gr.910003
Tyler, B. M., Tripathy, S., Zhang, X., Dehal, P., Jiang, R. H., Aerts, A., et al. (2006). Phytophthora genome sequences uncover evolutionary origins and mechanisms of pathogenesis. Science 313, 1261–1266. doi: 10.1126/science.1128796
Urban, M., Cuzick, A., Seager, J., Wood, V., Rutherford, K., Venkatesh, S. Y., et al. (2019). PHI-base: the pathogen–host interactions database. Nucleic Acids Res. 48, D613–D620. doi: 10.1093/nar/gkz904
van den Berg, N., Christie, J. B., Aveling, T. A. S., and Engelbrecht, J. (2018a). Callose and β-1,3-glucanase inhibit Phytophthora cinnamomi in a resistant avocado rootstock. Plant Pathol. 67, 1150–1160. doi: 10.1111/ppa.12819
van den Berg, N., Mahomed, W., Olivier, N. A., Swart, V., and Crampton, B. G. (2018b). Transcriptome analysis of an incompatible Persea americana-Phytophthora cinnamomi interaction reveals the involvement of SA-and JA-pathways in a successful defense response. PLoS One 13, e0205705. doi: 10.1371/journal.pone.0205705
van den Berg, N., Swart, V., Backer, R., Fick, A., Wienk, R., Engelbrecht, J., et al. (2021). Advances in Understanding Defense Mechanisms in Persea americana Against Phytophthora cinnamomi. Front. Plant Sci. 12, 123. doi: 10.3389/fpls.2021.636339
van der Werff, H., and Richter, v. d. H. (1996). Toward an improved classification of Lauraceae. Annal. Missouri Botanic. Gard. 6, 409–418. doi: 10.2307/2399870
van Esse, H. P., Van't Klooster, J. W., Bolton, M. D., Yadeta, K. A., van Baarlen, P., Boeren, S., et al. (2008). The Cladosporium fulvum virulence protein Avr2 inhibits host proteases required for basal defense. Plant Cell 20, 1948–1963. doi: 10.1105/tpc.108.059394
van Loon, L. C., Rep, M., and Pieterse, C. M. (2006). Significance of inducible defense-related proteins in infected plants. Ann. Rev. Phytopathol. 44, 135–162. doi: 10.1146/annurev.phyto.44.070505.143425
Vartapetian, A., Tuzhikov, A., Chichkova, N., Taliansky, M., and Wolpert, T. (2011). A plant alternative to animal caspases: subtilisin-like proteases. Cell Death Different. 18, 1289–1297. doi: 10.1038/cdd.2011.49
Veronese, P., Nakagami, H., Bluhm, B., AbuQamar, S., Chen, X., Salmeron, J., et al. (2005). The Membrane-Anchored BOTRYTIS-INDUCED KINASE1 plays distinct roles in arabidopsis resistance to necrotrophic and biotrophic pathogens. Plant Cell 18, 257–273. doi: 10.1105/tpc.105.035576
Vetukuri, R. R., Whisson, S. C., and Grenville-Briggs, L. J. (2017). Phytophthora infestans effector Pi14054 is a novel candidate suppressor of host silencing mechanisms. Euro. J. Plant Pathol. 149, 771–777. doi: 10.1007/s10658-017-1222-9
Voigt, C. A. (2014). Callose-mediated resistance to pathogenic intruders in plant defense-related papillae. Front. Plant Sci. 5, 168. doi: 10.3389/fpls.2014.00168
Warnes, G. R., Bolker, B., Bonebakker, L., Gentleman, R., Huber, W., Liaw, A., et al. (2009). gplots: Various R programming tools for plotting data [Online]. Available: https://cran.r-project.org/web/packages/gplots/index.html (accessed May 5, 2021).
Westermann, A. J., Gorski, S. A., and Vogel, J. (2012). Dual RNA-seq of pathogen and host. Nat. Rev. Microbiol. 10, 618–630. doi: 10.1038/nrmicro2852
Win, J., Morgan, W., Bos, J., Krasileva, K. V., Cano, L. M., Chaparro-Garcia, A., et al. (2007). Adaptive evolution has targeted the C-terminal domain of the RXLR effectors of plant pathogenic oomycetes. Plant Cell 19, 2349–2369. doi: 10.1105/tpc.107.051037
Wolf, T., Kämmer, P., Brunke, S., and Linde, J. (2018). Two's company: studying interspecies relationships with dual RNA-seq. Curr. Opin. Microbiol. 42, 7–12. doi: 10.1016/j.mib.2017.09.001
Wolstenholme, B. N., and Sheard, A. (2010). “Integrated management of phytophthora root rot the “pegg wheel,” in Updated. South African Avocado Growers' Association Avoinfo Newsletter [Online], 175. Available online at: https://www.researchgate.net/profile/Andrew_Sheard/publication/262494179_INTEGRATED_MANAGEMENT_OF_PHYTOPHTHORA_ROOT_ROT_THE_PEGG_WHEEL_UPDATED/links/0c960537dddb151fdc000000.pdf (accessed December 01, 2020).
Wu, M., and Gu, L. (2021). TCseq: Time course sequencing data analysis [Online]. Available online at: https://bioconductor.org/packages/release/bioc/html/TCseq.html (accessed August 21, 2021).
Xia, Y., Ma, Z., Qiu, M., Guo, B., Zhang, Q., Jiang, H., et al. (2020). N-glycosylation shields Phytophthora sojae apoplastic effector PsXEG1 from a specific host aspartic protease. Proceed. Nat. Acad. Sci. 117, 27685–27693. doi: 10.1073/pnas.2012149117
Xia, Y., Suzuki, H., Borevitz, J., Blount, J., Guo, Z., Patel, K., et al. (2004). An extracellular aspartic protease functions in Arabidopsis disease resistance signaling. EMBO J. 23, 980–988. doi: 10.1038/sj.emboj.7600086
Xiang, G., Yin, X., Niu, W., Chen, T., Liu, R., Shang, B., et al. (2021). Characterization of CRN-Like Genes From Plasmopara viticola: Searching for the Most Virulent Ones. Front. Microbiol. 12, 53207. doi: 10.3389/fmicb.2021.632047
Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11, 1–12. doi: 10.1186/gb-2010-11-2-r14
Yu, X., Feng, B., He, P., and Shan, L. (2017). From Chaos to harmony: responses and signaling upon microbial pattern recognition. Ann. Rev. Phytopathol. 55, 109–137. doi: 10.1146/annurev-phyto-080516-035649
Zhang, W., Fraiture, M., Kolb, D., Löffelhardt, B., Desaki, Y., Boutrot, F. F., et al. (2013a). Arabidopsis receptor-like protein30 and receptor-like kinase suppressor of BIR1-1/EVERSHED mediate innate immunity to necrotrophic fungi. Plant Cell 25, 4227–4241. doi: 10.1105/tpc.113.117010
Zhang, X., Han, X., Shi, R., Yang, G., Qi, L., Wang, R., et al. (2013b). Arabidopsis cysteine-rich receptor-like kinase 45 positively regulates disease resistance to Pseudomonas syringae. Plant Physiol. Biochemistr. 73, 383–391. doi: 10.1016/j.plaphy.2013.10.024
Zimmermann, D., Gomez-Barrera, J. A., Pasule, C., Brack-Frick, U. B., Sieferer, E., Nicholson, T. M., et al. (2016). Cell death control by matrix metalloproteinases. Plant Physiol. 171, 1456–1469. doi: 10.1104/pp.16.00513
Keywords: dual RNA-sequencing, Phytophthora cinnamomi, avocado, receptor-like kinase, protease, time-course
Citation: Backer R, Engelbrecht J and van den Berg N (2022) Differing Responses to Phytophthora cinnamomi Infection in Susceptible and Partially Resistant Persea americana (Mill.) Rootstocks: A Case for the Role of Receptor-Like Kinases and Apoplastic Proteases. Front. Plant Sci. 13:928176. doi: 10.3389/fpls.2022.928176
Received: 25 April 2022; Accepted: 25 May 2022;
Published: 28 June 2022.
Edited by:
Andreia Figueiredo, University of Lisbon, PortugalReviewed by:
Hugo Germain, Université du Québec à Trois-Rivières, CanadaSusana Traquete Serrazina, University of Lisbon, Portugal
Copyright © 2022 Backer, Engelbrecht and van den Berg. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Noëlani van den Berg, bm9lbGFuaS52ZGJlcmcmI3gwMDA0MDtmYWJpLnVwLmFjLnph