Erratum: Survival trade-offs in plant roots during colonization by closely related beneficial and pathogenic fungi

Nature Communications 7 Article number:11362 (2016); Published 6 May 2016; Updated 29 September 2016. The HTML version of this Article previously published omitted Supplementary Data 8. This has now been corrected in the HTML version of the Article; the PDF was correct from the time of publication.

F ungal endophytes are a ubiquitous and phylogenetically diverse group of organisms that establish stable associations with living plants, but in most cases their ecophysiological significance is poorly understood 1 . Species of the fungal genus Colletotrichum are best known as destructive pathogens on 43,000 species of dicot and monocot plants worldwide, causing anthracnose diseases and blights on leaves, stems, flowers and fruits 2 . However Colletotrichum species can also grow benignly as endophytes on symptomless plants 3 , and although only few pathogenic members of the genus attack plant roots 4 , Colletotrichum endophytes are frequently isolated from the roots of healthy plants 5,6 . Moreover, although the genome sequences and in planta transcriptomes were recently described for four species pathogenic on above-ground plant parts 2,7 , such information is not available for any root-associated Colletotrichum pathogens or endophytes.
We found recently that C. tofieldiae (Ct) is an endophyte in natural populations of Arabidopsis thaliana growing in central Spain 8 . The fungus initially penetrates the rhizoderm by means of undifferentiated hyphae, which then ramify through the root cortex both inter-and intracellularly, occasionally spreading systemically into shoots via the root central cylinder without causing visible symptoms. Under phosphate-deficient conditions (50 mM KH 2 PO 4 ), colonization by Ct promoted plant growth and fertility and mediated the translocation of phosphate into shoots, as shown by 33 P radiotracer experiments 8 . However, neither the plant growth promotion nor phosphate translocation activities were detectable under phosphate-sufficient conditions (625 mM KH 2 PO 4 ), indicating that plant fitness benefits conferred by Ct are strictly regulated by phosphate availability. In striking contrast, colonization of A. thaliana roots by the closely related pathogenic species C. incanum (Ci), which attacks members of the Brassicaceae, Fabaceae and Solanaceae, severely inhibited Arabidopsis growth and mediated only low levels of 33 P translocation into shoots 8 . These findings raise the possibility that in low-phosphate soils, root colonization by the Ct endophyte compensates for the absence of key genetic components required for mycorrhizal symbiosis in the Brassicaceae lineage, which is otherwise conserved in B80-90% of terrestrial plants 9 .
In the present study, we report the genomes of five isolates of beneficial Ct and one isolate of pathogenic Ci, and analyse the transcriptomes of each species during their colonization of Arabidopsis roots under phosphate-deficient and phosphatesufficient conditions. Comparison of the two species allows us to identify fungal adaptations to the endophytic lifestyle at the level of both gene repertoire and gene regulation, and provides insights into the evolutionary transition from parasitism to endophytism within a single fungal genus. On the host side, transcriptional responses of Arabidopsis roots to colonization by beneficial Ct are modulated by the phosphate status, providing evidence that trade-offs between defense and nutrition control the outcome of the interaction between Arabidopsis and Ct. Our findings also shed light on the ability of plants to maximize survival by prioritizing their responses to simultaneous biotic and abiotic stresses.

Results
Genome sequencing and evolution of Ct and Ci lifestyles. We sequenced the genome of the plant growth-promoting fungus Ct isolate 0861, a root endophyte isolated from natural populations of A. thaliana in Spain 8,10 , and those of four other Ct isolates isolated from diverse dicot and monocot hosts in Europe (Supplementary Note 1). We also sequenced the broad hostrange pathogen Ci, isolated from radish (Raphanus sativus) leaves in Japan, that strongly impairs plant growth when inoculated onto Arabidopsis roots 8,11 (Supplementary Fig. 1, Supplementary  Table 1 and Supplementary Note 1). Illumina short reads were used to build high-quality genome assemblies of similar size for all isolates, ranging from 52.8 to 53.6 Mb (Supplementary Table 2 and Supplementary Note 2). Molecular phylogeny, whole-genome alignment and divergence date estimates indicate that Ct and Ci are closely related taxa within the Colletotrichum spaethianum species complex and diverged only B8.8 million years ago (Fig. 1a, Supplementary Figs 2 Table 3 and Supplementary Note 3). Our phylogenetic analysis suggests that evolution from pathogenic ancestors towards the beneficial endophytic lifestyle in Ct is a recent adaptation in Colletotrichum fungi.

