Shifts in diversification rates and host jump frequencies shaped the diversity of host range among Sclerotiniaceae fungal plant pathogens.

The range of hosts that a parasite can infect in nature is a trait determined by its own evolutionary history and that of its potential hosts. However, knowledge on host range diversity and evolution at the family level is often lacking. Here, we investigate host range variation and diversification trends within the Sclerotiniaceae, a family of Ascomycete fungi. Using a phylogenetic framework, we associate diversification rates, the frequency of host jump events and host range variation during the evolution of this family. Variations in diversification rate during the evolution of the Sclerotiniaceae define three major macro-evolutionary regimes with contrasted proportions of species infecting a broad range of hosts. Host-parasite cophylogenetic analyses pointed towards parasite radiation on distant hosts long after host speciation (host jump or duplication events) as the dominant mode of association with plants in the Sclerotiniaceae. The intermediate macro-evolutionary regime showed a low diversification rate, high frequency of duplication events and the highest proportion of broad host range species. Our findings suggest that the emergence of broad host range fungal pathogens results largely from host jumps, as previously reported for oomycete parasites, probably combined with low speciation rates. These results have important implications for our understanding of fungal parasites evolution and are of particular relevance for the durable management of disease epidemics.

INTRODUCTION

The host range of a parasite has a central influence on the emergence and spread of disease (Woolhouse & Gowtage‐Sequeria, 2005). There is a clear demarcation between specialist parasites that can only infect one or a few closely related host species, and generalists that can infect more than a hundred unrelated host species (Barrett, Kniskern, Bodenhausen, Zhang, & Bergelson, 2009; Woolhouse, Taylor, & Haydon, 2001). Host specialization, when lineages evolve to infect a narrower range of hosts than related lineages, is a frequent occurrence in living systems and can be driven by a parasite sharing habitat with only a limited number of potential hosts. There are also clear examples of parasite adaptations that restrict the use of co‐occurring potential hosts. For instance, the red rust fungal pathogen Coleosporium ipomoeae infects fewer Ipomoea species from its native community than when inoculated to non‐native communities, implying that evolution within local communities has narrowed pathogen host range (Chappell & Rausher, 2016). Some isolates of the rice blast fungus Magnaporthe oryzae are able to infect Oryza sativa japonica varieties but not indica varieties co‐occurring in the Yuanyang area of China, because their genome harbours numerous avirulence effector genes (Liao et al., 2016). It has been proposed that specialization results from trade‐offs between traits needed to infect a wide range of hosts (Futuyma & Moreno, 1988; Joshi & Thompson, 1995). It is frequently associated with the loss of traits that are not required to infect a particular host, such as the loss of lipid synthesis in parasitoid wasps (Visser et al., 2010), and the loss of secondary metabolism and carbohydrate active enzymes in powdery mildew fungal pathogens (Spanu et al., 2010). Specialization is sometimes considered as an evolutionary dead end (Moran, 1988), as gene losses are often irreversible and may lead specialist lineages “down a blind alley” that limits transitions back to generalism (Day, Hua, & Bromham, 2016; Haldane, 1951). There is nevertheless evidence for transitions from specialist to generalist parasitism (Hu et al., 2014; Johnson, Malenke, & Clayton, 2009). For instance in plant pathogens, Pseudoperonospora cubensis differs from other downy mildew oomycete pathogens in that it is able to infect a wide range of Cucurbits (Thines & Choi, 2015). For many parasite lineages however, knowledge on host range diversity and evolution at the macro‐evolutionary level is lacking.

Sclerotiniaceae is a family of Ascomycete fungi of the class Leotiomycetes which includes numerous plant parasites. Among the most studied are the grey mould pathogen Botrytis cinerea, considered to be one of the 10 most devastating plant pathogens (Dean et al., 2012), and the white and stem mould pathogen Sclerotinia sclerotiorum. Both are economically important pathogens in agriculture that infect hundreds of host plant species (Bolton, Thomma, & Nelson, 2006; Mbengue et al., 2016). The Sclerotiniaceae family also includes host specialist parasites such as Ciborinia camelliae that causes flower blight on Camellia (Denton‐Giles, Bradshaw, & Dijkwel, 2013), Sclerotinia glacialis that specifically infects Ranunculus glacialis (Graf & Schumacher, 1995), and Monilinia oxycocci causing the cottonball disease on cranberry (McManus, Best, & Voland, 1999). Other species from the Sclerotiniaceae have intermediate host range (tens of plant species) such as Sclerotinia trifoliorum, S. subarctica and S. borealis (Clarkson, Carter, & Coventry, 2010; Farr & Rossman, 2016). While S. sclerotiorum and B. cinerea are considered as typical necrotrophic pathogens, rapidly killing host cells to cause disease, the Sclerotiniaceae include species with diverse lifestyles. For instance, the poplar pathogen Ciborinia whetzelii and several Myriosclerotinia species are biotrophs that can live as endophytes (Andrew, Barua, Short, & Kohn, 2012; Schumacher & Kohn, 1985), while Coprotinia minutula is coprophilous (Elliott, 1967). How this remarkable diversity evolved remains elusive. To gain insights into this question, knowledge of phylogenetic relationships and host range diversity at the macro‐evolutionary level is needed. Specifically, ancestral state reconstruction can provide insights into how fungal host range has changed over time. The relationship between host range and diversification rates is an active area of research in insect ecology (Hamm & Fordyce, 2015; Hardy & Otto, 2014). For instance, a positive correlation was reported between species richness in the butterfly family Nymphalidae and the diversity of plants they feed upon (Janz, Nylin, & Wahlberg, 2006). However in contrast, some studies found a negative relationship between host‐plant breadth and diversification rate (Hardy & Otto, 2014). Transitions between feeding strategies generally appeared to be associated with shifts in insect diversification rates (Hardy & Otto, 2014; Janz & Nylin, 2008). Estimating fungal species divergence times in the Sclerotiniaceae will allow testing of whether biological diversification is related to host range variation.

