Evolutionary dynamics of mycorrhizal symbiosis in land plant diversification

Mycorrhizal symbiosis between soil fungi and land plants is one of the most widespread and ecologically important mutualisms on earth. It has long been hypothesized that the Glomeromycotina, the mycorrhizal symbionts of the majority of plants, facilitated colonization of land by plants in the Ordovician. This view was recently challenged by the discovery of mycorrhiza-like associations with Mucoromycotina in several early diverging lineages of land plants. Utilizing a large, species-level database of plants’ mycorrhiza-like associations and a Bayesian approach to state transition dynamics we here show that the recruitment of Mucoromycotina is the best supported transition from a non-mycorrhizal state. We further found that transitions between different combinations of either or both of Mucoromycotina and Glomeromycotina occur at high rates, and found similar promiscuity among combinations that include either or both of Glomeromycotina and Ascomycota with a nearly fixed association with Basidiomycota. Our results portray an evolutionary scenario of evolution of mycorrhizal symbiosis with a prominent role for Mucoromycotina in the early stages of land plant diversification.

facilitated terrestrialisation 16 , or that early land plants formed dual Mucoromycotina-Glomeromycotina partnerships [18][19][20][21] . After this discovery, Rhynie Chert fossils where re-evaluated, revealing mycorrhiza-like colonizations resembling both Glomeromycotina and Mucoromycotina 11 . Moreover, mycorrhiza-formation genes from Mucoromycotina-associated liverworts recover the Glomeromycotina-associated phenotype in a transformed mutant of the angiosperm Medicago truncatula, which reveals that the genes required for symbiosis have been conserved among liverworts that associate exclusively with Mucoromycotina as well as higher plants that associate exclusively with Glomeromycotina 13,20 .
Given that Ascomycota, Basidiomycota, Glomeromycotina, and Mucoromycotina likely diverged prior to the divergence of land plants 22, 23 , it is possible to treat different combinations of symbiotic association with these phyla as categorical character states on the plant phylogeny and analyse transition dynamics between the states in a Bayesian phylogenetic comparative context. Considering the uncertainty of the evolutionary relationships of early Embryophytes 24,25 , we assessed the probability of all possible combinations of mycorrhizal associations for the most recent common ancestor of land plants.

Results
We obtained a dataset of 732 species of land plants for which the mycorrhizal fungi have been identified with molecular methods. 45 species were added to represent non-mycorrhizal lineages. We used the plant chloroplast DNA markers psbA, rbcL and rps4 to infer phylogenetic relationships between these species. Our estimates of phylogeny correspond well with the prevailing understanding of the systematics of the land plants at least so far as the monophyly of major groups and the relative branching order of these groups under the different rooting scenarios are concerned 26,27 .
Optimising the observed repertoires of mycorrhizal association as transitioning categorical states on our phylogenetic estimates resulted in a general pattern of phylogenetic conservatism: major plant groups associate quite uniformly with major fungal groups (Fig. 1). Our ancestral state reconstructions recover strong support for the presence of mycorrhiza-like association for the most recent common ancestor of the land plants. However, the particular state for the root was equivocal, showing comparable levels of support for i) an association just with Mucoromycotina, ii) a repertoire comprising both Mucoromycotina and Glomeromycotina; and iii) no mycorrhizal association at all ( Figure S1). The relative levels of support, and the inclusion of additionally supported root states, were influenced by different rooting scenarios ( Figure S2).
The pattern of transitions among different repertoires of symbiotic association suggests two main paths along which individual associations within a larger repertoire are gained and lost relatively promiscuously (Fig. 2). The first of these paths involves Mucoromycotina and Glomeromycotina: the association with Glomeromycotina is added to, and subtracted from, the association with Mucoromycotina at relatively high instantaneous transition rates. The association with Mucoromycotina within a repertoire that spans both is also lost at relatively high rates, but gained at much lower rates, suggesting that the association with Glomeromycotina is relatively more facultative within this repertoire. The second path includes gains and losses of Ascomycota, and losses of Glomeromycotina (but gains less so), at high rates within repertoires in which the association with Basiodiomycota appears more obligate.
Explicit hypothesis testing to quantify which transition away from a state of no mycorrhizal association is best supported prefers Mucoromycotina under all four rooting scenarios: in three out of four, the Bayes Factor (BF) was larger than 10, interpreted as strong support, in the fourth scenario (hornworts sister to all other land plants) the BF was ~8.35, which is generally interpreted as substantial support 28 (Table S1). Placing the evolution of mycorrhizal associations on a temporal axis in a sliding window analysis (Fig. 3) shows Mucoromycotina and Glomeromycotina dominating early associations, while associations that include Basidiomycota and Ascomycota become more pervasive later in land plant evolution.

