The fall of the summer truffle: Recurring hot, dry summers result in declining fruitbody production of Tuber aestivum in Central Europe

Abstract Global warming is pushing populations outside their range of physiological tolerance. According to the environmental envelope framework, the most vulnerable populations occur near the climatic edge of their species' distributions. In contrast, populations from the climatic center of the species range should be relatively buffered against climate warming. We tested this latter prediction using a combination of linear mixed effects and machine learning algorithms on an extensive, citizen‐scientist generated dataset on the fruitbody productivity of the Burgundy (aka summer) truffle (Tuber aestivum Vittad.), a keystone, ectomycorrhizal tree‐symbiont occurring on a wide range of temperate climates. T. aestivum's fruitbody productivity was monitored at 3‐week resolution over up to 8 continuous years at 20 sites distributed in the climatic center of its European distribution in southwest Germany and Switzerland. We found that T. aestivum fruitbody production is more sensitive to summer drought than would be expected from the breadth of its species' climatic niche. The monitored populations occurring nearly 5°C colder than the edge of their species' climatic distribution. However, interannual fruitbody productivity (truffle mass year−1) fell by a median loss of 22% for every 1°C increase in summer temperature over a site's 30‐year mean. Among the most productive monitored populations, the temperature sensitivity was even higher, with single summer temperature anomalies of 3°C sufficient to stop fruitbody production altogether. Interannual truffle productivity was also related to the phenology of host trees, with ~22 g less truffle mass for each 1‐day reduction in the length of the tree growing season. Increasing summer drought extremes are therefore likely to reduce fruiting among summer truffle populations throughout Central Europe. Our results suggest that European T. aestivum may be a mosaic of vulnerable populations, sensitive to climate‐driven declines at lower thresholds than implied by its species distribution model.

Fungal DNA from the soil cores was extracted using the PowerSoil MOBIO DNA 1

Truffle fruiting and tree growth season calculation 1
To calculate the lengths of the growing and fruiting seasons for host trees and truffles, S-shaped 2 Gompertz functions were fit to the relationship between the numerical day of the year vs tree 3 diameter and cumulative fruitbody production, respectively (Čufar et al., 2008). The Gompertz 4 function is expressed as: 5 6 eq (1) 7 Where a, b, and c are parameters that describe the maximum (a) and slope and timing (b and c) 8 of tree growth or fruitbody production. These three parameters were fit for each site x year x 9 tree combination (and site x year for truffle fruitbodies) using the least sum of squares method 10 in R (nls function in the base package). The start, end, and point of max growth and 11 fructification were calculated from these fit parameters using linear approximation of the 12 Gompertz function. Accordingly, 13 14 Start day: eq (2) End day: eq (3) Max growth rate: eq (4) Max growth day: eq (5)

15
For host trees, each individual tree was fit a different Gompertz function (for a maximum for 4 16 per site per year) and treated as a technical replicate. The average start and end days and max-17 growth-rates of all trees were computed as potential predictors of annual productivity. 18 1 + log

Variance Partitioning 1
For each site x year combination, we summed the annual total number of fruitbodies per unit 2 area (m -2 yr -1 ). We then calculated mean truffle mass as the ratio of total annual mass (g yr -1 ) 3 over the total number of truffles (# of truffles yr -1 ). This allowed us to express truffle yield (z, 4 in g m -2 yr -1 ) as the product of fruitbody yield (x, in # of fruitbodies m -2 yr -1 ) and mean truffle 5 mass (y, in g). In order to correct for the differences in scale between the variables x and y 6 (measured in different units with multiplicative effects on z), we log-transformed the product z 7 = x×y into the sum log(z) = log(x)+log(y). This allowed us to partition the variance in log(z) 8 into the sum of two covariance terms (for x and y), such that 9 When both covariance terms are positive (as they are in our analysis), then years with a greater 11 truffle yield (z) tend to have both a greater number of truffles (x) and a greater mean truffle 12 mass (y). This means that the ratio of the x and y covariance terms over var [ log (z) ] gives the 13 proportion of the variance accounted for by x (truffle numbers) and y (mean truffle mass), 14 respectively. 15

Variables and Model Selection 16
A large suite of predictors was assembled and exposed to model selection using the random 17 forest algorithm. These were plotted in Figure 4a. We provide a full list of these variables along 18 with the units and data source in Table S1. 19

Principal Component Analyses 20
For vegetative and soil meta-community variables principle component analyses (package 21 "prcomp" in R) were performed to collapse the high dimensional datasets into three orthogonal 22 axes (Supplemental Figures S1,S2). 23 For soil fungal meta-communities, there were 2,219 OTUs (after rarefaction and 24 removal of singletons). This makes visualizing loading values with vectors on top of the PCA labels. In Figure S2 (top row) we characterized the top 5 positive and negative loading values 1 on the first 3 PCA axes. 2 For the dendrometer data, we conducted a separate principal component analysis. We 3 included data for temperature and precipitation during the hottest quarter, mean tree growth 4 during the hottest quarter, and mean tree water deficit during the warmest quarter (Figure 6a). 5 . 6 Supplemental Table 1. All variables used in the random forest analysis of annual truffle yield (g m -2 yr -1 ). Variables are sorted according to 1 the rank of their Inc Node Purity importance scores (as in the x-axis of Figure 3)