Deciphering the Infectious Process of Colletotrichum lupini in Lupin through Transcriptomic and Proteomic Analysis

The fungal phytopathogen Colletotrichum lupini is responsible for lupin anthracnose, resulting in significant yield losses worldwide. The molecular mechanisms underlying this infectious process are yet to be elucidated. This study proposes to evaluate C. lupini gene expression and protein synthesis during lupin infection, using, respectively, an RNAseq-based transcriptomic approach and a mass spectrometry-based proteomic approach. Patterns of differentially-expressed genes in planta were evaluated from 24 to 84 hours post-inoculation, and compared to in vitro cultures. A total of 897 differentially-expressed genes were identified from C. lupini during interaction with white lupin, of which 520 genes were predicted to have a putative function, including carbohydrate active enzyme, effector, protease or transporter-encoding genes, commonly described as pathogenicity factors for other Colletotrichum species during plant infection, and 377 hypothetical proteins. Simultaneously, a total of 304 proteins produced during the interaction were identified and quantified by mass spectrometry. Taken together, the results highlight that the dynamics of symptoms, gene expression and protein synthesis shared similarities to those of hemibiotrophic pathogens. In addition, a few genes with unknown or poorly-described functions were found to be specifically associated with the early or late stages of infection, suggesting that they may be of importance for pathogenicity. This study, conducted for the first time on a species belonging to the Colletotrichum acutatum species complex, presents an opportunity to deepen functional analyses of the genes involved in the pathogenicity of Colletotrichum spp. during the onset of plant infection.


Introduction
The distribution of Colletotrichum spp. is global. This anthracnose agent has been classified in the top 10 fungal pathogens of scientific and economic importance [1]. Within this genus, Colletotrichum lupini, responsible for lupin anthracnose, belongs to the Colletotrichum acutatum species for three days until germination. Inoculated plants were grown at 25 • C day and night, with 16 h of light and 8 h of darkness, and 60% humidity. Plants were watered every two days throughout the experiment. Symptoms were evaluated on all inoculated plants by measurement of the lesion length whenever visible. An ANOVA was performed to study the influence of biological repetition on necrosis length, using the RStudio software (version 3.2.5). Differences between variables were determined with multiple comparisons tests with the post hoc Tukey HSD method. The level of significance was set at α = 0.05.

Inoculation and Incubation Conditions
Twenty microliters of the inoculum adjusted at 1 × 10 7 spores·mL −1 were dropped at the base of the radicle emerging from the seed. Inoculated seedlings were put in Jiffy ® pots (peat pellets, 41 mm diameter) (Jiffy Products International BV, Zwijndrecht, Netherlands) for 24, 36, 48, 60, 72, 84 and 96 hpi. Sampling times were selected to target the first stages of the infection, the putative switch between biotrophic and necrotrophic phases and the necrotrophic phase. At 24, 48, 60, 72 and 84 hpi, 10 inoculated-seedlings were used for RNAseq analysis, while at 36, 60, 72, 84 and 96 hpi, 10 additional ones were used for proteomic analysis. Four biological repetitions were performed (except at 96 hpi, when three biological repetitions were done). Seedlings inoculated with water were used as negative controls.
Control conditions for proteomic and transcriptomic analyses were performed in vitro in a Czapek Dox liquid medium with 0.5 g·L −1 of sucrose instead of 30 g·L −1 (adapted from BD Difco™ CzaPEK-Dox Broth composition, BD, Franklin Lakes, NJ, USA) to simulate the low availability of carbohydrates on plant surface during fungal infection and germination. For each of the four repetitions, five Erlenmeyer flasks containing 100 mL of liquid medium were inoculated with 2 × 10 6 spores of C. lupini, as described previously. Flasks were then incubated for 24 h under agitation at 120 rpm at 25 • C before filtration for RNA extraction.
Microscope observations of morphological structures of C. lupini infecting plants were obtained from stem sections of two week-old lupin plants grown in vitro on Farheus medium under the same photoperiod as described previously. Under sterile conditions, cut ends were sealed with molten paraffin wax and droplets of spore suspensions at 10 5 spores·mL −1 were placed at intervals along the hypocotyl. After 12, 24, 48 and 72 hpi, epidermal strips were carefully removed, colored by lactophenol cotton blue and then placed in a drop of water for microscope observation. The extraction of RNA from pure cultures of C. lupini was performed from mycelium grown in Czapek Dox media with 0.5 g·L −1 of sucrose as described above. Mycelium was filtered through 0.45 µm cellulose acetate filter by Büchner filtration method and frozen in liquid nitrogen. One hundred milligrams of frozen mycelium was ground in a Mixer Mill (Retsch, MM400, Éragny, France) four times at a frequency of 30 Hz for one minute in lysing matrix A tube (MP biomedicals, Santa Ana, CA, USA). RNA extraction was performed with the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions. Samples were shaken using Mixer Mill (Retsch, MM400, Éragny, France) three times at a frequency of 30 Hz for one minute, after adding lysing buffer RLT. A DNase digestion step with RNase-Free DNase Set (Qiagen, Hilden, Germany) was added after the nucleic acid precipitation step and following the manufacturer's instructions.

Extraction of Total RNA from Inoculated Plants
At each time-point, infected plants were removed from Jiffy ® pots. Approximately 100 mg of the radicle of ten plants inoculated in the same conditions were cut with a scalpel around the inoculation or necrosis area. The ten pieces of plants were pooled together and frozen in liquid nitrogen. Samples were first ground using mortar and pestle, then using gentleMACS™ Dissociator (Miltenyi Biotec GmbH, Bergisch Gladbach, Germany) with gentleMACS™ C tubes three times at 4000 rpm for 20 s. The RNA extraction was performed with the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions. The total volume of lysing buffer RLT was adjusted to 1 g of ground plant corresponding to the pool of ten plants. A DNase digestion step with RNase-Free DNase Set (Qiagen, Hilden, Germany) was added after nucleic acid precipitation. Total RNA quality and quantity control was validated using Nanodrop1000 (Thermo Fischer Scientific, Waltham, MA, USA) and Agilent 2100 Bioanalyzer system analysis (Agilent technologies, Santa Clara, CA, USA) with RNA 6000 Nano LabChip ® Kit according to the manufacturer's instructions. The samples that were examined were expected to have a RIN (RNA Integrate Number) above eight, a total area under 18S and 25S curve above 50% and an rRNA ratio (25S/18S) above 1.8.