and 3, Supplementary
SNP distribution and reproductive mode of Ct isolates. Although the five Ct isolates originate from widely separated geographical areas and distantly related plant hosts, they diverged only B0.29 million years ago and the aligned fractions (493%) of their genomes share 499% sequence identity ( Fig. 1a and Supplementary Tables 1,3 and 4). The overall frequency of single-nucleotide polymorphisms (SNPs) between isolates was similar (2.22-3.04 SNPs per kb) but the SNP distribution within each genome was uneven, with alternating tracts of low (0.22-0.32 SNPs per kb) and high (4.25-5.12 SNPs per kb) SNP density (Fig. 1b, Supplementary Fig. 4 and Supplementary Table 5). This peculiar SNP distribution, also visible in the genomes of other plant-interacting fungi 12,13 , is consistent with chromosome recombination events. However, the SNP density profiles are remarkably similar between isolates and large haplotype blocks are conserved between all (21%), four (19%), three (18%) or two (17%) of them, with only 22% being isolate specific (Fig. 1b,c, Supplementary Fig. 4, Supplementary Table 6 and Supplementary Note 4). These conserved SNP signatures in the genomes of geographically distant isolates were likely generated by rare or ancestral sexual/parasexual reproduction and maintained by frequent clonal propagation.
Evolutionary dynamics of multigene families in Colletotrichum. Similar numbers of protein-coding genes were predicted in Ct0861 and Ci (B13,000; Supplementary Table 2), with 411,300 orthologous genes shared between both species. By clustering protein-coding sequences into sets of orthologous genes using OrthoMCL, we identified 7,297 gene families conserved across all six analysed Colletotrichum species and 10,519 shared between Ct0861 and Ci (Fig. 2a,b and Supplementary Note 5). Using a maximum-likelihood approach, we also reconstructed ancestral genomes for each Colletotrichum lineage and predicted the number of gene families that were likely gained or lost in each species compared with its corresponding ancestor (Supplementary Fig. 5 and Supplementary Note 6). We found significantly more gene families gained (1,009) than lost (198) on the branch leading to Ct compared with other branches of the tree (Fisher's exact test, P ¼ 3.98 Â 10 À 136 ; Supplementary Fig. 5 and Supplementary Data 1). Functional enrichment analysis among the 1,009 gene families gained ( Supplementary Fig. 5) and the 1,486 Ct-specific gene families (Fig. 2b) revealed a significant enrichment for genes encoding secondary metabolite biosynthesis-related proteins in Ct (Fisher's exact test, P ¼ 5.89 Â 10 À 3 and 3.31 Â 10 À 8 , respectively). This result contrasts with the very low number of secondary metaboliterelated genes detected in the genomes of other root-associated fungal endophytes and mycorrhizal fungi 14   isolates with respect to the Ct0861 reference assembly. Tracks represent (from the outside) the five largest Ct0861 contigs (scale: kb); locations of predicted genes; locations of SNPs versus Ct0861 in CBS495, CBS130, CBS127, CBS168 (see Supplementary Table 1 for full culture IDs) and SNPs common to these four isolates; conserved regions with low SNP density between all the five isolates; mean read coverage (per 100 bases) for isolates CBS495, CBS130, CBS127 and CBS168. Coverage plot scales are 0 to 1,000 (CBS495) or 0 to 500 (CBS130, 127, 168). (c) SNP density (per 1 kb) in isolates CBS495, CBS130, CBS127 and CBS168 versus Ct0861, compared with gene density (per 10 kb) and GC content (%) on the five largest Ct0861 contigs.  Note 8). By comparing the Ct gene repertoires to those of five other plant-associated fungal endophytes from both ascomycete and basidiomycete lineages, we found no obvious common genomic signatures to indicate the convergent evolution of an endophyte 'toolkit' . Furthermore, the convergent loss of decay mechanisms characteristic of ectomycorrhizal fungi 15,16 is not a hallmark shared by the non-mycorrhizal root endophytes ( Supplementary  Fig. 12), suggesting that these fungi have followed different evolutionary trajectories to acquire the ability for intimate growth in living root tissues 14,17 .
Despite the overall similar secretome size of all analysed Colletotrichum species (13.3-15.9% of the total proteome), the proportion of genes encoding candidate secreted effector proteins (CSEPs), which may promote fungal infection 18 , varied considerably between species (6.6-15.8% of the total secretome; Fig. 3a Fig. 3b). Genomes from additional Ci isolates are now needed to determine whether there is differential host-selective pressure on the CSEP repertoires of endophytic Ct and pathogenic Ci that reflect their contrasting lifestyles. Similar to other Colletotrichum species 2 , CSEPs in Ct and Ci are not organized into large multigene families, possibly due to a low frequency of duplication events in their respective genomes (Fig. 3c,d and Supplementary Table 2).
Both Ct and Ci genomes encode a very broad range of CAZymes, including large arsenals of pectate lyases, carbohydrate esterases and glycoside hydrolases acting on all major plant cell wall constituents (Fig. 4a, Supplementary Fig. 12 and Supplementary Data 4). However, the number of predicted carbohydrate-binding modules is inflated in Ct compared with pathogenic Colletotrichum species, especially chitin-binding CBM18 (48 versus 28-40) and CBM50 (57 versus 30-54) modules (Fig. 4a, Supplementary Data 4), though few of the corresponding Ct genes were induced in planta ( Supplementary  Fig. 14). These two chitin-binding modules are similarly highly enriched in the genomes of two other non-mycorrhizal root symbionts 19,20 (Piriformospora indica and Harpophora oryzae; Supplementary Data 4), suggesting this is a genomic signature common to independently evolving root-associated fungal endophytes.
Dual RNAseq of Arabidopsis roots and fungal partners. We report elsewhere that Ct promotes Arabidopsis growth under phosphate-deficient ( À P) but not phosphate-sufficient ( þ P) conditions and that transfer of radioactive 33 P from Ct hyphae to host plants is strictly regulated by Pi (inorganic phosphate) availability 8 . To compare the transcriptional dynamics of beneficial Ct and pathogenic Ci during colonization of Arabidopsis roots and study the corresponding host responses, we extensively re-analysed the previously created RNA-seq data for the Ct-Arabidopsis interaction    in planta ( Table 10). We also found major differences between Ct and Ci in the expression of gene categories typically associated with fungal pathogenicity. In planta activation of CSEPs was more pronounced in Ci compared with Ct, with seven times more CSEPs highly expressed (top 1,000 expressed genes) and three times more upregulated in planta at 10 d.p.i. . Likewise genes encoding CAZymes and secondary metabolism enzymes displayed earlier and stronger transcriptional activation in planta and broader diversity in Ci (Fig. 4b,c and Supplementary Fig. 22). Consistent with this, we observed a reduced number of living cells and a depletion of beta-linked polysaccharides (including cellulose) from host cell walls in Ci-colonized roots at 10 d.p.i., but not in Ct-colonized roots ( Supplementary Fig. 1  of Ct might contribute to, or be a consequence of, the beneficial relationship. Overall, our results suggest that the recent transition from pathogenic to beneficial lifestyles might be partly controlled through transcriptional downregulation of pathogenicity-related genes in Ct. Host responses to Ct are phosphate-status dependent. To disentangle how Pi-starved and non-starved Arabidopsis roots respond to Ct colonization over time, we compared Ct-colonized and mock-inoculated roots under þ P and À P conditions. In total, 5,661 Arabidopsis genes were differentially expressed in at least one of the 16 pair-wise comparisons (moderated t-test, |log 2 FC|Z1, FDRo0.05) and grouped into 20 major gene expression clusters (Fig. 5a and Supplementary Data 8). GO term enrichment analysis among these clusters indicated that the phosphate level used in our study (50 mM) was sufficient to provoke a phosphate starvation response in Arabidopsis roots (clusters 2 and 4; Fig. 5b). Furthermore, our analysis indicates that 'response to stimulus', 'indole glucosinolate metabolic process', 'defense response' and 'ethylene metabolic process' are activated in Ct-colonized roots under þ P but not À P conditions (cluster 9) (Fig. 5b and Supplementary Data 9). In contrast, the genes related to 'root cell differentiation' (cluster 8, Fig. 5b) and phosphate uptake 8 were preferentially activated in Pi-starved Arabidopsis roots during Ct colonization, similar to mycorrhizal symbiont-host interactions 21 . To identify key regulatory genes (hub genes) that might orchestrate transcriptional reprogramming in the contrasting directions seen in clusters 8 and 9, we checked which of these genes are often co-regulated in other expression data sets using the ATTED-II gene coexpression database (Fig. 5c). Among the hub genes that showed high connectivity within cluster 8 (highlighted with black dots), many encode proteins involved in cell wall remodelling and root hair development. Particularly, genes encoding the root hair-specific proteins RHS8, RHS12, RHS13, RHS15 and RHS19 (ref. 22) are upregulated (moderated t-test, |log 2 FC|Z1, FDRo0.05) in Ct-colonized versus mock-treated roots under À P conditions, which was validated by RT-qPCR ( Fig. 5d and Supplementary Fig. 18). This expression pattern suggests that Ct-dependent remodelling of root architecture might play a key role to enhance phosphate uptake during starvation (Supplementary Note 9). Similarly, we identified 27 hub genes within cluster 9 (Fig. 5c, black dots)  patterns by receptors of the innate immune system and is needed for broad-spectrum defence to restrict the growth of fungal pathogens 26,27 . Notably, in Arabidopsis mutants that cannot activate PEN2-mediated antifungal defense, the promotion of plant growth by Ct is impaired, while the depletion of all Trpderived secondary metabolites renders Ct a pathogen on Arabidopsis 8 . These findings strongly suggest that the phosphate starvation response and Trp-derived indole glucosinolate metabolism are interconnected to control fungal colonization of Arabidopsis roots 28 . Phosphate status-dependent activation of defense responses was also observed among the 411 expressed Arabidopsis genes annotated as 'chitin-responsive' (Supplementary Fig. 23), based on GO term enrichment among all significantly regulated genes ( Supplementary Fig. 24) and this was validated by RT-qPCR ( Fig. 5d and Supplementary Fig. 18). These data reveal a remarkable capacity of Arabidopsis roots to prioritize different transcriptional outputs in response to Ct, favouring either defense responses under þ P conditions or root growth and phosphate metabolism under À P conditions.

Phosphate-starved roots activate defense responses to Ci.
To clarify whether the reduced activation of defense responses observed in Ct-colonized roots under À P conditions is not simply due to phosphate deficiency, we compared the transcriptomes of Pi-starved Arabidopsis roots in response to either Ci or Ct at 10 d.p.i. In total, 2,009 differentially expressed genes were identified (moderated t-test, |log 2 FC|Z1, FDRo0.05), including 988 genes induced in Ct-colonized roots (cluster 1) and 1,021 genes in Ci-colonized roots (cluster 2; Fig. 6a and Supplementary Data 10). GO term enrichment analysis revealed that ion transport and root cell differentiation mechanisms were activated in Ct-colonized roots, whereas strong defense responses were triggered in Ci-colonized roots (Fig. 6b). Thus, although Pi-starved Arabidopsis roots remain able to mount immune responses against pathogenic Ci, transport and root growth are instead prioritized during interaction with beneficial Ct.

Discussion
Deciphering the genetic basis of the transition from pathogenic to beneficial plant-fungal interactions is crucial for a better understanding of the evolutionary history of fungal lifestyles 20,29 . It was recently shown that the ectomycorrhizal lifestyle arose independently multiple times during evolution and that the transition was associated with (1) convergent loss of genes encoding PCWDEs present in their saprotrophic ancestors and (2) the repeated evolution of lineage-specific 'toolkits' of mycorrhiza-induced genes 15 . However in striking contrast with ectomycorrhizal fungi, this transition in Ct, P. indica and H. oryzae was not accompanied by contraction of their PCWDE repertoires 19,20 . In our study, the close phylogenetic relatedness of beneficial Ct and pathogenic Ci, and their ability to infect the same plant host, allowed us to resolve both genomic and transcriptomic signatures associated with this evolutionary transition. The overall high genomic similarity between Ct and Ci suggests that this transition involved only subtle remodelling of the gene repertoire (that is, a reduced set of CSEPs and expansion of chitin-binding and secondary metabolism-related protein families). The retention of abundant pathogenicity-or saprotrophy-related genes implies that they are still needed by Ct, perhaps for exploitation of other plant hosts or during plant senescence when Arabidopsis leaves are extensively colonized by Ct mycelium 8 . Our results also suggest that changes in fungal  Overrepresented (yellow to red) and underrepresented transcripts (yellow to blue) are shown as log 2 (fold changes) relative to the mean expression across all stages. E1 and E2 correspond to two fully independent experiments (see Supplementary Note 9). Gene expression fold changes (green: downregulated; violet: upregulated) were calculated between C. tofieldiae-colonized versus mock-treated roots, C. incanum-colonized versus mock-treated roots or C. incanum-colonized versus C. tofieldiaecolonized roots. (b) GO term enrichment analysis of Arabidopsis genes preferentially expressed in response to C. tofieldiae (Cluster 1) or in response to C. incanum (cluster 2). Each circle corresponds to a significantly enriched GO term (Po0.05, hypergeometric test, Bonferroni step-down correction). The colour code reflects P values and the circle size the number of genes associated to each GO term. Similar to Fig. 5b, the GO terms that are tightly connected are functionally linked and therefore only the major host-response outputs are indicated (dotted line). REI, Relative Expression Index.
gene expression patterns during host colonization, rather than extensive remodelling of the gene repertoire, provides an alternative and probably transient adaptation to a beneficial endophytic lifestyle. This may reflect the relatively recent transition from pathogenic to non-pathogenic lifestyles in Ct and, consequently, a latent capacity to revert to a pathogenic lifestyle. During the last decade, the molecular mechanisms by which plants respond to colonization by pathogenic or mutualistic fungi have been extensively studied 30 . However, it remains unclear how plants discriminate and respond appropriately to closely related fungal partners with different lifestyles. The sedentary nature of plants suggests they have evolved regulatory systems to integrate exposure to conflicting biotic and abiotic stresses and balance their resource allocation strategically to maximize growth and survival. A recent report showed that plant responses to multiple stresses are not cumulative and suggested that prioritization of stress responses does take place 31 . For plant-mycorrhizal associations, an inverse correlation was observed between phosphate levels and the number of arbuscules formed in roots 32 . Although the detailed molecular mechanism remains unclear, this suggests that the nutritional status of the plant impacts fungal colonization efficiency. Here, we show that host transcriptional responses to Ct are dependent on phosphate availability, with defense responses activated or suppressed under high-or low-phosphate conditions, respectively. The fact that immune responses are retained in phosphate-starved roots colonized by Ci makes it unlikely that metabolic competition between phosphate starvation and defense response systems attenuates defense gene activation during interactions with Ct under P-limiting conditions. Recently, a metabolic link between the phosphate starvation response and glucosinolate biosynthesis was described 28 and the functional relevance of this link is supported by our observation that Ct-mediated plant growth promotion is impaired in Arabidopsis mutants lacking regulatory components of indole glucosinolate metabolism or the phosphate starvation response 8 . Therefore, we hypothesize that connectivity between nutrient sensing and innate immunity systems in the host, combined with subtle genomic adaptations in Ct, has enabled the transition from pathogenic to beneficial Arabidopsis-Colletotrichum interactions ( Supplementary  Fig.  25). Consequently, the interaction with beneficial Ct, but not with pathogenic Ci, is tightly controlled in plant roots by trade-offs between nutrition and defense. Whether phosphate stressdependent defense attenuation renders Ct-colonized plants super-susceptible to other microbial pathogens remains to be tested. Our results are consistent with the fact that transfer of Pi from ramifying fungal hyphae to roots, and subsequent allocation to shoots for plant growth, occurs only under phosphate-deficient conditions 8 . Notably, where Ct naturally associates with Arabidopsis in central Spain, the level of bioavailable phosphate in soil at those locations is very low (5.5 to 17 p.p.m., Supplementary Table 11). Our findings suggest that both innate immune responses (that is, indole glucosinolate metabolism) and soil phosphate availability are important selective forces driving fungal adaptation and contributing to the evolutionary transition from parasitic to beneficial Arabidopsis-fungal associations.

Methods
Genome sequencing and assembly. C. incanum and the five C. tofiediae isolates were grown in liquid Mathur's medium (2.8 g glucose, 1.22 g MgSO 4 .7H 2 O, 2.72 g KH 2 PO 4 and 2.18 g Oxoid mycological peptone in 1 l deionized water) supplemented with 100 mg ml À 1 rifampicin and 125 mg ml À 1 streptomycin. Genomic DNA was isolated using the DNeasy Plant Mini Kit (Qiagen) from 100 mg of fungal mycelium. Library construction, quality control and DNA sequencing for 454 GFLX þ or Illumina Hiseq sequencing were performed at the Max Planck Genome Centre Cologne (http://mpgc.mpipz.mpg.de) using 1 mg genomic DNA. After the preparation of genomic DNA libraries, 454 reads (557 bp on average) and Illumina paired-end reads (100 bp) were obtained from Roche 454 FLX þ and Illumina HiSeq2500 sequencers, respectively. For the Ct0861 reference genome, a hybrid assembly strategy was used combining 454 and Illumina data. Unpaired 454 reads were first assembled using MIRA 4.0 (ref. 33) and filtered MIRA-contigs (45,000 bp) were further used for scaffolding of Illumina paired read assemblies from SPAdes 3.0 (ref. 34). The established SPAdes 3.0 pipeline was used in 'careful' mode providing 454 MIRA assemblies as untrusted-contigs for scaffolding only and a kmer scan using 21, 31, 41, 61, 75 and 81. All other assemblies were constructed only from Illumina data using a combination of VELVET 1. 2.1 (ref. 35) and SPAdes 34 . Using BLASTN searches, contigs were identified that were missing from combined SPAdes assemblies but present in VELVET assemblies. To integrate those contigs and extend further where possible, SPAdes was re-run as described above but in 'trusted-contigs' mode where trusted contigs were provided as fasta files with absent contigs only. All the assemblies were generated using 'careful' mode in SPAdes to avoid miss-pairing of contigs by scaffolding and for further analyses, contigs o100 bp were removed. To identify and remove potential contaminating sequences, assemblies were aligned to the genomes of A. thaliana, H. sapiens and PhiX (sequencing spike-in control) using MUMmer 36 with default parameter settings. Contigs that aligned with more than 50% of their sequence (coverage; 'COV') and at least 85% sequence identity ('IDY') to any of the tested contaminants were removed from the assemblies. In addition, contigs that aligned with 75-85% identity (and 450% coverage) or with 10-50% coverage (and 485% identity) were also removed, if the judgment of the sequence being non-fungal was confirmed through BLASTN searches in the NCBI nr database (with default settings). For the Ct0861 assembly, RNA-sequencing data were used for further clean-up. Finally, assembly quality was assessed on the basis of L50/75/90 and N50/75/90 values, percentage of error-free bases estimated with REAPR 37 (version 1.0.16, default settings) and gene space coverage estimated with CEGMA 38 (version 2.0, default settings).
Repetitive DNA analysis. We identified repetitive DNA in the genome assemblies using either de novo or homology approaches. For de novo searches, we used PILER and PALS 39 to identify repetitive sequences and classify them into families. The resulting libraries of consensus sequences were then used to scan the genome sequences using RepeatMasker 40 (version 4.0.3) to identify individual repetitive elements. For homology-based searches, we used RepeatMasker using a library of all fungal elements in the Repbase database 41 (version 20140131).
Phylogeny and divergence date estimation. All phylogenetic analyses performed in this study are described in the Supplementary Note 11. For evolutionary divergence date estimation, clustering, protein family selection and phylogenetic analyses were performed with scripts in the Mirlo package (https://github.com/ mthon/mirlo). The phylogeny was calibrated using the penalized-likelihood method implemented in r8s (ref. 42) using one primary and two secondary calibration points (Supplementary Note 11).
Short-read alignment and SNP analysis. To compare the genome sequences of Ct isolates, Illumina short reads of the four other isolates were mapped onto the genome assembly of Ct0861 using Bowtie2 (ref. 43) (default settings for paired-end data). Subsequently, duplicate reads were removed using the rmdups function from the SAMtools toolkit 44 (default settings). On the basis of the mapped genome sequencing reads, single-nucleotide polymorphisms (SNPs) were identified using the mpileup function in SAMtools 44 (version 0.1.18; with option -u) The obtained SNP sets were filtered by applying the bcftools script vcfutils.pl varFilter (SAMtools) with adjusted read depth settings according to the respective sequencing read coverage to -d 80 and -D 800 for CBS495 and to -d 40 and -D 400 for CBS130, CBS127 and CBS168. The SNP locations, read coverage for each isolate and locations of conserved regions were visualized using the Circos software package 45 (version 0.62.1). In addition, we also calculated SNP densities (SNPs per kb) relative to Ct0861 for each isolate as a function of the genomic location on all Ct0861 contigs larger than 50 kb, using a 10-kb sliding window that moved 1 kb at each step. For visualization of the SNP densities, these windows were sorted in the increasing order by contig number and position on the contig. To identify windows with a low SNP density, that is, a common haplogroup, between isolates we classified the SNP density in each window as either 'low' or 'high' using a two-state hidden Markov model (HMM). This HMM was created and fitted on the observed 10 kb SNP densities by the expectation-maximization algorithm using functions 'depmix' and 'fit' (R package depmixS4), and subsequently the posterior state sequence (with states 'low' and 'high'), computed via the Viterbi algorithm, was extracted with function 'posterior' (R package depmixS4).
Gene annotation. The prediction of Ct and Ci gene models was performed using the MAKER pipeline 46 (version 2.28) , which integrates different ab initio gene prediction tools together with evidence from EST and protein alignments. In a first step, for each genome, the pipeline was run using Augustus 47 (with species model Fusarium graminearum) and GeneMark-ES 48 for ab initio gene prediction together with transcript and protein alignment evidence. The resulting gene models from this first run were used as training set for a third ab initio prediction tool, SNAP 49 , NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11362 ARTICLE and subsequently the annotation pipeline was re-run, this time including all three ab initio prediction tools together with the transcript and protein alignment evidence to yield the final gene models. The alignment evidence was created from BLAST and Exonerate 50 alignments of both protein and transcript sequences of each respective fungus (Ct/Ci) and protein sequences of C. higginsianum and C. graminicola. Ct (isolate 0861) and Ci transcript and protein sequences were obtained from the corresponding RNA-seq data via a transcriptome de novo assembly. For this purpose, we extracted all RNA-seq read pairs that did not align to the host plant genome from four (Ct) to nine (Ci) in planta samples and combined these with the read pairs from one in vitro sample of the respective fungus. The combined RNA-seq reads were then used as input for Trinity 51 (with default parameter settings for paired-end reads) to assemble transcripts and extract peptide sequences of the best-scoring ORFs (using the Perl script 'transcripts_to_best_scoring_ORFs.pl' provided with the Trinity software). General functional annotations for the predicted gene models were obtained using Blast2GO (ref. 52). To perform Blast2GO searches and ensure stable databases over time for multiple genome annotations, the NCBI nr database was downloaded locally (version: 8 January 2015). In addition, a local b2gdb mysql database was generated (version 201402) and connected to the Blast2GO java tool. For each genome annotation, BLASTP was performed against the local NCBI nr database ( À e 1E À 3, -v 10 À b 10) and tabular BLAST output was loaded into Blast2GO using graphical java interface. Further analyses were performed according to the Blast2GO user manual.
MCL analysis. Gene families and clusters of orthologous genes were inferred using OrthoMCL 53 (version 2.0) with standard parameters and granularity 1.5 for the MCL clustering step. Functional enrichment and overrepresentation analyses were performed using a Fisher's exact test, adjusting for FDR. For each gene family inferred with orthoMCL, a multiple sequence alignment of the protein sequences was obtained using Clustal Omega 54 and an HMM model was generated with the hhmake program of the HHSuite toolkit 55 . Sequences from the fungal database fuNOG 56 were similarly aligned and HMM models generated. To annotate whole gene families, the hhsearch program was used to obtain matches between the gene family and the fuNOG HMMs and only hits with a probability equal to or higher than 0.99 were considered. To annotate whole gene families, the hhsearch program was used to obtain matches between the gene family and the fuNOG HMMs and only hits with a probability Z0.99 were considered.
Ancestral genome reconstruction. Gene families inferred with OrthoMCL were used to reconstruct the ancestral genomes of each Colletotrichum lineage. GLOOME 57 (maximum-likelihood approach) was used to infer ancestral gene gains and losses (GGLs) and to reconstruct the ancestral GGLs of gene families on the species tree of Ct0861 and the other five genomes available for this genus. Evolution of the GGLs along the branches of a phylogenetic tree was modelled as a continuous time Markov process using a binary character alphabet corresponding to gene family presence or absence. Default parameters were used, corresponding to a mixture model that allows varying GGL rates across gene families. We approximated the total number of gene families that were gained or lost on a branch by summing up the individual posterior probabilities for each gene family to be gained or lost on that branch and rounding this number to the closest integer. The number of genes either gained or lost (annotated with one specific category) was compared with the respective numbers detected for all other branches of the tree. The significance was assessed using Fisher's exact test and FDR corrected. d N /d S analysis. A multiple-sequence alignment (MSA) of orthologous groups of coding sequences (CDSs) was created with Clearcut 58 . Based on the MSA and the CDSs, a codon alignment was constructed for each protein family with pal2nal (ref. 59; version 14) using default parameters. Because of the data set size and the shorter runtime of neighbour joining algorithms compared with maximumlikelihood methods, Clearcut, a relaxed neighbour joining algorithm 58 , was chosen for reconstructing phylogenetic trees from the MSA of each protein family with slightly modified additive pairwise distances whereby gaps are not counted as mismatches. Gaps in this alignment were mostly of technical origin due to the alignment of short contigs to longer reference sequences. Using an in-house tool (phylorecon), CDSs and amino acid sequences were reconstructed for the internal nodes of each phylogenetic tree using maximum parsimony as a criterion 60 , and the synonymous and non-synonymous substitution rates per site were inferred with correction for multiple substitutions. The average d N /d S ratio was calculated for each protein family and a one-sided Fisher's test (FDR corrected) was performed to identify protein families with a significant enrichment of synonymous mutations per synonymous site versus non-synonymous mutations per non-synonymous site.
Annotation of specific gene categories. Secretomes of all species were predicted using WoLF-PSORT 61 with default settings. Colletotrichum CSEPs were defined as extracellular proteins with no significant BLAST homology (E-value o1 Â 10 À 3 ) to sequences outside the genus Colletotrichum in the UniProt database (SwissProt and TrEMBL components). To identify secreted proteases, sequences of predicted extracellular proteins were subjected to a MEROPS Batch BLAST analysis 62 .
Membrane transporters were identified and classified through BLAST searches against the Transporter Collection Database (http://www.tcdb.org/). To predict the repertoire of carbohydrate-active enzymes encoded by Colletotrichum species, we scanned their genomes using the CAZy annotation pipeline 63 (http:// www.cazy.org). For annotating genes encoding secondary metabolism key enzymes in Colletotrichum species, we used an in-house bioinformatics pipeline that was developed as described in Supplementary Note 12.
RNA sequencing. The RNA-seq samples presented in Hiruma et al. 8 and the new samples presented here were prepared as follows. Fungal cultures were maintained on Mathur's agar medium at 25°C, and conidia were harvested from 7-to 10-day-old cultures. For sample preparation, Arabidopsis Col-0 seeds were surface sterilized in 70% ethanol and subsequently in 2% hypochlorous acid (v/v) containing 0.05% (v/v) Triton. We inoculated A. thaliana Col-0 seeds with spores (5 Â 10 4 spores ml À 1 ) of Ct 0861 or Ci and transferred the inoculated seeds onto solid half-strength Murashige and Skoog medium (pH ¼ 5.1) either in normal [625 mM] or low phosphate [50 mM] conditions. For each biological replicate (n ¼ 3), the entire root system of at least 10 plants was collected at time intervals (6, 10, 16 or 24 d.p.i.) and pooled before RNA extraction. In addition, we grew Ct and Ci in liquid Mathur's medium (in vitro samples) for 2 days at 24°C with shaking at 50 r.p.m. and collected the hyphae by filtration. Total RNA was purified with the NucleoSpin RNA plant kit (Macherey-Nagel) according to the manufacturer's protocol. RNA-seq libraries were prepared from an input of 1 mg total RNA using the Illumina TruSeq stranded RNA sample preparation kit. Libraries were subjected to paired-end sequencing (100 bp reads) using the Illumina HiSeq2500 Sequencing System. To make sure the sequenced reads were of sufficiently high quality, an initial quality check was performed using the FastQC suite (http:// www.bioinformatics.babraham.ac.uk/projects/fastqc/). Subsequently, the RNA-seq reads were mapped to the assembled and annotated genomes of either Ct 0861 or Ci, and in parallel to the annotated genome of the host plant A. thaliana (TAIR10) using Tophat2 (ref. 64; a ¼ 10, g ¼ 10, r ¼ 100, mate-std-dev ¼ 40). The mapped RNA-seq reads were then transformed into a fragment count per gene per sample using the htseq-count script (s ¼ reverse, t ¼ exon) in the package HTSeq 65 . The complete RNA-Seq data presented by Hiruma et al. 8 and in this manuscript have been deposited under the GEO series accession number GSE70094. Statistical analysis of differential gene expression. All statistical analyses of plant and fungal gene expression were performed in R (codes are available upon request). For the analyses of plant gene expression, genes with less than 100 mapped fragments in total (that is, across all the analysed samples) were rated as 'not expressed' and therefore excluded. For analyses of fungal gene expression, we excluded genes that were not sufficiently expressed in the in planta samples, that is, genes with less than 100 (Ct, 24 samples) or less than 50 (Ci, 6 samples) mapped fragments across all the analysed samples. Subsequently, the count data for all expressed genes was TMM-normalized and log-transformed using the functions 'calcNormFactors' (R package EdgeR 66 ) and 'voom' (R package limma 67 ) to yield log 2 counts per million (log 2 cpm). To analyse the aspects of differential gene expression in Ct0861, Ci and their host plant Arabidopsis, we fitted for each analysis a distinct linear model to the respective log 2 -transformed count data using the function lmFit (R package limma 67 ) and subsequently performed moderated t-tests for specific comparisons of interest. Resulting P values were adjusted for false discoveries due to multiple hypotheses testing via the Benjamini-Hochberg procedure (FDR). To extract genes with significant expression differences, a cutoff of FDRo0.05 and |log 2 FC|Z1 was applied. Heatmaps of gene expression profiles were generated with the Genesis expression analysis package 68 and interactive Tree Of Life 69 was used to visualize CSEP gene expression data. To derive Arabidopsis, Ct and Ci gene expression profiles during the time-course experiment, log 2 expression ratios were calculated between the normalized number of reads detected for a given gene at a given developmental stage and the geometrical mean of the number of reads calculated across all developmental stages. This log 2 ratio is referred as the 'Relative Expression Index'. The Cytoscape plug-in ClueGO þ CluePedia 70 was used to construct GO term enrichment networks and to visualize functionally grouped terms among significantly regulated genes. Significant enrichments were determined using the hypergeometric test and Bonferroni step-down corrected P values are represented. Co-regulated genes that were also co-expressed in other Arabidopsis expression data sets were identified using ATTED-II (http://atted.jp/) and co-expression networks were generated using Cytoscape 71 (version 3.1.1).
RT-qPCR analysis. First-strand cDNA was synthesized from 1 mg DNase-treated total RNA using the iScript cDNA synthesis kit (Bio-Rad) and PCR amplification was performed using the iQ5 real-time PCR detection system (Bio-Rad). For each gene, specific primers were designed with the Primer 3 and AmplifX programs. BLASTN searches against the Ct and A. thaliana genomes were performed to rule out cross-annealing artefacts. Gene expression levels were normalized using the reference gene actin (ACT2, AT3G18780) for A. thaliana and the reference gene tubulin beta-1 chain (CT04_12898) for Ct. These genes were used to normalize gene expression levels using the Pfaffl calculation method 72 .
Microscopy methods. For cytology experiments, surface-sterilized A. thaliana Col-0 seeds were inoculated with either Ct or Ci conidia (5 Â 10 4 spores ml À 1 ). The seeds were then transferred to half-strength Murashige and Skoog agarose medium without sucrose and low-phosphate content (50 mM). Inoculated plants were grown at 22°C with a 10-h photoperiod (80 mE m À 2 s À 1 ) for 1 to 24 days. The roots were either mounted in water for viewing GFP or first stained with Calcofluor white (0.01 %, Sigma) or fluorescein diacetate (10 mg ml À 1 , Sigma). For visualizing GFP and FDA fluorescence, we used an Olympus FV1000 confocal microscope equipped with dry Â 20 and Â 40 objectives, using the 488-nm line of an Argon laser for excitation and fluorescence was collected at 490-520 nm. For imaging Calcofluor fluorescence, we used a Zeiss Axiophot epifluorescence microscope (filter set BP 365, FT 395, LP 397).