domain: sdv sde

text:
Benoît Perez-Lamarque

Maarja Öpik

Odile Maliet

Ana C Afonso

Marc- André Selosse

Florent Martos

Hélène Morlon

Analyzing diversification dynamics using barcoding data: the case of an obligate mycorrhizal symbiont

Keywords: microbial diversification, arbuscular mycorrhiza, obligate symbiosis, ecological niche, macroevolution, fungi

Analyzing diversification dynamics is key to understanding the past evolutionary history of clades that led to present-day biodiversity patterns. While such analyses are widespread in well-characterized groups of species, they are much more challenging in groups which diversity is mostly known through molecular techniques. Here, we use the largest global database on the small subunit (SSU) rRNA gene of Glomeromycotina, a subphylum of microscopic arbuscular mycorrhizal fungi that provide mineral nutrients to most land plants by forming one of the oldest terrestrial symbioses, to analyze the diversification dynamics of this clade in the past 500 million years (Myr). We perform a range of sensitivity analyses and simulations to control for potential biases linked to the nature of the data. We find that Glomeromycotina tend to have low speciation rates compared to other eukaryotes. After a peak of speciations between 200 and 100 Myr ago, they experienced an important decline in speciation rates toward the present. Such a decline could be at least partially related to a shrinking of their mycorrhizal niches and to their limited ability to colonize new niches. Our analyses identify patterns of diversification in a group of obligate symbionts of major ecological and evolutionary importance and illustrate that short molecular markers combined with intensive sensitivity analyses can be useful for studying diversification dynamics in microbial groups.

Introduction:

Understanding past dynamics of speciation and extinction, as well as the abiotic and biotic factors that modulate the frequency of speciation and extinction events (i.e. diversification rates) is key to understanding the historical processes that shaped presentday biodiversity patterns [START_REF] Barnosky | Distinguishing the effects of the red queen and court jester on miocene mammal evolution in the northern rocky mountains[END_REF][START_REF] Condamine | Testing the role of the Red Queen and Court Jester as drivers of the macroevolution of Apollo butterflies[END_REF][START_REF] Morlon | Phylogenetic approaches for studying diversification[END_REF][START_REF] Varga | Megaphylogeny resolves global patterns of mushroom evolution[END_REF] [START_REF] Barnosky | Distinguishing the effects of the red queen and court jester on miocene mammal evolution in the northern rocky mountains[END_REF][START_REF] Benton | The Red Queen and the Court Jester: Species diversity and the role of biotic and abiotic factors through time[END_REF][START_REF] Chomicki | The evolution of mutualistic dependence[END_REF][START_REF] Clarke | Climate, energy and diversity[END_REF]. While phylogenetic analyses of diversification are widespread in well-characterized groups of species, such as animals and plants [START_REF] Givnish | Orchid phylogenomics and multiple drivers of their extraordinary diversification[END_REF][START_REF] Magallón | Absolute diversification rates in angiosperm clades[END_REF][START_REF] Rolland | Faster speciation and reduced extinction in the tropics contribute to the mammalian latitudinal diversity gradient[END_REF][START_REF] Upham | Inferring the mammal tree: Species-level sets of phylogenies for questions in ecology, evolution, and conservation[END_REF], they are much more challenging in groups which diversity is mostly known through environmental DNA sequences and molecular techniques. In particular, the characterization of poorly cultivable microbial groups such as most bacteria and fungi is often limited to metabarcoding techniques, which consist in the specific amplification and sequencing of a short DNA region [START_REF] Taberlet | DNA amplification and multiplexing[END_REF]. One the one hand, these data often render species delineation, phylogenetic reconstruction, and the estimation of global scale diversity highly uncertain, which all affect the phylogenetic inference of diversification dynamics [START_REF] Lekberg | More bang for the buck? Can arbuscular mycorrhizal fungal communities be characterized adequately alongside other fungi using general fungal primers?[END_REF][START_REF] Moen | Why does diversification slow down?[END_REF]. On the other hand, it is possible to assess the robustness of phylogenetic diversification analyses to data uncertainty. Given the current limitations of sequencing technologies and the nature of the molecular data available for most microbial groups, using metabarcoding data and performing thorough robustness analyses is one of the only (if not the only) possible approach to analyze their diversification dynamics [START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF][START_REF] Lewitus | Clade-specific diversification dynamics of marine diatoms since the Jurassic[END_REF][START_REF] Louca | Bacterial diversification through geological time[END_REF].