Preparation of Libraries, RNA Sequencing and Analysis
Total RNA extracted from liquid culture and inoculated plants were sent to Genome Quebec Innovation Centre (McGill University, Montreal, Canada) for sequencing. Libraries were generated from 250 ng of total RNA as follows: mRNA enrichment was performed using the NEBNext Poly(A) Magnetic Isolation Module (New England BioLabs, Ipswich, MA, USA). cDNA synthesis was achieved with the NEBNext RNA First Strand Synthesis and NEBNext Ultra Directional RNA Second Strand Synthesis Modules (New England BioLabs, Ipswich, MA, USA). The remaining steps of library preparation were done using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs, Ipswich, MA, USA). Adapters and PCR primers were purchased from New England BioLabs. Libraries were quantified using the Quant-iT™ PicoGreen ® dsDNA Assay Kit (Life Technologies, Carlsbad, CA, USA) and the Kapa Illumina GA with Revised Primers-SYBR Fast Universal kit (Kapa Biosystems, Wilmington, MA, USA). Average fragment size was determined using a LabChip GX (PerkinElmer, Waltham, MA, USA) instrument.
The libraries were normalized and pooled before being denatured in 0.05N NaOH and neutralized using HT1 buffer. ExAMP was added to the mix following the manufacturer's instructions. The pool was loaded at 200 pM on an Illumina cBot, and the flowcell was run on a HiSeq 4000 for 2 × 100 cycles (paired-end mode). A phiX library was used as a control and mixed with libraries at 1% level. The Illumina control software was HCS HD 3.4.0.38, the real-time analysis program was RTA v. 2.7.7. The bcl2fastq2 v2.20 program was then used to demultiplex samples and generate fastq reads. These data have been submitted to the ENA (European Nucleotide Archive) with the bioproject accession number PRJEB40331.
Read quality check and cleaning were performed on ABiMS Galaxy platform (galaxy.sb-roscoff.fr). Low quality reads were cleaned using Trimmomatic with the default parameters [36]. Paired reads were pseudo-aligned and quantified onto the sequenced and annotated C. lupini RB221 genome using Salmon package (Version 0.8.2) of Galaxy platform [37]. Genes which had a minimum of 100 counts per million (cpm) values in at least three repetitions were kept for further analysis.
In planta, gene expression for each time post-inoculation was first compared to the liquid culture condition to identify differentially-expressed genes (DEGs) during host interaction, using DESeq2. Then, a pairwise comparison between two close time points (i.e., 84 hpi vs. 72 hpi, 72 hpi vs. 60 hpi, 60 hpi vs. 48 hpi and 48 hpi vs. 24 hpi) was performed to determine if DEGs expression was significantly different using DESeq2. DEGs were selected with a false discovery rate (FDR) below 0.05 and a Log2-Fold-Change (LFC) higher than 1 or lower than −1. It was also verified that most genes expressed in planta were also found in vitro, ensuring that important genes were not discarded by using the liquid culture condition as a point of comparison.
A heatmap of DEGs was generated using the heatmap.2 function for the gplots R package based on normalized reads counts transformed to Z-score. Data were scaled by row, and Pearson correlation distances between genes were used for hierarchical clustering based on the expression profiles of DEGs.

Extraction and Digestion of Total Proteins
A pool of ten inoculated stems was ground in liquid nitrogen using a mortar and pestle. Proteins were extracted from 200 µL of stem powder using the TCA acetone protocol described in Méchin et al. (2007) [38].
Proteins were solubilized in 30 µL per mg of extract with a buffer containing 6 M urea, 2 M thiourea, 10 mM dithiothreitol (DTT), 30 mM Tris-HCl pH 8.8, and 0.1% zwitterionic acid labile surfactant (ZALS I, Proteabio, Morgantown, WV, USA) in water. Protein powders were mixed in the buffer before vortexing the tubes for 3 min. Samples were centrifuged (14,000 g, 25 min, 25 • C) and supernatants were transferred into new tubes. Protein concentrations were estimated using the plusOne 2DQuant Kit (GE Healthcare, Little Chalfont, UK), and adjusted to 3 µg·µL −1 prior digestion.
Digestion was performed in 0.2 mL strip tubes from 10 µL of diluted proteins representing 30 µg of proteins. Tubes were incubated for 30 min at room temperature for protein reduction by the 10 mM DTT present in the solubilization buffer. Then, 2 µL of a solution containing 300 mM of iodoacetamide (IAA) in 50 mM ammonium bicarbonate (BICA) was added, and proteins were alkylated by 1 h incubation at room temperature in the dark. After adding 90 µL of 50 mM of BICA, proteins were digested overnight with 800 ng of trypsin (Promega V5111, Promega, Fitchburg, WI, USA) at 37 • C. Trypsin digestion was stopped by adding 5.5 µL of 18.6% trifluoroacetic acid (TFA) representing 1% of the final concentration. Samples were incubated for 1 h at room temperature to allow TFA to cleave ZALS I. Peptides were desalted using C18 solid phase extraction (SPE) cartridges (strata XL 100 µm ref 8E-S043-TGB, Phenomenex, Torrance, CA, USA) as follows. The samples were first diluted in 2% acetonitrile (ACN), and 0.6% acetic acid in water (washing buffer), to a final volume of 500 µL. Then, the cartridges were conditioned with 500 µL of ACN and rinsed three times with 500 µL of washing buffer before loading the samples. Peptides were rinsed four times with 500 µL of washing buffer, and eluted two times by adding 300 µL of 40% ACN and 0.6% acetic acid in water (final pH around 2). Finally, eluted peptides were dried in a speedvac.

Orbitrap Mass Spectrometry and Analysis
A total of 19 protein digests were analysed with an Eksigent nlc425 nanoHPLC (SCIEX) coupled with a Q exactive+ mass spectrometer (Thermo Fisher, Waltham, MA, USA). Peptides were solubilized in 150 µL of a loading buffer containing 2% ACN and 0.1% formic acid (FA) in water.
For each injection, 4 µL (400 ng) of solubilized peptides were loaded onto a Biosphere C18 precolumn (particle size: 5 µm, pore size: 12 nm, inner/outer diameters: 360/100 µm, length: 20 mm; NanoSeparations, Nieuwkoop, Netherlands) and desalted for 4 min with the loading buffer at 7.5 µL·min −1 . Peptides were then separated on a Biosphere C18 column (particle size: 3 µm, pore size: 12 nm, inner/outer diameters: 360/75 µm, length: 300 mm; NanoSeparations) using buffer A (0.1% FA in water) and buffer B (0.1% FA in ACN) at 300 nL·min −1 as follows: (i) the column was equilibrated for 9 min with 95% of buffer A and 5% of buffer B; (ii) a linear gradient from 95% of buffer A and 5% of buffer B to 65% of buffer A and 35% of buffer B was applied for 110 min; (iii) the column was regenerated with 5% of buffer A and 95% of buffer B for 5 min. Electrospray ionization was performed at 1.8 kV with an uncoated capillary probe (noncoated capillary silica tips, 360/20-10, New Objective Inc., Woburn, MA, USA). S-lens RF level was set to 50.
Step 2 was repeated for the eight most intense ions detected in step 1 with the following criteria: intensity threshold superior as 8.3 × 103, precursor charge step of 2 and 3, dynamic exclusion of 50 s, peptide match on and exclusion of isotopes. Raw data were transformed to mzXML format using msconvert (proteowizard 3.0.7069, [39]).
Data were searched with X!Tandem (version 2015.04.01.1, [40]) against the RB221 strain C. lupini genome, Lupinus angustifolius cv. Tanjil genome (AOCW00000000.1) and a homemade database containing standard contaminants. Trypsin digestion was set in strict mode with one authorized missed cleavage. Cysteine carbamidomethylation was set as a fixed modification. Methionine oxidation, protein Nter acetylation with or without excision of methionine, Nter glutamine deamidation, Nter carbamidomethyl cysteine deamidation, and Nter glutamic acid dehydration were set as potential modifications. To identify additional peptides for proteins identified after this first pass, all samples were submitted to a second pass (refine mode of X!Tandem), in which the excision of Nter peptide signal, tryptophane oxidation, glutamine and asparagine deamidation were added as potential modifications. Additional missed cleavages were authorized at this step. Protein inference was performed by using X!TandemPipeline (cpp) v0.2.18 [41] with the following parameters: peptide e-value less than 0.01, protein e-value less than 10 −5 , two identified peptide by protein. Inference was performed using all samples together. Using a reverse version of the Colletotrichum lupini and Lupinus angustifolius database as a decoy, the FDR was estimated by X!Tandem to 0.05% and 0.01% for peptide-spectrum match and protein identification, respectively.
Peptide quantification was performed on extracted ion currents (XIC) by using Masschroq 2.2.2 [42]. Only the isotope theoretically the most intense was considered for peptide quantification. Identification and raw data are publicly available at http://moulon.inra.fr/protic. To quantify proteins, we used two complementary approaches. We used XIC to quantify proteins based on peptide intensities. Only protein-specific peptides present in at least 85% of the injections and showing a significant correlation (r > 0.5) with the other peptides of the same protein were kept for further data analysis. Normalization was performed, taking into account peptide retention time as described in Lyutvinskiy et al. (2013) [43]. Only proteins quantified with at least two peptides were considered. Protein relative abundances were computed as the sum of the normalized values of the selected peptides. One-way analyses of variance were performed on log10-transformed abundance values. The resulting p-values were adjusted for multiple comparisons [44]. Proteins with adjusted p-value < 0.01 were considered as showing significant abundance variations. All data analyses were performed using the RStudio-software (version 3.2.5).