Discussion
For each evaluated scenario of land plant evolution, our results support the hypothesis that the most recent common ancestor of land plants was involved in symbiotic interactions with fungi. This result is in accordance with evidence from the fossil record 11 and genomics [12][13][14] . For the small, rootless, leafless plants with rhizoid-based absorbing systems that started colonizing the land, the alliance with fungi is hypothesized to have been essential in overcoming major issues of nutrient and water limitation in the absence of existing soils 29,30 . Our analyses suggest that the fungal associates of these earliest land plants most likely included Mucoromycotina, and not exclusively Glomeromycotina, as commonly assumed 1,31,32 . An exclusive association with Mucoromycotina for the root of the land plants received the highest support of all possible mycorrhizal repertoires, for all hypotheses of the relationships between the main land plant lineages. Furthermore, our hypothesis tests supported Mucoromycotina over Glomeromycotina as the initial gain for the most recent common ancestor of the land plants. However, our reconstructions also suggest that a repertoire comprising both Mucoromycotina and Glomeromycotina cannot be ruled out, and we find high rates for transitions in which Glomeromycotina are gained and lost in combination with Mucoromycotina (Fig. 2), suggesting a versatile scenario for the evolution of association with both groups. Mucoromycotina have been recorded in the rhizoids and roots of extant liverworts 16 , hornworts 18 , lycophytes 21 , ferns 21 , gymnosperms 33,34 , and angiosperms 35 , but within early diverging land plant lineages (except for the liverwort lineage Haplomitriopsida) 16 they were mostly found simultaneously with Glomeromycotina 16,18,21 . The association with both fungal lineages was likely also present in the Devonian fossil plant Horneophyton ligneri 11 , and Field et al. 19 speculated that the ability to associate with more than one fungal partner was an ancient strategy that allowed the earliest land plants to occupy highly heterogeneous and dynamic environments. However, this plasticity appears to not be maintained: once association with Mucoromycotina is lost, reversals occur at a low rate (Fig. 2)  From the prevalent association with strictly Glomeromycotina, there have been multiple independent evolutionary shifts towards Ascomycota and Basidiomycota, leading to increasingly prevalent reconstruction of these interactions over the course of plant diversification (Fig. 3). Our results suggest that these transitions started with a gain of Basidiomycota, rather than Ascomycota (Fig. 2). Subsequent gains of Ascomycota and losses of Glomeromycotina occur at high rates, leading to various association repertoires that include either or both Ascomycota and Basidiomycota. These repertoires are present in several extant land plant lineages and represent the ectomycorrhizal, orchid mycorrhiza, and ericoid mycorrhizal types 9,17 . The ability to recruit saprotrophic lineages of wood and litter decaying fungi from among Ascomycota and Basidiomycota into novel symbioses was likely instrumental for plant adaptation to various ecological challenges 5 . For example, for Orchidaceae, the most species-rich lineage of non-arbuscular mycorrhizal plants, the transition from associations with Glomeromycotina to Ascomycota and Basidiomycota is linked to niche expansions and radiations, which in synchrony with the development of specialized pollination syndromes has promoted speciation in the largest family of plants on earth 38,39 . Similarly, the independent evolution of ericoid mycorrhiza in Diapensiaceae and Ericaceae, estimated to date back to the Cretaceous 40,41 , is a potential adaptation to nutrient poor, acidic soils 31 . Also, transitions to ectomycorrhiza independently evolved in various gymnosperm (e.g. Pinaceae, Gnetum, Taxus) and angiosperm lineages (e.g. Nyctaginaceae, Polygonaceae, Myrtaceae, Malvales, Malpighiales, Fabaceae, Fagales; Fig. 1). Parallel to the latter, a shift towards fungi involved in the ectomycorrhizal and ericoid symbiosis has also occurred in liverworts (Fig. 1). Although relatively few plant species -mostly trees and shrubs -are ectomycorrhizal, the worldwide importance of the ectomycorrhizal association is considerable, due to its dominance in temperate and boreal forests, and in tropical rainforests in Southeast Asia 17 . Ectomycorrhizal symbioses likely emerged in semi-arid forests dominated by conifers under tropical to subtropical climates and diversified in angiosperms and conifer forests driven by a change to cooler climate during the Cenozoic 42,43 . Loss of mycorrhizal symbiosis has occurred from all single association states, mostly at relatively low transition rates (Fig. 3). These transitions are explained by plant adaptations to either nutrient-rich or extremely nutrient-poor soils, for which the benefits of the symbiosis do not outweigh its costs 44 . However, transition rates towards the non-mycorrhizal state may have been underestimated here, since several non-mycorrhizal angiosperm lineages (all with a recent evolutionary origin 45 ) have not been included. A notable increase in the proportion of non-mycorrhzial lineages around 450-400 mya is caused by the origin of mosses and the diversification of non-mycorrhizal liverworts. Similar to Maherali et al. 45 we reconstructed a regain of symbiosis from a non-mycorrhizal ancestor for a few lineages. Because it is not known whether the mycorrhizal symbiosis can be recovered after loss, it is possible that this pathway may not occur in nature.
Our results portray an evolutionary scenario of evolution of mycorrhizal symbiosis with a prominent role for Mucoromycotina in the early stages of land plant diversification. In most plant lineages, Glomeromycotina, the dominant mycorrhizal symbionts of extant land plants, subsequently replaced Mucoromycotina. Later on, several  transitions from Glomeromycotina to various Ascomycota and Basidiomycota lineages have occurred, establishing novel mycorrhizal syndromes, such as orchid, ericoid, and ectomycorrhizas. Our findings demonstrate the importance of Mucoromycotina fungi for our understanding of the early evolution of the mycorrhizal symbiosis. We still know very little about the biology of symbiotic Mucoromycotina 36 , but experimental evidence suggests they form mycorrhizas that are physiologically and functionally different from symbioses with Glomeromycota 20 . Their presence as symbiotic fungi in land plants has been overlooked until recently 16,20 , and it is likely that further screening of land plants will reveal that many more plant taxa are associated with Mucoromycotina.