Here we analyze the diversification dynamics of arbuscular mycorrhizal fungi from the subphylum Glomeromycotina. These fungi are obligate symbionts that have been referred to as an "evolutionary cul-de-sac, albeit an enormously successful one" [START_REF] Malloch | The evolution of mycorrhizae[END_REF][START_REF] Morton | Species and clones of arbuscular mycorrhizal fungi (Glomales, Zygomycetes): their role in macro-and microevolutionary processes[END_REF]. This alludes to their ecological success despite limited morphological and species diversities: they associate with the roots of >80% of land plants, where they provide mineral resources in exchange for photosynthates [START_REF] Smith | Mycorrhizal Symbiosis[END_REF]. Present in most terrestrial ecosystems, Glomeromycotina play key roles in plant protection, nutrient cycling, and ecosystem processes (van der Heijden, Martin, [START_REF] Van Der Heijden | Mycorrhizal ecology and evolution: the past, the present, and the future[END_REF]. Fossil evidence and molecular phylogenies suggest that Glomeromycotina contributed to the emergence of land plants [START_REF] Feijen | Evolutionary dynamics of mycorrhizal symbiosis in land plant diversification[END_REF][START_REF] Field | Symbiotic options for the conquest of land[END_REF][START_REF] Selosse | The land flora: a phototroph-fungus partnership?[END_REF][START_REF] Strullu-Derrien | The origin and evolution of mycorrhizal symbioses: from palaeomycology to phylogenomics[END_REF] and coevolved with them for more than 400 million years (Myr) [START_REF] Lutzoni | Contemporaneous radiations of fungi and plants linked to symbiosis[END_REF][START_REF] Simon | Origin and diversification of endomycorrhizal fungi and coincidence with vascular land plants[END_REF][START_REF] Strullu-Derrien | The origin and evolution of mycorrhizal symbioses: from palaeomycology to phylogenomics[END_REF].

Glomeromycotina are microscopic soil-and root-dwelling fungi that are hard to differentiate based on morphology and difficult to cultivate without host plant. Although their classical taxonomy is mostly based on the characters of spores and root colonization [START_REF] Smith | Mycorrhizal Symbiosis[END_REF][START_REF] Stürmer | A history of the taxonomy and systematics of arbuscular mycorrhizal fungi belonging to the phylum Glomeromycota[END_REF], Glomeromycotina species delineation has greatly benefited from DNA sequencing [START_REF] Krüger | Phylogenetic reference data for systematics and phylotaxonomy of arbuscular mycorrhizal fungi from phylum to species level[END_REF].

Experts have defined "virtual taxa" (VT) based on a minimal 97% similarity of a region of the 18S small subunit (SSU) rRNA gene and monophyly criteria [START_REF] Öpik | DNA-based detection and identification of Glomeromycota: the virtual taxonomy of environmental sequences[END_REF][START_REF] Öpik | The online database MaarjAM reveals global and ecosystemic distribution patterns in arbuscular mycorrhizal fungi (Glomeromycota)[END_REF]. As for many other pragmatic species delineation criteria, VT have rarely been tested for their biological relevance [START_REF] Powell | Evolutionary criteria outperform operational approaches in producing ecologically relevant fungal species inventories[END_REF], and a consensual system of Glomeromycotina classification is still lacking [START_REF] Bruns | Glomeromycotina: what is a species and why should we care?[END_REF]. Besides the rDNA region, Glomeromycotina remain poorly known genetically: other gene sequences are available for only a few species [START_REF] James | Reconstructing the early evolution of Fungi using a six-gene phylogeny[END_REF][START_REF] Lutzoni | Contemporaneous radiations of fungi and plants linked to symbiosis[END_REF] and less than 30 complete genomes are currently available [START_REF] Venice | At the nexus of three kingdoms: the genome of the mycorrhizal fungus Gigaspora margarita provides insights into plant, endobacterial and fungal interactions[END_REF].

Hence, despite the ecological ubiquity and evolutionary importance of Glomeromycotina, large-scale patterns of their diversification dynamics, as well as the factors that correlate with these dynamics, remain poorly known. A previous dated phylogenetic tree of VT found that many speciation events occurred after the last major continental reconfiguration around 100 Myr ago [START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF], suggesting the radiation of Glomeromycotina is not linked to vicariant speciation during this geological event. Indeed, vicariant speciation might only play a minor role in Glomeromycotina diversification, as these organisms have spores that disperse efficiently, promoting gene flow [START_REF] Bueno | How do arbuscular mycorrhizal fungi travel?[END_REF][START_REF] Correia | First evidence for the joint dispersal of mycorrhizal fungi and plant diaspores by birds[END_REF][START_REF] Egan | Detection of arbuscular mycorrhizal fungal spores in the air across different biomes and ecoregions[END_REF]. Based on the diversity and abundance of Glomeromycotina in tropical grasslands [START_REF] Read | Mycorrhizas in ecosystems[END_REF], it has been suggested (but never tested) that these habitats are diversification hotspots for Glomeromycotina [START_REF] Pärtel | Historical biome distribution and recent human disturbance shape the diversity of arbuscular mycorrhizal fungi[END_REF]. In this case, the pace of Glomeromycotina diversification through time could be tightly linked to changes in the total area of tropical grasslands. Finally, Glomeromycotina are currently obligate symbionts and their evolutionary history could thus have been largely influenced by their interactions with their host plants [START_REF] Lutzoni | Contemporaneous radiations of fungi and plants linked to symbiosis[END_REF][START_REF] Sauquet | Key questions and challenges in angiosperm macroevolution[END_REF][START_REF] Zanne | Three keys to the radiation of angiosperms into freezing environments[END_REF]. Over the last 400 Myr, land plants have experienced massive extinctions and radiations [START_REF] Cleal | Composition and dynamics of the great Phanerozoic Evolutionary Floras[END_REF][START_REF] Zanne | Three keys to the radiation of angiosperms into freezing environments[END_REF], adaptations to various ecosystems [START_REF] Bredenkamp | On the origin of northern and southern hemisphere grasslands[END_REF][START_REF] Brundrett | Evolutionary history of mycorrhizal symbioses and global host plant diversity[END_REF], and associations with different soil microorganisms [START_REF] Werner | Symbiont switching and alternative resource acquisition strategies drive mutualism breakdown[END_REF][START_REF] Werner | A single evolutionary innovation drives the deep evolution of symbiotic N2-fixation in angiosperms[END_REF]. All these events could have influenced the diversification dynamics of Glomeromycotina, although their relative generalism (Perez-Lamarque, [START_REF] Perez-Lamarque | Cheating in arbuscular mycorrhizal mutualism: a network and phylogenetic analysis of mycoheterotrophy[END_REF][START_REF] Sanders | Preference, specificity and cheating in the arbuscular mycorrhizal symbiosis[END_REF][START_REF] Van Der Heijden | Mycorrhizal ecology and evolution: the past, the present, and the future[END_REF] could buffer this influence.

We aim to characterize the pace of Glomeromycotina diversification in the last 500 Myr and to test the association between diversification rates and a variety of biotic and abiotic factors. We begin by reconstructing several thoroughly sampled phylogenetic trees of Glomeromycotina, considering several criteria of species delineations and uncertainty in phylogenetic reconstructions. We combine this phylogenetic data with paleoenvironmental data and data of current Glomeromycotina geographic distributions, ecological traits, interaction with host plants, and genetic diversity. Finally, we apply a series of birth-death models of cladogenesis to answer specific questions and test hypotheses related to Glomeromycotina diversification: (i) how often do speciation events occur? (ii) were speciation rates relatively constant, or were they higher during specific periods of evolutionary history? and do speciation rates decline through time, as observed for many macroorganisms [START_REF] Moen | Why does diversification slow down?[END_REF]? (iii) are speciation rates positively correlated with past temperature, CO2 concentration, and/or land plant diversity? (iv) are present-day speciation rates correlated with geographic distribution, spore size (itself often inversely related to dispersal capacity, [START_REF] Nathan | Mechanisms of long-distance seed dispersal[END_REF], degree of specialization toward plant species, and genetic diversity? For each of these questions, we thoroughly assess the robustness of our results to uncertainty in the data.

Material & methods:

Virtual taxa phylogenetic reconstruction:

We downloaded the Glomeromycotina SSU rRNA gene sequences from MaarjAM, the largest global database of Glomeromycotina gene sequences updated in June 2019 (Öpik et al., 2010). We reconstructed several Bayesian phylogenetic trees of the 384 virtual taxa (VT) from the corresponding representative sequences available in the MaarjAM database (Supplementary Methods 1). We used the full length (1,700 base pairs) SSU rRNA gene sequences from [START_REF] Rimington | Ancient plants with ancient fungi: liverworts associate with early-diverging arbuscular mycorrhizal fungi[END_REF] to better align the VT sequences using MAFFT [START_REF] Katoh | MAFFT Multiple sequence alignment software version 7: Improvements in performance and usability[END_REF]. We selected the 520 base pair central variable region of the VT aligned sequences and performed a Bayesian phylogenetic reconstruction using BEAST2 [START_REF] Bouckaert | BEAST 2: A software platform for Bayesian evolutionary analysis[END_REF]. We set the crown root age at 505 Myr [START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF], which is coherent with fossil data and previous dated molecular phylogenies [START_REF] Lutzoni | Contemporaneous radiations of fungi and plants linked to symbiosis[END_REF][START_REF] Strullu-Derrien | The origin and evolution of mycorrhizal symbioses: from palaeomycology to phylogenomics[END_REF]. We also used the youngest (437 Myr) and oldest (530 Myr) crown age estimates from [START_REF] Lutzoni | Contemporaneous radiations of fungi and plants linked to symbiosis[END_REF] in diversification analyses that may be particularly sensitive to absolute dates.

Delineation into Evolutionary Units (EUs):

We considered several ways to delineate Glomeromycotina species based on the SSU rRNA gene. In addition to the VT species proxy, we delineated Glomeromycotina de novo into evolutionary units (EUs) using a monophyly criterion and 5 different thresholds of sequence similarity ranging from 97 to 99%. We gathered Glomeromycotina sequences of the SSU rRNA gene from MaarjAM, mainly amplified by the primer pair NS31-AML2 (variable region) [START_REF] Lee | Improved PCR primers for the detection and identification of arbuscular mycorrhizal fungi[END_REF][START_REF] Simon | Specific amplification of 18S fungal ribosomal genes from vesicular-arbuscular endomycorrhizal fungi colonizing roots[END_REF] (dataset 1, Supplementary Table 1). There were 36,411 sequences corresponding to 27,728 haplotypes.

We first built a phylogenetic tree of these haplotypes and then applied to this tree our own algorithm of EU delineation (R-package RPANDA [START_REF] Morlon | RPANDA: An R package for macroevolutionary analyses on phylogenetic trees[END_REF][START_REF] Core | R: A language and environment for statistical computing[END_REF]) that traverses the tree from the root to the tips, at every node computes the average similarity of all sequences descending from the node, and collapses the sequences into a single EU if their sequence dissimilarity is lower than a given threshold (Supplementary Methods 2). In other words, Glomeromycotina sequences are merged into the same EU if they form a monophyletic clade and if they are on average more similar than the sequence similarity threshold. Finally, we performed Bayesian phylogenetic reconstructions of the EUs using BEAST2, using the same crown ages as above (Supplementary Methods 1).

Coalescent-based species delineation analyses:

Finally, we considered the Generalized Mixed Yule Coalescent method (GMYC) [START_REF] Fujisawa | Delimiting species using single-locus data and the generalized mixed yule coalescent approach: A revised method and evaluation on simulated data sets[END_REF][START_REF] Pons | Sequence-based species delimitation for the DNA taxonomy of undescribed insects[END_REF], a species delineation approach that does not require specifying an arbitrary similarity threshold. GMYC estimates the time t in a reconstructed calibrated tree that separates species diversification (Yule process -before t) and intraspecific differentiation (coalescent process -after t). GMYC is too computationally intensive to be applied to the 36,411 SSU sequences; we used it here on three clades of manageable size (the family Claroideoglomeraceae; the order Diversisporales; and an early-diverging clade composed of the orders Archaeosporales and Paraglomerales) to (i) investigate whether the SSU gene evolves fast enough to accumulate substitutions between Glomeromycotina speciation events [START_REF] Bruns | Glomeromycotina: what is a species and why should we care?[END_REF] and (ii) evaluate the biological relevance of the VT and various EUs delineations. For each clade, we reconstructed Bayesian phylogenetic trees of haplotypes (Supplementary Methods 1). We then ran GMYC analyses (splits R-package [START_REF] Ezard | SPLITS: SPecies' LImits by Threshold Statistics[END_REF]) on each of these trees and evaluated the support of the GMYC model compared to a null model in which all tips are assumed to be different species, using a likelihood ratio test (LRT). If the LRT supports the GMYC model, different SSU haplotypes belong to the same Glomeromycotina species, i.e. the SSU rRNA gene has time to accumulate substitutions between Glomeromycotina speciation events.

Total diversity estimates:

We evaluated how thoroughly sampled our species-level Glomeromycotina phylogenetic trees are by estimating the total number of VT and EUs using rarefaction curves and the Bayesian Diversity Estimation Software (BDES [START_REF] Quince | The rational exploration of microbial diversity[END_REF])) (Supplementary Methods 3). The BDES estimates the total number of species by extrapolating a sampled taxa abundance distribution at global scale [START_REF] Quince | The rational exploration of microbial diversity[END_REF].

Additional molecular markers:

We explored the possibility to carry some of our analyses using two other molecular markers: the large subunit (LSU) rRNA gene and the ITS2 region. We downloaded the Glomeromycotina LSU database of [START_REF] Delavaux | Utility of large subunit for environmental sequencing of arbuscular mycorrhizal fungi: a new reference database and pipeline[END_REF] as well as the LSU sequences available in MaarjAM. We obtained a total 2,044 sequences that we aligned using MAFFT and TrimAl. We retained the 1,760 unique haplotypes, reconstructed the phylogenetic tree of the LSU sequences using BEAST2 and used the resulting calibrated tree to delineate Glomeromycotina LSU units with the GMYC model (same pipeline as above). We similarly downloaded the Glomeromycotina ITS dataset of [START_REF] Lekberg | More bang for the buck? Can arbuscular mycorrhizal fungal communities be characterized adequately alongside other fungi using general fungal primers?[END_REF]. We tried to align them but confirmed that the ITS sequences of Glomeromycotina are very difficult to align, making them unsuitable for phylogenetic reconstruction and subsequent diversification analyses (Supplementary Fig. 1).

Diversification analyses:

Unless specified differently, our diversification analyses were performed using the SSU rRNA gene. In order to account for various sources of uncertainties in the SSU rRNA data, we replicated all our diversification analyses across different species delineations, phylogenetic reconstructions and dating, and total diversity estimates. For each species delineation criterion, we obtained a consensus tree and selected 12 trees equally spaced in 4 independent Bayesian chains, hereafter referred to as the replicate trees. When the 12 trees were not sufficient to conclude, we used 100 replicate trees.

We estimated lineage-specific speciation rates using ClaDS, a Bayesian diversification model that accounts for speciation rate heterogeneity by modeling small rate shifts at speciation events [START_REF] Maliet | A model with many small shifts for estimating species-specific diversification rates[END_REF]. At each speciation event, the descending lineages inherit new speciation rates sampled from a log-normal distribution with an expected value log [α×λ] (where λ represents the parental speciation rate and α is a trend parameter) and a standard deviation σ. We considered the model with constant turnover ε (i.e. constant ratio between extinction and speciation rates; ClaDS2) and ran a newly-developed ClaDS algorithm based on data augmentation techniques which enables us to estimate mean rates through time [START_REF] Maliet | Fast and accurate estimation of species-specific diversification rates using data augmentation[END_REF]. We ran ClaDS with 3 independent chains, checked their convergence using a Gelman-Rubin diagnostic criterion [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF], and recorded lineage-specific speciation rates. We also recorded the estimated hyperparameters (α, σ, ε) and the value m=α×exp(σ 2 /2), which indicates the general trend of the rate through time [START_REF] Maliet | A model with many small shifts for estimating species-specific diversification rates[END_REF]. We replicated these analyses using the LSU gene.

In addition, we applied CoMET (TESS R-package (Höhna, May, & Moore, 2016;May, Höhna, & Moore, 2016)), another diversification approach that does not consider rate variation across lineages, but models temporal shifts in speciation and extinction rates affecting all lineages simultaneously. CoMET is a piecewise-constant model in a Bayesian framework. We chose the Bayesian priors according to maximum likelihood estimates from TreePar [START_REF] Stadler | Mammalian phylogeny reveals recent diversification rate shifts[END_REF], disallowed or not mass extinction events, and ran the MCMC chains until convergence (minimum effective sample sizes of 500).

We also fitted a series of time-dependent and environment-dependent birth-death diversification models using RPANDA [START_REF] Condamine | Macroevolutionary perspectives to environmental change[END_REF][START_REF] Morlon | RPANDA: An R package for macroevolutionary analyses on phylogenetic trees[END_REF] to confirm the observed temporal trends and test the influence of temperature, pCO2, and land plant fossil diversity on rates of Glomeromycotina speciation. For the timedependent models, we considered models with constant or exponential variation of speciation rates through time and null or constant extinction rates (fit_bd function). As extinction is notoriously hard to estimate from reconstructed phylogenies [START_REF] Rabosky | Challenges in the estimation of extinction from molecular phylogenies: A response to Beaulieu and O'Meara[END_REF], we tested the robustness of the inferred temporal trend in speciation when fixing arbitrarily high levels of extinction (Supplementary Methods 4). For the environmentdependent models, we considered an exponential dependency of the speciation rates with the environmental variable (env), i.e. speciation rate=b*exp(a*env), where a and b are two parameters estimated by maximum likelihood (fit_env function). Best-fit models were selected based on the corrected Akaike information criterion (AICc), considering that a difference of 2 in AICc indicates that the model with the lowest AICc is better. We replicated these analyses using the LSU gene.

The influence of temperature was tested on the complete Glomeromycotina phylogenetic trees, using estimates of past global temperature [START_REF] Royer | CO2 as a primary driver of Phanerozoic climate[END_REF]. We also carried a series of simulation analyses to test the robustness of our temperature-dependent results (Supplementary Methods 5). The influence of pCO2 [START_REF] Foster | Future climate forcing potentially without precedent in the last 420 million years[END_REF] and of land plant fossil diversity was tested starting from 400 Myr ago, as these environmental data are not available for more ancient times. For these analyses we sliced the phylogenies at 400 and 200 Myr ago, and applied the diversification models to the sliced sub-trees larger than 50 tips. Estimates of land plant diversity were obtained using all available Embryophyta fossils from the Paleobiology database (https://paleobiodb.org) and using the shareholder quorum subsampling method (Supplementary Methods 6; [START_REF] Alroy | Geographical, environmental and intrinsic biotic controls on Phanerozoic marine diversification[END_REF]).

We considered missing species in all our diversification analyses by imputing sampling fractions, computed as the number of observed VT or EUs divided by the corresponding BDES estimates of global Glomeromycotina diversity (Table 1). We used a global sampling fraction for all Glomeromycotina, as the main Glomeromycotina clades had a similar sampling fraction (Supplementary Table 2). To assess the robustness of our results to global diversity estimates, we replicated all diversification analyses using a range of lower sampling fractions (from 90% to 50%, i.e. assuming that only that percentage of the global Glomeromycotina species diversity is in fact represented in our dataset).

Testing for correlates of present-day Glomeromycotina speciation rates:

To further investigate the potential factors correlating with Glomeromycotina speciation rates, we assessed the relationship between lineage-specific estimates of present-day speciation rates (obtained with the ClaDS analyses) and characteristics of each Glomeromycotina taxonomic unit, i.e. VT or EUs.

First, to assess the effect of specialization on speciation rates, we characterized Glomeromycotina relative niche width using a set of 10 abiotic and biotic variables recorded in MaarjAM database for each Glomeromycotina unit. In short, among a curated dataset containing Glomeromycotina sequences occurring only in natural ecosystems (dataset 2; Supplementary Table 2; [START_REF] Perez-Lamarque | Cheating in arbuscular mycorrhizal mutualism: a network and phylogenetic analysis of mycoheterotrophy[END_REF], for each Glomeromycotina unit, we reported the number of continents, ecosystems, climatic zones, biogeographic realms, habitats, and biomes where it was sampled, as well as its number of plant partners, their phylogenetic diversity, and its centrality in the plant-fungus bipartite network, and performed a principal component analysis (PCA; Supplementary Methods 7). For Glomeromycotina units represented by at least 10 sequences, we tested whether these PCA coordinates reflecting Glomeromycotina niche widths were correlated with the presentday speciation rates using both linear mixed-models (not accounting for phylogeny) or MCMCglmm models [START_REF] Hadfield | MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R package[END_REF]. For MCMCglmm, we assumed a Gaussian residual distribution, included the fungal phylogenetic tree as a random effect, and ran the MCMC chains for 1,300,000 iterations with a burn-in of 300,000 and a thinning interval of 500.

Next, we tested the relationship between speciation rates and geographic characteristics of Glomeromycotina units. To evaluate the effect of latitude on speciation rates, we associated each Glomeromycotina unit with its set of latitudes and used similar MCMCglmm with an additional random effect corresponding to the Glomeromycotina unit. To account for inhomogeneous sampling along the latitudinal gradient, we re-ran the model on jackknifed datasets (we re-sampled 1,000 interactions per slice of latitude of twenty degrees). Similarly, we tested the effect of climatic zone and habitat on speciation rates.

Then, to assess the effect of dispersal capacity on speciation rates, we evaluated the relationship between spore size and speciation rate for the few (n=32) VT that contain sequences of morphologically characterized Glomeromycotina isolates [START_REF] Davison | Microbial island biogeography: isolation shapes the life history characteristics but not diversity of root-symbiotic fungal communities[END_REF]. We gathered measures of their average spore length [START_REF] Davison | Microbial island biogeography: isolation shapes the life history characteristics but not diversity of root-symbiotic fungal communities[END_REF] and tested their relationship with speciation rate by using a phylogenetic generalized least square regression (PGLS).

Finally, as a first attempt at connecting Glomeromycotina macroevolutionary diversification to microevolutionary processes, we measured intraspecific genetic diversities across Glomeromycotina units. For each Glomeromycotina unit containing at least 10 sequences, we computed genetic diversity using Tajima's estimator [START_REF] Tajima | Evolutionary relationship of DNA sequences in finite populations[END_REF](θπ; Supplementary Methods 8). Using similar statistical tests as above, we investigated the correlation of Glomeromycotina genetic diversity with speciation rate, niche width, geographic characteristics, and spore size. We tested the robustness of the results to the minimal number of sequences per Glomeromycotina unit (10, 15, or 20) used to compute genetic diversity and to perform the PCA.

These statistical models were replicated on the different phylogenetic trees (consensus or replicates) for each delineation and we reported p-values (P) corresponding to two-sided tests.

Simulation analyses:

The use of a short and slowly evolving gene such as the central region of the SSU rRNA gene to delineate species may lead to an artificial lumping of species into the same unit that would reduce the number of phylogenetic branching events toward the present and result in a biased inference of temporal diversification dynamics, including an artifactual detection of a diversification slowdown [START_REF] Moen | Why does diversification slow down?[END_REF]. We used simulations mimicking the evolution of the SSU rRNA gene as Glomeromycotina diversified to quantify this potential bias.

We simulated the diversification of a clade of species in the last 505 Myr, according to two scenarios: (i) constant speciation rate and no extinction and (ii) constant speciation and extinction rates (Supplementary Figure 2a). To model intraspecific differentiation, we added intraspecific splits on these simulated species trees by grafting coalescent events at each tip: for each species, we uniformly sampled between 2 and 15 individuals and we considered that all these individuals had to coalesce before the last speciation event; the age of the coalescent tree within each species was uniformly sampled between 0 and the age of the last speciation event (with a maximum of 30 Myr). We used the functions pbtree and rcoal from the R-packages phytools and ape [START_REF] Paradis | APE: Analyses of phylogenetics and evolution in R language[END_REF][START_REF] Revell | phytools: An R package for phylogenetic comparative biology (and other things)[END_REF] to simulate the species phylogenies and the intraspecific coalescences respectively. We used two net diversification rates (r=0.010 and r=0.015) for simulating the species phylogenies, in order to reach a total number of species similar to that obtained with our empirical data when using the VT and EU99 delineations, respectively. Next, we simulated the evolution of short 520 bp DNA sequences on the obtained trees, using the function simulate_alignment (R-package HOME; Perez-Lamarque & Morlon, 2019). We used a substitution rate of 0.001 event per Myr and only 25% of variable sites, which resulted in an alignment that mimicked the Glomeromycotina SSU rDNA alignment. We performed 10 simulations per scenario. For each of these simulations we kept the unique haplotypes at present and applied the same pipelines as above, using the EU99 species delineation criteria: after delineating the EU99 units, we reconstructed the EU99 phylogenetic trees, ran the ClaDS analyses on these trees, and recorded mean estimated speciation rates at present and 50, 100, and 150 Myr ago.

Results:

Glomeromycotina species delineations & phylogenetic reconstructions:

We automatically delineated Glomeromycotina into evolutionary units (EU) using a monophyly criterion and several thresholds of SSU rRNA sequence similarity (from 97% to 99%). The EU97.5 and EU98 delineations (obtained using a threshold of 97.5% and 98% respectively) provided a number of Glomeromycotina units (340 and 641) relatively comparable to the 384 currently recognized virtual taxa (VT), while the EU97 delineation had much less units (182). Conversely, the EU98.5 and EU99 delineations yielded a much larger number of Glomeromycotina units (1,190 and 2,647). These numbers obtained with the EU98.5 and EU99 delineations were consistent with the numbers obtained using GMYC analyses, which delineate species-like units based on detecting when splitting events in the haplotype tree start to follow branching patterns consistent with intra-specific differentiations (i.e. coalescent patterns) instead of speciation events (i.e. birth-death patterns; Supplementary Tables 3, 4, &5). The GMYC results therefore support the idea that some VT might lump together several cryptic species [START_REF] Bruns | Glomeromycotina: what is a species and why should we care?[END_REF](Supplementary Note 1), and that a 98.5 or 99% similarity threshold is more relevant for Glomeromycotina species delineation. In addition, the GMYC model is significantly supported over the model where all SSU rRNA haplotypes correspond to a different species (GMYC LRT: P<0.05; Supplementary Fig. 3), with on average 10 SSU haplotypes per species-like unit, and a mean intraspecific sequence similarity of 99% (Supplementary Table 5 & Supplementary Fig. 3). This indicates that the region of the SSU marker used to characterized Glomeromycotina evolves fast enough to accumulate substitutions between Glomeromycotina speciation events, meaning that it is an informative (although not perfect) marker for delineating Glomeromycotina species-like units. In comparison, the same pipeline carried on the LSU database delineated only 181 GMYC units, suggesting that it was much less complete than the SSU database. We replicated the subsequent diversification analyses using the LSU region, even though we put more trust in our results using the SSU database given the incompleteness of the LSU database.

Rarefaction curves as well as BDES (Bayesian Diversity Estimation Software) and Chao2 estimates of diversity suggested that more than 90% of the total Glomeromycotina diversity is represented in our SSU dataset regardless of the delineation threshold (Fig. 1, Table 1, Supplementary Tables 5 &6), which is consistent with the proportion of new Glomeromycotina units detected in recent studies [START_REF] Sepp | Non-random association patterns in a plant-mycorrhizal fungal network reveal host-symbiont specificity[END_REF].

The reconstructed Bayesian phylogenetic trees based on VT and EU delineations did not yield high support for the nodes separating the main Glomeromycotina orders; yet, the trees had no significantly-supported conflicts either, and similar branching times of the internal nodes overall (Fig. 2, Supplementary Fig. 4). As expected, finer delineations resulted in an increase in the number of nodes close to the present (Supplementary Fig. 5). However, we observed a slowdown in the accumulation of new lineages close to the present in all lineage through time plots (LTTs), including those with the finest delineations (EU98.5 and EU99; Supplementary Fig. 6).

Temporal diversification dynamics:

We found that speciation rates for Glomeromycotina ranged from 0.005 to 0.03 events per lineage per Myr, using both the VT and EU SSU rRNA delineations (Fig. 2; Supplementary Fig. 7). Speciation rates varied both within and among Glomeromycotina orders, with Glomerales and Diversisporales having the highest present-day speciation rates (Supplementary Fig. 8). As expected we observed higher present-day speciation rates for finer delineations, but at the haplotype level (i.e. at the level of the individual SSU rRNA sequences within each unit) we found a significant correlation of present-day speciation rates computed with ClaDS using different delineations (Supplementary Fig. 9). Whatever the delineations, Glomeromycotina experienced their highest speciation rates between 200 and 100 Myr ago according to estimates obtained with ClaDS (Fig. 2; Supplementary Fig. 10) and between 150 and 50 Myr ago according to CoMET (Fig. 2; Supplementary Fig. 11).

ClaDS estimates of speciation rates at 150 Myr ago were 26% (± s.d. 17) higher than those at 300 Myr with the EU99 delineation. With the VT delineation, the increase was of 3% (±8).

The peak was even stronger using CoMET: 30% ± 20 higher at 150 Myr in comparison to 300 Myr with the EU99 delineation (71% ± 40 with the VT delineation; Fig. 2).

The peak of speciation rates was followed by a decline in the recent past (Fig. 2; Supplementary Fig. 10), as suggested by the plateauing of the LTTs. A global decline of the speciation rates through time was independently supported by ClaDS and CoMET analyses, as well as time-dependent models in RPANDA (Morlon, Parsons, & Plotkin, 2011)(Supplementary Figs. 11, 12, & 13). This speciation rate decline was robust to all species delineations, the branching process prior (Supplementary Table 7), phylogenetic uncertainty, and assumed sampling fractions as low as 50%, except in ClaDS analyses where the trend disappeared in some EU99 trees and for sampling fractions lower than 70% (Supplementary Figs. 14 & 15). We also found a period of high speciation rates between 200 and 100 Myr ago followed by a decline in our analyses with the LSU region, for assumed sampling fractions as low as 60% (Supplementary Figure 16).

We did not find a strong signal of extinction in our analyses: the turnover rate estimated from ClaDS was generally close to zero (Supplementary Fig. 12b), and models including extinctions were never selected in RPANDA (Supplementary Fig. 13). Similarly, the extinction rates estimated in piecewise-constant models (CoMET) were not significantly different from 0 and we did not find significant support for mass extinction events (Supplementary Fig. 17). Yet, forcing the extinction rate to high positive values did not modify the general trend of speciation rate slowdown (Supplementary Figs. 18 & 19).

Correlates of Glomeromycotina diversification:

When fitting environment-dependent models of diversification, we found that temperature-dependent models better fit Glomeromycotina diversification than timedependent models, with higher speciation rates during warm climatic periods (Fig. 3; Supplementary Fig. 20). This was true for all Glomeromycotina delineations, sampling fractions, and crown ages (Supplementary Figs. 21,22,23,& 24), with the exception of some EU99 trees with a 50% sampling fraction (Supplementary Fig. 24). It was also true in our analyses using the LSU region, for sampling fractions down to 50% (Supplementary Figure 29). This signal of temperature dependency was not due to a temporal trend (Supplementary Figs. 25 & 26) nor to an artefact caused by rate heterogeneities (Supplementary Fig. 27). Evidence for temperature dependency, however, decreased in some clades closer to the present, as small trees tend to be best fit by constant or timedepend models (Supplementary Fig. 28). We detected a significant positive dependency of the speciation rates on CO2 concentrations in some sub-trees, but rarely found a significant effect of plant fossil diversity (Supplementary Fig. 28).

The PCA of Glomeromycotina relative niche width characteristics had a first principal component (PC1) that indicated the propensity of each Glomeromycotina unit (VT or EUs) to be vastly distributed among continents, ecosystems and/or associated with many plant species and lineages (i.e. high generalism), whereas the second principal component (PC2) indicated the propensity of a given Glomeromycotina unit to associate with few plant species on many continents (i.e. high specialism toward plants; Supplementary Figs. 30, 31, & 32). Hence, PC1 reflects Glomeromycotina niche width, whereas PC2 discriminates the width of the abiotic relatively to the biotic niche (Fig. 4a-b).

We found a positive correlation between lineage-specific speciation rates and PC1 in the majority of the VT and EU99 trees, but no significant correlation with PC2 (Fig. 4c-d; Supplementary Fig. 33a). However, these results were no longer significant when controlling for phylogenetic non-independence between Glomeromycotina units (Supplementary Fig. 33b), likely because a single Glomeraceae clade, including the abundant and widespread morphospecies Rhizophagus irregularis and R. clarus (high PC1 values), had both the highest speciation rates and the largest niche widths among Glomeromycotina (Supplementary Fig. 34).

Although Glomeromycotina diversity is currently higher in the (sub)tropics (Supplementary Fig. 35), we found no effect of latitude on speciation rates, regardless of the Glomeromycotina delineation or the minimum number of sequences per Glomeromycotina unit (MCMCglmm: P>0.05). In addition, we actually did not detect a higher total number of Glomeromycotina species in grasslands compared to forests (Supplementary Figure 36; confirming the results of [START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF], and it is thus not surprising that we reported no effect of habitat or climatic zone on speciation rates (Supplementary Fig. 37), suggesting that tropical grasslands are not particular diversification hotspots for Glomeromycotina. Similarly, we recovered no significant correlation between spore size and speciation rate (Supplementary Fig. 38), nor between spore size and level of endemism (Supplementary Fig. 39).

Finally, Tajima's estimator of Glomeromycotina genetic diversity was significantly and positively correlated with niche width (PC1) for all Glomeromycotina delineations and minimal number of sequences per Glomeromycotina unit considered, and in particular with abiotic aspects of the niche (PC2) in many cases (Fig. 4e-h; Supplementary Fig. 33).

Genetic diversity was not correlated with speciation rate (Supplementary Fig. 33), latitude, habitat, climatic zone (MCMCglmm: P>0.05), or spore size (PGLS: P>0.05).

Simulation results:

When we simulated the evolution of a short DNA gene mimicking the SSU rRNA marker and used it to delineate species, we found that the number of EU99 delineated units was generally lower than the number of simulated species (~10% to 20% lower; Supplementary Figure 2b). Hence, even the EU99 delineation tends to lump together some closely related species. As expected, this lumping resulted in an artefactual inference of a decline of speciation rates toward the present, but this artifactual decline was significantly smaller in magnitude than that observed in Glomeromycotina (Figure 5). Hence, these analyses suggest that the lumping of species resulting from the use of a small, slowly evolving marker is unlikely to fully explain the strong temporal decline in speciation rate we found in Glomeromycotina.

Discussion

Glomeromycotina species delineations, diversity, and phylogeny:

It is difficult to delineate species in Glomeromycotina, which are poorly differentiated morphologically and mainly characterized by environmental sequences [START_REF] Bruns | Glomeromycotina: what is a species and why should we care?[END_REF]. Our GMYC analyses suggest that Glomeromycotina species-like units correspond to SSU rRNA haplotypes with a sequence similarity between 98.5 and 99%.

With this criterion of species delineation, we estimate that there are between 1,300 and 2,900 Glomeromycotina 'species'. These estimates are largely above the number of currently described morphospecies or VT (Supplementary Note 1) but remain low in comparison with other fungal groups, like the Agaricomycetes that include taxa forming ectomycorrhiza [START_REF] Varga | Megaphylogeny resolves global patterns of mushroom evolution[END_REF].

Our phylogenies based on the SSU rRNA gene did not resolve the branching of the Glomeromycotina orders, with node supports similar to those of previous studies [START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF][START_REF] Krüger | Phylogenetic reference data for systematics and phylotaxonomy of arbuscular mycorrhizal fungi from phylum to species level[END_REF][START_REF] Rimington | Ancient plants with ancient fungi: liverworts associate with early-diverging arbuscular mycorrhizal fungi[END_REF](Supplementary Note 2). These findings confirm that additional genomic evidence is required to reach consensus. We considered this uncertainty in species delineation and phylogenetic reconstruction by repeating our diversification analyses across species delineation criteria and on a set of trees spanning the likely tree space. We found effects of species delineation consistent with a priori expectations: criteria that lump together more dissimilar sequences (e.g. those that use a lower percentage of similarity cut-off) result in lower diversity estimates, lower estimates of speciation rates, and patterns of diversification through time that reflect longer terminal branch-lengths, such as peaks of diversification that occur earlier. Despite this variability, we found that general patterns, such as the observed temporal decline in speciation rates and the significant association between temperature and speciation rates, were consistent across species delineations and trees. Therefore, our study based on a short SSU (or LSU) rRNA region should encourage both efforts to obtain more genetic data, including longer reads [START_REF] Krehenwinkel | Nanopore sequencing of long ribosomal DNA amplicons enables portable and simple biodiversity assessments with high phylogenetic resolution across broad taxonomic scale[END_REF][START_REF] Tedersoo | Perspectives and benefits of high-throughput long-read sequencing in microbial ecology[END_REF] and additional genomic information, with the aim of reconstructing better supported, comprehensive phylogenies and efforts to conduct diversification analyses despite uncertainty in the data for groups where better data is not yet available.

Glomeromycotina diversify slowly:

We found speciation rates for Glomeromycotina an order of magnitude lower than rates typically found for macro-eukaryotes [START_REF] Maliet | A model with many small shifts for estimating species-specific diversification rates[END_REF][START_REF] Upham | Inferring the mammal tree: Species-level sets of phylogenies for questions in ecology, evolution, and conservation[END_REF], like plants [START_REF] Zanne | Three keys to the radiation of angiosperms into freezing environments[END_REF], or Agaricomycetes [START_REF] Varga | Megaphylogeny resolves global patterns of mushroom evolution[END_REF]. Low speciation rates in Glomeromycotina may be linked to their multinucleate hyphal state [START_REF] Yildirir | Parasexual and sexual reproduction in arbuscular mycorrhizal fungi: Room for both[END_REF], to their occasional long-distance dispersal that homogenizes populations globally over evolutionary timescales [START_REF] Savary | A population genomics approach shows widespread geographical distribution of cryptic genomic forms of the symbiotic fungus Rhizophagus irregularis[END_REF], and/or to the fact that they are generalist obligate symbionts [START_REF] Morlon | Explosive radiation of a bacterial species group[END_REF].

Regardless of the proximal cause, and contrary to Agaricomycetes for example, which present a large diversity of species, morphologies, and ecologies, Glomeromycotina have poorly diversified in the last 500 Myr despite their ubiquity; their niche space is restricted to plant roots and the surrounding soil because of their obligate dependence on plants for more than 400 Myr [START_REF] Rich | Diet of arbuscular mycorrhizal fungi: Bread and butter?[END_REF][START_REF] Tisserant | Genome of an arbuscular mycorrhizal fungus provides insight into the oldest plant symbiosis[END_REF].

Our estimates of speciation rates were highly variable across lineages. We reported the highest speciation rates in Glomeraceae and Diversisporaceae. Speciation rates in Paraglomeraceae and Archaeosporaceae, which are thought to be less beneficial for the plants than the fast diversifying Glomeraceae and Diversisporaceae [START_REF] Säle | Ancient lineages of arbuscular mycorrhizal fungi provide little plant benefit[END_REF], were an order of magnitude lower. We can therefore speculate that good symbiotic abilities may favor diversification, although this remains to be tested in further investigations.

We found little evidence for species extinction in Glomeromycotina, including at mass extinction events. Because Glomeromycotina are relatively widespread and have a an ancient tendency toward generalism, they might therefore be quite resilient to land plant mass extinctions and low extinction rates have been predicted before based on their ecology [START_REF] Morton | Species and clones of arbuscular mycorrhizal fungi (Glomales, Zygomycetes): their role in macro-and microevolutionary processes[END_REF]). Yet, these low extinction rate estimates could also come from the difficulty of estimating extinction from molecular phylogenies [START_REF] Rabosky | Challenges in the estimation of extinction from molecular phylogenies: A response to Beaulieu and O'Meara[END_REF], one of the limitations of phylogeny-based diversification analyses (Supplementary Note 3). Fossils of Glomeromycotina that can be ascribed to species or genera are too scarce to support or conflict with this finding.

Glomeromycotina diversification through time:

The observed peak of Glomeromycotina speciations detected between 200 and 100 Myr (or 150-50 Myr depending on the models) was mainly linked to the frequent speciations in the largest family Glomeraceae. This peak was concomitant with the radiation of flowering plants [START_REF] Sauquet | Key questions and challenges in angiosperm macroevolution[END_REF], but also with a major continental reconfiguration, including the breakdown of Pangea and the formation of climatically contrasted landmasses [START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF]. This period was also characterized by a warm climate potentially directly or indirectly favorable to Glomeromycotina diversification, such that disentangling the impact of these various factors on Glomeromycotina diversification rates is not straightforward. Interestingly, a peak of speciations at this period was also found in the Agaricomycetes, a clade of fungi including lineages forming ectomycorrhizae [START_REF] Varga | Megaphylogeny resolves global patterns of mushroom evolution[END_REF]. This peak in the occurrence of speciation events was followed by a decline in speciation rates. The detection of temporal declines in speciation rates in phylogenetic diversification analyses can sometimes be artifactual, for example if some species are incorrectly lumped together during species delineation or if the proportion of species not represented in the phylogeny is under-estimated [START_REF] Moen | Why does diversification slow down?[END_REF]). We considered these potential biases, conducted sensitivity analyses, and found that the observed slowdown was robust, and even amplified under scenarios of high extinction. Some Glomeromycotina species are likely lumped together into the same SSU haplotypes [START_REF] Krüger | Phylogenetic reference data for systematics and phylotaxonomy of arbuscular mycorrhizal fungi from phylum to species level[END_REF], but both our use of an overly small assumed sampling intensity (50%) and our simulation analyses demonstrated that this lumping is not sufficient to explain the observed slowdown. In addition, we also detected a temporal decline in speciation rates when using another marker (the LSU rRNA gene).

Temporal declines in speciation rates have been observed in many clades, including microorganisms [START_REF] Condamine | Assessing the causes of diversification slowdowns: temperature-dependent and diversity-dependent models receive equivalent support[END_REF][START_REF] Morlon | Explosive radiation of a bacterial species group[END_REF][START_REF] Rabosky | Density-dependent diversification in North American wood warblers[END_REF]. They have often been interpreted as a progressive reduction of the number of available niches as species diversify and accumulate [START_REF] Moen | Why does diversification slow down?[END_REF][START_REF] Rabosky | Ecological limits and diversification rate: Alternative paradigms to explain the variation in species richness among clades and regions[END_REF]. In Glomeromycotina, this potential effect of niche saturation could be exacerbated by a reduction of their niches linked to both repetitive breakdowns of their symbiosis with plants and climatic changes. Indeed, since the Cretaceous, many plant lineages evolved alternative root symbioses or became non-symbiotic [START_REF] Brundrett | Evolutionary history of mycorrhizal symbioses and global host plant diversity[END_REF][START_REF] Maherali | Mutualism persistence and abandonment during the evolution of the mycorrhizal symbiosis[END_REF][START_REF] Selosse | The land flora: a phototroph-fungus partnership?[END_REF][START_REF] Werner | Symbiont switching and alternative resource acquisition strategies drive mutualism breakdown[END_REF]: approximately 20% of extant plants do not interact with Glomeromycotina anymore [START_REF] Van Der Heijden | Mycorrhizal ecology and evolution: the past, the present, and the future[END_REF]. Additionally, the cooling of the Earth during the Cenozoic reduced the surface of tropical regions [START_REF] Meseguer | Ancient tropical extinctions at high latitudes contributed to the latitudinal diversity gradient[END_REF][START_REF] Ziegler | Tracing the tropics across land and sea: Permian to present[END_REF], which tend to be a reservoir of ecological niches for Glomeromycotina [START_REF] Brundrett | Evolutionary history of mycorrhizal symbioses and global host plant diversity[END_REF][START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF][START_REF] Read | Mycorrhizas in ecosystems[END_REF].

The difficulty of reconstructing past symbiotic associations prevents direct testing the hypothesis that the emergence of new root symbioses in plants led to a decline in speciation rates in Glomeromycotina. However, we tested the hypothesis that global temperature changes affected speciation rates and found a strong relationship. Such associations between temperature and speciation rates have been observed before in eukaryotes and have several potential causes [START_REF] Condamine | Assessing the causes of diversification slowdowns: temperature-dependent and diversity-dependent models receive equivalent support[END_REF]. In particular, the productivity hypothesis states that resources and associated ecological niches are more numerous in warm and productive environments, especially when the tropics are large, which entail higher speciation rates [START_REF] Clarke | Climate, energy and diversity[END_REF]. This hypothesis is particularly relevant for Glomeromycotina, which have many host plant niches in the tropics, as shown by their latitudinal diversity gradient, and potentially relatively less in temperate and polar regions [START_REF] Toussaint | Asymmetric patterns of global diversity among plants and mycorrhizal fungi[END_REF], where a higher proportion of plants are non-mycorrhizal [START_REF] Bueno | Plant mycorrhizal status, but not type, shifts with latitude and elevation in Europe[END_REF] or ectomycorrhizal [START_REF] Brundrett | Evolutionary history of mycorrhizal symbioses and global host plant diversity[END_REF][START_REF] Varga | Megaphylogeny resolves global patterns of mushroom evolution[END_REF].

Hence, the observed effect of past global temperatures could reflect the shrinkage of tropical areas and the associated decrease of the relative proportion of arbuscular mycorrhizal plants. Future developments of diversification models incorporating interspecific interactions would allow us to better test these hypotheses.

A few Glomeromycotina clades displayed a significant support for diversification models with a positive dependency of speciation rates on CO2 concentrations, which reinforces the idea that for the corresponding Glomeromycotina, benefits retrieved from plants could have been amplified by high CO2 concentrations and fostered diversification [START_REF] Field | Functional analysis of liverworts in dual symbiosis with Glomeromycota and Mucoromycotina fungi under a simulated Palaeozoic CO2 decline[END_REF][START_REF] Humphreys | Mutualistic mycorrhiza-like symbiosis in the most ancient group of land plants[END_REF]. Conversely, we found a limited effect of land plant fossil diversity, which indicates that variations in the tempo of Glomeromycotina diversification did not systematically follow those of land plants. Still, the possible concordance of the peak of Glomeromycotina speciations with the radiation of the Angiosperms is noteworthy, in particular in Glomeraceae that frequently interact with present-day Angiosperms [START_REF] Rimington | Ancient plants with ancient fungi: liverworts associate with early-diverging arbuscular mycorrhizal fungi[END_REF]. Plant diversification might have fostered the diversification of Glomeromycotina from the emergence of land plants until the Mesozoic [START_REF] Lutzoni | Contemporaneous radiations of fungi and plants linked to symbiosis[END_REF][START_REF] Morton | Species and clones of arbuscular mycorrhizal fungi (Glomales, Zygomycetes): their role in macro-and microevolutionary processes[END_REF], but less so thereafter, when Glomeromycotina diversification declined while some flowering plants radiated, including Glomeromycotina-associated groups, like the Poaceae, but also Glomeromycotina-free groups such as the extraordinary radiation of Orchidaceae [START_REF] Givnish | Orchid phylogenomics and multiple drivers of their extraordinary diversification[END_REF], blurring co-diversification patterns (Supplementary Fig. 40) [START_REF] Cleal | Composition and dynamics of the great Phanerozoic Evolutionary Floras[END_REF][START_REF] Ramírez-Barahona | The delayed and geographically heterogeneous diversification of flowering plant families[END_REF].

Correlates of Glomeromycotina recent speciation rates:

Looking at the correlates of Glomeromycotina present-day speciation rates, we found no effect of habitat or climatic zone, even though Glomeromycotina are more frequent and diverse in the tropics [START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF][START_REF] Pärtel | Historical biome distribution and recent human disturbance shape the diversity of arbuscular mycorrhizal fungi[END_REF][START_REF] Toussaint | Asymmetric patterns of global diversity among plants and mycorrhizal fungi[END_REF] and a positive correlation with global temperature. Further work, including a more thorough sampling of the distribution of Glomeromycotina species across latitudes and habitats, would be required to confirm these patterns and to distinguish whether speciation events are indeed no more frequent in the tropics or, if they are, whether long-distance dispersal redistributes the new lineages at different latitudes over long time scales. Contrary to previous predictions [START_REF] Pärtel | Historical biome distribution and recent human disturbance shape the diversity of arbuscular mycorrhizal fungi[END_REF], we did not find that tropical grasslands are diversification hotspots for Glomeromycotina; we actually did not even find higher Glomeromycotina species richness in (tropical) grasslands versus forests at global scale, in agreement with [START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF].

Similarly, although the temporal changes in the availability of Glomeromycotina niches likely influenced the diversification of the group, we found little support for Glomeromycotina species with larger niche width having higher lineage-specific speciation rates. We also note that there are important aspects of the niche that we do not (and yet cannot) account for in our characterization of Glomeromycotina niche width: it is thought that some Glomeromycotina species may mainly provide mineral nutrients extracted from the soil, whereas others may be more specialized in protecting plants from biotic or abiotic stresses [START_REF] Chagnon | A trait-based framework to understand life history of mycorrhizal fungi[END_REF] and such (interor intra-specific) functional variations may have evolutionary significance. Finally, although spore size is often inversely related to dispersal capacity [START_REF] Nathan | Mechanisms of long-distance seed dispersal[END_REF], which can limit speciation by increasing gene flow, we found no significant correlation between spore size and speciation rates, which may be explained either by a weak or absent effect or by the low number of species for which this data is available. In addition, the absence of correlation between spore size and level of endemism suggests that even Glomeromycotina with large spores experience long-distance dispersal [START_REF] Davison | Microbial island biogeography: isolation shapes the life history characteristics but not diversity of root-symbiotic fungal communities[END_REF][START_REF] Kivlin | Global mycorrhizal fungal range sizes vary within and among mycorrhizal guilds but are not correlated with dispersal traits[END_REF]. Thus, if large spores might limit dispersal at smaller (e.g. intracontinental) scales in Glomeromycotina [START_REF] Bueno | How do arbuscular mycorrhizal fungi travel?[END_REF][START_REF] Chaudhary | Trait-based aerial dispersal of arbuscular mycorrhizal fungi[END_REF], this does not seem to affect speciation rates.

In Glomeromycotina, intraspecific variability is an important source of functional diversity [START_REF] Munkvold | High functional diversity within species of arbuscular mycorrhizal fungi[END_REF][START_REF] Savary | A population genomics approach shows widespread geographical distribution of cryptic genomic forms of the symbiotic fungus Rhizophagus irregularis[END_REF] and their genetic diversity may indicate the intraspecific variability on which selection can act, potentially leading to speciation. Here, geographically widespread Glomeromycotina species appear to be more genetically diverse, as previously suggested by population genomics [START_REF] Savary | A population genomics approach shows widespread geographical distribution of cryptic genomic forms of the symbiotic fungus Rhizophagus irregularis[END_REF], but do not necessarily speciate more frequently. Along with a decoupling between genetic diversity and lineage-specific speciation rate, this suggests that the accumulation of genetic diversity in the SSU region among distant subpopulations is not enough to spur Glomeromycotina speciation.

Analyzing diversification dynamics using a short marker gene:

Short DNA regions, like those used in metabarcoding surveys, typically do not allow to robustly delineate species, estimate global-scale diversity, and reconstruct phylogenetic trees. As these three aspects can all affect results of diversification analyses [START_REF] Moen | Why does diversification slow down?[END_REF], such analyses are rarely performed with these types of data. Yet, for many species-rich groups of organisms, in particular microorganisms, no other data currently provide a thorough representation of diversity at the "species" level. Hence, these data, although far from ideal, are the only one that can be used to study the past diversification of such groups (see [START_REF] Lewitus | Clade-specific diversification dynamics of marine diatoms since the Jurassic[END_REF]Louca et al., 2018 for antecedents).

The approach we took here is to recognize all these potential sources of uncertainty and biases and to test the robustness of our results. We demonstrated the usefulness of this approach: while some results inevitably depend on the choices made for species delineation, phylogenetic reconstruction, and the estimation of global scale diversity, others are sufficiently strong to hold despite uncertainty in the data. Our results therefore illustrate that using a short DNA marker (e.g. a metabarcode) combined with intensive sensitivity analyses can be useful for studying the diversification dynamics of poorlyknown groups.

Conclusion:

Our findings that Glomeromycotina have low speciation rates, likely constrained by the availability of suitable niches, reinforce the vision of Glomeromycotina as an "evolutionary cul-de-sac" [START_REF] Malloch | The evolution of mycorrhizae[END_REF]. We interpret the significant decline in speciation rates toward the present as the conjunction of the emergence of plant lineages not associated with Glomeromycotina and the reduction of tropical areas induced by climate cooling, in the context of obligate dependence of Glomeromycotina on plants.

Temporal declines in speciation rates have often been interpreted as the signal of adaptive radiations [START_REF] Harmon | Tempo and mode of evolutionary radiation in iguanian lizards[END_REF][START_REF] Moen | Why does diversification slow down?[END_REF], that is clades that experienced a rapid accumulation of morphological, ecological, and species diversity [START_REF] Simpson | The major features of evolution[END_REF]. Conversely, Glomeromycotina provide here a striking example of a clade with slow morphological, ecological, and species diversification that features a pattern of temporal decline in speciation rates, that might reflect the reduction of the global availability of their mycorrhizal niches.

Tables:

Table 1: Estimation of the total diversity of Glomeromycotina: Estimated sampling fraction using the Bayesian Diversity Estimation Software (BDES; [START_REF] Quince | The rational exploration of microbial diversity[END_REF] for the different species delineations (VT or EU) assuming a Sichel species abundance distribution. The estimated number of units corresponds to the median value and we indicated the 95% confidence interval. We indicated the sampling fractions for each delineation, computed as the number of observed VT or EUs divided by the corresponding BDES estimates of global Glomeromycotina diversity.

The Sichel distribution was selected compared to other distributions (log-normal, log-Student, and inverse gaussian) based on lowest deviance information criterion (DIC). We compared the virtual taxa (VT) delineation from (Öpik et al., 2010) because some of the 12 replicate trees showed different trends, we replicated ClaDS inferences using 100 replicate trees. Unlike most replicate trees, the EU99 consensus tree tends to present a limited decline in speciation rates, which reinforces the idea that consensus trees can be misleading [START_REF] Janzen | Inferring the role of habitat dynamics in driving diversification: evidence for a species pump in Lake Tanganyika cichlids[END_REF]. For both delineations, the boxplots represent the results obtained for the consensus tree and the 12 independent replicate trees. Boxplots indicate the median surrounded by the first and third quartiles, and whiskers extend to the extreme values but no further than 1.5 of the inter-quartile range. The horizontal dotted lines highlighted the values estimated for the consensus trees. Compared to the replicate trees, the consensus trees tends to present extreme values (stronger support for temperature-dependent model), which reinforces the idea that consensus trees can be a misleading representation [START_REF] Janzen | Inferring the role of habitat dynamics in driving diversification: evidence for a species pump in Lake Tanganyika cichlids[END_REF]. Average temperature

• • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -500 -
• • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -500 -
• • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
-500 -400 -300 -200 -100 0

Time (in Myr)

Average temperature

• • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
-500 -400 -300 -200 -100 0

Time (in Myr)

Average temperature Tested variables were: the numbers of continents on which the Glomeromycotina unit occurs (nb_continent), of realms (nb_realm), of ecosystems (nb_ecosystems), of habitats (nb_habitats), of biomes (nb_biomes), and climatic zones (nb_climatic) (Öpik et al., 2010), as well as information about the associated plant species of each unit, such as the number of plant partners (nb_plants), the phylogenetic diversity of these plants (PD), and the betweenness and closeness measurement of each fungal unit in the plant-fungus interaction network (see Methods).

• • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -500 -
(c-d): Speciation rates as a function of the PC1 coordinates for each VT (c) or EU99 (d) unit.

Only the Glomeromycotina consensus tree is represented here (other replicate trees are presented in Supplementary Fig. 33).

(e-h): Genetic diversity (Tajima's θπ estimator) as a function of the PC1 (e-f) or PC2 (g-h) coordinates for each VT (e-g) or EU99 (f-h) unit. Only the Glomeromycotina consensus tree is represented here (other replicate trees are presented in Supplementary Fig. 33). The grey lines indicate the statistically significant linear regression between the two variables inferred using MCMCglmm.

1213 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -6 -4 -2 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0 5 0.006 0.010 0.016 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -6 -4 -2 0 0.005 0.020 0.100 PCA2 Genetic diversity (Nei) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -10 -8 -6 -4 -2 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • -5 0 5 0.01 0.05 PCA1 Genetic diversity (Nei) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -6 -4 -2 0 0.005 0.020 0.100 PCA2 Genetic diversity (Nei) • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0 5 0.005 0.020 PCA1 Genetic diversity (Pi) θ π θ π d f h • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -5 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • •• • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • •• • • • • -10 -8 -6 -4 -2 0
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • •• • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • •• • • • • -10 -8 -6 -4 -2 0

  discussions and C. Fruciano and E. Lewitus for providing codes. This work was supported by a doctoral fellowship from the École Normale Supérieure de Paris attributed to BPL and the École Doctorale FIRE -Programme Bettencourt. MÖ was supported by the European Regional Development Fund (Centre of Excellence EcolChange) and University of Tartu (PLTOM20903). Funding of the research of FM was from the Agence Nationale de la Recherche (ANR-19-CE02-0002). HM acknowledges support from the European Research Council (grant CoG-PANDA).

Figures Figure 1 :

 1 Figures

  Figure 2: The speciation dynamic of Glomeromycotina (arbuscular mycorrhizal fungi) varies significantly through time and between lineages.

  Figure 3: Temperature-dependent diversification models reveal that global temperature positively associates with the speciation rates of Glomeromycotina (arbuscular mycorrhizal fungi) in the last 500 million years.

  Figure 4: Abiotic and biotic correlates of speciation rates and genetic differentiation in Glomeromycotina (arbuscular mycorrhizal fungi) (a-b): Projection of 10 abiotic and biotic variables on the two principal coordinates according to the VT (a) or EU99 (b) delineations. Principal coordinate analysis (PCA) was performed for the Glomeromycotina units represented by at least 10 sequences. Colors represent the contribution of the variable to the principal coordinates. The percentage for each principal coordinate (PC) indicates its amount of explained variance.

Acknowledgment:

The authors acknowledge C. Strullu-Derrien, M. Elias, D. de Vienne, A. Vogler, J.-Y. Dubuisson, C. Quince, S.-K. Sepp, and M. Chase for helpful discussions. They also thank L. Aristide, S. Lambert, J. Clavel, I. Quintero, I. Overcast, and G. Sommeria for comments on the manuscript, D. Marsh for English editing, and the Editor and three reviewers for improvements of an earlier version of this manuscript. BPL acknowledges B. Robira, F. Foutel-Rodier, F. Duchenne, E. Faure, E. Kerdoncuff, R. Petrolli, and G. Collobert for useful

Data Accessibility:

All of the data used in this study are available in the open-access MaarjAM database (https://maarjam.botany.ut.ee). Spore lengths of Glomeromycotina were collected in the supplementary data of [START_REF] Davison | Microbial island biogeography: isolation shapes the life history characteristics but not diversity of root-symbiotic fungal communities[END_REF]. New scripts for delineating evolutionary units (EU) are available in the R-package RPANDA (through this GitHub branch for the moment: https://github.com/hmorlon/PANDA/tree/Benoit_phylosignal). The alignment of all the SSU rRNA sequences and their associated metadata publicly accessible through the Open Science Framework (osf) portal: osf.io/y2ts5.

Author contributions:

All the authors designed the study. MÖ gathered the data and BPL performed the analyses.

OM and ACAS provided some codes. BPL and HM wrote the first version of the manuscript and all authors contributed substantially to the revisions.

Competing Interests statement:

The authors declare that there is no conflict of interest. Comparison of the magnitude of the decline in speciation rates observed in Glomeromycotina (in orange) and on simulated data (in blue). The intensity of the slowdown is measured as the difference between the mean speciation rate estimated at present and the mean speciation rate estimated 150 Myr ago (a), 100 Myr ago (b), or 50 Myr ago (c). Negative differences indicate a speciation rate decline. Sequence alignments were simulated on phylogenetic trees obtained under a scenario of constant speciation rate and no extinction (i.e. pure birth "PB") and constant speciation and extinction rates (i.e. birth death "BD"), with characteristics mimicking the slow evolution of the SSU rRNA marker.

PCA1 Speciation rates (/Myr)

We simulated phylogenies with two different net diversification rates, such that we obtained simulations with total numbers of species similar to the total numbers of VT or EU99 units (Supplementary Fig. 2). Boxplots indicate the median surrounded by the first and third quartiles, and whiskers extend to the extreme values but no further than 1.5 of the inter-quartile range. Each boxplot represents results for the consensus trees and the 12 independent replicate trees for each of the 10 simulations, and for the consensus trees and the 100 independent replicate trees for the Glomeromycotina. Differences between the magnitude of the decline measured in Glomeromycotina (VT or EU99) and in the corresponding simulations were tested using linear models (reported at the top of the plots).

Supplementary Tables Legends:

Supplementary Figures Legends:

Supplementary