Molecular evidence of the avocado defense response to Fusarium kuroshium infection: a deep transcriptome analysis using RNA-Seq

Strain, media, and culture conditions

The F. kuroshium HFEW-16-IV-019 strain (first identified as F. euwallaceae and later as Fusarium spp. associated with KSHB) (Ibarra-Laclette et al., 2017; Sanchez-Rangel et al., 2018) was provided by the Centro Nacional de Referencia Fitosanitaria (CNRF). All experiments were conducted following biosecurity protocols. Conidia were stored in 25% glycerol at −80 °C and propagated on potato dextrose agar (PDA) when necessary, and the plates were incubated for 5–7 days at 28 °C in darkness. Fungal spores were collected by gently shaking the plate with 3–5 mL of sterile water at room temperature. The suspension was collected and transferred to 2-mL Eppendorf tubes and centrifuged at 13,000 rpm for 15 min. The supernatant was decanted, and the conidia were washed twice with sterile water. Total conidia were counted under an optical microscope in a hemocytometer, and a suspension containing 1 × 10⁸ conidia/mL was prepared in sterile water.

Avocado stem inoculation

One-year-old healthy avocado trees (P. americana cv. Hass) were acquired from the Salas certified plant nursery at Uruapan, Michoacán, Mexico. The avocado trees were acclimated for 4–6 months in a greenhouse until used. Selected branches of approximately 30–40 cm in length and 1 cm in diameter were cut with garden scissors and cleaned with sterile water, and each stem was then cut into fragments of 3.5–4 cm in length. A hole was drilled into the center of each avocado stem with a 1/16″ Dremel^®. The injured stems were placed into humid chambers prepared by placing sterile filter paper discs in the bottom of Petri dishes (150 × 20 mm) and adding 4–5 mL of sterile distilled water. Each stem was inoculated with 40 µL of conidial suspension (1 × 10⁸ conidia/mL) or H₂O for the control. The humid chambers were maintained at 27 °C and 16 h light/8 h dark for 1, 4, 7 and 14 days post-inoculation (dpi).

RNA isolation, library preparation, and sequencing

Infected and uninfected stems (at 1, 4, 7 and 14 dpi) were frozen with liquid nitrogen and then pulverized using a mortar and pestle. Total RNA was isolated from each sample using PureLink^®\ Plant RNA Reagent (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s protocol. RNA concentration was measured on a NanoDrop 2000c spectrophotometer (Thermo Scientific, Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity (RIN ≥ 8.5) was evaluated on an Agilent 2100 Bioanalyzer using an RNA 6000 Nanokit (Agilent Technologies, Santa Clara, CA, USA). Then, 200 ng of total RNA from each selected sample was used to generate the corresponding sequencing libraries. We used the RNA-seq method based on poly(A) selection, which enriches eukaryotic mRNA and other polyadenylated RNA molecules. Three independent biological replicates from each selected sample were processed. Using the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) and index codes to identify each sample independently, RNA sequencing libraries were generated in the Massive Sequencing Unit of the Ecology Institute (INECOL, Veracruz, México) located in Xalapa, Veracruz, Mexico, using Illumina sequencing instrumentation (NextSeq500). Multiplexed libraries were sequenced on a single flow cell of the NexSeq500 platform (Illumina, San Diego, CA, USA) to generate 150 bp paired-end reads. The RNA-seq data have been deposited in the Short Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) with accession number PRJNA655460.

Assembly and identification of protein-coding regions in RNA transcripts

