Comparative transcriptome analysis of R. necatrix growing on avocado roots vs PDA medium.
A transcriptome analysis was carried out to capture genes expressed during necatrix growth on susceptible `Dusa´ avocado roots and on PDA medium, in order to compare their expression profiles (Fig. 1). The RNA-Seq data including the raw reads from three biological replicates of R. necatrix CH53 virulent strain colonizing avocado roots (RGA1; RGA2 and RGA3) and growing on culture medium (RGPDA1; RGPDA2 and RGPDA3) were processed. A total of 12,104 transcripts were obtained, among which 11,807 were present in both conditions, while 137 and 160 transcripts were exclusively expressed in either RGA or RGPDA, respectively (Fig. 2). Total transcripts were subjected to statistical analysis to evaluate differential gene expression between RGA vs RGPDA test situations. Analyses resulted in 1,937 differentially expressed genes (DEG), 61.9% induced and 38.1% repressed (-2 > fold change (FC) > 2; P-value < 0.05) (Fig. 3). A heat map of DEGs showed consistence in expression patterns among RGA1, RGA2 and RGA3 and among RGPDA1, RGPDA2 and RGPDA3, supporting the reliability of the RNA-Seq data (Fig. 4).
Validation of the RNA-Seq analysis
Differences found in gene expression profiles between RGA vs RGPDA were further verified through a quantitative real time PCR (qRT-PCR) assay on total cDNA samples from mycelia of three biological replicates. For this, five randomly selected genes over-expressed in RGA vs RGPDA and with different FC, were analyzed. Actin gene was used as reference gene for data normalization. The expression levels of these genes amplified by qRT-PCR are shown in Table 1. Although higher expression values were obtained by qRT-PCR than those observed on the RNA-Seq, results corroborated the overall differences found between the two samples (RGA and RGPDA) in the RNA-Seq analysis.
Functional annotation and pathways analysis of differentially expressed genes (DEGs)
To better understand the infection process of necatrix colonizing susceptible avocado roots, all differentially expressed genes were functionally enriched and categorized based on blast sequence homologies and gene ontology (GO) annotations using Blast2GO software [18] (P< 05), selecting the NCBI blast Fungi as taxonomy filter and default parameters. DEGs were significantly grouped into the regulation of eight molecular function (MF), such as heme binding (GO:0020037), iron ion binding (GO:0005506), oxidoreductase activity acting on CH-OH group of donors (GO:0016614), flavin adenine dinucleotide binding (GO:0050660), cellulose binding (GO:0030248), NADP binding (GO:0050661), peroxidase activity (GO:0004601) and N,N-dimethylaniline monooxygenase activity (GO:0004499), and three biological process (BP), such as carbohydrate transport (GO:0008643), cellular oxidant detoxification (GO:0098869) and mycotoxin biosynthesis (GO:0043386) (Fig. 5A). To identify processes and functions over-represented in R. necatrix during infection, GO term enrichment analysis was also applied to the Top 100 over-expressed genes (Fig. 5B). The functions of these DEGs were significantly enriched in the regulation of five BP, such as oxido-reduction process (GO:0055114), cellulose catabolic process (GO:0030245), mycotoxin biosynthesis (GO:0043386), glucose import (GO:0046323) and response to hydrogen peroxide (GO:0042542), and 13 MF (Fig. 5B) among which activities related to plant cell wall degradation, including glucosidase activity (GO:0015926); endo-1,4-β-xylanase activity (G0:0031176); cellulose 1,4-beta-cellobiosidase activity (GO:0016162); xyloglucan-specific exo-β-1,4-glucanase activity (GO:0033950) and arabinogalactan endo-1,4-β-galactosidase activity (GO:0031218) were found.
To investigate the metabolic pathways affected in R. necatrix during avocado root infection, a KEGG pathway analysis was performed with Blast2go [18]. For the total of 1,937 DEGs, 100 metabolic pathways that involved 208 genes were identified (P-value <0.05). The metabolic pathways were reorganized into eleven categories (Table 2) being the nucleotides metabolism the one with the highest number of genes (n=64). Interestingly, metabolic pathways involved in antibiotic and drug metabolism were also affected, in accordance with GO enrichment analysis results, where mycotoxin biosynthetic process was one of the molecular functions over-represented.
Candidate genes involved in the pathogenesis of R. necatrix
At least 69 transcripts showing homology to genes previously reported to be involved in fungal infection were identified among the 1,937 DEGs. These include homologs to genes involved in the production of CWDE (Table 3), proteases, fungal toxins, detoxification and transport of toxic compounds, gibberellin biosynthesis and gene silencing (Table 4) as well as gene effectors (Table 5). Out of the 69 selected genes, 30 were associated with cell wall hydrolysis, among which 16 showed fold change (FC) values above 50, with three of them (SAMD00023353_0503130, SAMD00023353_6500680 and SAMD00023353_4001240) allocated in the top20 over-expressed genes in R. necatrix during avocado root-colonization (Table 3 and Additional file 1: Table S1). Five genes were identified as proteases, two aspartic proteases and three serine proteases, with the contig SAMD00023353_1500930 expressed over 411 times in RGA vs RGPDA (Table 4). Five contigs showed homology to genes encoding fungal toxins, among which the contig SAMD00023353_5500610 encoding the putative aflatoxin B1 aldehyde reductase member 2 showed the higher transcript abundance with a FC value of 18.65 (Table 4).
Nineteen genes were related to degradation of toxic compounds such as reactive oxygen species (SAMD00023353_5200870), aflatoxins (SAMD00023353_0902760, SAMD00023353_12800020, SAMD00023353_3200110), and antibiotics (SAMD00023353_3600430, SAMD00023353_6600160, SAMD00023353_0702510, SAMD00023353_0100280, SAMD00023353_2201610), among other drugs. R. necatrix also over-expressed genes related to transport of toxic compounds, in particular, four (SAMD00023353_2601150, SAMD00023353_2501030, SAMD00023353_3000620 and SAMD00023353_6200040) and two contigs (SAMD00023353_10000080 and SAMD00023353_2200710) showed homology with genes encoding ATP-binding cassette (ABC) transporters and major facilitator superfamily (MFS) transporters, respectively. Expression values of genes homologous to ABC transporters were higher (FC values ranging from 5 to 7) than those observed for MFS transporters (ranging from 2-3) (Table 4).
Two genes were selected for being associated with hormone biosynthesis (GA4 desaturase family protein SAMD00023353_10100030 and gibberellin 20-oxidase SAMD00023353_1901120) showing FC values of 38.2 and 2.39 respectively and one gene, the argonaute siRNA chaperone complex subunit Arb1 (SAMD00023353_0801000), postulated to play a role in RNA induced transcriptional silencing (Table 4).
The RNAseq analysis also revealed 137 genes only expressed in R. necatrix during its growth on avocado roots. From those contigs, 24 were predicted as candidate effector proteins (CEP) by the CSIRO tool EffectorP2 (a machine learning method for fungal effector prediction in secretomes) [19] with a probability above 60% (Table 5). All CEPs, except for SAMD00023353_2100110, SAMD00023353_2801560, SAMD00023353_3900800, SAMD00023353_11900020 and SAMD00023353_1700590, showed no similarity with proteins in the public database. Out of the 24 CEP, 13 were predicted to be secreted by SignalP3 server and ten were determined to have an apoplastic localization by the CSIRO tool ApoplastP (a machine learning method for predicting localization of proteins) [20] (Table 5).
To test any existing relationship within the candidate effectors proteins identified in this study with previously described effectors proteins, the PHI (Pathogen Host Interaction) database was used; i.e., PHI-base is a database of virulence and effector genes that have been experimentally proven via pathogen-host interaction [21]. Blastp was used to match PHI-base with an e-value cutoff of 1E-03 and 30% identity. As result, 3 R. necatrix candidate effectors were annotated, SAMD00023353_11900020 encoding a putative glycoside hydrolase, showed the higher percentage of identity with the effector Lysm from Penicillium expansum (Identity 44.58%, E-value 9.94 E-53). SAMD00023353_2100110 and SAMD00023353_1700590 showed identity with effectors BEC1040 and Mocapn7 from Blumeria graminis (Identity 32.76%, E-value 1.32 E-05) and Magnaporthe oryzae (Identity 35.82%, E-value 1.32 E-03), respectively.