Host range is a trait determined not only by the evolutionary history of a parasite, but also by that of its potential hosts (Poulin & Keeney, 2008). Accounting for host association patterns should therefore prove useful to understand host range evolution in the Sclerotiniaceae. The Leotiomycete class diverged less than 200 million years ago (Beimforde et al., 2014; Prieto & Wedin, 2013) and likely radiated with the diversification of flowering plants (Smith, Beaulieu, & Donoghue, 2010). Molecular phylogenetic studies distinguished Myriosclerotinia, Sclerotinia sensu stricto, Botrytis and Botryotinia, and Monilinia sensu stricto as monophyletic clades within the Sclerotiniaceae family (Holst‐Jensen, Vaage, & Schumacher, 1998; Holst‐Jensen, Vrålstad, & Schumacher, 2004). A phylogenetic analysis of Monilinia species suggested that cospeciation with host plants was the dominant pattern in this clade (Holst‐Jensen, Kohn, Jakobsen, & Schumacher, 1997). By contrast, there was no evidence for cospeciation with host plants in a phylogenetic analysis of Botrytis species (Staats, Van Baarlen, & Van Kan, 2005). Botrytis species were thus proposed to have evolved through host jumps to unrelated host plants followed by adaptation to their new hosts (Dong, Raffaele, & Kamoun, 2015; Staats et al., 2005). Host–parasite cophylogenetic analyses are required to test whether variations in the frequency of host jumps may have impacted on host range variation in the Sclerotinaceae.

In this study, we performed a phylogenetic analysis on 105 Sclerotiniaceae species to reveal multiple independent shifts and expansions of host range. We show that three macro‐evolutionary regimes with distinct diversification rates and dominant host association patterns have shaped the diversity of the Sclerotiniaceae and lead to contrasted proportions of broad host range species. Specifically, we highlight an increased emergence of broad host range parasites during the transition between macro‐evolutionary regimes dominated by distinct patterns of host–pathogen association. These results suggest that reduced diversification rates and high host jump frequency could associate with the emergence of generalist pathogens.

MATERIALS AND METHODS

Taxon and host range data selection