Before the assembly process, low-quality reads (those with an average Phred score of less than 30 and in which only 10% of the bases that belonged to the reads had a Phred score below the quality permitted (20)) were removed with a Phyton-based script (https://github.com/Czh3/NGSTools/blob/master/qualityControl.py). Prior to de novo assembly, screening was performed to identify fungus-like sequences using the Hisat2 v2.1.0 program (Kim, Langmead & Salzberg, 2015). The sequences used as references were a unique transcript collection (Sanchez-Rangel et al., 2018) and the complete genome sequence (Ibarra-Laclette et al., 2017) available in the GenBank public database.

High-quality (HQ) reads were divided into two data sets of avocado sequences and fungus-like sequences. HQ paired-end reads obtained from the libraries of the uninfected samples plus sequences from infected samples that did not match F. kuroshium were combined and assembled using Trinity v2.4.0 software (Grabherr et al., 2011) to produce contigs (UniGenes). The UniGenes were assayed with SeqClean software (https://sourceforge.net/projects/seqclean/files/) to remove poly A/T tails, ends rich in Ns (undetermined bases) and/or low-complexity regions from the sequences. Only the resulting UniGenes with a minimum length of 200 bp were selected and processed with AlignWise (Evans & Loose, 2015). After correcting the frameshifts, a data set for the avocado transcriptome with coding sequences (CDS) and their corresponding proteins was created. Redundant sequences were eliminated using the BlastClust algorithm (https://www.ncbi.nlm.nih.gov/Web/Newsltr/Spring04/blastlab.html). A CDS was considered redundant if it showed an identity of at least 95% over 90% or more of the length of the sequences being compared. Only CDSs that, once translated, generated peptides or proteins longer than 75 nucleotides were retained and used for further analyses.

Functional annotation of avocado UniGenes

To annotate the avocado transcriptome, we searched homologous protein sequences and performed a functional classification based on Gene Ontology (GO) terms (Ashburner et al., 2000). First, the resulting proteins translated from the assembled UniGenes were compared with a reference database using the BLASTp algorithm (e-value 10⁻⁵) (Altschul et al., 1990). The list of plant species used as references includes representative, well-annotated species that were recently used (Rendon-Anaya et al., 2019) to inform ancient evolutionary relationships between avocado and some of the major land plant lineages eudicotyledons: Amborella trichocarpa (basal angiosperms), Solanum lycopersicum (asterids, Solanales), Vitis vinifera (rosids incertae sedis, Vitales), Prunus persica (rosids, Rosales), and Arabidopsis thaliana (rosids, Brassicales); monocotyledons: Musa acuminata (commelinids, Zingiberales) and Zea mays (commelinids, Poales). In addition, the putative protein domains contained within the translated UniGenes were identified using Hidden Markov model (HMM)-based searches against the Pfam database (e-value 10⁻²) (Finn et al., 2016; Punta et al., 2012). A PfamScan Perl-based script was used for this purpose (https://github.com/SMRUCC/GCModeller/tree/master/src/interops/scripts/PfamScan/PfamScan).

After homolog identification, avocado UniGenes were classified according to Gene Ontology terms into at least one of the three major categories (biological processes, molecular functions, and cellular components). These GO terms were attributed to avocado UniGenes based mainly on the identified A. thaliana homologs (https://www.arabidopsis.org/index.jsp). Finally, the avocado UniGenes were assigned to specific metabolic pathways using the Kyoto Encyclopedia of Genes and Genomes database (KEGG; Kanehisa et al., 2016). The KEGG Automatic Annotation Server (KAAS; https://www.genome.jp/tools/kaas/) was used for this purpose with the single-directional best hit (SBH) method to assign orthologues.

Orthologous genes were identified in silico using the OrthMCL v.2.0.9 pipeline (Li, Stoeckert & Roos, 2003). Then, considering an inflation value (Enright, Van Dongen & Ouzounis, 2002) (1.5 in this study) and using the Markov cluster (MCL) algorithm (Enright, Van Dongen & Ouzounis, 2002), groups of orthologues (and paralogs) can be inferred across multiple taxa. We considered a minimum input length of 30 amino acids in all compared proteins, as well as a threshold e-value of 10⁻¹⁰ in the BLAST step. This stringent cutoff value was chosen to avoid false-positive results. The list of species used as references was the same as that used in the BLAST step performed during the annotation process.

Differential expression analysis

To identify differentially expressed genes during the infection process in both avocado and F. kuroshium, HQ reads were sequentially mapped against a specific data set that was used as a reference. In the case of avocado, HQ reads (avocado-like sequences) from infected and uninfected stems (at 1, 4, 7 and 14 dpi) were independently mapped against UniGenes obtained from the de novo assembled transcriptome. Fungus-like HQ reads identified in infected samples (at 1, 4, 7 and 14 dpi) were independently mapped to the gene models predicted for the F. kuroshium genome. The aligner tool Bowtie2 v2.3.5.1 (Langmead & Salzberg, 2012) was used to map the HQ reads to the corresponding data set used as references. Subsequently, the RSEM v1.3.1 package (Li & Dewey, 2011) with default parameter settings was used to calculate, for each sample, the posterior mean estimates, 95% confidence intervals, and maximum likelihood abundance estimates or expected counts (EC) for each of the UniGenes/genes used as the reference. For the identification of differentially expressed UniGenes/genes, the Bioconductor tool with the DESeq v1.24.0 package (Anders & Huber, 2010) was used. The values from the RSEM package (EC) were input into DESeq, which models the over-dispersed Poisson count data using a negative binomial model to generate the base-mean of each UniGene/gene and of each of the samples, which were compared in pairs. DESeq also computes the “fold-change” (FC) as the ratio of the base-mean calculated in the control samples divided by the base-mean estimated from the treatment sample. The data for the uninfected stem segments in the avocado sample data set were compared against those for the infected stem segments at each of the sampling times (1, 4, 7 and 14 dpi). The FC was then log₂ transformed [log₂(FC)]. Upregulated genes had positive log₂ values (≥1), while downregulated genes had negative log₂ values (≤−1). As part of its statistical tests, DESeq provides adjusted p-values that are calculated from all UniGene p-values using the Benjamini-Hochberg adjustment procedure (Benjamini & Hochberg, 1995). In the present study, an adjusted p-value of ≤ 0.05 was considered the threshold for identifying differentially expressed UniGenes.

Although the steps for analyzing the avocado and F. kuroshium data sets were similar, comparisons of the samples in the F. kuroshium data set at 4, 7 or 14 dpi with those at 1 dpi were made independently. Again, both a log₂FC value ± 1.0 and an adjusted p-value of ≤0.05 were the criteria for identifying genes as differentially expressed genes.

Gene clustering analysis: hierarchical clustering and k-means

Gene expression differences and their resulting profiles over time were calculated based on FC increases or decreases at the different sampling times. The infected and uninfected samples were compared at each time point, and hierarchical clustering of the log₂-transformed data sets was performed across the different genes (and time points) using Euclidean distance and average linkage clustering. The clustering results were represented visually as a heat map dendrogram. In addition, the log₂(FC) values calculated for each gene and throughout the time course were also grouped by k-means analysis. A maximum of 1,000 iterations and a conservative number of clusters that allowed the inclusion of the largest number of differentially expressed genes (n = 6) were the options selected. Both hierarchical clustering and k-means analysis were performed using Genesis software version 1.8.1 (Sturn, Quackenbush & Trajanoski, 2002).

Functional enrichment analysis

g:Profiler (http://biit.cs.ut.ee/gprofiler/; (Raudvere et al., 2019)), a widely used web toolset for identifying functional categories enriched in gene lists, was used to identify the enriched GO terms in the list of avocado UniGenes identified as differentially expressed in response to the interaction with F. kuroshium. The g:SCS method was applied to perform multiple testing correction using the p-adjusted values with a ≤ 0.05 threshold.

Identification of avocado immune receptors

The coding of avocado UniGenes to pattern recognition receptors (PRRs) was performed based on their A. thaliana homologs (Nürnberger & Kemmerling, 2006) and the presence of the LysM motif in some of them (Bozsoki et al., 2020) was confirmed through the Pfam annotation.

To identify UniGenes that coded NLR proteins, the characteristic domains, namely, a variable N-terminal domain (TIR, CC or PWR8), an NBS domain and/or a C-terminal LRR domain, were recognized (Collier, Hamel & Moffett, 2011; Xiao et al., 2001; Zhong & Cheng, 2016; Monteiro & Nishimura, 2018). First, we identified all UniGenes that were grouped as orthologues of the well-characterized A. thaliana NLR receptors (Meyers et al., 2003). Then, using the information from the Pfam annotation, each of the UniGene sequences containing at least one of the previously mentioned domains was individually checked using the NCBI BLAST interface and the generated graphical summary that illustrates the domains contained in the sequences; these domains are identifiable in the Conserved Domain Database (CDD) (Marchler-Bauer et al., 2011).

Avocado NLR proteins were aligned in the SeaView v4.6.1 phylogeny package (Gouy, Guindon & Gascuel, 2010) using the MUSCLE v3.8.31 program (Edgar, 2004). Phylogenies were generated in a maximum likelihood (ML) framework using PhyML v3.0 software (Guindon et al., 2010; Guindon & Gascuel, 2003) and the following parameters: model: Le Gascuel (LG); branch support: approximate likelihood-ratio test (aLRT: (Anisimova & Gascuel, 2006)); amino acid equilibrium frequencies: model-given; invariable sites: optimized; across-site rate variations: optimized; tree searching operations: best algorithm of neighbor interchange (NNI) and subtree pruning and regrafting (SPR); and starting tree: Bio Neighbor-joining (BioNJ; (Gascuel, 1997)).

Finally, the 3D structures of avocado NLR proteins were modeled by the rigid body grouping method using the SWISS-MODEL workspace (http://swissmodel.expasy.org/; Arnold et al., 2006). The structures of the A. thaliana ZAR1 protein in its inactive and active forms (Protein Data Bank entries 6j5t.1.A.pdb and 6j5t.1.C.pdb, respectively) were used as templates. Once the proteins were independently modeled, all of them were superimposed and visualized using the PyMOL molecular graphics system, v 2.3.4 Schrödinger, LLC. In addition, two of these avocado NLR proteins were also modeled as part of the ZAR1-RKS1-PBL2 resistosome (PDB entry 6j5t; (Wang et al., 2019a)). To assess its accuracy, each model generated was checked by its global model quality estimation (GMQE) and quaternary structure quality estimation (QSQE) scores.

Fungal infection and disease progression

The avocado stems infected with F. kuroshium displayed deposits of a powdery, white-colored material along their surface after 7 dpi. This sugary exudate is recognized as a typical symptom of FD (Eskalen et al., 2013). To assess vascular damage, longitudinal cuts were made on all stems (infected and uninfected). As shown in Fig. 1, the uninfected stems displayed a response associated with the injury, showing a brownish zone around the mechanically damaged area, but the zone did not extend more than 0.2–0.4 mm, and the integrity of the vascular tissue was maintained throughout the entire stem. In contrast, the infected stems clearly showed symptoms such as brownish discoloration starting in the early infection stage (1 dpi). This symptom remained pronounced as the infection progressed, and the loss of tissue integrity was also observed at 7 and 14 dpi.

Construction of the set of UniGenes for avocado

Of the total paired reads generated (454,682,877), approximately 30.6% were identified as poor-quality reads and were removed. In addition, 7.82% (8,516,500, 6,163,902, 5,488,581 and 15,419,814 at 1, 4, 7 and 14 dpi, respectively) were identified as Fusarium-like sequences and were discarded. As expected, the highest number of Fusarium-like sequences were identified in the samples corresponding to the late stages of infection (Fig. S1; Table S1).

After filtering out low-quality reads and removing contaminant sequences, a total of 279,667,960 paired-end reads were included in the assembly process. The same sequences were used for mapping and the differential expression analysis. As a result of the assembly process, a total of 87,343 UniGenes were generated. The corrected, nonredundant sequence data set comprised a collection of CDSs and corresponding proteins from a total of 51,379 UniGenes whose sequence lengths ranged from 75 to 15,321 bp (average length of 740.17 bp) (File S1). Only this sequence data set was considered representative of the transcriptome of P. americana stems and was used in the downstream analyses.

Homolog search and functional annotation of avocado UniGenes

From the total number of nonredundant avocado UniGenes, 47,833 UniGenes (93.09%), showed high similarity (e-value ≤ 10⁻⁵) with proteins from at least one of the reference species (Table S2). As expected, the number of avocado UniGenes that could be annotated with respect to a reference species was distinct and depended on the quality of the predicted gene models of the reference genomes as well as their phylogenetic relationships (Fig. 2A; Table S2).

The unique accession numbers (or loci) that resulted from counting the homologs identified in each species that was compared with the avocado UniGenes could be considered an estimation of the minimal number of genes codified in the avocado genome or, at least, of the number of genes that represent the stem transcriptome. The number varied between 13,809 and 16,799 proteins (Table S2), depending on the species considered as a reference. According to the recently published avocado genome (Rendon-Anaya et al., 2019), this number of UniGenes represents approximately 62.17% of the total protein-coding genes identified in avocado. However, a better estimation of the completeness of the assembled transcriptome was obtained from the BUSCO v3.0.1 pipeline (Simao et al., 2015); 88.0% completeness was achieved considering a predefined set of 1,375 single-copy genes present in the Embryophyte clade.

Then, avocado UniGenes were classified and grouped into at least one of the three major GO categories. Each of the subcategories assigned was inherited based mainly on the GO annotations available (ftp://ftp.arabidopsis.org/home/tair/Ontologies/) for each A. thaliana protein, which were identified as homologs, and by the BLAST algorithm used in the first step of the annotation process. At least one GO term was assigned to a total of 43,493 avocado UniGenes, and the assignments included 2,393 unique GO terms from the “biological processes” category, 3,516 from “molecular functions” and 782 from “cellular components” (Table S3). Considering that more than one GO term can be assigned to a single gene, we estimated that after the annotation process, an average of four GO terms were allocated to each of the avocado UniGenes (Table S3). The most prominent GO terms in the biological process, molecular function and cellular component categories were “cellular metabolic process”, “catalytic activity”, and “intracellular”, respectively (Fig. 2B; Table S3). Finally, the avocado UniGenes were assigned to some metabolic pathways, resulting in the assignment of 3,763 unique KEGG orthology (KO) identifiers (KEGG numbers) and 2005 distinct enzyme commission numbers to a total of 20,588 avocado UniGenes (Table S4).

The last step in the annotation process was orthologue (and paralog) identification between avocado and the other angiosperm plant species used as references in this study. A total of 281,493 proteins (29,399 from A. trichopoda; 36,334 from P. americana; 44,070 from M. acuminata; 48,669 from Z. mays; 33247 from S. lycopersicum; 38,554 from V. vinifera; 30,498 from P. persica; and 22,474 from A. thaliana) were grouped in 27,751 orthogroups or OrthoMCL-defined protein families (Fig. S2; Table S5). Among them, 7,600 orthogroups were shared among all eight reference species; 257 orthogroups contained proteins from eudicotyledon plants (S. lycopersicum, V. vinifera, P. persica, and A. thaliana) but no monocots (M. acuminate and Z. mays) or proteins from ancient clades of flowering plants (A. trichopoda and P. americana) (Fig. S2; Table S5).

Gene expression profile changes in avocado stem segments in response to infection with F. kuroshium

In total, 2,403 differentially expressed UniGenes were identified (fold change = 2, Log₂FC = ±1), with an adjusted significant p-value of ≤ 0.05 for at least one of the sampling points (Fig. 3A; Table S6). The comparison of all the data sets indicated that the most differentially expressed UniGenes were identified at 1 and 14 dpi. At 14 dpi, a significant number of UniGenes showed significant changes in their expression profiles (Fig. 3B). As expected, principal component analysis showed the overall profile of differentially expressed UniGenes and demonstrated a clear separation between noncorrelated groups (explaining 50.4% of the variance); one group included the early stages of infection (1 and 4 dpi), and the other group represented later stages (7 and 14 dpi). There were pronounced differences between the groups, especially at 14 dpi (Fig. 3C).

Finally, to identify expression patterns throughout the infection process (1–14 dpi), we carried out a k-means clustering analysis. The preestablished number of groups to obtain was six (n = 6) (Fig. S3; Table S6). Cluster 1 included 310 UniGenes that showed a general trend of downregulated expression throughout the infection process, displaying a maximum level at 1 dpi and reaching a minimum level at 14 dpi. In contrast, the expression profiles of 416 UniGenes (cluster 4) showed an upregulated trend from 1 to 14 dpi, suggesting that these genes play an important role during pathogenesis. UniGenes belonging to clusters 5 and 6 (373 and 391, respectively) showed no significant changes in their expression profiles at 1, 4 and 7 dpi; however, at 14 dpi, they displayed opposite trends, with a downward trend in cluster 5 and an upward trend in cluster 6. Cluster 2 had the highest number of UniGenes (630) and displayed a slight decrease at 1 dpi, reached minimum levels at 7 dpi, and then increased at 14 dpi, returning to the expression levels displayed at 1 dpi. Finally, 284 UniGenes belonging to cluster 3 showed similar expression profiles at 1, 4 and 14 dpi, with a maximum peak of expression at 7 dpi (Fig. S3).

Functional enrichment analysis of differentially expressed avocado genes

The GO enrichment analysis of the differentially expressed UniGenes displayed interesting data, which are shown in a Manhattan plot (Fig. 4A; Table S7). In the biological processes category, 83 GO terms were identified as enriched functional categories, and we selected eight primary or highly representative categories. We suggest that these categories not only effectively explain the biotic processes related to the global changes in transcriptomic events but also summarize some redundant terms (at different GO levels). The eight categories were as follows: response to stress (GO:0006950, 536 UniGenes), response to hormones (GO:0006950; 263 UniGenes), oxidation-reduction process (GO:0055114; 275 UniGenes), organization and biogenesis of cell wall processes (GO:0071554; 104 UniGenes), response to abscisic acid (GO:0009737; 100 UniGenes), polysaccharide metabolic process (GO:0005976; 82 UniGenes), response to water deprivation (GO:0009414; 68 UniGenes), and photosynthesis (GO:0015979; 47 UniGenes). The UniGenes contained in each of these categories are listed in Table S8. In addition to these eight categories, Table S7 shows all GO terms enriched in biological processes, and the five categories with the highest number of differentially expressed UniGenes can be identified. Those five categories were the following: cellular process (GO:0009987; 1,371 UniGenes), metabolic process (GO:0008152; 1,133 UniGenes), organic substance metabolic process (GO:0071704; 1,160 UniGenes), cellular metabolic process (GO:0044237; 1,133 UniGenes), and primary metabolic process (GO:0044238; 1,061 UniGenes). Additionally, enriched KEGG pathways were identified for the avocado UniGenes (Table S9). These data revealed that biosynthesis of secondary metabolites (KEGG:01110; 1,008 UniGenes), phenylpropanoid biosynthesis (KEGG:00940; 279 UniGenes), flavonoid biosynthesis (KEGG:00941; 144 UniGenes), phenylalanine metabolism (KEGG:00360; 72 UniGenes), and photosynthesis (KEGG:00195; 26 UniGenes) were the most significantly enriched pathways (Table S9).

Moreover, in Figs. 4B–4E, the circle plots represent the genes included in the eight enriched categories that were selected as relevant in the avocado stem-F. kuroshium infection process. The radar chart shows the distribution of individual GO terms at 1, 4, 7 and 14 dpi, and no significant differences between the number of genes that were up- and down-regulated were observed. The GO terms that changed the most over all the infection stages were response to stress (GO:0006950), oxidation-reduction process (GO:0055114), and response to hormones (GO:0009725); 226, 112 and 101 UniGenes were upregulated and 198, 117 and 109 UniGenes were downregulated for these GO terms, respectively.

The UniGenes included in the eight selected enriched GO terms were also identified in the different clusters. The data showed 1,475 UniGenes in total, and we detected that 235 belonged to cluster 1, 404 to cluster 2, 138 to cluster 3, 328 to cluster 4, 188 to cluster 5 and 182 to cluster 6. To analyze the infection process in more detail, we surveyed the 12 UniGenes that showed the greatest change in expression level (6 UniGenes that were upregulated and 6 UniGenes that were downregulated) in each of the six distinct clusters of differentially expressed UniGenes with different expression patterns (Fig. S3). The UniGenes selected are listed in Table S10, and their relevance is described in the discussion section.

Avocado immune PRRs and NLRs

To contribute to knowledge related to immune receptors that may participate in the avocado defense response against F. kuroshium, we identified in the avocado data two types of receptors: the pattern recognition receptors (PRRs) and nucleotide-binding domain and leucine-rich repeat-containing proteins (NLRs), which active PAMP-triggered immunity (PTI) and effector-triggered immunity (ETI), respectively.

Regarding PRRs, either by the annotation of their corresponding homologs or based on the presence of some characteristic motifs, a total of 223 UniGenes were identified as transcripts codifying for PRRs (Table S11). Interestingly, only 10 of these UniGenes (4.48%) were differentially expressed (Table S12). The UniGenes UN009421, UN019209, UN021431 and UN025891 belong to cluster 2. UN009421 is homolog to CoAc2 (AT4G26180) and encodes a mitochondrial CoA transporter. UN025891 encodes a plasma membrane-localized leucine-rich repeat receptor kinase involved in brassinosteroid signal transduction. UN019209 and UN021431 UniGenes are homologs to a LysM motif-containing F-box protein (AT1G55000). The six remaining differentially expressed UniGenes showed distinct expression patterns and encoding to homologs of RKL1, namely, AT1G48480 (UN004766); PR5K; AT5G38280 (UN010960); and AT2G19130 (UN004495, UN020702, UN021172, and UN061464), which is a member of a multigene family that encode to no well characterized G-type lectin S-receptor-like serine/threonine-protein kinase.

Regarding NLRs, key players in plant defense which recognize pathogen effectors and trigger immunity (Van Wersch et al., 2020), a total of 716 UniGenes contained at least one of the NLR characteristic domains in their coding region (Fig. 5A; Table S13). Interestingly, one unique UniGene (UN000626) encoded a TIR-type NLR; this protein seems to be complete since the conserved domains were identified in the structure. Another unique UniGene (UN036979) was identified as a putative PWR8-type NLR; in this case, the protein sequence seems to be partial, and only the N-terminal domain PWR8-type (CCr) and the central NBS domain were identified. The rest of the sequences were classified as follows: full CC-NLR proteins, 46 UniGenes; CC-NBS- and CC- truncated versions, 44 and 66 UniGenes, respectively; NLR proteins in which none of the canonical variable N-terminal domains (TIR, CC or CCr) could be identified, 26 UniGenes; and partial sequences in which the translated coding region contained only an NBS domain or an LRR domain, 226 and 311 UniGenes, respectively (Fig. 5A; Table S13).

Considering that mainly CC-type NLR receptors were identified in our avocado transcriptome, we resolved the phylogenetic relationships within these families (Fig. 5B). Only UniGenes that encode full CC-NLR proteins were considered. In addition, specific functional motifs were also identified, such as the MADA motif, which is evolutionarily conserved and is present in a certain percentage (~20%) of the total CC-NLR protein-coding genes in plants. Interestingly, in avocado, the MADA motif was clearly identified in a high percentage of the full CC-NLR proteins (67.39%, 31 proteins; Fig. S4). In the remaining proteins, the MADA motif was absent or it was unclear whether it was present because short sequences were missing in the N-terminal domain and near the start codon.

We also identified the EDVID motif or the EDVID-like motif (Fig. S4). Based on the results presented here, avocado CC-NLR proteins, like those in other plant species (Lukasik & Takken, 2009; Meyers et al., 2003; Wróblewski et al., 2018), can be classified into at least three major groups based on this EDVID-like motif (Fig. S4).

Regarding the NBS domains, we observed conserved motifs (p-loop, kinase, RNBS, GLPL, MHD; Fig. S5) that are crucial for ADP/ATP binding and shifting (Danot et al., 2009; Van Ghelder et al., 2019) and are involved in intra- and extramolecular interactions (Moffett et al., 2002; Van Ooijen et al., 2008). It has been consistently proven that, despite some variations, these motifs are conserved in all NLR subfamilies (including those present in gymnosperms) and in the proteins of prokaryotes with an NBS domain (NB-ARC), which do not display upstream N-terminal domains (Yue et al., 2012). Interestingly, only four UniGenes that encoded full CC-NLR proteins were identified as differentially expressed in response to F. kuroshium infection (Fig. 5B; Tables S8 and S13). Based on their expression profiles (Fig. S3), two of them, UN003976 and UN001791, belong to cluster 3, which contains differentially expressed UniGenes with maximum upregulation at 7 dpi. UN003288 belongs to cluster 4, a group of genes whose expression pattern increased over time, and UN003220 belongs to cluster 1, corresponding to UniGenes that show a rapid increase in expression level at 1 dpi and then a gradual decrease in expression.

Then, we performed modeling of both their inactive and active forms based on the A. thaliana orthologue ZAR1 (Fig. 6). We noticed that the CC-type N-terminal domain was complete only in UN001791 and UN003288, and the MADA motif was clearly distinguishable in these UniGenes (Fig. S4). These two CC-NLR proteins were also independently modeled in a structure consisting of five monomers (Fig. S6), mimicking the core region of the ZAR1-RKS1-PBL21 reconstituted resistosome (Wang et al., 2019a).

Interestingly, an avocado orthologue to A. thaliana RKS1 protein was not identified in any of the OrthoMCL-defined protein families (Table S5). RKS1 is a member of an OrthoMCL-defined protein family (OrthoGroup02788) containing a total of 22 proteins (7 proteins from A. thaliana (including RKS1), 5 from P. persica, 7 from V. vinifera, and 3 from S. lycopersicum). Notably, orthologues were absent not only in avocado but also in monocots (Z. mays and M. acuminata) and A. trichopoda, a member of the most ancient known clade of flowering plants. The PBL21 protein is a member of OrthoGroup02770, another family that also contains a total of 22 orthologous proteins but is apparently present in all reference plants used in this study (4 from A. thaliana, 2 from P. persica, 1 from V. vinifera, 3 from S. lycopersicum, 2 from Z. mays, 6 from Musa acuminata, 1 from A. trichopoda, and 3 from avocado (P. americana)). Despite its presence in the avocado transcriptome, none of the orthologue UniGenes were identified as differentially expressed in response to F. kuroshium infection.

Identification of fungal transcripts present in avocado stems

According to the criteria (a log₂FC value ± 1.0 and an adjusted p-value of ≤ 0.05), a total of 57 F. kuroshium genes were identified as differentially expressed during the plant-fungi interaction and the development of infection (Table S14). Heat map and Venn diagram analyses at 4, 7 and 14 dpi showed that more than 50% of genes maintained a similar expression profile over time; moreover, 74% of the genes were upregulated and 26% were downregulated (Fig. 7A and 7B, respectively). Orthologues (and paralogs) of these F. kuroshium genes were identified in the genomes of Nectria haematococca (commonly referred to by its asexual name, F. solani), F. fujikuroi, F. oxysporum, F. verticillioides, F. pseudograminearum and F. graminearum. Interestingly, the copy numbers of these genes were variable. According to the annotation, fourteen differentially expressed F. kuroshium genes had no functional or descriptive assignment (uncharacterized proteins). For the rest of the genes, the enriched GO terms within the biological process category were summarized using the REVIGO tool (Supek et al., 2011). The summary showed 13 genes for “alcohol metabolism” and 15 genes for “endoplasmic reticulum-plasma membrane tethering” (Fig. 7C). Additionally, “negative regulation of catalytic activity”, “metabolism”, “carbohydrate metabolism”, “cellular aromatic compound metabolism” and “proteolysis” were represented by fewer than 4 genes each (Fig. 7C). To identify genes related to pathogenicity, each protein was compared by the BLASTp algorithm against the pathogen-host interaction database. The analysis showed that 21 expressed F. kuroshium genes (37%) had significant similarity (e-value ≥ 10-6) to some PHI-base sequences, and the mutant phenotype observed in other pathogenic fungi showed reduced virulence, suggesting that these genes play an important role during pathogenesis (Table S14).

The comparison of the three data sets (4, 7 and 14 dpi) in the Venn diagram indicated that 39 up- and 10 downregulated genes were shared among the three stages of infection, 4 downregulated genes were shared between 4 and 14 dpi, 1 upregulated gene was shared between 4 and 7 dpi and one gene was only upregulated at 4 dpi. The only two genes upregulated at 4 dpi were Scf14_Gene-0.24 and Scf41_Gene-21.18, which were annotated as “hypothetical protein like phosphatidyl synthase” and “aldehyde dehydrogenase”, respectively. Scf27_Gene-18.41 was downregulated at both 4 and 7 dpi and was annotated as “probable endoglucanase type K”; this gene showed the greatest decrease in expression on both days. Scf201_Gene-0.1, Scf38_Gene-0.53, Scf71_Gene-7.29 and Scf135_Gene-0.11 were downregulated at 4 and 14 dpi and were identified as “uncharacterized protein”, “ammonium transporter MEPa”, “CAP20-virulence factor” and “urea transporter DUR3”, respectively. Among the genes upregulated at all three stages of infection, Scf30_Gene-0.42 exhibited the greatest change in expression and was annotated as “alcohol oxidase”, and Scf38_Gene-1.54 and Scf9_Gene-3.23 displayed the highest expression at 4 dpi and were both annotated as uncharacterized proteins. Scf15_Gene-17.40, Scf63_Gene-21.71 and Scf25_Gene-1.60 were annotated as “hypothetical protein similar to extracellular invertase”, “serine-type carboxypeptidase F precursor” and “related to tetracycline resistance proteins”, respectively, and displayed the highest expression at 7 dpi. Finally, at 14 dpi, Scf2_Gene-8.62, “related to allantoate permease”; Scf61_Gene-6.97, “related to nonribosomal peptide synthetase”; and Scf14_Gene-8.44, “uncharacterized protein” exhibited the highest expression. The genes that were the most downregulated at all infection stages were Scf15_Gene-2.40, annotated as “hypothetical uncharacterized protein”; Scf7_Gene-2.51, annotated as “related to glycosyl hydrolase family protein”; and Scf13_Gene-3.51, annotated as “related to uracil permease”.

The avocado- Fusarium kuroshium pathosystem

FD is a disease that is independently caused by two distinct complexes formed by ambrosia beetles that are morphologically related to Euwallacea fornicatus (Eichhoff) (Coleoptera: Curculionidae: Scolytinae), the polyphagous shot hole borer (PSHB) and KSHB, and their respective fungal symbionts (Carrillo et al., 2016; Eskalen et al., 2012, 2013; Mendel et al., 2012; Na et al., 2017). The complexes have a broad host range, including landscape and agricultural trees, particularly avocado trees (Carrillo et al., 2016; Eskalen et al., 2013; Na et al., 2017). The relevance of FD has led to the identification of the fungal symbionts associated with each insect; for PSHB, the three pathogenic fungal symbionts are F. euwallaceae, G. euwallaceae and Paracremonium pembeum (Freeman et al., 2013; Lynch et al., 2016), and for KSHB, the associated fungi are F. kuroshium and G. kuroshium.

Here, we developed a reproducible, practical and low-cost model system consisting of Hass avocado stems inoculated with F. kuroshium. This system allowed us to control some biotic and abiotic variables to have confidence in the molecular data, which was the main goal of this research. F. kuroshium pathogenicity under in vitro or greenhouse conditions has been analyzed only by Na et al. (2017). They tested the pathogenicity of the fungus in 2-year-old avocado cv. Zutano and reported the mean lesion length of xylem discoloration; there was a clear difference from the control but less pathogenicity than was observed with F. euwallaceae. For our pathosystem, we selected the Hass avocado variety since it is the most important cultivated and marketed avocado variety in Mexico. We injured avocado stems to mimic the ambrosia beetle and to inoculate the fungus directly into the vascular tissue. Similar to a previous report (Na et al., 2017), we detected that the avocado stems displayed brownish discoloration starting in the early infection stages (1 and 4 dpi), and the discoloration increased over time (Fig. 1). In contrast to the data obtained by the Eskalen research group in whole trees, we detected a white exudate exclusively in the infected avocado stems at 7 and 14 dpi (Fig. 1). The white exudate has been recognized as a typical symptom of FD in some but not all host trees. This phenomenon may be associated with vessel occlusion, resulting in the formation of tyloses or with the deposition of gums composed of phenols, pectin, and lipids that are already recognized as a response of some trees to various stimuli, including mechanical damage and vascular pathogen infection (Inch et al., 2012; Veronica De et al., 2016). It will be interesting to characterize the precise chemical nature and properties of this white exudate associated with the response of Hass avocado stems to F. kuroshium infection in future studies and to analyze the implications of the anatomical changes during the production of this exudate. Additionally, it would be relevant to coinfect avocado stems with G. kuroshium and F. kuroshium to understand whether this coexistence potentiates the pathogenicity capabilities of both fungal pathogens.

Major gene categories involved in avocado defense responses to infection by F. kuroshium

The avocado transcriptomics data were generated mainly to investigate molecular changes related to the avocado defense response after inoculation with the F. kuroshium. Both groups of differentially expressed UniGenes (up- and downregulated) showed, in general, oscillatory behavior (Fig. 3). Overall, our analysis indicated that plant defenses were activated and that the defense response persisted in the avocado stems at later stages but was not sufficient to prevent infection from progressing. These increases and decreases in differential UniGene expression may fluctuate, similar to the changes in their expression profiles over time (Fig. 3; Fig. S3). Moreover, the fluctuations were also evident by the six clusters generated.

Our analysis was focused on the UniGenes within the different clusters that were well identified based on the A. thaliana annotation and showed their major expression changes. The UniGenes that displayed upregulation are as follows. In Cluster 1, an ABC transporter F family member 1-like protein was identified that belongs to one of the largest protein families and is implicated in various processes, such as pathogen response and phytohormone transport (Kang et al., 2011). Information regarding this type of transporter in bacteria related to ABC-F transporter function indicated that they play an important role in antibiotic resistance and ribosome protection (Ero et al., 2019); this result combined with our data suggests that in avocado, they play an important role in ribosome protection and in toxic compound resistance during infection with F. kuroshium. Another family of proteins that was likely involved in avocado infection was the peptidyl-prolyl cis-trans isomerase family (PPiase). These proteins are involved in protein folding processes and are recognized in at least three families, namely, cyclophilins, FK506-binding proteins (FKBPs) and parvulins (He, Li & Luan, 2004). In the A. thaliana genome, there are 29 PPiases (Romano, Horton & Gray, 2004), and their expression is induced by fungal infection (Godoy et al., 2000) and hormones (Romano, Horton & Gray, 2004). Furthermore, A. thaliana knock-out mutations in three immunophilins (AtCYP19, AtCYP57, and one FKBP) result in an increased susceptibility to P. syringae (Pogorelko et al., 2014). Similar results were observed when cotton cyclophilin was overexpressed in tobacco, which conferred higher tolerance to both salt stress and P. syringae infection (Zhu et al., 2011).

In cluster 2, UN028398 (not identified) was the most highly upregulated, followed by UN023639, which encodes the mediator complex subunit Med28. The mediator is a protein complex that functions as a moderator between transcription factors and transcription activation by RNA polymerase II (Björklund & Gustafsson, 2005; Kornberg, 2005). Mediator studies in A. thaliana plant Med12, Med25 and Med8 mutants showed that the mutation enhanced susceptibility to Botrytis cinerea and A. brassicicola. Additionally, the Med16 mutant was defective in cold acclimation to freezing temperatures and was compromised in salicylic acid-induced PR gene expression as well as resistance to P. syringae (Wathugala et al., 2012).

Moreover, Med16 plants had reduced PR-1 levels and exhibited impaired systemic acquired resistance (Fu & Dong, 2013). In cluster 3, UniGene UN003714, encoding glutamine-dependent asparagine synthase 1 (ASN1), was the most upregulated. ASN1 is involved in the synthesis of L-asparagine from glutamine and has an important function in nitrogen metabolism and defense responses to microbial pathogens; in transgenic A. thaliana plants, overexpressing pepper asparagine synthetase resulted in enhanced resistance to P. syringae. Moreover, ASN1 induction can have opposite effects against the hemibiotrophic bacterial pathogen Xanthomonas campestris depending on the time post-infection (Hwang, An & Hwang, 2011). Finally, it has been proposed that ASN1 is involved in metabolic changes that facilitate host cell death during plant-pathogen interactions (Seifi et al., 2014). In cluster 4, the UniGene encoding multidrug resistance-associated protein 10 was the most upregulated. Multidrug resistance-associated proteins (MRPs) are a subfamily of full-molecule ABC transporters (Rea, 2007) that have an integral role in the plant detoxification mechanism (Klein, Burla & Martinoia, 2006). They work via the sequestration vacuole (“excretion storage”) to remove both endogenous and xenobiotic compounds such as glutathione conjugates (Klein, Burla & Martinoia, 2006; Rea, 2007).

Moreover, two maize MRPs are important in the transport of the anthocyanin pigment into the vacuole (Goodman, Casati & Walbot, 2004). Additionally, MRPs can transport other amphipathic anions, including glucuronate conjugates, linearized tetrapyrrole catabolites, folate and its derivatives (Rea, 2007). This evidence suggests the important role of MRPs in the plant-fungus interaction. The UniGene in cluster 5 with the highest expression level encodes a pentatricopeptide repeat (PPR) superfamily protein and has been recognized for its role in plant growth and development (Xing et al., 2018). Additionally, significant functions for several PPRs in response to biotic and abiotic stresses have been identified, such as functions in photooxidative stress responses (Karpinski et al., 1999) and ROS stress responses (Zsigmond et al., 2008) as well as in isoprenoid biosynthesis. Finally, the most highly upregulated UniGene in cluster 6 encodes a kinesin-13A (a microtubule-based motor protein) that participates in regulating the branching pattern of leaf trichomes (Li et al., 2017; Lu et al., 2005) in A. thaliana and Gossypium hirsutum. Trichomes serve as physical barriers against insect attack and fungal infection (Levin, 1973). Taken together, these data on the upregulated UniGenes suggest that in avocado stems during F. kuroshium infection, part of the defense response is related to both detoxification and protection mechanisms.

In addition, we analyzed the genes that displayed the greatest downregulation in each of the six clusters identified. In cluster 1, a nucleolar protein, gar2 (UN023639), named NSR1 in yeast and exhibiting high similarity to parallel1 (Par1; (Petricka & Nelson, 2007)), was found; in A. thaliana, parl1 mutants exhibit auxin-dependent growth defects and aberrant leaf venation. No direct information regarding its function during interactions with fungal pathogens has been obtained. In cluster 2, UN083934 and UN004830 encoded ubiquitin 11. Ubiquitin 11 is part of the ubiquitin–proteasome system, and its important function in plant-pathogen interactions is well documented (Miricescu, Goslin & Graciet, 2018). UN007206 was the most downregulated in cluster 3 and encodes light-regulated zinc finger protein 1, a transcriptional factor called LZF1, STH3, DBB3 and BBX22 (Khanna et al., 2009), and it modulates early chloroplast development and anthocyanin accumulation (Chang et al., 2008). These observations suggest that BBX22 could participate in the response to pathogens, since anthocyanins play an important role in protecting plant cells against damage by reactive oxygen species and are the first sign of programmed cell death (Dauphinee et al., 2017). In cluster 4, the UniGene for PHE ammonia lyase 1 (PAL) was the most downregulated. PAL catalyzes the first step in the phenylpropanoid pathway and is involved in the synthesis of many important metabolites, flavonoids, phenylpropanoids, lignin, and glyceollins (soybean-derived phytoalexins; Zhang et al., 2017). Many studies have demonstrated the important function of PAL during the defense response in different pathosystems (Abbasi, Graham & Graham, 2001; Cass et al., 2015; Kim & Hwang, 2014). In cluster 5, the most downregulated UniGene encodes a member of the amino acid permeases (AAPs) within the subfamily of amino acid/auxin permeases. AAPs are tissue-specific proteins that predominantly transport neutral and acidic amino acids with moderate affinity (Wan et al., 2017). No direct evidence relates AAPs to fungal biotic stress. However, it is obviously involved in important events related to the nitrogen mobilization pathway, in which amino acid transporters are associated with plant responses to pathogens (Sonawala et al., 2018). However, a protein with unknown function was the most downregulated in cluster 6, followed by MLP-like protein 423, a member of the major latex protein family; these proteins belong to the Bet v1 family, also known as the pathogenesis-related 10 (PR10)-like protein family (Wang et al., 2015). Recently, He and collaborators described the biological role of MLP-like protein 423, finding that its expression was inhibited by Botryosphaeria berengeriana f. sp. piricola and A. alternata apple pathotype infection (He et al., 2020). Additionally, the overexpression of MLP-like protein 423 in apple calli resulted in lower expression of resistance-related genes, and RNA-seq analysis indicated that MLP-like protein 423 negatively regulates apple resistance to infections by inhibiting the expression of defense and stress-related genes, as well as transcription factors. The identity of each of the UniGenes identified as downregulated suggests that the number of certain proteins is modulated through general regulation pathways, such as proteasome and transcription pathways. Finally, the combined data for the up- and downregulated UniGenes showed changes in avocado nitrogen metabolism.

Avocado immune receptors expression during F. kuroshium interactions and in the resistosome

Recognition of pathogen invasion through immune receptors is performed mainly by two structurally different proteins located on different subcellular compartments. One is the plasma membrane-localized PRRs that perceive extracellular pathogen-associated molecular patterns (PAMPs) activating PAMP-triggered immunity (PTI) (Boutrot & Zipfel, 2017). The other type of proteins is composed of intracellular receptors of the nucleotide-binding domain leucine-rich repeat protein superfamily (NLRs) (Jones, Vance & Dangl, 2016), which induces the effector-triggered immunity (ETI), which is often accompanied by hypersensitive response (HR) cell death. Together, these immune receptors act as a network of surveillance machines that recognize extracellular and intracellular pathogen invasion-derived molecules, ranging from conserved structural epitopes to virulence-promoting effectors (Kim & Castroverde, 2020).

Membrane-bound plant PRRs include receptor-like kinases (RLKs) (Shiu & Bleecker, 2003) that have an extracellular domain, such as leucine rich repeats (LRRs), lectin, lysine motif (LysM) or wall-associated kinases (WAK), with a single transmembrane spanning region and a cytoplasmic kinase domain;receptor- like proteins (RLPs) (Wang et al., 2008), which possess an extracellular LRR domain and a C-terminal membrane anchor but lack the cytoplasmic kinase domain; and polygalacturanase-inhibiting proteins (PGIP) (Di Matteo et al., 2003), which have only an extracellular LRR domain.

During F. kuroshium interaction in avocado, a total of 223 UniGenes were identified as PRRs, but only ten were selected as differentially expressed. Interestingly, UN025891 encodes a plasma membrane-localized leucine-rich repeat receptor kinase involved in brassinosteroid signal transduction (ATBRI1; AT4G39400). Additionally, it has been proven that disease resistance caused by some bacterial and fungus species is augmented in barley lines modified in the brassinosteroid receptor BRI1 (Ali et al., 2014). UN019209 and UN021431 UniGenes are homologs to a LysM motif-containing F-box protein (AT1G55000). These proteins (with F-box+LysM combination) are found not only in angiosperm plants but also in mosses and in green algae (Zhang, Cannon & Stacey, 2009). With a few exceptions, these proteins are single-copy genes and highly conserved in the plant kingdom, which implicates a conservation of their biological function, which is still unknown. Due to the presence of the LysM domain, it has been speculated that these proteins can recognize glycoproteins and presumably participate in their degradation (Zhang, Cannon & Stacey, 2009). Moreover, recently, in Lotus japonicas, receptors associated with nodulation factors and chitin were found to have a very similar structure but contain two diverging motifs in the LysM1 domain that are necessary for discriminating between immunity and symbiotic functions (Bozsoki et al., 2020).

The six remaining differentially expressed UniGenes associated with PRRs show distinct expression patterns and encoding to homologs of RKL1 (UN004766), PR5K (UN010960) and AT2G19130 (A. thaliana G-type lectin receptor kinases; UN004495, UN020702, UN021172, and UN061464). RKL1 is a salicylic acid-responsive gene that has been suggested to mediate defense responses (Ohtake, Takahashi & Komeda, 2000). In A. thaliana, the PR5K is a polypeptide with the structural and biochemical features of a transmembrane receptor protein that is possibly involved in the perception of microbial signals (including fungus) (Wang et al., 1996). Finally, the G-type lectin S-receptor-like serine/threonine-protein kinase (G-LecRKs) includes 38 members that are represented in A. thaliana (Teixeira et al., 2018) and contain an α-mannose binding bulb lectin domain. LORE is a member of G-LecRKs that can directly recognize bacterial lipopolysaccharide (Ranf et al., 2015). Therefore, the participation of the G-LecRks in the avocado defense processes may be relevant.

With respect to NLRs, four UniGenes (UN003976, UN001791, UN003288, and UN003220) encoding CC-type NLR proteins were identified as differentially expressed in response to F. kuroshium infection (Figs. 5B and 6; Table S13). Interestingly, these proteins are transcriptionally active mainly at early (1 dpi) and later (7 and 14 dpi) stages of infection. This activity pattern reflects their role during pathogen recognition, defense response activation and disease establishment, which provoked symptoms such as discoloration and vascular tissue damage (Fig. 1). These results are consistent with those of previous reports showing that genes encoding NLR proteins show distinct expression profiles and even tissue-specific expression patterns (Adachi, Derevnina & Kamoun, 2019); therefore, it is argued that plants have evolved tissue-specific NLR networks in a species-dependent manner, possibly to match specific pathogen attacks on different organs and tissues (Deng et al., 2017).

At least two of these four CC-type NLR proteins identified in avocado that were differentially expressed in response to F. kuroshium infection were modeled (Fig. S6) based on the structure of the ZAR1-RKS1-PBL21 A. thaliana resistosome (Wang et al., 2019a). The ZAR1-RKS1-PBL2 resistosome structure reveals how NLR immune receptors are activated. The active intermediate state of ZAR1 (Wang et al., 2019b) forms a wheel-like pentamer and resides in a 1:1 heterodimeric complex with RKS1. When a PBL2 protein (modified by the AvrAC effector of Xanthomonas campestris) is added to the ZAR1:RKS1 complex, it can be found in a different PBL2-induced conformation (Wang et al., 2019a).

Interestingly, these four avocado NLR receptors, which are orthologues to ZAR1, can be properly modeled with their active and inactive forms (Fig. 6) and even form a wheel-like pentamer (Fig. S6); no orthologues from RKS1 were detected in the avocado transcriptome. The results suggest that in avocado, as in other plant species (e.g., monocots (Z. mays and M. acuminata) and A. trichopoda), orthologues to RKS1 are absent. These observations are consistent with previous reports showing that ZAR1 can indirectly recognize other unrelated bacterial effector proteins (e.g., PBL2 modified by the AvrAC effector of X. campestris), all through an association with closely related pseudokinases that belong to receptor-like cytoplasmic kinase subfamily XII-2 (RKS1, ZED1 and ZRK3), respectively (Lewis et al., 2013; Seto et al., 2017; Wang et al., 2019a). Together, these results suggest that in avocado, even when some CC-type NLR proteins form a similar structure to the ZAR1-RKS1-PBL2 resistosome and even perform a similar function, specific proteins may be required to form a complex that represents the structure of the resistosome in its active intermediate state. Regarding avocado orthologues to PBL2, some of them (only partial sequences) can be identified in the avocado transcriptome; however, all of them were differentially expressed. This is expected considering that PBL2 needs to be modified by AvrAC effector recognition; until now, no plant serine/threonine-protein kinase (similar to PBL2) structure has been resolved in the presence of fungal effectors.

Regarding the domains and motifs present in CC-type NLR proteins, not only the highly conserved domains but also the motifs present on each of the domains play important roles in the structural conformation of resistosomes, and these motifs and domains are required for resistosome function. For example, RKS1 interacts with the ZAR1 LRR domain; it has been predicted that the LRR domain is the autoinhibitory region in the inactive state and is the sensor during the activation process (Burdett et al., 2019). The LRR domain was identified in all UniGenes encoding full CC-NLR proteins (Fig. 5A). Regarding the NBS domain, although the only known 3D structure other than ZAR1 is from the tomato protein NRC1 (Steele, Hughes & Banfield, 2019), we demonstrate that at least the differentially expressed avocado UniGenes that encode CC-type NLR proteins can be modeled in both their active and inactive forms (Fig. 6) when the resolved 3D structures for ZAR1 (PDB: 6j5t.1. A and 6j5t.1. C) are used as templates. In addition, the superposition of the modeled structure of these proteins shows that, as in the available structures from NRC1 and ZAR1, the P-loop (kinase a), kinase 2, and RNBS-A, -B, -C and -D motifs adopt analogous positions, and the orientation of the ADP nucleotide is virtually identical (Fig. 6) (Wang et al., 2019a). The similarities in the structures also extend to other conserved regions, such as GLPL and MHD motifs (Fig. S5; Fig. 6) (Van Ghelder et al., 2019).

Regarding the CC-type N-terminal domain, we note the presence of the EDVID motif (or the EDVID-like motif corresponding area) in all full CC-NLR proteins (Fig. S4). This motif has been shown to be important for intramolecular interactions between the NBS and LRR domains (Bentham et al., 2018), and it is necessary for both the autoinhibition and activation of the protein (Burdett et al., 2019). Finally, after confirming the presence of the MADA motif in a high percentage of the full CC-NLR proteins (Fig. S4), including the differentially expressed UN003976 and UN001791 (Fig. S6), we suggest that at least these two CC-type NLR proteins in avocado are involved in the defense mechanism induced against F. kuroshium infection.

Candidate genes involved in F. kuroshium pathogenesis

F. kuroshium is vectored by the ambrosia beetle and gains access to the plant to colonize the host xylem. According to Hulcr & Stelinski (2017), FD is the consequence of mass F. kuroshium accumulation on stressed trees. However, it remains unknown whether F. kuroshium contributes to disease development through an active pathogenic mechanism. In answering this question, it is important to consider that F. kuroshium belongs to a widely recognized genus of necrotrophic phytopathogens that can produce effectors, toxins and hydrolytic enzymes that, in combination, provoke host cell death (Ma et al., 2013). The notion of F. kuroshium as an active pathogen is supported by the fact that F. kuroshium infects healthy avocado tissue (Na et al., 2017) and was reaffirmed by the present study, in which the transcriptome coverage allowed us to identify 57 genes that are differentially expressed in F. kuroshium during the infection process. Interestingly, the analysis revealed that alcohol metabolism is strongly involved in the F. kuroshium infection process in avocado stems, as the most upregulated genes were associated with this GO term (Fig. 7). The enzymes encoded by these genes are alcohol oxidases, 5-diketo-D-gluconic acid reductase, aldehyde dehydrogenase, protocatechuate 3,4-dioxygenase beta subunit, choline dehydrogenase, and cellobiose dehydrogenase. There are few reports of these enzymes in filamentous fungi, which is additional evidence regarding their roles during pathogenesis. For example, alcohol oxidases are flavoenzymes that catalyze the oxidation of alcohols to carbonyl compounds, producing hydrogen peroxide (H₂O₂) (Couderc & Baratti, 1980; Cregg et al., 1989; Koch et al., 2016). However, the lack of an alcohol oxidase negatively impacts Cladosporium fulvum pathogenicity in tomato plants (Lycopersicon esculentum L.) (Segers et al., 2001). The role of this type of enzyme is not yet clear, and we hypothesize that it detoxifies the ethanol content that is produced in avocado tissue during infection and that H₂O₂ release may impact the redox environment. In our study, Scf30_Gene-042, which encodes an alcohol oxidase, was strongly upregulated during all three infection stages. Another interesting enzyme was aldehyde dehydrogenase (Scf41_Gene-21.18), which catalyzes the irreversible oxidation of endogenous and exogenous aldehydes to nontoxic carboxylic acids. In Magnaporthe grisea, two putative family-four aldehyde dehydrogenase genes were silenced by RNAi, severely compromising the pathogenesis of the rice blast fungus. Additionally, the mutant strains were highly sensitive to membrane stress (Abdul et al., 2018), suggesting that aldehyde dehydrogenases play a conserved role in sustaining membrane integrity by scavenging reactive aldehydes, fatty acid radicals, and other alcohol derivatives. All these compounds may be formed during membrane lipid peroxidation, which is triggered mainly by lipoxygenases and reactive oxygen species and is a hallmark of plant pathogen responses (Shah & Chaturvedi, 2008).

The protocatechuate 3,4-dioxygenase beta subunit is another enzyme identified here that participates in the catabolic pathway of protocatechuate (3,4-dihydroxybenzene, PCA); PCA is an aromatic compound and an intermediate product of the degradation of plant biopolymers such as lignin and other aromatic compounds (Brown et al., 2004). Finally, cellobiose dehydrogenase (Scf74_Gene-3.18) is an extracellular flavocytochrome produced by several wood-degrading fungi and appears to be involved in the degradation of both lignin and cellulose by fungi (Ma et al., 2017; Zhang, Fan & Kasuga, 2011). Together, these data highlight an active alcohol metabolic process that implies both lignin degradation in the avocado tissue and detoxification mechanisms to ensure fungal survival. To our knowledge, there is no information regarding the enzymes described above in the Fusarium genus.

Another interesting group of proteins that is well known to be involved in fungal pathogenicity is the hydrolytic enzymes. In the case of glycoside hydrolases, an endoglucanase type K (Scf27_Gene-18.41) showed the strongest downregulation at early stages of infection. The UniProt database describes that the catalytic activity of endoglucanase type K is the endohydrolysis of (1-4)-beta-D-glucosidic linkages in cellulose, lichenin and cereal beta-D-glucans; homologs in F. oxysporum (OMG_07426) and F. oxysporum f. sp. cubense tropical race 4 (FOIG_08964) were identified, but no further information was available. Interestingly, during the interaction between avocado and F. kuroshium, this protein was repressed; this is not surprising since some plants are capable of secreting inhibitors of certain pathogen hydrolytic enzymes, including glucanases, as part of their defense response (Misas-Villamil & Van der Hoorn, 2008). Therefore, it is possible that this mechanism occurs in avocado. A pectate lyase was also upregulated during the interaction (Scf61_Gene-6.25), and previous evidence clearly demonstrated the importance of this class of enzymes during plant-pathogen interactions since pectin is one of the major components of the plant cell wall. Pectin hydrolysis is a useful invasion or nutrition strategy for phytopathogens. In Alternaria brassicicola, a strain with a deletion of the PL1332 gene that encodes a pectate lyase was approximately 30% less virulent than the wild-type in Brassica oleracea (Cho et al., 2015). Moreover, purified VdPEL1, a pectate lyase of Verticillium dahliae, conferred resistance to Botrytis cinerea and V. dahliae in tobacco and cotton plants, and the mutant VdPEL1 lacked the ability to induce either cell death or plant resistance, revealing new insights into the role of a pectate lyase during host-pathogen interactions (Yang et al., 2018).

Additionally, proteases, other hydrolytic enzymes, are recognized as powerful tools that are produced by necrotrophic phytopathogens during infection processes; in this context, serine-type carboxypeptidase F (Scf18_Gene-0.42) was upregulated in F. kuroshium during the infection process. Serine proteases are hydrolytic enzymes that utilize serine to cleave peptide bonds in proteins. In general, fungi have an extensive repertoire of these enzymes, which are classified into more than 50 families and participate in several cellular activities, including nutrient degradation for subsequent assimilation and protection from the host’s immune system. The type of enzyme we detected in F. kuroshium belongs to the S10 superfamily that catalyzes extracellular degradation (Muszewska et al., 2017).

Like previous work in which RNA-seq was performed to document the response of F. kuroshium to growth at different pH levels (Sanchez-Rangel et al., 2018), transcripts for proteases, ABC transport and the genes that encode the enzymes involved in FA biosynthesis, a well-documented virulence factor present in other Fusarium species (Ding et al., 2018; Liu et al., 2020; López-Díaz et al., 2018), were identified. Interestingly, genes related to FA biosynthesis and another already characterized Fusarium toxin were not identified here, possibly because most of the coverage was related to the avocado transcriptome.

However, related to virulence factors, in the fungal gene set, a probable CAP20-virulence factor (Scf71_Gene-7.29) was recorded, although it was completely repressed in all the infection stages in our experiment. It is relevant that F. kuroshium is capable of encoding this virulence factor, since studies conducted in Colletotrichum gloeosporioides suggested that the perilipin homolog protein involved in the functional appressoria development affects virulence by reducing the penetration of the immature appressoria into the host cuticle (Lin et al., 2018). Moreover, in F. proliferatum, secretome profiling was performed during its interaction with banana, and a homolog of the CAP20 virulence factor was identified (CZR38705.1). Because the CAP20 virulence factor may be exclusive to fungi, it is a good candidate for a more detailed study to investigate its participation in F. kuroshium infection.

Finally, based on the annotation, 24% of differentially expressed F. kuroshium genes were described as uncharacterized proteins, and most of them showed strong upregulation in some of the three infection stages analyzed. It will be interesting to characterize these genes in greater detail in future studies with respect to their function during the infection process.