Data collection.
To compile a dataset of plants and their symbiotic fungi, we searched the NCBI Nucleotide databank for records of Glomeromycotina (at the time of the search 'Glomeromycota'), Mucoromycotina, Ascomycota and Basidiomycota that had annotations recording the plant host species ( Figure S3). Subsequently, for each of these plant host species we conducted a GenBank search and reduced our dataset to all records with an rbcL sequence available for the plant host. For the remaining records, we verified mycorrhizal status through literature study and discarded all unconfirmed records from the dataset. We then performed a literature search for plant orders that were not in the dataset as well as for early diverging lineages of land plants. Because it is difficult to discriminate among symbioses formed by Glomeromycotina and Mucoromycotina by morphological observations, we only included mycorrhizal associations based on DNA identification for these fungi. For lycopods, polypod ferns, hornworts and liverworts, species that were not found to harbour mycorrhiza-like associations during literature surveys were classified as non-mycorrhizal, although this could be a sampling artefact for some species. Furthermore, mosses and Nymphaea alba were included to represent major non-mycorrhizal lineages. The final dataset covers 732 plant species distributed over 78 plant orders. The dataset includes 24 hornworts (11% of extant species diversity 46 ), 7 mosses (0.06%) 46 , 76 liverworts (0.84%) 46 , 518 angiosperms (0.18%) 46 , 73 gymnosperms (6.77%) 46 , 16 lycopods (1.24%) 46 , and 18 polypod ferns (0.17%) 46 . For these plants species, we found associations with 150 Ascomycota, 305 Basidiomycota, 385 Glomeromycotina, 28 Mucoromycotina and 45 non-mycorrhizal species (Table S2).
DNA sequence data of the plants, including members of the mosses, were obtained from GenBank to reconstruct phylogenetic relationships. For liverworts, hornworts, polypod ferns, and lycopods, we added several species to the dataset to increase taxon sampling, resulting in a total of 759 species for phylogenetic analysis. For 146 species, full or partial chloroplast genomes were available, which we used to extract sequences for psbA, rbcL and rps4. For other species, rps4 and psbA sequences were downloaded, where possible, to supplement the rbcL dataset. Accession numbers are listed in the supplementary data (Table S3).

Phylogenetic analysis and divergence dating.
For each marker, we aligned the sequences with MAFFT v.7 47 using the FFT-NS-i Iterative refinement method, and then selected the substitution model with jModel-Test 2.1.10 48,49 . For each marker, 3 substitution schemes where tested on a neighbor-joining topology, including models with unequal base frequencies, rate heterogeneity, and a proportion of invariant sites. The GTR + I + γ model was selected for all partitions using the AIC. We performed divergence dating with BEAST2 v2.3.2 50 using four fossil calibration points and one age estimate from literature for the crown node of liverworts to date the phylogeny. We selected a uniform distribution for each of the calibration points using the minimum and maximum estimates for these nodes from literature (Table S4). We chose a Yule prior with a uniform birth rate for the analysis, a lognormal relaxed clock model, and estimated the clock rate. We applied the GTR substitution model with a Gamma category count of 4 and estimated shape parameter value of 1.0. The proportion of invariant sites was estimated (initial value 0.01) and the mean substitution rate fixed. We selected an exponential distribution for the prior on the mean substitution rate. To test the effect of different phylogenetic hypotheses 24 Figure S1). The first three hypotheses have recently been proposed as the best-supported explanations using a large transcriptomic dataset 27 . The fourth hypothesis has been traditionally regarded as the most accurate representation of early land plant diversification (e.g. Field et al. 19 ). During the MCMC analyses, trace files were updated every 1000 generations, and trees sampled every 10,000 generations, until the effective sample size of major traced parameters exceeded 200 (and all others exceeded 100) using a burn-in of 100 * 10 6 generations. We thus terminated the runs after, respectively, 374,735,000 generations for ABasal; 354,497,000 generations for ATxMB; 345,720,000 for MBasal; and 374,254,000 for TBasal. We then constructed the maximum clade credibility tree using Tree Annotator v2.2.1.