Transcriptome and Proteome Prediction and Analysis Pipeline
The genome of C. lupini RB221 was used as a reference sequence (https://mycocosm.jgi.doe.gov/ Collup1/Collup1.home.html). Predicted protein sequences of DEGs and identified proteins were used to (i) predict the function, (ii) classify them into functional classes and (iii) explore the metabolism induced during the infectious process.
The prediction of other protein sequences was performed using Blast2GO, which attributed Gene Ontology (GO) term, Interpro number and PFAM number to genes that showed conserved domains. Different families of transcription factors were identified based on conserved domains.
In addition, Necrosis and Ethylene-inducing Proteins (NEPs) were predicted using conserved domains specific to NEP-like proteins (IPR008701/PF05630) [54]. Cytochrome P450 were identified via the family conserved domain (PF00067) [55]. Nudix proteins were identified by a family conserved domain (PF00293; IPR015797) [56]. The remainders were manually classified in putative functional groups according to Hmmscan, BLASTP and based on their predicted function. It is noteworthy that a pathogenicity-related class was created which included genes involved in stress responses and other pathogenicity factors such as genes involved in toxin synthesis. Nudix and NEP were further included in this class while cytochrome P450 genes were included in the oxidoreduction class. Proteins without predicted function, without predicted conserved domain, or with contradictory prediction were characterized as hypothetical proteins.
All predicted functions were further validated using Hmmscan program and Pfam database (version 33.0) with a cut-off e-value of 1.0 × 10 −3 , and a BLASTP searches to the GenBank nr database (release 233) using an e-valuethreshold of 1.0 × 10 −5 .
The agriGO v2.0 platform [57] was used to perform a customized Singular Enrichment Analysis (SEA) aiming to identify enriched GO terms in the up-and down-regulated DEGs, and in the detected proteins. Statistical analysis concerning RNAseq was performed using the predicted C. lupini proteome as reference on a query gene set which included genes that were up-or down-regulated in infected plants compared to the liquid culture condition. The percentage of DEGs in each GO category was calculated compared to the total of genes containing the same GO term in RB221 genome. Enriched GO terms were detected using the Fisher's exact test (multi-test adjustment method: Hochberg, FDR = 0.05).

Symptom Evaluation and Morphological Structures
The inoculation of lupin radicle with C. lupini first induced stem and petiole twisting as early as 60 hpi. Necrosis symptoms around the inoculation point began to appear at 60 and 72 hpi on a few plants, before being clearly visible at 84 hpi on 65% of plants ( Figure 1A,C). Necrosis size measured at 84 hpi significantly varied according to biological repetition (p = 1.36 × 10 −5 ). Appressoria were visible on epidermis of lupin hypocotyls as early as 12 hpi, and infection structures, like primary and secondary hyphae, were further observed at 48 hpi ( Figure 1B).    lupini DEGs compared to liquid culture (LC) (LFC < −1 and > 1, FDR < 0.05), grouped by putative functional class. Pathogenicity-related genes included stress response genes, Nudix, NEP and other functions described as associated with pathogenicity. Colored bars indicated gene expression by Z-score calculated from normalized reads of DEGs, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red). The green or black rectangles at the right of each gene represent genes respectively predicted as encoding secreted proteins or small secreted proteins.

Transcriptomic and Proteomic Data and Selection
Transcriptomic analysis was performed on RNA extracted from a liquid culture of C. lupini RB221, and from lupin hypocotyl inoculated with the same strain at 24, 48, 60, 72 and 84 hpi. For each sample, the total number of reads varied between 50.6 M and 85.5 M. Between 94.83% and 98.41% of total read pairs were kept after the cleaning step performed with Trimmomatic. Of these, 83.19% to 93.14% were paired and kept for further analysis (Table S1). The average mapping rate for the reads from RB221 liquid culture growth was 84.91%. The average mapping rate from infected plants ranged from 0.18% at 24 hpi to 3.54% at 84 hpi (Table S1).
With the aim of identifying C. lupini proteins produced during the onset of lupin infection, total proteins were extracted from lupin plants at 36, 60, 72, 84 and 96 hpi with RB221. Proteins were then identified and quantified by nLC Q exactive orbitrap MS-MS. Peptides quantified based on XICs led to the identification of 4377 putative proteins of both C. lupini and L. albus. Of these, 304 were identified and associated with C. lupini after filtering (below named QPs for quantified proteins) (Table S2).

Overview of Gene Expression and Protein Synthesis
Compared to in vitro liquid cultures, 897 genes were differentially expressed by C. lupini RB221 in planta for at least one sampling time. Among these genes, 593 were upregulated ( Figure 2A and Table S3) and 304 were downregulated ( Figure 2B and Table S3) compared to the liquid culture condition. At 24, 48, 60, 72 and 84 hpi, 63, 103, 223, 384 and 483 genes were respectively upregulated. In addition, 7, 1, 7, 54, and 170 genes were uniquely upregulated at 24, 48, 60, 72 and 84 hpi, respectively ( Figure 2A). In total, only 16 genes were uniquely expressed in planta (i.e., there were no counts in the in vitro condition, and more than one read for at least three out of the four repetitions in at least one sampling time), most of which encoding hypothetical proteins. By Blast2GO analysis, a total of 7132 genes with at least one GO term were found in the RB221 predicted proteome. Among the lists of DEGs and QPs, at least one GO term could be attributed to 331 upregulated genes, 163 downregulated genes and 245 proteins. Of these, respectively 18 and 26 GO terms were significantly enriched in the up-and down-regulated gene lists, as well as and 142 among QPs, respectively ( Figure 3).
Among the upregulated genes, the most significant GO terms related to biological processes, molecular functions and cellular components were respectively associated with the carbohydrate catabolic process, motor activity and cytoskeleton modifications. Among downregulated genes, they Upregulated genes compared to liquid cultures were detected from 24 to 84 hpi, while downregulated genes were detected later, from 60 to 84 hpi (Figure 2A,B).
For protein analysis, samples were collected 12 h after those for RNA analysis, to take into account the duration of translation steps (i.e., 36, 60, 72, 84 and 96 hpi). In total, 272 proteins from the 304 QPs were detected from 36 to 96 hpi. Among them, 16 were specific to the later stages, from 72 to 96 hpi ( Figure 2C). For 99.3% of the QPs, the highest abundance level was reached at the last sampling time, i.e., 96 hpi.
Among the 304 identified proteins, only 93 corresponded to DEGs identified in the transcriptomic approach ( Figure 2D).
By Blast2GO analysis, a total of 7132 genes with at least one GO term were found in the RB221 predicted proteome. Among the lists of DEGs and QPs, at least one GO term could be attributed to 331 upregulated genes, 163 downregulated genes and 245 proteins. Of these, respectively 18 and 26 GO terms were significantly enriched in the up-and down-regulated gene lists, as well as and 142 among QPs, respectively ( Figure 3).

