European mushroom assemblages are phylogenetically structured by temperature

Our results show that macrofungal assemblages are phylogenetically structured by temperature across Europe, suggesting phylogenetically constrained specialization towards temperature extremes. Predicted anthropogenic warming is likely to affect species composition and phylogenetic diversity with additional consequences for the carbon- and nutrient cycles.


Introduction
There is evidence that climate change affects biodiversity and key ecosystem processes (Lindner et al. 2010, Bellard et al. 2012. To predict the effects of climate change on biodiversity, we need to better understand the processes structuring species assemblages along climate gradients (Urban et al. 2016). Despite the importance of fungi for ecosystem processes (Cavicchioli et al. 2019), studies of fungi across climate gradients at large spatial scales are scarce compared to those of plants and animals (Davies et al. 2011, Tedersoo et al. 2014. Recent studies have shown that climatic variables are the main factors that correlate with fungal diversity at large spatial scales. Tedersoo et al. (2014) revealed mean annual precipitation as the main factor explaining global soilrelated fungal diversity based on metabarcoding. In contrast, Větrovský et al. (2019) found that temperature-related variables are more important at a similar spatial scale. At subcontinental (central and northern Europe) scale, temperature mean and variability are the main drivers affecting fungal fruit body diversity (Abrego et al. 2017, Andrew et al. 2018. However, besides temperature-related variables, soilrelated (e.g. pH, soil organic content, bulk density) variables also showed significant relationships with fungal diversity but to a lesser extent than temperature-related factors (Tedersoo et al. 2014, Andrew et al. 2018, Větrovský et al. 2019. There is some evidence that macrofungi are affected by climate through their fruit body traits, e.g. a relationship between moisture and spore size (Kauserud et al. 2011). More specifically, macrofungal communities seem to assemble via thermally-relevant morphological fruit body traits such as colour and size at a large spatial scale (Krah et al. 2019, Bässler et al. 2021. These studies suggest evolutionary adaptations to average temperatures and, therefore, ecological selection by temperature as a major assembly filter of macrofungal communities at large scales. However, it is currently not known whether fungi are phylogenetically structured along with temperature in Europe. Phylogenetic alpha diversity indicates whether extreme temperature conditions filter local assemblages towards specific fungal lineages, which would lead to phylogenetic clustering within assemblages. Extreme temperature conditions can occur at the warm and cold edge of the mean annual temperature gradient. Further, extreme temperature conditions can occur at the temperature seasonality edge, characterized by a high intra-annual variability (i.e. continental in contrast to maritime temperature). Phylogenetic beta diversity indicates whether assemblages in extreme temperature conditions (warm versus cold, maritime versus continental) are recruited from similar or different lineages.
A previous study on macrofungal fruiting diversity, at the same spatial scale, revealed that species richness decreased with extreme mean temperatures (cold or warm) and increased with temperature seasonality (Andrew et al. 2019). Temperature seasonality varies across continents, with an increase of temperature ranges and the occurrence of extreme temperatures from maritime towards continental environments. We hypothesize that the observed species richness decrease with extreme mean temperatures is mirrored by a decrease in phylogenetic alpha diversity (phylogenetic clustering within communities, when controlling for species richness, Fig. 1B). Phylogenetic clustering in extreme climate conditions has been shown for bacteria (Wang et al. 2012) and plants (Li et al. 2014), supporting this expectation. We also hypothesize a phylogenetic clustering with temperature seasonality because we expect adaptations to tolerate large temperature ranges and extremes to be phylogenetically constrained (Fig. 1C). Studies using temperature seasonality as a predictor for phylogenetic alpha diversity are scarce. However, bat communities from the Atlantic Forests of South America were characterized by phylogenetic clustering in areas with the greatest seasonality (Stevens and Gavilanez 2015). Phylogenetic clustering in extreme mean and seasonal temperature environments would indicate environmental filtering to be the dominant assembly process (Gotelli 2000, Tucker et al. 2017, via the selection of lineages with extreme temperature-adapted abilities. Andrew et al. (2018) revealed a strong turnover of macrofungal fruiting assemblages related to temperature mean and seasonality across Europe. However, it is currently unclear whether this turnover is caused by entire lineages or species across different lineages. Analyses of phylogenetic beta diversity can inform about the evolutionary mechanisms behind observed community change with temperature variables across large spatial scales (Graham and Fine 2008). For example, the phylogenetic beta diversity of tropical tree species between two plots increases with an increasing difference in temperature at a regional scale, implying that habitat specialization by entire lineages has an important role in structuring local communities (González-Caro et al. 2014). If disparate fungal lineages are adapted to extreme temperature environments, we would observe a strong correlation between phylogenetic beta diversity and mean and seasonal temperature differences (distance, Fig. 1C). If, however, many clades across the phylogenetic tree include species with adaptations to extreme temperature conditions, we would observe no or only a weak relationship between fungal phylogenetic beta diversity and temperature distances.
To address whether macrofungal fruiting assemblages are phylogenetically structured along temperature gradients, we used a dataset of 2882 macrofungal species, assembled from citizen science and fungarium data from eight European countries and covering 40 years from 1970 to 2010 (Andrew et al. 2017). We used a previously established mega-phylogeny and applied a null model approach to test the following four hypotheses ( Fig. 1): 1) phylogenetic alpha diversity shows a hump-shaped relationship with annual mean temperature and 2) a linear decrease with temperature seasonality. 3) Phylogenetic beta diversity increases with increasing annual mean temperature differences (distance) and 4) increasing seasonal temperature differences (distance).

Species observational data
This study utilized data from the ClimFun meta-database, a source of unified, quality-controlled, multi-source data which itself originated from many independent data repositories of fungal fruit body records across Europe (Andrew et al. 2017). The data are comprehensive in temporal and spatial coverage, extending back decades and across a large geographic range in northern, western and central Europe (Andrew et al. 2017(Andrew et al. , 2018. For this study, we selected national-scale data from the ClimFun meta-database from those countries with substantial data (i.e. Austria, Denmark, Germany, Netherlands, Norway, Slovenia, Switzerland and the UK) across a climatically-relevant timespan of 1970-2010. Species were restricted to the macroscopic fruit body forming Agaricomycotina (the classes Agaricomycetes, Tremellomycetes and Dacrymycetes, and removing the Cystofilobasidiales and Trichosporales). Other taxonomic groups comprised very low proportions of the dataset (Andrew et al. 2017), and are biased at this spatiotemporal scale, hence they were omitted. Fungal species records were summarized across 50 × 50 km grids covering the countries mentioned above, and utilizing the UTM coordinate system (zone 32).

Environmental data
To address our hypotheses, we compiled data on annual mean temperature and temperature seasonality. From previous studies, we know that mean and seasonal precipitation, soil organic carbon, nitrogen deposition and tree Figure 1. (A) Spatial distribution of mean temperature (°C) and temperature seasonality (amount of temperature variation within the year based on the standard deviation of monthly temperature averages) across our study area in Europe. (B) Expected response of phylogenetic alpha diversity to mean and seasonal temperature. If mean and seasonal temperature acts as a habitat filter, we hypothesize a lower phylogenetic diversity than expected in extreme temperature conditions within fungal assemblages. (C) Expected response of phylogenetic beta diversity to mean and seasonal temperature. If specialization to temperature extremes (warm, cold, pronounced seasonality) is evolutionarily conserved, we hypothesize a linear relationship between phylogenetic diversity and temperature. diversity correlate with the spatial distribution of fungal species (Andrew et al. 2018, Yu et al. 2021. To account for possible offset effects, we considered these variables as covariates. The environmental data were averaged for each 50 × 50 km grid cell, based on values linked directly to each species record. Thus, the reported means for the 50 × 50 km grids are based on the actual fungal data points found within each grid and therefore reflect the abundance of a species per location. This approach captures more precise environmental conditions for each grid as it is weighted by the presence of fungal occurrences within the grid (Andrew et al. 2017(Andrew et al. , 2018. Collinearity among all environmental variables was low (|r| < 0.70), and no variables were removed from the dataset for this reason (Dormann et al. 2013).
Temperature and precipitation data were obtained from WorldClim (representative of 1960-1990 bioclimatic variables, 2.5-0.5 min original resolution, <www.worldclim.org/ current>). We used the variables bio1 (mean annual temperature, °C), bio4 (temperature seasonality, amount of temperature variation within the year based on the standard deviation of monthly temperature averages, SD), bio12 (mean sum of yearly precipitation, mm) and bio15 (precipitation seasonality, variation in monthly precipitation totals over the course of the year calculated as the ratio of the standard deviation of the monthly total precipitation to the mean monthly total precipitation, i.e. coefficient of variation, %). Soil organic carbon (%) of the topsoil was obtained from the European Soil Data Centre (ESDAC) (1 km original resolution, <http://esdac.jrc. ec.europa.eu/content/octop-topsoil-organic-carbon-content-europe>). Nitrogen deposition data (NHx per hectare and year) were accessed through the greenhouse gas management database (GHG Europe, FP7, Reduced and oxidized downward nitrogen deposition velocities, 1850-2010, CTM data: <www.europe-fluxdata.eu/ghg-europe/data/others-data>). NHx was selected for further analyses, being a proxy for NHx and NOy due to high collinearity (Andrew et al. 2018). Host tree data were obtained from EU-Forest (Mauri et al. 2017). This source provides occurrence records at 1 × 1 km resolution for over 200 tree species across Europe. We calculated the number of host tree genera per grid from this data set and used the host tree genera community matrix for further analyses. We used host tree data at genera level since possible specialization of fungi to their hosts has been suggested to take place at higher taxonomic ranks (Peay et al. 2015, Krah et al. 2018. Note, however, that the number of host tree genera and the number of host tree species per grid is closely correlated (r = 0.99 based on a log-log scale). Besides using the number of host tree genera and the host tree genera community composition, we also calculated the mean pairwise distance (MPD) of the host tree genera within (alpha analyses) and among (beta analyses) grids based on the host tree genera phylogeny. The host tree genera phylogeny was built based on the function 'tol_induced_subtree' in the add-on R package 'rotl' (Michonneau et al. 2016, OpenTreeOfLife et al. 2019. Exploring these phylogenetic-related measures as predictors did not change ecological inferences of the models, and we, therefore, present only the results based on the number of host tree genera and the host tree community composition. To account for spatial autocorrelation, we additionally considered longitude and latitude for each grid as covariates in the model (below).

Phylogeny
To calculate phylogenetic alpha and beta diversity, we used a previously published phylogeny based on the fungal dataset used in Krah et al. (2019), which followed a mega-phylogeny construction protocol (Krah et al. 2018). In brief, we used five gene regions (28S and 5.8S rRNA, rpb1, rpb2, tef1), with gene partitioning, a comprehensive backbone guide tree, and a column reliability score and used RAxML tree inference software. We further conducted 1000 approximate Shimodaira-Hasegawa likelihood ratio tests to assess branching support (SH-aLRT branch support) using RAxML (flag -f J). The final phylogeny used for analyses consisted of 2882 fungal species. We produced an ultrametric tree of the phylogeny using penalized likelihood estimation of divergence times as implemented in the R function 'chronos' in the addon R package 'vegan' (Oksanen et al. 2020).

Phylogenetic alpha and beta diversity metrics and null model approach
To make inferences about the relationship between phylogenetic alpha and beta diversity and temperature variables, adequate diversity measures and a rigorous null model approach have been recommended (Ulrich and Gotelli 2013). We, therefore, calculated the mean pairwise distance (MPD) and mean nearest taxon distance (MNTD) between each species within a sample (alpha) and among samples (beta). While MPD uses the full ancestral range of species across the phylogenetic tree, MNTD is sensitive to the changes of lineages close to the phylogenetic tips (Webb et al. 2002). In the case of phylogenetic beta diversity, MPD is comparable to other metrics like the phylogenetic Sørensen index (Bryant et al. 2008). We calculated the standardized effects sizes (SES) for MPD and MNTD on both alpha (function 'ses.mpd' and 'ses.mntd', add-on R package 'picante', Kembel et al. 2010) and beta level (function 'ses.comdist' and 'ses.comdistnt', Stegen et al. 2012). To calculate the SES, species were randomized across the tips of the phylogeny (null model 'taxa. labels'). SES quantifies the number of standard deviations by which the observed MPD or MNTD differs from the mean of the null distribution (100 randomizations, note that using a test-wise higher number of randomizations produced the same results). SES values > 2 indicate that coexisting species within (alpha) or among (beta) samples are more distantly related than expected by chance (phylogenetic overdispersion). In contrast, values < −2 indicate that coexisting species within (alpha) or among (beta) samples are more similar than expected by chance (phylogenetic clustering). As we were particularly interested to interpret effects more or less different from a random expectation, we used the SES metrics for all analyses.

Models
For the analyses, we considered only grid cells with more than 50 species. This number reflects a minimum sampling effort per grid cell and guarantees robust communities (Andrew et al. 2018, Krah et al. 2019. Furthermore, analyses were carried out based on presence-absence data of species per grid cell. Our final data set consisted of 2882 species and 582 grid cells. To approach our study questions, we fitted phylogenetic diversity measures versus the mean temperature and temperature seasonality or their distance matrices for alpha and beta diversity, respectively. As phylogenetic alpha level response, we used the standardized effect size (SES) of the mean pairwise (MPD) and mean nearest taxon (MNTD) distance. According to our expectation (Fig. 1B), the response of the phylogenetic alpha measure to the temperature variables can either follow monotonous or non-linear curves. At the alpha level, we, therefore, used the framework provided by the generalized additive models (GAMs, functions 'gamm', package 'mgcv', Wood 2017). The predictors of interest (annual mean temperature, temperature seasonality) and the covariates (mean sum precipitation, precipitation seasonality, soil organic carbon, nitrogen deposition and number of host tree genera) were considered as smooth terms in the models. Temperature seasonality, mean sum of precipitation and soil organic content were log 10 -transformed to reach a Gaussian distribution. We also used the spatial surface based on the coordinates (latitude and longitude) of the grids as an additional term ('s(x,y)').
At the beta level, we used the phylogenetic dissimilarity matrix based on the standardized effects sizes (SES) of the mean pairwise (MPD) and mean nearest taxon (MNTD) distance as responses. We furthermore subjected each environmental predictor and the spatial coordinates (x, y) to distance matrices (Euclidean distances) calculations. For the host tree genera community matrix, we used Bray-Curtis as the distance measure. We finally standardized the distance matrices to zero mean and unit variance (function 'standardize' in the package 'robustHD') to make the coefficients comparable. We then ran multiple regressions on distance matrices for each response matrix using the function 'MRM' within the 'ecodist' package (Goslee and Urban 2007). Note that within the figures, we display the raw (untransformed) data.
For illustrative purposes, we calculated the percentages of taxonomic ranks (species, genera, families, order) to extreme thermal conditions. We used a priori cut-offs based on visualization of the histograms (mean annual temperature: cold < 5°C, warm > 9°C; temperature seasonality: continental seasonality > 66 SD, maritime seasonality < 56 SD). For this purpose, we used 2870 out of 2882 species with information on the full taxonomy (genera, families, order), i.e. excluding those characterized by insertae sedis at the family or order level.

Phylogenetic alpha diversity
Mean temperature and temperature seasonality had the strongest effects on phylogenetic alpha diversity based on MPD and MNTD, compared to all other covariates (Table 1). As hypothesized, models based on MPD and MNTD revealed a hump-shaped relationship with mean temperature indicating phylogenetic clustering, i.e. more phylogenetically similar lineages within assemblages than expected by chance in warm and cold environments. However, contrary to hypothesis 2, models revealed a positive relationship with temperature seasonality, with overdispersion in strongly seasonal environments, i.e. more phylogenetically dissimilar lineages than expected within assemblages (Fig. 2). Relationships between our alpha diversity measures and the other covariates are displayed in the Supporting information.

Phylogenetic beta diversity
Phylogenetic beta diversity based on MPD was mainly explained by mean temperature distance (Table 1), showing a significant positive trend (Fig. 2). The slope of the relationship was weak (Fig. 2), indicating a share of lineages in cold and warm environments throughout Europe. The model of phylogenetic beta diversity based on MNTD revealed the largest effect size for temperature seasonality distance (Table 1). As for the model based on MPD, the slope of the relationship was weak (Fig. 2), indicating a share of lineages in environments with different extents of temperature seasonality. Relationships between our beta diversity measures and the other covariates are shown in the Supporting information.

Percent of species, genera, families and orders related to extreme temperature conditions
Eight percent of all species occurred exclusively in cold (< 5°C) or warm (> 9°C) grid cells, respectively, and 4% in strong seasonal (continental) environments (Fig. 3, Supporting information). Associations with extreme temperature environments were found mainly at genus level but rarely at family and absent at order level. The majority of species in extreme temperature environments belonged to families and genera with a broad temperature range, i.e. more than 85% of the species related to extreme temperature conditions belonged to families and genera that included species occurring outside the extreme environments. However, 5% of all genera occurred exclusively in cold, 6% exclusively in warm and 9% exclusively in strong seasonal (continental) environments (Fig. 3). Finally, 2% of all families were exclusively related to warm and 4% exclusively to strongly seasonal environments. No family occurred exclusively in cold environments (Fig. 3, Supporting information).

Discussion
By considering phylogenetic relatedness among species, we found that mean and seasonal temperature are important drivers explaining the structuring of macrofungal fruiting assemblages across Europe. In support of our first hypothesis, we found that both warm and cold environments act as filters constraining the phylogenetic diversity of fungal communities across temperate and boreal Europe. In contrast, we found no support for hypothesis 2; temperature seasonality did not reduce phylogenetic alpha diversity. A positive relationship between phylogenetic beta diversity with mean and seasonal temperature distance supports hypotheses 3 and 4. However, the weak relationship indicates that most coexisting species in extreme temperature habitats are recruited from phylogenetically diverse lineages across the phylogeny. Nevertheless, some species-poor families and genera appear exclusively related to extreme temperature environments (Supporting information, Fig. 3).

Phylogenetic alpha diversity
We observed phylogenetic clustering within assemblages at both ends of the mean temperature gradient (cold and warm environments), indicating niche selection for lineages adapted to such conditions. The effect seems slightly more pronounced at the cold than the warm end of the gradient, and this is in line with studies of butterflies (Pellissier et al. 2013) and plants (Li et al. 2014). Cold temperatures, therefore, could act as an environmental filter that allows a limited number of macrofungal species from similar lineages to coexist within assemblages (see Supporting information for results related to the species diversity).
We expected that strongly seasonal temperatures might also act as a filter similarly to temperature extremes at both ends of the mean temperature gradient. Seasonal and hence continental habitats are characterized by more extreme temperatures (cold and hot) and a higher temperature variability than mild temperature conditions in maritime areas (Xu et al. 2013). However, our results showed a strong positive response of phylogenetic alpha diversity with increasing seasonality (Fig. 2). This effect appeared to be even more pronounced than effects of mean temperature, as indicated by the effect sizes (Table 1). Thus, within areas characterized by less pronounced seasonality (maritime conditions), phylogenetic diversity within assemblages was lower than expected (clustered), but grid cells characterized by a high level of seasonality (continental conditions) had larger phylogenetic diversity than expected (overdispersed). Furthermore, species richness was also positively related to temperature seasonality (Supporting information), in line with findings of Andrew et al. (2019). We can only speculate about the underlying mechanisms. One explanation might be that habitats with a more pronounced seasonality are characterized by a higher number of temperature niches summed across the year that allow more species to coexist from different phylogenetic lineages (MacArthur and Wilson 1967). In support of this explanation, experimental studies on fungal community assembly have shown species richness to increase under variable rather than stable microclimatic (moisture) conditions (McLean and Huhta 2000). However, there is also broad support for climate variability (at small temporal scale, e.g. intraannual) causing decreases in species richness and community stability at different spatial scales (Zhang et al. 2018). More generally, intra-annual climate variability (seasonality) has Table 1. Model output from the generalized additive models (GAMs) of phylogenetic alpha diversity (standardized effect size of the phylogenetic mean pairwise -MPD-distance and mean nearest taxon distance -MNTD within grids) as response. For these models, we present the estimated degrees of freedom (edf), the F-values and the adjusted R 2 s. Model outputs of the multiple regressions on distance matrices (MRMs) using phylogenetic beta diversity (standardized effect size of the phylogenetic mean pairwise -MPD-and mean nearest taxon distance -MNTD among grids). For MRM models, we present the coefficients (coef) and the R 2 s. The predictors of interest, according to our research questions, are grey shaded.  Effects of mean temperature and temperature seasonality on phylogenetic alpha diversity. Alpha diversity plots represent partial effects from generalized additive models (GAM). Shaded areas represent the fitted smooth ± 2 standard errors. (E-H) Effects of mean temperature and temperature seasonality on phylogenetic beta diversity. Beta diversity plots are based on raw scatterplots and linear regressions of the distance matrices. Note that the 95% confidence level is not visible.

Model
been suggested as a hypothesis explaining large scale diversity pattern (e.g. latitudinal-diversity gradient, LDG) (Pielou 1975). However, this hypothesis was rejected by Currie (1991) and Rohde (1992). Instead, climate instability at a geological time scale is among the hypotheses which are still under discussion to explain the LDG (Jansson and Dynesius 2002). However, the spatial scale of our study seems too small to contribute to the knowledge explaining the LDG. Further studies are needed to shed more light on the mechanistic response of fungal diversity to climate variability.

Phylogenetic beta diversity
We found a weakly significant positive response of phylogenetic beta diversity based on MPD with mean temperature distance but not based on MNTD. MNTD is more sensitive towards shallow nodes of the phylogenetic tree, while MPD is more sensitive towards deeper nodes (Webb et al. 2002). Our finding, therefore, indicates an early adaptation of fungal lineages to mean temperature extremes. However, it is important to note that the strength of the relationship between species beta diversity versus temperature distance far exceeds the strength for phylogenetic beta diversity, indicating a stronger species-than phylogenetic turnover (Andrew et al. 2018, Supporting information). This finding, therefore, indicates that a broad range of lineages can occur at a wide range of temperatures, and contrasts with several plant studies showing a high phylogenetic turnover in beta diversity with temperature distances (González-Caro et al. 2014, Qian et al. 2014. A weak relationship of phylogenetic beta diversity with mean temperature distance is further supported by our additional exploration of the distribution of species, genera, families and orders in extreme thermal environments. We found that no order was exclusively related to the coldest or warmest conditions (Fig. 3). Moreover, most species associated with cold or warm environments belong to families and genera with broad temperature niches. However, the significant relationship between phylogenetic beta diversity and mean temperature distance was driven by a few families and genera, exclusively reported from the coldest or warmest environments. It is important to focus more closely on these taxa with psychrophilic capacities in future studies so as to learn more about climate adaption mechanisms. Adaptions might act at the genetic (e.g. genome or gene duplication) and phenotypic level (e.g. morphology of mycelia, fruit bodies and spores) (Yusof et al. 2021). Recent studies at the morphological level have shown that the size and the color of the fruit bodies of assemblages are correlated with the thermal environment. For example, mean fruit body size of assemblages are small in cold environments supporting the heatup-cool-down hypothesis (Bässler et al. 2021). Further, mean fruit body color lightness of assemblages are darker in cold environments supporting the thermal-melanism hypothesis (Krah et al. 2019). Finally, fungi isolated from the Arctic and Antarctic showed increased intracellular trehalose and polyol concentrations and unsaturated membrane lipids as well as secretion of antifreeze proteins and enzymes active at low temperatures suggesting physiological adaption to cold conditions (Robinson 2001). However, our knowledge of the breadth of genetic and phenotypic adaptations to cold environments including all fungal modules (mycelia, fruit bodies, spores) is rather limited. Our models revealed phylogenetic beta diversity based on MPD to be significantly related to mean annual temperature distance. In contrast, phylogenetic beta diversity based on MNTD was significantly related to temperature seasonality Below the bold frame, is the percent of species exclusively related to moderate mean temperatures (5-9°C). (B) Percent of species exclusively related to extreme continental (> 66 SD) temperature seasonality (bold frame). Above the bold frame, is the percent of species exclusively related to mild maritime seasonality (< 56 SD). Below the bold frame, is the percent of species exclusively related to moderate seasonality (56-66 SD). We used 2870 out of 2882 species for this exploration with information on the full taxonomy (genera, families, order). differences. Therefore, deeper and thus larger clades might reflect niche conservatism and basic biochemical or morphological adaptations (Cordero and Casadevall 2017, Krah et al. 2019, Luo et al. 2021 to be able to grow and fruit at different mean temperature environments, while shallow clades seem characterized by species with phenotypic plasticity to cope with varying temperature environments. For example, for maritime coastal angiosperm tree assemblages, phylogenetic beta diversity based on MNTD showed a significant relationship with the environment but not based on MPD, suggesting niche conservatism and dispersal limitation of more derived (recent) lineages explaining species assembly (Massante and Gerhold 2020). This underpins the importance of considering both MPD and MNTD in phylogenetic beta diversity studies to illuminate the role of evolutionary history in explaining the diversity pattern observed today.

Climate change implications
Our aim was to increase our understanding of macrofungal community processes along temperature gradients (based on fruit bodies). Therefore, we cautiously suggest implications for future assemblages and ecosystem processes. 1) With an increase in mean annual temperature, the observed environmental filter at the cold edge of the gradient might, at least partly, diminish and fungal assemblages shift towards more dissimilar lineages. This process might result in a loss of cold-adapted species in more northern latitudes or at higher elevations. Examples of species related to cold environments, which might be vulnerable to further warming, are Amylocystis lapponica, Phellopilus nigrolimitatus and Lentaria byssiseda (Supporting information). Alternatively, these species may be able to track climate warming by shifting northwards and to higher elevations if suitable habitat conditions exist. Predictions of the biological response at the warm end of our gradient are speculative since our study does not cover the complete climate gradient available within Europe due to data limitations in southern Europe. More generally, we expect species that have narrow temperature niches to be particularly vulnerable to rising temperatures. 2) It has been suggested that with increasing temperature in winter, seasonality will diminish in cold regions (Xu et al. 2013). Changed seasonality may lead to a loss of cold-adapted species and consequently to lower phylogenetic alpha diversity.
3) The phylogenetic beta analysis revealed that few genera are exclusively related to cold environmental conditions and thus may be expected to be particularly vulnerable to climate change, similar to what is expected for tree species (Shooner et al. 2018). The clades observed in our study outside cold extremes might expand further as climate change continues. This has already been demonstrated for plant and animal species (Scheffers et al. 2016). 4) Fungal communities consist of members characterized by different lifestyles related to important ecosystem processes (e.g. primary production, carbon and nutrient cycling via mutualist and decay fungi). Several studies have demonstrated significant relationships between fungal diversity and ecosystem processes (Fukami et al. 2010, Kahl et al. 2017, Averill et al. 2019. These studies suggest that changes in fungal assemblages caused by climate change might also affect ecosystem processes at a large scale. Even though speculative, a substitution of cold-adapted by warm-tolerant species may cause increased mycelial growth (Damialis et al. 2015) and thus decomposition rates (A' Bear et al. 2014), leading to increased CO 2 emission via enhanced fungal respiration, potentially further increasing the global greenhouse effect.

Cautionary notes
1) Expectations of a change in fungal assembly processes caused by climate change might fall short when based only on environmental envelopes ignoring the ability of macrofungi to adjust their physiology to the changing conditions, i.e. to acclimate either via physiological adaption or by rapid evolution. Here, we suggest further experimental studies and the use of macroecological frameworks to account for evolution in predictive models (Diniz-Filho et al. 2019). 2) Our analyses were based on fruit body inventories but their presence does not perfectly reflect the presence and extend of mycelium. Thus, although our data set is extensive in space and time (Andrew et al. 2017), the mycelial community might not respond in exactly the same way to temperature gradients as indicated by our current results. Distribution, extend and species composition could now be determined by extracting DNA from soils. However, gathering molecular data, particularly at the same grain size (and extent), is not easy to achieve. Nonetheless, using both approaches in the future would allow deeper insights into fungal community dynamics in times of climate change. 3) Finally, the patterns we observed may not only be attributed to thermal selection but also dispersal limitation. Dispersal ability is expected to differ considerably among species (Golan and Pringle 2017), but it is unclear which species in our data set might be dispersal limited. Peay et al. (2012) suggested community-wide dispersal limitation at the kilometre scale for ectomycorrhizal fungi, while Komonen and Müller (2018) concluded that wood-inhabiting fungi are not dispersal limited at landscape scale. Whether dispersal limitation contributes to explaining the diversity pattern at our scale is unknown. However, if we consider long time scales, dispersal limitation might become less important, as demonstrated by Geml et al. (2012). We, therefore, suggest that the fungal assemblage structure is related to the thermal environment.

Conclusions
Significant responses of the effect sizes of phylogenetic alpha and beta diversity based on null models to temperature mean and seasonality suggest that fungal assembly processes are deterministically structured by the thermal environment at the European scale. Our data further suggest that many species across lineages have evolved strategies to cope with extreme temperatures. In contrast, only some speciespoor genera seem to occur exclusively in extreme thermal environments. If temperature regimes (mean and seasonal temperature) change further as predicted by climate models, we thus expect a re-organization of macrofungal assemblages. To further understand the response of macrofungal diversity to climate gradients in space and time, and the effects of this on ecosystem functioning, the following two questions particularly need to be addressed in future studies: 1) Which physiological (biochemical and morphological) adjustments and adaptations allow a species to cope with climate extremes, focusing on both fruit bodies and the mycelium? 2) Which ecosystem consequences, e.g. carbon and nutrient cycling, could be related to a loss of cold-adapted species in times of climate change?