Comparative analysis and hypothesis tests.
In our analysis we assume that the four major fungal groups of which members participate in symbiotic associations were already in existence prior to the diversification of land plants 23 . Therefore, we treat each distinct repertoire of associations that land plants form with members of these groups as a discrete state whose evolutionary transition dynamics we modelled subsequent to two additional assumptions. First, because there are qualitative differences between the types of symbiotic associations that are formed with some of the different fungal groups (e.g. intracellular versus ectomycorrhizal association), we assumed that the evolutionary adaptations required to enable such associations are not gained (or lost) instantaneously. Hence, we disallowed state shifts that implied multiple, simultaneous gains and losses such that, for example, a change from a state representing a repertoire confined to Glomeromycotina to one confined to Mucoromycotina has to pass through an intermediate state where the repertoire is broadened to include both groups. Second, because the respective adaptations that enable different types of mycorrhiza-like association are likely subject to evolutionary trade-offs such that repertoires of associations cannot expand infinitely we limited any intermediate states to those we observe in nature. For example, simultaneous association with both Glomeromycotina and Mucoromycotina does occur in our dataset of extant taxa, but complete generalism that includes all fungal groups in a single repertoire does not, which is why we allowed the former, but not the latter, as possible ancestral states.
A convenient side effect of these assumptions was that this limited the number of free parameters in the state transition (Q) matrix, which otherwise would have undergone a combinatorial explosion had we included all possible permutations in the repertoires of mycorrhiza-like association as distinct states, which would have impeded convergence in our analyses. To mitigate such proliferation of potentially unneeded, free parameters further, we performed our analyses using Reversible-Jump MCMC, as implemented in BayesTraits's 'multistate' analysis mode. We ran each of our analyses in triplicate for 10 6 generations, as initial experimentation had demonstrated reasonable convergence in our data under these settings. In cases where we required estimates of marginal likelihoods, i.e. for hypothesis testing by Bayes factor analysis, we approximated these using a stepping stone sampler that we ran for 100 stones, with 200,000 iterations per stone.
Using this approach, we reconstructed the ancestral states for the four different rootings of our phylogeny. However, although such analyses result in estimates for the posterior distribution of states at any given node (such as the root), they do not necessarily provide the false certainty on which to base a single, unambiguous scenario for the order in which symbiotic associations are acquired, especially not when multiple states are reconstructed with similarly large posterior probabilities at deep nodes (as was the case). Given the number of fungal groups and the differences and similarities among these with respect to the types of mycorrhizal associations they participate in, we expected there to be distinct paths along which repertoires of association have evolved. Interrogation and visualisation of the Q matrix showed that, broadly, two such paths appear to exist: one where various permutations of association with Glomeromycotina and Mucoromycotina are gained and lost, and another that traverses Ascomycota, Basidiomycota in addition to Glomeromycotina. However, which of these paths was taken first was not yet evident.
We therefore constructed explicit hypothesis tests to distinguish between various plausible scenarios. To do so, in addition to the assumptions affecting the Q matrix outlined above, we further constrained our analyses to require the absence of any mycorrhizal association on the root node, and then tested which initial gain was best supported by the data. To quantify this, we estimated the marginal likelihood of the model where the root is constrained to have no association but without any additional constraints on the order in which subsequent associations are acquired (beyond the general assumptions already discussed), and compared this with models where, respectively, each of the initial gains of a single fungal group is disallowed. The logic here is that disallowing the initial shift that best fits the data will result in the marginal likelihood that differs most significantly from the less-constrained model.
Lastly, to place the expansion of repertoires of symbiotic association on a temporal axis, we placed the ancestral state reconstructions for the scenario where the root node has no mycorrhizal association in bins of 50 Myr to visualise these in a states-through-time plot (Fig. 3). All data and scripts are available (https://doi.org/10.5281/zenodo.1037586).
Data availability. The datasets generated during and/or analysed during the current study are available in the GitHub repository, https://doi.org/10.5281/zenodo.1037586.