Expression of the Predicted Functional Class of DEGs and QPs in Lupin Infection
Secreted proteins constitute mobile molecules involved in pathogenicity which can directly interact with the host. The transcriptomic and proteomic analysis respectively showed that the secretion of 103 DEGs and 27 QPs was predicted. Twenty six out of these secreted proteins were The three main GO categories represent biological processes (orange), molecular functions (green), and cellular components (yellow). For each GO term, the percentage of up-and down-regulated genes compared to the number of genes in the genome associated with the same GO term is represented. Significant GO terms associated with DEGs were ordered by increasing FDR value.
Among the upregulated genes, the most significant GO terms related to biological processes, molecular functions and cellular components were respectively associated with the carbohydrate catabolic process, motor activity and cytoskeleton modifications. Among downregulated genes, they were respectively associated with oxidative phosphorylation, transporter activities and membrane modifications.

Secreted Proteins, Small-Secreted Peptides (SSP)
Secreted proteins constitute mobile molecules involved in pathogenicity which can directly interact with the host. The transcriptomic and proteomic analysis respectively showed that the secretion of 103 DEGs and 27 QPs was predicted. Twenty six out of these secreted proteins were identified by both analyses. They contained a signal peptide, were localized in the extracellular space and had neither a TMH domain nor a GPI anchor; these three criteria were used as a definition of secreted proteins by   [12] ( Table 1). The number of secreted proteins from the DEG list was narrowed down to 46, which had less than 300 amino acid residues; of these, eight were also found in the QP list (Table 1).

Candidate Effectors
Effectors are SSP that enable the pathogen to manipulate host cell metabolism and overcome defense responses [58]. Among SSPs, 27 DEGs were predicted as effectors, among which, five were also found in the QP list ( Figure 4A and Table 2). The analysis of C. lupini genome led to the identification of 973 genes encoding secreted proteins, among which 294 were identified as candidate effectors, indicating that 9.1% of total candidate genes encoding effectors were differentially expressed from control conditions and 1.7% were proteins identified during the first 84 hpi. Among the candidate effectors predicted from the DEG list, five contained a putative conserved hypothetical protein domain also detected by mass spectrometry (CLUP02_14495), a LysM domain (CLUP02_16647), a ToxB N-terminal domain (CLUP02_01667), a fungal hydrophobin domain (CLUP02_07198) or a meiotically upregulated gene family domain (CLUP02_12726). One candidate effector (CLUP02_16164) matched with a Colletotrichum spp. hypothetical protein (from 95.6% to 97.3% of identity, e-value = from 6.0 × 10 −60 to 2.0 × 10 −75 ) but also with a putative secreted-in-xylem 11 (SIX 11) effector of Fusarium oxysporum in the second best result (79.28% identity, e-value = 8.0 × 10 −60 ) [59]. This candidate effector was also identified in the QP list. Most candidate effectors were upregulated compared to liquid cultures from 60 hpi to 84 hpi ( Figure 1D) and the number of upregulated candidate effectors increase strongly between 24 hpi and 48 hpi (Table 2). Five candidate effectors were identified by mass spectrometry, of which three were detected as early as 36 hpi, and their abundance increased steadily to 96 hpi. CLUP02_16164 and CLUP02_14495 were detected only after 72 hpi ( Figure 4B).  1  1  3  10  8  13  Methyltransferase  1  1  3  6  10  13  Mitochondria  1  1  2  2  1  2  Protein kinase  1  1  2  3  6  7  Transcription  4  3  3  6  4  8  Translation  3  5  9  10  5  12  Ubiquitinilation  1  1  1  2  4  4  Metabolic  process   Alpha/beta hydrolase  0  0  1  1  2  2  AMP-binding enzyme  0  0  1  1  2  2  Carbohydrate metabolism  1  1  2  2 Figure 4. Overview of C. lupini expression and synthesis of candidate effectors. Venn diagram showing the overlap between candidate effector-encoding genes upregulated compared to liquid cultures in the RNAseq study and the candidate effector proteins revealed by this study (A). Heatmap of the five quantified proteins associated with candidate effectors function, regarding their abundance and the corresponding transcript expression profile. Hatching areas indicate that data are not available (B). Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red).

Carbohydrate Active Enzymes (CAZymes)
CAZymes play an important role in many carbohydrate metabolism processes, some of which are related to host pathogenicity as degrading or binding enzymes to the host plant cell wall. A total of 63 candidate CAZymes were detected from the DEG list, 14 from the QP list and 12 were shared by both analyses ( Figure 5A). These CAZymes belonged to six classes: 13 associated with auxiliary activities, 38 to glycosyl hydrolases including one associated with a starch binding domain (CBM20) and another associated with a polysaccharide deacetylase (CE4), seven to glycosyl transferases, one polysaccharide lyases, five to carbohydrate esterases and one to carbohydrate-binding modules ( Table 2). The analysis of C. lupini genome led to the identification of 798 CAZyme-encoding genes, indicating that 7.9% of total CAZyme-encoding genes were differentially expressed compared to control conditions. Among the 65 predicted CAZymes, 28 were secreted proteins (including six from the QP list) and two were small-secreted proteins. The DEGs associated with CAZymes were mostly upregulated in the late stage of the infection, including a set of 23 genes specifically upregulated 84 hpi. Among significantly downregulated genes, downregulation was also more intense after 72 hpi Figure 4. Overview of C. lupini expression and synthesis of candidate effectors. Venn diagram showing the overlap between candidate effector-encoding genes upregulated compared to liquid cultures in the RNAseq study and the candidate effector proteins revealed by this study (A). Heatmap of the five quantified proteins associated with candidate effectors function, regarding their abundance and the corresponding transcript expression profile. Hatching areas indicate that data are not available (B). Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red).    CAZymes play an important role in many carbohydrate metabolism processes, some of which are related to host pathogenicity as degrading or binding enzymes to the host plant cell wall. A total of 63 candidate CAZymes were detected from the DEG list, 14 from the QP list and 12 were shared by both analyses ( Figure 5A). These CAZymes belonged to six classes: 13 associated with auxiliary activities, 38 to glycosyl hydrolases including one associated with a starch binding domain (CBM20) and another associated with a polysaccharide deacetylase (CE4), seven to glycosyl transferases, one polysaccharide lyases, five to carbohydrate esterases and one to carbohydrate-binding modules ( Table 2). The analysis of C. lupini genome led to the identification of 798 CAZyme-encoding genes, indicating that 7.9% of total CAZyme-encoding genes were differentially expressed compared to control conditions. Among the 65 predicted CAZymes, 28 were secreted proteins (including six from the QP list) and two were small-secreted proteins. The DEGs associated with CAZymes were mostly upregulated in the late stage of the infection, including a set of 23 genes specifically upregulated 84 hpi. Among significantly downregulated genes, downregulation was also more intense after 72 hpi ( Figure 1D). Regarding CAZymes identified by mass spectrometry, all, except for CLUP02_07510 (which was detected later from 72 hpi to 96 hpi), were detected from 36 hpi to 96 hpi, and their abundance increased throughout the infection process ( Figure 5B). Identified CAZymes were mainly involved in degradation of fungal or plant cell wall, glycosylation of proteins, energy, storage and recovery or synthesis of fungal cell wall (Table 3). ( Figure 1D). Regarding CAZymes identified by mass spectrometry, all, except for CLUP02_07510 (which was detected later from 72 hpi to 96 hpi), were detected from 36 hpi to 96 hpi, and their abundance increased throughout the infection process ( Figure 5B). Identified CAZymes were mainly involved in degradation of fungal or plant cell wall, glycosylation of proteins, energy, storage and recovery or synthesis of fungal cell wall (Table 3). Hatching areas indicate that data are not available (B). Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red). Hatching areas indicate that data are not available (B). Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red).