We used all 105 Sclerotiniaceae species for which at least ITS marker sequence data were available in the GenBank database. As outgroups we selected 56 Rutstroemiaceae species and 39 representative species of Leotiomycetes and Sordariomycetes for a total of 200 species. Host range data were obtained from (Boland & Hall, 1994; Melzer, Smith, & Boland, 1997), Index Fungorum (http://www.indexfungorum.org), the SMML fungus‐host distribution database (Farr & Rossman, 2016) and references therein. In total, we retrieved 7,101 fungus‐host association records that did not show strong geographic or crop/wild species bias (Figure S1). Fungal host range was not correlated to the number of database records (Figure S1c). Calibrated trees deciphering the relationships between the host families based on seven gene regions (18S rDNA, 26S rDNA, ITS, matK, rbcL, atpB and trnL‐F) were extracted from (Qian & Zhang, 2014) and updated with (Hedges, Marin, Suleski, Paymer, & Kumar, 2015). The tree of host families used for cophylogenetic analyses is provided as File S1.

Initial phylogenetic analysis

ITS sequences were aligned using mafft version 7 (Katoh & Standley, 2013). The alignment was manually adjusted to minimize possible homoplasic positions using Seaview 4 (Gouy, Guindon, & Gascuel, 2010). We retained gaps shorter than 41 positions present in a maximum of 39 sequences, with ungapped blocks being at least five characters long with a maximum of 5% ambiguous nucleotides per position. Unalignable and autapomorphic regions were excluded from the analysis, yielding an alignment with 797 informative sites (File S2). Maximum‐likelihood phylogeny was inferred with PhyML 3 (Guindon et al., 2010) with Smart Model Selection which permits an automatic substitution model selection supporting the general time‐reversible model with gamma distribution using six substitution rate categories (GTR+G6) model as the best fit. Statistical branch support was inferred with the SH‐like approximate likelihood ratio test (SH‐aLRT) (Anisimova, Gil, Dufayard, Dessimoz, & Gascuel, 2011) and bootstrap analysis with 100 replicates (Files 3 and 4). The topology of the tree was confirmed by three additional methods. First, we used a neighbour‐joining approach in FastME 2.0 (Lefort, Desper, & Gascuel, 2015) with the LogDet substitution model and tree refinement by subtree pruning and regrafting (File S5). Second, we used the parsimony ratchet approach implemented in the phangorn package in r (Schliep, 2011; File S6). Third, we used Bayesian analysis in MrBayes 3.2.6 (Ronquist et al., 2012) using the GTR substitution model with gamma‐distributed rate variation across sites and a proportion of invariable sites in 100 million MCMC generations, sampling parameters every 1,000, removing 10% of the tree files as burn‐in. The resulting trees, available in Dryad Digital Repository (https://doi.org/10.5061/dryad.7cs3g), were edited with Figtree (available at http://tree.bio.ed.ac.uk/software/figtree/).

Ancestral state reconstruction

To infer possible ancestral hosts, we used Reconstruct Ancestral State in Phylogenies 3.1 (RASP) (Yu, Harris, Blair, & He, 2015). We used the S‐DIVA (Statistical Dispersal‐Vicariance Analysis), S‐DEC (Statistical Dispersal–Extinction–Cladogenesis model), BBM (Bayesian Binary MCMC) and BayArea methods to verify congruence between the methods and assess the robustness of their output, as recommended by (Yu et al., 2015). We report the results of the S‐DIVA analysis which is considered the best adapted for host–parasite association analyses (Razo‐Mendivil & De Leon, 2011; Yu et al., 2015). To implement S‐DIVA, we delimited host groups as follows: Vitales (A), Asterids (B), Campanuliids (C), Commelinids (D), coprophilous (E), core Eudicots (F), Eudicots (G), Fabids (H), Polypodiidae (I), Lamiids (J), Magnoloidae (K), Malvids (L), Monocotyledones (M), Pinidae (N), allowing a maximum of four groups at each node. Among these groups, the latest divergence is that of Malvids and Fabids, estimated around 82.8 and 127.2 Mya (Clarke, Warnock, & Donoghue, 2011) and predates the emergence of Sclerotiniaceae estimated at ~69.7 Mya in this work; we therefore did not restrict associations in the ancestral reconstruction analysis. To account for incomplete sampling in this analysis, ancestral state reconstruction was computed for every plant group by the re‐rooting method (Yang, Kumar, & Nei, 1995) under an entity‐relationship (ER) model available in phytools (Revell, 2012) and derived from the ape library in R (Paradis, Claude, & Strimmer, 2004). This method, based on maximum likelihood, computed the probability to infect a given group of hosts for every nodes of the Sclerotiniaceae phylogeny. We considered plant groups as hosts when the probability of the ancestral state was >50%. For every plant group, up to 10% of terminal nodes were pruned randomly 100 times and the ancestral state reconstructed on pruned trees. For each plant group, we then extracted the variation in the age of the most ancestral inclusion into the Sclerotiniaceae host range (Figure S2). For all plant groups, this variation was not significantly different from 0, indicating that ancestral state reconstruction was robust to tree pruning.

Divergence dating analyses

We used a Bayesian approach to construct a chronogram with absolute times with the program beast 1.8.2. As no fossil is known in this fungal group, we use Sordariomycetes–Leotiomycetes divergence time to calibrate the tree (~300 Mya, Beimforde et al., 2014). Using our initial phylogeny as a constraint, the partition file was prepared with the BEAUti application of the beast package. Considering that our data set includes distantly related taxa, subsequent analyses were carried out using an uncorrelated relaxed clock model with lognormal distribution of rates. We used a uniform birth–death model with incomplete species sampling as prior on node age, following (Beimforde et al., 2014) to account for incomplete sampling in our phylogeny. Analyses were run three times for 100 million generations, sampling parameters every 5,000 generations, assessing convergence and sufficient chain mixing using Tracer 1.6 (Drummond, Suchard, Xie, & Rambaut, 2012). We removed 20% of trees as burn‐in, and the remaining trees were combined using LogCombiner (beast package), and summarized as maximum clade credibility (MCC) trees using TreeAnnotator within the beast package. The trees were edited in FigTree. The r packages ape (Paradis et al., 2004) and (Revell, 2012) were used for subsequent analyses.

Diversification analysis

We used the section containing all 105 Sclerotiniaceae species of the chronogram with absolute times generated by beast to perform diversification analyses with three different methods. We first used Bayesian Analysis of Macroevolutionary Mixtures (BAMM) version 2.5 (Rabosky et al., 2013, 2014), an application designed to account for variation in evolutionary rates over time and among lineages. The priors were set using the setBAMMPriors command in the bammtools r‐package (Rabosky et al., 2014). We ran four parallel Markov chains for 5,000,000 generations and sampled for every 1,000 trees. The output and subsequent analyses were conducted with bammtools. We discarded the first 10% of the results and checked for convergence and the effective sample size (ESS) using the coda r‐package. Lineage‐through‐time plots with extant and extinct lineages were computed using the phytools package in r (Revell, 2012) with the drop.extinct=FALSE parameter. The exact number of Sclerotiniaceae species is currently unknown, and no intentional bias was introduced in data collection (Figure S1). To detect shifts in diversification rates taking into account incomplete sampling of the Sclerotiniaceae diversity and tree uncertainty, we used MEDUSA, a stepwise approach based upon the Akaike information criterion (AIC) (Alfaro et al., 2009; Drummond, Eastwood, Miotto, & Hughes, 2012). For this analysis, we used the birth and death (bd) model, allowed rate shifts both at stem and nodes, used the AIC as a statistical criterion with initial r = .05 and ε = 0.5. Sampling of the fungal diversity is typically estimated around 5% (Hawksworth & Luecking, 2017). We used random sampling rates between 5% and 100% for each individual species in 100 MEDUSA bootstrap replicates to estimate the impact of sampling richness. Next, we randomly pruned 1% to 10% of species from the tree in 100 MEDUSA bootstrap replicates to estimate the impact of tree completeness. To control for the impact of divergence time estimates on the diversification analysis, we altered all branching times in the tree randomly by −15% to +15% in 100 MEDUSA bootstrap replicates (Figure S3). Finally, we used rpanda which includes model‐free comparative methods for evolutionary analyses (Morlon et al., 2016). The model‐free approach in rpanda compares phylogenetic tree shapes based on spectral graph theory (Lewitus & Morlon, 2015). It constructs the modified Laplacian graph of a phylogenetic tree, a matrix with eigenvalues reflecting the connectivity of the tree. The density profile of eigenvalues (Figure S4a) provides information on the entire tree structure. The algorithm next used k‐medioids clustering to identify eight modalities in the phylogenetic tree of the Sclerotiniaceae (Figure S4b,c). A post hoc test comparing Bayesian information criterion (BIC) values for randomly bifurcated trees (BICrandom) with that of the tree of the Sclerotiniaceae (BICSclerotiniaceae) supported the eight modalities (BICrandom/BICSclerotiniaceae = 10.34, well above the significance threshold of 4.0). We obtained significant BIC ratios with two or more modalities, supporting the existence of at least two macro‐evolutionary regimes in the Sclerotiniaceae.

Cophylogeny analysis

To detect potential codivergence patterns, we tested for cophylogeny between fungal species and host family trees with CoRe‐PA 0.5.1 (Merkle, Middendorf, & Wieseke, 2010), PACo (Balbuena, Míguez‐Lozano, & Blasco‐Costa, 2013) and Jane 4 (Conow, Fielder, Ovadia, & Libeskind‐Hadas, 2010) as these programs can accommodate pathogens with multiple hosts. The only three species which were not able to interact with an angiosperm host (saprotroph or only described on gymnosperms) were excluded from the analysis (namely: Elliottinia kerneri, Coprotinia minutula and Stromatinia cryptomeriae). We used a full set of 263 host–pathogen associations and a simplified set of 121 associations minimizing the number of host families involved to control for the impact of sampling bias (Table 1). We also performed the analyses independently for each macro‐evolutionary regime. To estimate cost parameters in CoRe‐PA, we ran a first cophylogeny reconstruction with cost values calculated automatically using the simplex method on 1,000 random cycles, with all host switches permitted, direct host switch and full host switch permitted (other parameters set to CoRe‐PA default). The best reconstruction had a quality score of 1.348 with total costs 37.076 and calculated costs of 0.0135 for cospeciation, 0.1446 for sorting, 0.2501 for duplication and 0.5918 for host switch. With these costs and parameters, we then computed reconstruction with 1,000 random associations to test for the robustness of the reconstructions. Next, we used 100,000 permutations of the host–pathogen association matrix in PACo to test for the overall congruence between host and pathogen trees. To test for the contribution of each association on the global fit, we performed taxon jackknifing. We used one‐tailed z tests to compare squared residuals distribution for each host–pathogen link with the median of squared residuals for the whole tree and infer likely co‐evolutionary and host shift links at p < .01. We used optimal costs calculated by CoRe‐PA to run cophylogeny reconstruction in Jane 4 Solve Mode with 50 generations and a population size of 1,000. Statistical tests were performed using 100 random tip mapping with 20 generations and population size of 100. For the analysis of each macro‐evolutionary regime in Jane 4, we used simplified host–pathogen association set to reduce the amount of Losses and Failure to Diverge predicted and facilitate comparisons with CoRe‐PA and PACo results. We visualized the interactions with TreeMap3 (Charleston & Robertson, 2002), using the untangling function to improve the layout.

Table 1

Results of the cophylogeny analyses under CoRe‐PA optimized cost settings using various methods and host association sets

		Event frequency (% of host associations)
Method	Host associationsa	Cospeciation	Sorting/Loss	Duplication	Host switch	FD	Cost/p‐valb
CoRe‐PAc	Full set (263)	15.82 ± 2.34	15.90 ± 4.80	34.35 ± 3.35	33.94 ± 3.49		38.21 ± 0.79
	Simplified (121)	9.63 ± 2.63	11.09 ± 4.97	40.20 ± 3.25	39.07 ± 4.00		40.44 ± 0.81
	G1 only (68)	15.75 ± 3.19	16.32 ± 7.35	37.78 ± 4.67	30.14 ± 5.23		13.95 ± 0.44
	G2 only (130)	17.03 ± 3.57	16.21 ± 7.58	33.16 ± 5.27	33.60 ± 5.63		14.57 ± 0.47
	G3 only (65)	11.78 ± 5.17	14.38 ± 8.69	32.86 ± 6.75	40.98 ± 6.87		8.95 ± 0.34
PACod	Full set (263)	20.91	48.29		30.79		p < 1.0−05
	Simplified (121)	15.70	57.02		27.27		p < 1.0−05
	G1 only (68)	27.94	55.88		16.18		p = 2.0−05
	G2 only (130)	6.15	71.54		22.31		p = .00045
	G3 only (65)	0.00	66.15		33.85		p = .00607
Jane 4	Full set (263)	0.86	67.73	8.00	3.57	19.83	1,291 (p < .01)
	Simplified (121)	2.90	65.22	14.78	11.59	5.51	719.47 (p < .01)
	G1 simple (43)	1.39	41.67	25.00	26.39	5.56	206.98 (p < .01)
	G2 simple (54)	2.28	75.80	11.87	3.65	6.39	374.08 (p < .01)
	G3 simple (24)	9.68	25.81	19.35	41.94	3.23	105.41 (p = .13)

Refers to the set of host association tested for cophylogeny: “Full set” corresponds to the complete list of all plant‐Sclerotiniaceae associations; “Simplified” corresponds to a reduced set covering the whole Sclerotiniaceae family; “G1,” “G2” and “G3” corresponds to associations involving Sclerotiniaceae species from macro‐evolutionary regime G1, G2 or G3 only.

Values correspond to CoRe‐Pa total reconstruction costs, p‐value of the observed host–parasite association matrix m² in 100,000 permutations with PACo, or the p‐value of the observed reconstruction cost in 100 random tip mappings in Jane 4. FD, “failure to diverge,” corresponding to events when a host speciates and the parasite remains on both new host species.

Standard deviation corresponds to frequencies calculated for 1,000 reconstruction with randomized host–parasite associations. Associations are classified as “cospeciation” when speciation of host and pathogen occurs simultaneously; “duplication” when pathogen speciation occurs independently of host speciation; “sorting or loss” when a pathogen remains associated with a single descendant host species after host speciation; and “host switch” when a pathogen changes host independently of speciation events.

In PACo taxon jackknifing, associations that contributed significantly and positively to cophylogeny were classified as “cospeciation,” significantly and negatively as “host switch” and associations with no significant contribution to cophylogeny are classified as either sorting/loss or duplication.

RESULTS

Multiple independent shifts and expansions of host range in the evolution of the Sclerotiniaceae

To document the extant diversity in the Sclerotiniaceae fungi, we collected information on the host range of 105 species in this family. For comparison purposes, we also analysed the host range of 56 species from the Rutstroemiaceae family, a sister group of the Sclerotiniaceae (Figure 1a, Table S1). To reduce biases that may arise due to missing infection reports, we analysed host range at the family level instead of genus or species level. We found one species in the Sclerotiniaceae (Coprotinia minutula) and two species in the Rutstroemiaceae (Rutstroemia cunicularia and Rutstroemia cuniculi) reported as nonpathogenic to plants that are coprophilous (Elliott, 1967). Most species in the Sclerotiniaceae are reported necrotrophic pathogens, such as Botrytis, Monilinia and Sclerotinia species. Exceptions include Ciborinia whetzelii which was reported as an obligate biotroph of poplar (Andrew et al., 2012) and Myriosclerotinia species that specifically infect monocot families as facultative biotrophs or symptomless endophytes (Andrew et al., 2012). The most frequently parasitized plant group was Fabids (by 47 Sclerotiniaceae species and 26 Rutstroemiaceae species), a group of the Rosids (Eudicots) including notably cultivated plants from the Fabales (legumes) and Rosales (rose, strawberry, apple…) orders. Plants from the order Vitales (e.g., grapevine), from the Magnoliids and from the Polypodiidae (ferns) were colonized by Sclerotiniaceae but not Rutstroemiaceae fungi.

Figure 1

Multiple independent shifts and expansions of host range in the evolution of the Sclerotiniaceae. (a) Distribution of plant hosts of parasites from the Sclerotiniaceae and Rutstroemiaceae fungi. (b) Maximum‐likelihood phylogeny of 105 Sclerotiniaceae and 56 Rutstroemiaceae species showing host range information and ancestral host reconstruction. Host range is shown as circles at the tips of branches, sized according to the number of host families and coloured as in (a) according to the earliest diverging plant group in host range. Numbers at the tips of branches refer to species listed in Table S1. Branch support indicated in light red for major clades corresponds to SH‐aLRT (regular), bootstrap (bold) and Bayesian posterior probabilities (italics). Reconstructed ancestral host is shown as triangles at intermediate nodes when a change compared to the previous node is predicted. Endophytes and biotrophic parasites are shown with empty circles. (c) Distribution of Sclerotiniaceae and Rutstroemiaceae species according to their number of host families [Colour figure can be viewed at http://wileyonlinelibrary.com]

We next used the internal transcribed spacer (ITS) region of rDNA sequences to construct a phylogenetic tree of the 105 Sclerotiniaceae and 56 Rutstroemiaceae species (Figure 1b, Files S1–S7). This marker is the only one available for all species and is currently considered as the universal marker for taxonomic use (Schoch et al., 2012). Maximum‐likelihood, neighbour‐joining, parsimony and Bayesian analyses yielded convergent tree topologies. The tree topologies confirmed the familial classification recognized in previous phylogenies (Andrew et al., 2012; Holst‐Jensen et al., 1997, 1998), identifying a clearly supported Botrytis genus including 23 species; a Myriosclerotinia genus including eight species; a Sclerotinia genus sensu stricto including four Sclerotinia; and a Monilinia genus including 16 Monilinia species and three other species. Botrytis calthae and Amphobotrys ricini were the only two species with imprecise positions in maximum‐likelihood reconstructions, due to the ITS marker being insufficient to infer complete phylogenetic placement (Lorenzini & Zapparoli, 2016; Staats et al., 2005).

We used RASP (Yu et al., 2015) and phytools (Revell, 2012) to reconstruct the ancestral host range across the phylogeny (Figures 1b, S1 and S5). This identified Fabids as the most likely ancestral hosts of the Sclerotiniaceae family (relative probability of 100% in S‐DIVA and S‐DEC analyses). Random tree pruning in phytools indicated that this result is robust to sampling biases. Over 48% of the Sclerotiniaceae species are pathogens of host plants that evolved prior to the divergence of the Fabids, suggesting numerous host jumps in this family of parasites. Notably, a jump to Malvids and then to Monocots was identify at the base of the Botrytis genus (87% and 85% probability in S‐DIVA, respectively), a jump to Commelinids occurred in the Myriosclerotinia genus (91% probability in S‐DIVA), a jump to the Ranunculales occurred at the base of the Sclerotinia genus (76% probability in S‐DIVA); a jump to Asterids was found at the base of a major group of Monilinia (89% probability in S‐DIVA).

A total of 73 Sclerotiniaceae species (69.5%) and 41 (73.2%) Rutstroemiaceae infected a single host family (Figure 1c, Table S1). Moellerodiscus lentus was the only Rutstroemiaceae species infecting hosts from more than five plant families whereas eleven Sclerotiniaceae species (10.4%) exhibited this trait, including Botrytis cinerea, Sclerotinia sclerotiorum, Sclerotinia minor and Grovesinia pyramidalis, each of which colonizes plants from more than 30 families. Each of these species belongs to a clearly distinct phylogenetic group with a majority of species infecting a single host family. This may result from radiation following host jumps (Choi & Thines, 2015) and suggests that the ability to colonize a broad range of plant was acquired multiple times independently through the evolution of the Sclerotiniaceae.

Two major diversification rate shifts in the evolution of the Sclerotiniaceae

To test for a relationship between host range variation and biological diversification in the Sclerotiniaceae, we estimated divergence times for the fungal family Sclerotiniaceae using the ITS marker in a Bayesian framework. A calibrated maximum clade credibility chronogram from these analyses is shown in Figure 2a (and Figure S6). This showed that Sclerotiniaceae fungi shared a most recent common ancestor between 33.9 and 103.5 million years ago (Mya), with a mean age of 69.7 Mya. The divergence of Botrytis pseudocinerea, estimated from the ITS data set, occurred ca. 3.35–17.8 Mya (mean age 9.8 Mya) and is similar to a previous estimate of 7–18 Mya (Walker et al., 2011).

Figure 2

Two major diversification rate shifts in the evolution of the Sclerotiniaceae. (a) Dated ‐based species tree for the Sclerotiniaceae with diversification rate estimates. The divergence times correspond to the mean posterior estimate of their age in millions of years calculated with beast. Mean age is shown for selected nodes with bars showing 95% confidence interval of the highest posterior density (HPD). Branches of the tree are colour‐coded according to diversification rates determined with BAMM. Major rate shifts identified in BAMM are shown as red circles, noted s 1, , s 2 and and labelled with the posterior distribution in the 95% credible set of macro‐evolutionary shift configurations. Species names are shown in black if host range includes less than five plant families, in yellow for five to nine plant families and in red for 10 or more plant families. Diversification rate shifts define three macro‐evolutionary regimes noted G1, G2 and G3 and boxed in blue, grey and brown, respectively. (b) Distribution of broad host range (five or more host families) parasites in Rutstroemiaceae, Sclerotiniaceae and under each macro‐evolutionary regime of the Sclerotiniaceae. p‐Values calculated by random permutations of host ranges along the tree are indicated above bars. Holo., Holocene; Mya, Million years ago; Plei., Pleistocene; Plio., Pliocene; Quat. Quaternary; Rut., Rutstroemiaceae [Colour figure can be viewed at http://wileyonlinelibrary.com]

Using this time‐calibrated phylogeny, we calculated speciation and diversification rates using three different methods, to control for the limits of individual methods (Rabosky, Mitchell, & Chang, 2017). First, we used the BAMM framework, which implements a Metropolis Coupled Markov chain Monte Carlo method to calculate diversification rates along lineages (Rabosky, 2014) (Figures 2a and 7a–c). The analysis identified two significant rate shifts (s 1 and s 2, accounting for >30% of the posterior distribution in the 95% credible set of macro‐evolutionary shift configurations) each immediately adjacent to a minor rate shift ( and , respectively; ≤10% of the posterior distribution) within two different clades (Figure 2a). Each shift affects a specific clade and is not detected coincidentally across the whole phylogeny. The earliest shift occurred between 34.2 and 42.7 Mya. It resulted in an increase in instantaneous diversification rate from ~0.05 lineage/million year in Monilinia and related clades to ~0.08 lineage/million year in Sclerotinia and derived clades. The second shift occurred between 9.8 and 21.7 Mya. It resulted in a further increase in instantaneous diversification to >0.15 lineage/million year in Botrytis lineage. These two shifts defined three distinct macro‐evolutionary regimes in the Sclerotiniaceae. Second, to determine how speciation rates changed over time, we computed lineage‐through‐time plots for the Sclerotiniaceae tree and 1,000 trees in which we altered all branching times randomly by −15% to +15% to control for the sensitivity to divergence time estimates (Figure 3a). We observed an exponential trend of species accumulation consistent with a late burst of cladogenesis or early extinction (Pybus γ = 0.284). Third, we used MEDUSA (Alfaro et al., 2009; Drummond, Eastwood, et al., 2012) to detect shifts in diversification rate in the Sclerotiniaceae phylogeny. MEDUSA reported two shifts, corresponding to shifts s 1 and identified in BAMM (Figure S3d). We used a bootstrap approach to test for the sensitivity of these shifts to sampling richness, tree completeness and dating accuracy (Figure 3b, Figure S3). Shift was detected in 100% of the bootstraps. Shift s 1 was the second most frequent shift detected in our bootstrap analysis although it appeared sensitive to sampling richness in particular. Estimates of speciation rates in MEDUSA were robust to sampling richness, tree completeness and dating accuracy (r = .038 in G1, r = .074 in G2 and r = .31 in G3). These values were consistent with diversification rate estimated in BAMM for each regime (Figure 3c). Fourth, we used model‐free analysis of branching patterns in phylogenetic trees with R‐PANDA (Morlon et al., 2016). The spectrum of eigenvalues suggested eight modes of division in the phylogeny of the Sclerotiniaceae (Figure S4a,b). Comparison to randomly bifurcating trees suggested that these modalities are significant (BIC ratio 0.09). The modalities identified by rpanda k‐means clustering algorithm overlapped with the three macro‐evolutionary regimes identified previously (Figure S4c). Overall, these analyses converged towards the identification of three macro‐evolutionary regimes in the Sclerotiniaceae: species in regime G1 are early diverging and showed low diversification rates, and they encompassed notably Monilinia, Encoelia and some Ciboria species. Regime G2 had intermediate diversification rates and encompassed notably Myriosclerotinia species and most Sclerotinia species. Regime G3 corresponding to Botrytis genus was the most recently diverged and presented the highest diversification rates.

Figure 3

Robustness of diversification rate shifts identification in the Sclerotiniaceae phylogeny. (a) Lineage‐through‐time plots for the Sclerotiniaceae tree and 1,000 trees in which branching times were altered randomly by −15% to +15% to control for the sensitivity to divergence time estimates. Pybus γ for the Sclerotiniaceae tree is provided. (b) Frequency (n) of diversification rate shift detection and diversification rate estimates (r) in a 100 MEDUSA bootstrap replicates in which sampling richness, tree completeness and divergence times were randomly altered. Labels indicate average diversification rate estimates (r) for each macro‐evolutionary regime (blue for G1, grey for G2, brown for G3), with standard deviation of the mean for a 100 replicates. (c) Net diversification rates over time estimated by BAMM for each macro‐evolutionary regime (blue for G1, grey for G2, brown for G3) [Colour figure can be viewed at http://wileyonlinelibrary.com]

We compared the proportion of broad host range (≥5 plant families) fungal species that emerged in the Rutstroemiaceae, in the Sclerotiniaceae, and under each the three macro‐evolutionary regimes in the Sclerotiniaceae (Figure S2b). We randomly permuted host range values across the tree 10,000 times to estimate the p‐values of these proportions occurring by chance. We counted one (1.8%, p = .9959) broad host range species in the Rutstroemiaceae (Moellerodiscus lentus) and eleven (10.4%, p = .0359) in the Sclerotiniaceae. In the Sclerotiniaceae, regime G2 with intermediate diversification rates showed the highest proportion of broad host range species (16.2%, p = .0163), followed by regime G3 (8.7% of broad host range species, p = .537) and regime G1 (5.1% of broad host range species, p = .8407). These analyses suggest that the evolutionary history of regime G2 could have favoured the emergence of broad host range parasites in the Sclerotiniaceae.

Diversification rate shifts associate with variations in rates of cospeciation, duplication and host jump in the Sclerotiniaceae

Within the last 50 Ma, the world has experienced an overall decrease in mean temperatures but with important fluctuations that dramatically modified the global distribution of land plants (Donoghue & Edwards, 2014; Nürk, Uribe‐Convers, Gehrke, Tank, & Blattner, 2015; Zachos, Pagani, Sloan, Thomas, & Billups, 2001). We hypothesized that modifications in the distribution of plants could have an impact on host association patterns in Sclerotiniaceae parasites. To test this hypothesis, we performed host–parasite cophylogeny reconstructions using CoRe‐PA (Merkle et al., 2010), PACo (Balbuena et al., 2013) and Jane 4 (Conow et al., 2010) (Table 1). We found that hosts of Sclerotiniaceae fungi include ~1,800 plant species. To consider this overall host diversity in a cophylogenetic analysis, we used a tree including 56 plant families, a tree of 102 Sclerotiniaceae parasitic species, and 263 host–parasite interactions (File S3). The resulting tanglegram (Figure 4a) highlighted clear parasite duplications notably on Ericaceae and Rosaceae in the G1 group of the Sclerotiniaceae (Monilinia species). There was no obvious topological congruence between the plant tree and groups G2 and G3 of the Sclerotiniaceae. To take into account sampling bias, we performed cophylogeny reconstructions using the full set of 263 host–parasite associations and a simplified set of 121 associations. Reconstructions were also performed independently on each of the three macro‐evolutionary regimes. CoRe‐PA classifies host–pathogen associations into (i) cospeciation, when speciation of host and pathogen occurs simultaneously, (ii) duplication, when pathogen speciation occurs independently of host speciation, (iii) sorting or loss, when a pathogen remains associated with a single descendant host species after host speciation, and (iv) host switch (also designated as host jump) when a pathogen changes host independently of speciation events (Merkle et al., 2010). The CoRe‐PA reconstructions indicated that duplications and host switch each represented ~34% of host associations from the full set of host–Sclerotiniaceae associations. Cospeciation and sorting each represented ~16% of host associations. The analysis of each macro‐evolutionary regime indicated that G1 is characterized by high duplication and low host jump frequencies and G3 by low cospeciation and sorting frequencies and high host jump frequencies (Figure 4b). We next used Procrustean superimposition in PACo (Balbuena et al., 2013) to identify host–pathogen associations that contribute significantly to cophylogeny between Sclerotiniaceae and host plant families. The PACo analysis suggested that overall, the Sclerotiniaceae lineages are not randomly associated with their host families (p < .01). The PACo analysis also includes taxon jackknifing to test for the relative contribution of each host–pathogen association in the cophylogeny pattern. In the full set of host–Sclerotiniaceae associations, ~21% contributed positively and significantly to cophylogeny, likely representing cospeciation events, while ~31% contributed negatively and significantly to cophylogeny, therefore likely representing host jump events (Table 1). Taxon jackknifing on individual macro‐evolutionary regime revealed a high frequency of associations with positive contribution to cophylogeny (likely cospeciation) in G1, and a high frequency of associations with negative contribution to cophylogeny (likely host switch) in G3, in good agreement with the CoRe‐PA analysis. In addition to the types of host–parasite association mentioned before, Jane 4 identifies “failure to diverge” associations, corresponding to events when a host speciates and the parasite remains on both new host species (Conow et al., 2010). Jane 4 analysis for the whole Sclerotiniaceae family (simplified set of associations) identifies losses as the dominant form of host association (~65%), followed by duplications (~15%) and host jumps (~12%) while failure to diverge (~5.5%) and cospeciation (~3%) was rare event (Table 1). Consistent with CoRe‐PA and PACo analyses, Jane 4 found the highest rate of host jumps in regime G3 (~42%). Unlike previous analyses, Jane 4 found the highest rate of cospeciation in G3 (~9.7%). G1 was characterized by a high duplication rate, while G2 was characterized by a high rate of losses and failure to diverge events but the lowest host jump rate. Overall, our cophylogeny analyses converged towards the conclusion that cospeciations represent a minor proportion of host associations in the Sclerotiniaceae and that the proportion of host jumps varied markedly between the three macro‐evolutionary regimes, with G3 showing the highest host jump rate (between 33% and 42%).

Figure 4

Diversification rate shifts associate with variations in rates of cospeciation, duplication and host switch in Sclerotiniaceae fungi. (a) Tanglegram depicting the associations between 102 Sclerotiniaceae species and 59 plant families. The three macro‐evolutionary regimes are indicated by coloured boxes on the Sclerotiniaceae tree. Fungal species labels are colour‐coded as in Figure 2. (b) Proportion of cospeciation, sorting/loss, duplication and host switches in host–Sclerotiniaceae associations as predicted by CoRe‐PA in 1,000 cophylogeny reconstructions. ** indicate large effect size in a macro‐evolutionary compared to the complete set of associations as assessed by Cohen's d test. (c) Proportion of host–Sclerotiniaceae associations contributing significantly and positively (likely cospeciation), non significantly and significantly and negatively (likely host switch) to cophylogeny in PACo analysis. The black dotted line indicates the percentage of broad host range species (five or more host families) in each group [Colour figure can be viewed at http://wileyonlinelibrary.com]

DISCUSSION

Our analyses lead to a model in which the extant diversity of Sclerotiniaceae fungi is the result of three macro‐evolutionary regimes characterized by distinct diversification rates and host association patterns. Patterns of cophylogeny decreased from regime G1 to G3, while the frequency of inferred host jump events increased from G1 to G3. Regime G2, which includes the highest proportion of broad host range parasites, showed a frequency of host jumps intermediate between G1 and G3. Our cophylogeny analyses are consistent with the view that long‐term plant‐pathogen cospeciation is rare (Devienne et al., 2013). The decrease in host‐pathogen cophylogeny signal from G1 to G3 regime (Figure 4) could indicate more frequent true cospeciation events in early diverging Sclerotiniaceae species or could result from host jumps being restricted to closely related hosts for G1 species whereas G2 and G3 species progressively gained the ability to jump to more divergent hosts (Devienne, Giraud, & Shykoff, 2007; Devienne et al., 2013). Low diversification rates in regimes G1 and G2 may have resulted from increased extinction rates. It is conceivable that a high extinction rate of specialist parasites during regime G2 is responsible for the reduced clade diversification rates, while increasing the apparent frequency of emergence of generalist species. This phenomenon, consistent with the view of specialization as “an evolutionary dead end” (leading to a reduced capacity to diversify), is notably supported in Tachinidae parasitic flies and hawkmoths pollinating Ruellia plants (Day et al., 2016). An analysis of Papilionoidea and Heliconii revealed lower rates of diversification for butterfly species feeding on a broad range of plants compared to specialist butterflies (Day et al., 2016; Hardy & Otto, 2014). Consistently, theoretical models of sympatric speciation predict that competition for a narrow range of resources (specialization) to be a strong driver of speciation (Dieckmann & Doebeli, 1999). Similar to plant‐feeding insects, the diversity of fungal and oomycete pathogens is considered largely driven by host jumps rather than host specialization that may follow (Choi & Thines, 2015; Hardy & Otto, 2014). Low diversification rates in G1 and G2 may also result from a strong increase in diversification rate during the transition from regime G2 to G3. Our ancestral state reconstruction analysis inferred a jump to monocots at the base of regime G3 (Figure 1). A recent study on Hesperiidae butterflies reported a strong increase in diversification coincident with a switch from dicot‐feeders to monocot‐feeders (Sahoo, Warren, Collins, & Kodandaramaiah, 2017). As postulated for Hesperiidae, the emergence of open grasslands and global temperature decrease may have affected diversification in the Sclerotiniaceae.

Competition for resources is likely to lead to speciation if host range expansion is costly (Ackermann, Doebeli, & Gomulkiewicz, 2004). Notably, if different host families are scattered in space and phylogenetically related hosts are clustered, the cost of host range expansion is expected to increase, due to higher costs for dispersal or trade‐offs with other traits (Ackermann et al., 2004). In agreement with this theory, genomic signatures associated with metabolic cost optimization were stronger in generalist than specialist fungal parasites (Badet et al., 2017). Strong climate oscillations during the Cenozoic Era, leading to population isolation through range fragmentations and dispersal events, likely contributed to the radiation of plant lineages (Nyman, Linder, Peña, Malm, & Wahlberg, 2012). For instance, major events in the diversification of Brassicaceae plant family coincide with glaciations and arid conditions of the Eocene–Oligocene and the Oligocene–Miocene transitions (Hohmann, Wolf, Lysak, & Koch, 2015). In addition, the mid‐Miocene (20–10 Mya) corresponds to the emergence of the first open grassland habitats in northern Eurasia (Strömberg, 2011). The fragmentation and diversification of plant host populations may have increased the cost for host range expansion in fungal pathogens. A jump to monocots at the base of Sclerotiniaceae regime G3 ~10 Mya may have favoured the conquest of grassland habitats and the rapid diversification of specialist species over the emergence of generalists in this group.

Similar to herbivore diet breadth (Forister et al., 2015), host range in the Sclerotiniaceae shows a continuous distribution from specialists to broad host range generalists, with the majority of species being specialists. Mathematical analyses and studies of the host range of insect herbivores suggest that host range expansion could involve the emergence of a “pre‐adaptation” followed by the colonization of new hosts (Janz & Nylin, 2008). The existence of such pre‐existing enabler traits has been proposed as a facilitator for several shifts in plant species distribution (Donoghue & Edwards, 2014). Notably, the pre‐existence of a symbiotic signalling pathway in algae is thought to have facilitated the association of land plants with symbiotic fungi (Delaux et al., 2015). In the plant genus Hypericum, increased diversification rates in the Miocene Epoch were likely facilitated by adaptation to colder climates (Nürk et al., 2015). In the bacterial pathogen Serratia marcescens and in yeast, adaptation to temperature change was associated with improved tolerance to other stresses (Caspeta & Nielsen, 2015; Ketola et al., 2013). Analysis of the complete predicted proteomes of S. borealis, S. sclerotiorum and Botrytis cinerea has revealed protein signatures often associated with cold adaptation in the secreted protein of all three species (Badet, Peyraud, & Raffaele, 2015). Cold tolerance might have been a pre‐adaptation that facilitated the emergence of generalist parasites under low diversification rates (regime G2) and a rapid diversification following host jumps on fragmented host populations (regime G3). Indeed, Sclerotinia borealis, S. glacialis, S. subarctica and S. nivalis, that diverged during regime G2, are largely restricted to hemiboreal climates (circumboreal region) and can have a lower optimal growth temperature than their sister species (Hoshino, Terami, Tkachenko, Tojo, & Matsumoto, 2010; Saito, 1997).

These findings suggest that global climate instability and host diversification in the Cenozoic might have impacted on the diversity of fungal parasites within the Sclerotiniaceae. This effect could have been direct, through the emergence of cold adaptation as an enabling trait, or indirect through changes in host population structures and host–parasite association patterns. Knowledge on the dynamics of pathogen evolution increases the understanding of the complex interplay between host, pathogen and environment governing the dynamics of disease epidemics. These evolutionary principles are useful for the design of disease management strategies (Vander Wal et al., 2014) and provide new insights into the factors that influenced the diversity of extant fungal parasite species.

DATA ACCESSIBILITY

Species identifiers: GenBank accessions provided in Table S1

DNA sequences and alignments: provided as File S2

Phylogenetic trees: provided as Files S1, S3, S4 and S5

Host range data, DNA sequences and alignments and Phylogenetic trees are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.7cs3g

AUTHOR CONTRIBUTIONS

O.N. collected data, performed analyses, wrote the original draft and revised the manuscript; A.B. performed analyses, wrote and revised the manuscript; A.T. collected data and revised the manuscript; J.P.C. involved in funding acquisition, data collection, revision of manuscript; S.R. involved in supervision, funding acquisition, project administration, writing original draft and revision of manuscript.

Click here for additional data file.

Click here for additional data file.