Peptidases
Peptidases also contribute to the degradation of host proteins and are involved in pathogenicity. The analysis of C. lupini genome led to the identification of 504 peptidases and 25 peptidase inhibitor-encoding genes. During the interaction with lupin, 20 DEGs and 13 QPs, among which four were identified by both analyses, were predicted as peptidases ( Figure 6B). The DEGs associated with peptidase activity included one peptidase inhibitor as well as 11 serine peptidases, three cysteine peptidases, 11 metallopeptidases, two threonine peptidases and one ubiquitin protease, among which seven were predicted as secreted ( Table 2). The upregulated putative peptidases and the peptidase inhibitor were mostly upregulated after 72 hpi, except for CLUP02_12375, a putative cysteine peptidase mainly upregulated as early as 24 hpi compared to liquid cultures ( Figure 1D), and for CLUP02_18115, a putative metallopeptidase quantified as early as 36 hpi ( Figure 6C) ( Table 2). All proteins were detected as early as 36 hpi and their abundance increased to 96 hpi except for CLUP02_18115, a putative metallopeptidase family M3. The abundance of the latter protein decreased significantly between 36 and 60 hpi before increasing again from 72 to 96 hpi ( Figure 6A,C).

Peptidases
Peptidases also contribute to the degradation of host proteins and are involved in pathogenicity. The analysis of C. lupini genome led to the identification of 504 peptidases and 25 peptidase inhibitorencoding genes. During the interaction with lupin, 20 DEGs and 13 QPs, among which four were identified by both analyses, were predicted as peptidases ( Figure 6B). The DEGs associated with peptidase activity included one peptidase inhibitor as well as 11 serine peptidases, three cysteine peptidases, 11 metallopeptidases, two threonine peptidases and one ubiquitin protease, among which seven were predicted as secreted ( Table 2). The upregulated putative peptidases and the peptidase inhibitor were mostly upregulated after 72 hpi, except for CLUP02_12375, a putative cysteine peptidase mainly upregulated as early as 24 hpi compared to liquid cultures ( Figure 1D), and for CLUP02_18115, a putative metallopeptidase quantified as early as 36 hpi ( Figure 6C) ( Table  2). All proteins were detected as early as 36 hpi and their abundance increased to 96 hpi except for CLUP02_18115, a putative metallopeptidase family M3. The abundance of the latter protein decreased significantly between 36 and 60 hpi before increasing again from 72 to 96 hpi (Figure 6 A,C). Heatmap of the 13 candidate peptidases regarding their abundance and their corresponding transcript expression profile. Hatching areas indicate that data are not available (C). Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red). Asterisks between histograms indicate significant changes over time.

Transmembrane Transporters
During plant-pathogen interaction, transmembrane transporters are actively deployed for the uptake of carbohydrates, oligopeptides, ions and other molecules released by plant cell lysing. The analysis of C. lupini genome led to the identification of 884 transmembrane transporter-encoding genes. By transcriptomic and mass spectrometry analyses, respectively 41 (4.6% of total transmembrane transporters in the genome) and seven putative transmembrane transporters were identified, among which only one was detected by both analyses ( Figure 7A). The transcripts belonged to 16 transporter families, including 14 belonging to the Major Facilitator superfamily (MFS), seven to the ATP-binding Cassette (ABC) superfamily and five to the P-type ATPase (P-ATPase) superfamily (Table 2). From 11 to 31 DEGs encoding transmembrane transporters were Figure 6. Overview of C. lupini expression and synthesis of candidate peptidases. Boxplot showing the mean abundance of the candidate peptidase from four biological repetitions and quantified by mass spectrometry (A). Venn diagram showing the overlap between peptidase-encoding genes upregulated in the RNAseq analysis and the peptidase proteins revealed by proteomic analysis (B). Heatmap of the 13 candidate peptidases regarding their abundance and their corresponding transcript expression profile. Hatching areas indicate that data are not available (C). Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red). Asterisks between histograms indicate significant changes over time.

Transmembrane Transporters
During plant-pathogen interaction, transmembrane transporters are actively deployed for the uptake of carbohydrates, oligopeptides, ions and other molecules released by plant cell lysing. The analysis of C. lupini genome led to the identification of 884 transmembrane transporter-encoding genes. By transcriptomic and mass spectrometry analyses, respectively 41 (4.6% of total transmembrane transporters in the genome) and seven putative transmembrane transporters were identified, among which only one was detected by both analyses ( Figure 7A). The transcripts belonged to 16 transporter families, including 14 belonging to the Major Facilitator superfamily (MFS), seven to the ATP-binding Cassette (ABC) superfamily and five to the P-type ATPase (P-ATPase) superfamily (Table 2). From 11 to 31 DEGs encoding transmembrane transporters were upregulated compared to liquid cultures from 60 hpi onwards while detected from 48 hpi. Downregulation of genes encoding transmembrane transporters was more intense in the late stages of the infection ( Figure 1D). All seven proteins were detected as early as 36 hpi and their abundance increased until 96 hpi ( Figure 7B).  Figure 1D). All seven proteins were detected as early as 36 hpi and their abundance increased until 96 hpi ( Figure 7B). Hatching areas indicate that data are not available (B). Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red).

Genes associated with Secondary Metabolism
The analysis of C. lupini genome led to the identification of 56 genes associated with secondary metabolism. Within the DEG and QP lists, only two secondary metabolism associated genes were predicted, being 3.6% of total genes associated with secondary metabolism in the genome ( Figure  8B). CLUP02_01848 and CLUP02_05705 were characterized by multiple conserved domains such as polyketide synthase dehydratase, condensation domain, phosphopantetheine attachment site and ketoreductase domain (KR; Table 2) and were both identified as hybrid (NRPS/PKS1). CLUP02_05705 shared 44.3% identity with a Beauveria bassiana Tenellin synthetase (e-value = 0), a protein involved in the pathogenicity of this fungus against insects [67] and 37.7% identity with an Aspergillus flavus hybrid PKS-NRPS synthetase (e-value = 0) that mediates biosynthesis of leporins, a precursor of toxins such as aflatoxins [68]. Likewise, CLUP02_01848 shared 39.7% of identity with Aspergillus clavatus polyketide synthase-nonribosomal peptide synthetase (e-value = 0) involved in mycotoxin biosynthesis [69]. Neither of these two genes was specific to any infection stage, but the maximum level of expression was reached at the end of the kinetics ( Table 2). The same proteins were detected from 36 to 96 hpi, but no significant increase was noticed over time ( Figure 8A). Hatching areas indicate that data are not available (B). Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red).

Genes Associated with Secondary Metabolism
The analysis of C. lupini genome led to the identification of 56 genes associated with secondary metabolism. Within the DEG and QP lists, only two secondary metabolism associated genes were predicted, being 3.6% of total genes associated with secondary metabolism in the genome ( Figure 8B). CLUP02_01848 and CLUP02_05705 were characterized by multiple conserved domains such as polyketide synthase dehydratase, condensation domain, phosphopantetheine attachment site and ketoreductase domain (KR; Table 2) and were both identified as hybrid (NRPS/PKS1). CLUP02_05705 shared 44.3% identity with a Beauveria bassiana Tenellin synthetase (e-value = 0), a protein involved in the pathogenicity of this fungus against insects [67] and 37.7% identity with an Aspergillus flavus hybrid PKS-NRPS synthetase (e-value = 0) that mediates biosynthesis of leporins, a precursor of toxins such as aflatoxins [68]. Likewise, CLUP02_01848 shared 39.7% of identity with Aspergillus clavatus polyketide synthase-nonribosomal peptide synthetase (e-value = 0) involved in mycotoxin biosynthesis [69]. Neither of these two genes was specific to any infection stage, but the maximum level of expression was reached at the end of the kinetics ( Table 2). The same proteins were detected from 36 to 96 hpi, but no significant increase was noticed over time ( Figure 8A). . Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red). Asterisks between histograms indicate significant changes over time.

Transcription Factors
Based on the analysis of conserved domains, 489 genes encoding transcription factors were predicted in the entire genome. A total of 16 DEGs (3.3% of total genes encoding transcription factors) and one QP were identified as transcription factors distributed in eight families (Table 2). Four putative transcription factors contained conserved domains including a fungal specific transcription factor with a Zinc finger C2H2 type or a fungal Zn(2)-Cys(6) binuclear cluster domain. One was predicted as secreted (CLUP02_11984). Most DEGs encoding candidate transcription factors were upregulated compared to liquid culture starting at 60 hpi, except two genes which were mainly overexpressed at 24 hpi followed by a significant decrease of expression between 60 and 72 hpi (CLUP02_01379 and CLUP02_02013; Figure 1D). The only protein identified as a transcription factor (CLUP02_17433) belonged to the TFIID complex and was detected from 36 to 96 hpi and its abundance increased after 60 hpi ( Figure 8C).

Cytochromes P450
Within C. lupini genome, 228 genes contained a cytochrome P450 conserved domain, and 11 were identified in the DEG list suggesting that 4.8% of the 228 genes were differentially expressed during the first 84 hpi in planta (Table 2). Putative functions were attributed to four of them, being L-ascorbate oxidase, ent-kaurene oxidase, benzoate 4-monooxygenase and N-alkane-inducible cytochrome P450. The 11 genes were mainly upregulated between 72 and 84 hpi, but one gene (CLUP02_02376) was significantly downregulated between 72 and 84 hpi. No cytochrome P450 was detected by mass spectrometry analysis. Hatching areas indicate that data are not available (B). Boxplots showing the mean abundance of the transcription factor from four biological repetitions quantified by mass spectrometry. (C). Colored bars indicated gene expression by Z-score calculated from normalized reads, and ranging from −2 (downregulated gene, blue) to 2 (upregulated gene, red), or protein abundance by Z-score calculated from protein log10 abundance, and ranging from −2 (lowest abundance, blue) to 2 (highest abundance, red). Asterisks between histograms indicate significant changes over time.

Transcription Factors
Based on the analysis of conserved domains, 489 genes encoding transcription factors were predicted in the entire genome. A total of 16 DEGs (3.3% of total genes encoding transcription factors) and one QP were identified as transcription factors distributed in eight families (Table 2). Four putative transcription factors contained conserved domains including a fungal specific transcription factor with a Zinc finger C2H2 type or a fungal Zn(2)-Cys(6) binuclear cluster domain. One was predicted as secreted (CLUP02_11984). Most DEGs encoding candidate transcription factors were upregulated compared to liquid culture starting at 60 hpi, except two genes which were mainly over-expressed at 24 hpi followed by a significant decrease of expression between 60 and 72 hpi (CLUP02_01379 and CLUP02_02013; Figure 1D). The only protein identified as a transcription factor (CLUP02_17433) belonged to the TFIID complex and was detected from 36 to 96 hpi and its abundance increased after 60 hpi ( Figure 8C).

Cytochromes P450
Within C. lupini genome, 228 genes contained a cytochrome P450 conserved domain, and 11 were identified in the DEG list suggesting that 4.8% of the 228 genes were differentially expressed during the first 84 hpi in planta (Table 2). Putative functions were attributed to four of them, being L-ascorbate oxidase, ent-kaurene oxidase, benzoate 4-monooxygenase and N-alkane-inducible cytochrome P450. The 11 genes were mainly upregulated between 72 and 84 hpi, but one gene (CLUP02_02376) was significantly downregulated between 72 and 84 hpi. No cytochrome P450 was detected by mass spectrometry analysis.

Factors Associated with Pathogenicity and Stress Response
Additional factors were detected by both RNA sequencing and mass spectrometry analyses and identified as involved in pathogenicity but did not belong to the previous functional class described earlier, including putative Necrosis and ethylene-inducing proteins (NEP), Nudix and stress responses proteins. NEPs, poorly studied in Colletotrichum spp., are known to elicit plant defenses, especially triggering hypersensitive response in dicots [70]. In C. lupini genome, 18 NEPs were identified of which six were upregulated during the infection and of these, two were detected by mass spectrometry (CLUP02_00595 and CLUP02_08026). They were all secreted and two of them (CLUP02_00596 and CLUP02_06520) were small-secreted peptides. The upregulation of NEPs was specific to medium infection stages since the maximum number of upregulated genes encoding NEPs was reached at 60 hpi (Table 2). Synthetized NEPs were detected as early as 36 hpi and their abundance increased until 96 hpi ( Figure 9). The six NEPs we identified showed mutations in the pattern "GHRHDWE" especially "D" and "E" which are crucial for the NEP activity suggesting that the NEP activity of these proteins was abolished [14] ( Figure 9B). CLUP02_13373 was identified as a putative intracellular hyphae protein. It contained a LysM domain probably involved in host interaction for establishment and maintenance of biotrophy, prevention of host recognition of the fungus and a barrier to host defense molecules [14,71,72].
Nudix were previously described as proteins belonging to effectors and potential markers of the end of the biotrophic phase and the start of the necrotrophic phase for hemibiotrophic pathogens [56]. In the entire genome of C. lupini, 20 genes were predicted as Nudix. Two genes encoding Nudix proteins were upregulated between 72 and 84 hpi (Table 2), among which one was also detected by mass spectrometry (CLUP02_07504) at 36 hpi and the highest amount of this protein was reached 96 hpi ( Figure 9A). CLUP02_07504 was predicted as a secreted protein. Among fungal stress response proteins, 16 genes were detected in the DEG list and 13 stress response proteins in the QP list, four were detected during infection by both analyses, all predicted as heat shock protein. The abundance of most stress response QPs was significantly higher at the end of the kinetics (i.e., 11/13 with significantly higher amounts at 96 hpi; Figure 9A). Then, other factors were predicted with function associated with pathogenicity, like protein with CFEM domain, Alternaria alternata allergen, cerato platanin. Among them, nine genes encoding proteins related to pathogenicity were upregulated and whose one was detected by mass spectrometry (CLUP02_08708) ( Figure 9A). Respectively, two, one and two genes were significantly upregulated between 24 and 48 hpi, between 48 and 60 hpi and between 72 and 84 hpi. The only protein detected was quantified from 72 to 96 hpi and its abundance significantly increased between 72 and 84 hpi ( Figure 9A).

Other Functions
Among the 593 upregulated DEGs, 213 were classified in the families described above, i.e., CAZymes, candidate effectors, transcription factors, secondary metabolism enzymes, transporters, peptidases, NEPs, Nudix, stress response, pathogenicity related and cytochrome P450. Most of the remainder genes (230) were upregulated between 60 and 84 hpi while a smaller part (40) was upregulated earlier between 24 and 48 hpi. These early-expressed genes were mainly involved in binding, cell structure, transcription, translation and oxidoreduction activities. Furthermore, 144 putative hypothetical proteins were identified, among which 17 were predicted as secreted proteins and seven were small-secreted proteins. Most hypothetical proteins were also upregulated between 60 and 84 hpi, except 26 which were mainly expressed between 24 and 48 hpi ( Figure 1D and Table 2). In the QP list, 215 proteins were identified with other putative functions especially proteins related to oxidoreduction activities (42 proteins), translation (48 proteins) mainly represented by ribosomal proteins and binding (28 proteins). Thirteen hypothetical proteins were found, among which eight were assigned as hypothetical by blast and 22 resulted from contradictory or incompatible predictions of functions depending on protein function assignment tools.

Discussion
C. lupini is a major threat for lupin crops and control solutions need to be found in order to use this leguminous crop in rotations as a promising alternative for protein sources. A few studies have been performed to identify its taxonomy through phylogenetic analysis [74][75][76][77][78], describe its morphology [75,77,78] and symptoms in the field [4] or evaluate the effect of temperature on its in vitro and in planta growth [78][79][80]. In this study, Illumina RNA sequencing and nLC mass spectrometry were used to identify candidate genes and proteins associated with C. lupini pathogenicity over the course of infection in lupin plants. Analyses were performed on living lupins inoculated with the pathogen and associated with symptom evaluation. We chose to work on plants inoculated at a pre-seedling stage, and not on detached organs, since detachment may also represent stressful conditions for the plant, ultimately altering fungal responses during host colonization. In addition, the inoculation was performed on the radicle emerging from the seed in order to mimic the primary inoculum source in fields where infected seeds are planted. The deciphering of the first steps of the infection by other Colletotrichum species has already been performed, e.g., for C. higginsianum on Arabidopsis [11], C. falcatum on sugarcane [13], C. fruticola on strawberry [12] or C. lentis on lentil [56]. Among these, both shared and species-specific molecular patterns of infection were reported. In this study, 897 DEGs were found in planta at 24, 48, 60, 72 and 84 hpi compared to C. lupini liquid cultures (filter of LFC > 1 and LFC < −1) while 304 proteins could be quantified in planta at 36, 60, 72, 84 and 96 hpi. Candidate functions were associated with some of these genes such as CAZymes, transmembrane transporters, peptidases, candidate effectors, transcription factors, secondary metabolism and cytochrome P450 which have already been associated with pathogenicity of Colletotrichum spp. in past transcriptomic analyses. For the first time, Nudix and NEPs known as pathogenicity markers were also reported here as one of the molecular mechanisms deployed during the infection, by both transcriptomic and proteomic analyses. Taken together, both analyses evidenced that only 93 genes and their corresponding proteins were differentially expressed during lupin infection. This low number of genes could be explained by the difference in kinetics applied for these analyses, i.e., protein analysis was performed 12 h after RNAseq analysis, and another limit is that the sampling for protein quantification ended at 96 hpi, when their abundance was still increasing. Futhermore, for both approaches, the drastic selection applied on data, and the different pipelines used may partially explain the low number of shared candidates, indicating also the reliability of this list of genes that may be involved in the pathogenicity of C. lupini. Below are described the putative molecular mechanisms used by C. lupini during the different stages of lupin infection based on transcriptomic and proteomic analyses.
During the early stage of infection, Colletotrichum spp. are known to develop an appressorium which enables to penetrate host cuticle and epidermal cells [81], the formation of appressoria is induced by molecular markers. In C. orbiculare, a signal transduction pathway led by cAMP dependent protein kinases induced spore germination, appressorium development and infection hyphae formation [81]. Similarly, putative transcripts encoding cAMP-dependent kinases were also identified as candidates for activation of melanin biosynthesis in C. gloeosporioides [9]. In this study, after 24 hpi, appressoria and then primary hyphae of C. lupini were observed by light microscopy. In addition, several DEGs associated with appressorium development were significantly upregulated at the early stages or downregulated at the end of the infection. Within our DEG list, a gene encoding a putative cAMP dependent protein kinase was identified with 88.22% similarity to a cAMP dependant candidate gene of C. higginsianum involved in appressorium production (e-value = 0) [11,19]. This gene was upregulated from 24 hpi to 60 hpi. In addition, a putative CAP protein, which is known to be involved in appressorium development, was identified as perilipin homolog [82] and was significantly downregulated at 72 and 84 hpi [82]. In addition, CLUP02_15638 identified as a putative CAS1 appressorium-specific protein (81% similarity to C. higginsianum; e-value = 2.0 × 10 −126 [83]) was also found to be downregulated from 60 hpi to 84 hpi. Early stages are also associated with the production of adhesion proteins and hydrophobins which are extracellular membrane proteins necessary for the interaction with the environment and as such, play an important role in pathogenicity. They are also shown to act as a developmental sensor of hydrophobic surface and to be involved in appressorium formation, attachment of hyphae to hydrophobic surfaces and to influence hyphal wall composition (reviewed in [84]). During the infection of lupin by C. lupini, a putative hydrophobin was significantly downregulated from 60 to 84 hpi. Likewise, CFEM-domain proteins are extracellular membrane proteins which were only described in fungi, with a higher occurrence in pathogenic fungi, and are suggested to play an important role in promoting pathogenicity although the exact function is still unknown [85]. For instance, CFEM-domain containing proteins from Magnaporthe grisea, namely adenylate cyclase-interacting protein ACI1 and PTH11 proteins, were shown to be involved in appressorium development [86,87]. In C. lupini, two genes encoding a protein carrying a CFEM domain were uniquely upregulated either at 60 hpi or 84 hpi, with the latter sharing 79% identity with a putative adhesin Mad1 of Colletotrichum spp. (CLUP02_08454; e-value = 0). Taken together, these results also suggest that appressorium formation may still be at stake later during infection, notably for spores which may germinate later. While transcriptomic analysis provided clear evidence that the genes involved in appressorium formation were upregulated early during infection onset, we did not detect any protein to support this observation. Interestingly, at 24 hpi, a set of 27 DEGs encoding hypothetical proteins was upregulated, evidencing the need to further explore the functions of these genes which are likely to be necessary for the onset of infection. Among them, we found a putative ribosomal protein SL 7 (CLUP02_04935) and a putative conserved hypothetical protein (PF11309, CLUP02_02353).
At 48 and 60 hpi, the observation of secondary hyphae and the upregulation of genes encoding CAZymes, effectors and pathogenicity-related proteins in C. lupini indicates the beginning of host penetration and the start of the setup of the necrotrophic phase. This resulted in a general increase of abundance of the corresponding proteins until 96 hpi, as observed by proteomic dynamics. At these stages, the fungus was supposedly penetrating host cells while almost no necrosis symptoms were visible. Penetration inside host cells is supported by the expression of CAZymes that degrade cell walls. The detection of the corresponding genes began at 48 hpi (9 genes) with mainly glycosyl hydrolase families GH1, GH18, GH28, and GH13 and carbohydrate esterases including pectinesterases while CAZymes were generally detected and quantified as early as 36 hpi. In C. fructicola, genes encoding CAZymes were upregulated as early as 24 hpi [12]. In C. lupini, the number of upregulated CAZymes continued to increase at 60 hpi with a total number of 23 genes out of the 63 identified in the DEG list. Most of them belonged to the families GH28, GH30, GH31, GH37 and GH131, but also AA9 (formerly GH61), CE like pectinesterases and pectin lyases. Except for the upregulated putative chitin synthase (CLUP02_17349), all other CAZymes were associated with plant cell wall degradation. One CAZyme was identified as a putative glucose-methanol-choline (GMC) oxidoreductase, described previously in Glomerella cingulata, C. higginsianum and C. graminicola and other phytopathogenic ascomycetes as reductases of plant-produced anti-fungal quinones and phenoxy-radicals [88,89]. Early stages of plant penetration are also characterized by the induction of genes encoding effectors. For instance, in C. higginsianum, candidate effectors were mostly upregulated during the biotrophic phase at 40 hpi [11]. In C. lupini, we found only one DEG encoding a candidate effector (CLUP02_08501) upregulated at 24 hpi versus 11 at 48 hpi and 19 at 60 hpi. Among the genes encoding candidate effectors, some contained conserved domains associated with pathogenicity functions such as a lysin-motif domain (IPR018392/PF01476; CLUP02_16647) described as involved in the sequestration of chitin oligosaccharides to avoid recognition by the host and plant defense induction [90], or a fungal hydrophobin domain (PF06766) generally found on the outer surface of conidia and hyphal wall, and involved in mediating contact and communication between the fungus and its environment [91]. Among pathogenicity-related genes, 3 DEGs encoding NEPs were detected here and were expressed as early as 48 hpi and 6 were upregulated from 60 hpi onwards. This result is further supported by the presence of NEP proteins that were found to be increasingly synthesized from 36 hpi onwards. The switch between the biotrophic and necrotrophic phase of C. higginsianum in A. thaliana appeared between 40 hpi and 60 hpi [11]. Likewise, C. fructicola developed necrotrophic secondary hyphae on strawberry before 72 hpi [12]. Like NEPs, Nudix are considered as other putative markers of the biotrophy-necrotrophy switch and may induce hypersensitive responses of the host plant. This protein, notably described in C. lentis (previously identified as C. truncatum), contains a signal peptide and a Nudix hydrolase domain that may be unique to hemibiotrophic fungal plant pathogens [56]. Of the two Nudix-encoding genes reported in this study (CLUP02_07504), one was found to be upregulated at 60 hpi while the other was upregulated at 84 hpi, surprisingly the candidate Nudix protein encoded by CLUP02_07504 was detected as early as 36 hpi but was significantly produced in higher amounts 96 hpi. Taken together, these results may indicate a transition from biotrophy to necrotrophy around 48 hpi and are consistent with the observation of primary hyphae at that time.
Transmembrane transporters were also highly expressed at 60 hpi among which MFS, APC and ABC transporters were mainly represented while these proteins were detected between 72 and 96 hpi based on proteomic analysis.
Overall, the setup of the necrotrophic phase seemed to occur around 60 hpi, based on our gene expression analysis. This was further corroborated by proteomic analysis, and also supported by observations on infected plants where necrotic lesions were almost absent at that time (only three hypocotyls out of 120 were found to display a 1 mm necrotic lesion) before clearly starting to appear at 72 hpi and keeping increasing at 84 hpi. Based on symptom observations, the necrotrophic phase was in progress between 72 hpi and 84 hpi. Among the DEGs, 54 genes were uniquely upregulated at 72 hpi and 168 at 84 hpi. These stages of infection, and even more at 84 hpi, were characterized by a general increase of the expression of genes encoding almost all families associated with the infectious process such as CAZymes, especially glycosyl hydrolases, peptidases, and transmembrane transporters, including sugar transporters (MFS), together with the general increase of the amount of corresponding QPs reaching a maximum 96 hpi. In contrast, the expression of most of the DEGs encoding candidate effectors remained stable or decreased significantly from 48 hpi onwards, except for 2 genes which expression was significantly increased from 72 and 84 hpi (CLUP02_01667 and CLUP02_12726). The former was identified as a putative ToxB-encoding gene, a toxin also produced by the rice blast fungus Pyrenophora tritici-repentis [92]. The latter gene contained two previously-described domains, the IPR015131/PF09044 domain identified as a toxin killer and the IPR029167/PF15474 involved in meiosis mechanism. These domains were already described among the predicted effectors of Fusarium graminearum [93,94]. Concerning CAZymes, 38 and 59 genes out of 63 were respectively upregulated at 72 and 84 hpi, mostly represented by glycosyl hydrolases and cellobiose dehydrogenases involved in plant cell degradation. Peptidases were also mainly expressed during the necrotrophic phase although detected in the transcriptomic analysis as early as 24 hpi. Likewise, proteomic analysis revealed that the highest levels of detected candidate CAZymes and proteases were reached 96 hpi. At the same time, genes encoding transmembrane transporters can be upregulated as they are required for assimilating the products of the degradative activity, as observed in C. higginsianum [11,12]. In particular, ABC and MFS transporters are known in other species to detoxify and eliminate antifungal secondary metabolites and various toxins produced by plants [95,96]. This function, which was characterized in C. acutatum CaABC1 transporter, may contribute to resistance against some fungicides, and was also necessary for conidiation and abiotic stress resistance [23], whereas some MFS transporters like Mfs1 in C. lindemuthianum are induced during the necrotrophic phase for mediating hexose uptake [97]. At this stage, genes encoding candidate enzymes involved in the response to oxidative stress were also upregulated, such as one catalase (CLUP02_03311) and two peroxidases (CLUP02_09271, CLUP02_09357). The expression of these genes in the late stages of the infection could indicate that C. lupini set up an alleviation of the oxidative stress, further contributing to detoxify the cell in response to the production of reactive oxygen species (ROS) such as hydrogen peroxide (H 2 O 2 ) by the host [98,99]. Likewise, in C. gloeosporioides f. sp. malvae, a catalase-encoding gene was highly expressed in the necrotrophic phase of an infection of round-leaved mallow, Malva pusilla [100]. Last, another gene encoding a cerato-platanin protein putatively involved in pathogenicity was uniquely upregulated at 84 hpi, previously described as both an effector and an elicitor protein abundant in filamentous fungi. This dual role was evidenced in C. truncatum infecting many leguminous species, where the cerato-platanin protein was detected in the biotrophic phase as an effector suppressing plant defense responses, and may also be involved as an elicitor of plant defense responses in the necrotrophic phase [101,102]. In C. falcatum, one cerato platanin EPL1 was also found as a secreted protein eliciting systemic resistance in sugarcane and HR response in tobacco [103]. Likewise, in our study, the cerato-platanin-encoding gene (CLUP02_11344) seemed to be associated with the necrotrophic stage only.
This study proposes that the evaluation of the pathogenicity associated with a transcript and the proteome analysis approach can provide a better understanding of the infectious process of C. lupini. Taken together, these results suggest that C. lupini may be a hemibiotrophic fungal pathogen with a switch between the biotrophic and necrotophic phase probably occurring around 48 hpi under our experimental conditions. The first stages of the infection began with an appressorium development during the first 24 h post-inoculation as well as an increase of expression of genes involved in cellular processes. Some genes encoding hypothetical proteins were uniquely upregulated at the early stage of infection and their future identification will probably contribute to better understand the infectious process of the biotrophic phase. During the biotrophic-necrotrophic transition, expression of genes encoding proteins like NEP, Nudix, CAZymes, peptidases, transmembrane transporters and candidate effectors increased suddenly, this increase was consistent with the quantification of proteins from the same family and was associated with the apparition of necrotic lesions on infected plants. A few toxin-encoding genes were also upregulated in the necrotrophic phase. Overall, transcriptomic and proteomic analyses of C. lupini infecting lupins shed a new light on the putative mechanisms involved in the pathogenicity of Colletotrichum spp. over the course of the infection. However, the findings presented in this study need to be confirmed by functional experiments to validate the roles played in the infectious process.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-2607/8/10/1621/s1, Table S1: RNAseq read counts and percentage mapping statistics to C. lupini (RB221 strain) genome. Table S2: Log10 abundance of the 368 proteins. Abundance of proteins written in red have a p-value > 0.01 (64 proteins) and those written in green have a p-value < 0.01 (304 proteins). R1-4 corresponds to the number of the replicate. Table S3: Expression analysis for all C. lupini DEGs in comparison with liquid culture and measured during the infectious process using RNA-Seq. LC: Liquid Culture; IP: In planta; the number after LC or IP corresponds to time post-inoculation in hour.