Forest management and former land use have no effect on soil fungal diversity in uneven-aged mountain high forests

Key message Metabarcoding analysis of soil fungal communities in French mountain forests revealed that harvest‑ ing intensity, time since last harvest and former land use had no effect on fungal community composition compared to key abiotic factors. Low‑intensity management in these uneven‑aged mountain forests therefore has limited effects on soil fungal community composition which is mainly driven by elevation and edaphic properties. Context Past and current human activities are known to affect forest biodiversity. However, the effects of former land use and forest management have been studied much more extensively on higher plants than on fungi. Aims Our objectives were to assess the effects of harvesting intensity, duration since last harvest and former land use on soil fungal communities in uneven‑aged mountain high forests. Methods On the basis of historical land‑use maps drawn between 1862 and 1864 and on historical forest manage‑ ment archives, we selected 62 sites in the French Alps with contrasting land‑use histories (ancient forests, which were already forested on historical maps vs recent forests


Introduction
Both former human activities and current forest management are known to affect forest biodiversity because of dispersal and recruitment limitation or habitat modifications (Nordén et al. 2014;Paillet et al. 2010).With the extinction of many species worldwide (Ceballos et al. 2017), we need to better understand the effects of past and current human activities on forest biodiversity (Janssen et al. 2018).
Fungi are recognized as fundamental components of biodiversity and ecosystem functioning (Mueller and Bills 2004), and biologists generally consider that unmanaged old-growth forests are key habitats for this taxonomic group (Kjučukov et al. 2022).Indeed, forest management can influence fungal communities by altering stand structure or soil conditions.For instance, harvesting may cause changes in vegetation, nutrient availability, soil microclimate and structure and litter quantity and quality (Hartmann et al. 2012).Forest cover and tree basal area (two parameters strongly influenced by forest management) can also affect the composition and richness of soil fungal communities by modifying microclimate and soil conditions (Santos-Silva et al. 2011;Spake et al. 2016).In addition, depending on harvesting conditions, logging can produce a variable amounts of residual ground woody debris, which can locally modify soil conditions and below-ground microbial communities (Peršoh and Borken 2017;Perreault et al. 2020) and may mitigate the effects of tree removal (Mayer et al. 2022).On the other hand, if forest management is stopped for a long time, forest succession is no longer interrupted, and the stands gradually evolve into old-growth forests, which are characterised by large amounts of standing and lying deadwood and high diversity of tree dimensions and ages (Fuhr et al. 2022).Spake et al. (2015) showed that ectomycorrhizal fungal species richness recovered to its old-growth level 90 years after management abandonment, but other authors found contradictory results, i.e. that ectomycorrhizal richness was higher in managed than in unmanaged stands, though saprotrophic fungi were favoured in unmanaged forests (Dvořák et al. 2017).
While the effects of forest management on fungal communities are well documented, few studies have focused on the effects of former land use (Nordén et al. 2014).Indeed, the effects of former agricultural activities on soil properties can persist for a long time and remain noticeable even after forest recolonization and canopy closure (Dupouey et al. 2002).These effects may induce recruitment limitation in recent forests.For instance, Diedhiou et al. (2009) showed that ancient Roman occupation partly explained the current composition of ectomycorrhizal communities in an oak forest.In addition, some soil fungi may have a low dispersal capacity, which could make them dependent on the temporal continuity of the forest, i.e. the age of the forested state of the land (Bergès and Dupouey 2021).In this regard, Boeraeve et al. (2018) showed that communities in recently restored forest stands were more similar to those in ancient forest stands when they were bordered by ancient forest than when they were isolated.In addition, Mennicken et al. (2020) and Flensted et al. (2016) showed that historic forest cover was a better predictor of soil fungal diversity and red-listed species richness than current forest cover at the landscape scale, indicating a time lag effect of more than 150 years.There is also evidence that macrofungal diversity is correlated with the diversity of plant species associated with long forest continuity (Hofmeister et al. 2014).
Soil fungal diversity is also strongly influenced by soil properties, climate and stand conditions (Leclerc et al. 2023), and any analysis must take these effects into account (Janssen et al. 2018).Furthermore, differentiation of fungal lifestyles is necessary, as they may respond differently to environmental predictors (Finlay and Thorn 2019; Leclerc et al. 2023).For instance, saprotrophic fungi might be more affected by the abundance of lying woody debris, while symbiotrophic fungi, which depend on the presence of trees, might be more affected by harvesting and forest continuity (Dvořák et al. 2017;Collado et al. 2020).
In their review of the effects of forest management on fungal communities, Tomao et al. (2020) identified four knowledge gaps that we have attempted to fill in our study.They point out that few studies have focused on the following: (1) forests besides temperate or boreal forests, (2) soil fungi using environmental DNA, (3) the duration of the effects of anthropogenic disturbances and (4) the effect of close-to-nature silvicultural practices on soil fungal communities in mountain areas.Specifically, we assumed that forests with different land-use histories would have different soil fungal community structures, due to dispersal and recruitment limitations, but that logging might interact with former land use: the shorter the time since the last harvest, the weaker the signal from former land use.Based on a metabarcoding analysis of soils sampled in the French Alps, we assess the respective effects of harvesting intensity, time since last harvest and former land use on the taxonomic and functional composition of below-ground fungal communities in uneven-aged high mountain forests while taking into account soil conditions, climate, stand structure and deadwood volume.

Study area
The study area is located in the public forests of the Tarentaise valley in the French Alps (Fig. 1) and has a mountain climate with a continental influence characterised by cold winters and cool summers (for the period 1981 to 2000, the average temperature in January is −1.8 °C and in July 16.6 °C; the annual rainfall is 984 mm).The geology of the Tarentaise valley is characterised by eroded sediments of the Hercynian orogenic belt (black shales and micaceous sandstones), partly covered by more recent glacial deposits.The soils are dominated by umbric leptosols and hypereutric cambisols (IUSS Working Group WRB 2022) according to the regional soil map of Savoy (Pary 2016).The current forests are predominantly coniferous, dominated by Picea abies (L.) H. Karst.( 1881) and Abies alba Mill.( 1768), and are being managed for lumber production and hazard protection as uneven-aged high forests.In the mid-nineteenth century, the forests covered 10% of the landscape, but the abandonment of agriculture and reforestation policies have led to some forest recovery.Today, forests cover 19% of the study area (Thomas et al. 2017).

Sampling design
We stratified our sample according to former land use and time since last harvest (divided into four classes at 20-year intervals, Fig. 1).Current forest cover (according to the FAO definition of forest, FAO 2012) was extracted from the BD FORET ® V2 map (scale 1:25 000), based on aerial photographs taken in 2014 by the French National Geographic Institute.Former land use was identified using the Ordnance Survey maps (scale 1:40 000, surveyed in the study area between 1862 and 1864), which are used as a reference in France to identify forests with long temporal continuity (Bergès and Dupouey 2021).These ordinance maps were vectorised and georeferenced by the Vanoise National Park (Thomas et al. 2017).They allowed us to identify which present-day forests were already existed in 1864 (hereafter called "ancient forests", AF) and which present-day forests developed on former agricultural land since the mid-nineteenth century (hereafter called "recent forests", RF).We selected our sample sites in public forests located on the north-facing slopes of the valley and managed by the National Forest Service.This institution archives all the silvicultural operations Fig. 1 Study area and distribution of sampling sites among ancient and recent forests and dates of the last harvest in the Tarentaise valley (French Alps).Biogeographical regions defined according to the EEA (2016).Points in black are located in ancient forests; points in white are located in recent forests undertaken, and the Vanoise National Park had previously digitised all these data for our study area.We were therefore able to select sites with contrasting times since last harvest, from 1 to 75 years ago.
We used QGIS 3.16 to identify potential sampling areas.We first selected conifer-dominated stands according to BD FORET ® V2 and then used aerial photographs from 1958 to restrict our sample to sites with closed canopies.This allowed us to exclude very recent forests as well as ancient forests that may have undergone a landuse change between 1850 and the time of the study.To eliminate areas of uncertain land use (errors in historical mapping or vectorisation), a 20-m buffer was removed from the edge of each polygon before final selection of the sampling sites.Within the selected polygons, 151 points were randomly placed, according to the stratification described in Table 1.As we focused on public forests, all recent forest sites were located on former common land used as pasture in the nineteenth century, and there were no former cropland or hay meadow in our sample.In addition, all sites were located on leptosols.
These sites were then visited during June 2021 to select sites where (1) the state of stump decomposition corresponded to the date of the last harvest mentioned in the management plan archives, (2) the canopy cover was above 40% and (3) the soil was homogeneous (i.e. the entire plot had the same conditions in terms of stoniness, soil depth and moisture).This resulted in a final selection of 62 sites.

Stand description
In a circular plot with a radius of 20 m, we recorded the species and the diameter at 1.30 m of all living trees with a diameter greater than 17.5 cm.The diameter, length and decay stage of each piece of lying woody debris over 20 cm in diameter were measured throughout the plot.Five stages of wood decay were estimated based on resistance to knife penetration 1: hard or unaltered; 2: rot < ¼ of the diameter; 3: rot between ¼ and ½ of the diameter; 4: rot between ½ and ¾ of the diameter; and 5: rot greater than ¾ of the diameter.
In addition, we ran three radial transects starting from the centre of the plot and separated from each other by an angle of 120°, and we measured the diameter of each piece of small lying woody debris with a diameter between 7.5 and 20 cm that intersected the transect.Warren and Olsen (1964) developed this method to estimate the volume of small lying woody debris using the formula as follows: where V SWD is the volume (in m 3 /ha) of the small woody debris, L the total length (in m) of the three transects and d i the diameter (in m) of each piece of woody debris intersecting the transect.
These measurements were then used to calculate the stand basal area (G), the mean square diameter of living trees (Dg), the volume of lying coarse woody debris (V CWD ), the volume of well decayed lying woody debris (decay stage above 3, V CWD.WD ) and the volume of lying small woody debris (V SWD ).
We also recorded the circumference of each stump (c 0 ) and converted this to a circumference at 1.30 m (c 1.30 ) using regression Eq. 1 below.We characterised harvesting intensity by the basal area of the previously harvested trees (G harv ) and by their mean square diameter (Dg harv ).
The regression coefficient used in Eq. 1 was estimated using c 0 and c 1.30 extracted from data collected by the French National Forest Inventory on Picea abies in the forests of the French Northern Inner Alps (sylvoecoregion H22, IFN 2022) and by fixing the intercept at zero (n = 1582 trees, R 2 = 0.99).

Soil fungal communities assessment
We sampled soils from the 62 selected sites and surveyed soil fungal communities using environmental DNA (eDNA) metabarcoding (Thomsen and Willerslev 2015).This method is sensitive to spatial and temporal variation, so it is recommended to adopt broad spatial coverage in eDNA sampling, either by sampling several points within sites or by combining all the samples into one large sample representing the entire site (Grey et al. 2018).As our sites were located on leptosols with high stoniness and low thickness, we focused only on the topsoil fungal communities.Therefore, we collected 30 soil samples (at a depth of 10 cm after humus removal) taken every 2 m along two orthogonal diameters of the plot and mixed them together to make them homogeneous.Humus was excluded from the sample because it can be acidic and very rich in organic matter (MOR type) in this region, which can be a problem for DNA extraction in the laboratory, as DNA can be degraded or adsorbed onto organic particles and organic matter can clog analyser filters.Then, a subsample of 200 cm 3 of (Eq. 1) c 1.30 = 0.73 c 0 soil was taken and stored in a sterilized and closed bottle along with 40 g of desiccant silica gel sachets, which were changed regularly to dry out the soil as quickly as possible.Samples were then stored in a dry, dark place for 3 months prior to metabarcoding analyses.To avoid cross-contamination, the corer sampler was sterilized by a flame, and the gloves worn during all the manipulations were changed before changing sites.In addition, in order to assess within-site variability, 10 of the sites were randomly chosen to replicate soil sampling.In these 10 sites, the 2 orthogonal diameters previously materialized were rotated by 45°, the corer was sterilized and 15 new soil samples were collected every 4 m on the new orthogonal transects of the plot.Again, 200 cm 3 of soil were subsampled and stored in a second sterilized and closed bottle with silica gel.We assessed variability in soil fungal composition due to soil sampling by comparing between-and within-site dissimilarities using a Sørensen-type dissimilarity index of order 1 as recommended by Alberdi and Gilbert (2019) (see details below Appendix Fig. 4).Soil extracellular DNA was isolated in October 2021 following a published protocol based on the Nucle-oSpin ® Soil kit (Macherey-Nagel) (Taberlet et al. 2012).Fungal communities were assessed with the Fung02 primer pair (forward 5′-GAA GTA AAA GTC GTA ACA AGG-3′ and reverse 3′-CAA GAG ATC CGT TGY TGA AAGTK-5′), which targets the nuclear ribosomal internal transcribed spacer region 1 (ITS1) (Epp et al. 2012;Taberlet et al. 2018).All samples were homogenized and amplified by PCR in quadruplicates as recommended by Alberdi et al. (2018).Each reaction included a total volume of 20 μL containing 1 U of AmpliTaq Gold Polymerase (Applied Biosystems), 2 mM of MgCl 2 , 0.2 mM of each dNTPs, 0.5 μM of each forward and reverse primer, 0.2 mg/mL of bovine serum albumin (BSA, Roche Diagnostics) and 2 μL of extracted DNA.The associated PCR programme involved an initial denaturation step of 3 min at 95 °C followed by 40 cycles of 30 s at 95 °C, an annealing of 30 s at 56 °C, an elongation step of 30 s at 72 °C and a final elongation step of 7 min at 72 °C.Library construction and sequencing for these samples were performed at Fasteris (http:// faste ris.com, Geneva, Switzerland).Libraries were prepared following the Metafast PCR-free protocol (www.faste ris.com/ en-us/ NGS/ DNA-seque ncing/ Metab arcod ing/ Metag enomi cs-16S-18S-ITS-or-custom-PCR-ampli cons) so as to limit common sequencing artefacts such as tag jumps (Schnell et al. 2015).This proprietary protocol is analogous to the PCR-free library protocol published by Carøe and Bohmann (2020) but replaces some consumables by those of the Illumina library preparation kits.The amplified samples were then sequenced in forward and reverse on a MiSeq-v3 device (2 × 250, 15 m reads, Illumina, San Diego, USA).The DNA sequence reads were then filtered and clustered into molecular operational taxonomic units (MOTU) by the eDNA-AnaEE platform, hosted at the Laboratoire d'Ecologie Alpine (Grenoble, France), and following an established pipeline using the dedicated OBITools software suite and the metabaR R package (Boyer et al. 2016;Zinger et al. 2021).For the following biodiversity analyses, we opted for a "relaxed restrictive" PCR strategy (Alberdi et al. 2018).Only MOTUs detected in more than one of the four PCR replicates were retained, and the number of reads in each PCR replicate was averaged.To ensure that the diversity recorded was not influenced by sequencing depth, the sampling completeness of each site was assessed using rarefaction curves (Appendix Fig. 5).
We used Hill numbers of order 1 to assess fungal diversity, which is similar to the exponential Shannon index ( 1 D ≈ e − R i=1 p i ln (p i ) , with p i = number of reads of the MOTU i divided by the total number of reads and R = the total number of MOTUs retained in the sample) (Hill 1973).This index was calculated for all retained MOTUs (total diversity) and for six fungal lifestyles (ectomycorrhizae, wood saprotrophs, litter saprotrophs, soil saprotrophs, other saprotrophs and pathotrophs).These fungal lifestyles were assessed for MOTUs whose genus could be assigned to a lifestyle using the FungalTrait database (Põlme et al. 2020), which is reported to be superior in quality and coverage to the FunGuild trait database (Tanunchai et al. 2022).
Even if the number of reads cannot be related to true abundance or biomass (Lamb et al. 2019), using the number of reads to estimate p i makes the metabarcoding analysis less sensitive to PCR or identification errors (Alberdi and Gilbert 2019).Calderón-Sanou et al. (2020) recommend using a diversity order between 1 and 2 when studying associations between diversity and environmental variables, as these indices yield more robust results than did 0 D (= richness).Here, we set q = 1 to study the relationship between environmental variables and fungal diversity, as 1 D and 2 D were strongly correlated (Appendix Table 5).
In addition, chemical analyses were carried out on the remainder of each soil sample at the INRAE soil analysis laboratory (Arras, France).These analyses included the available phosphorus content extracted by the Duchaufour method (hereafter called P, Duchaufour and Bonneau 1959), the total organic carbon content extracted by sulfochromic oxidation (hereafter called C org , ISO standard 14235) and the pH using a glass electrode in a 1:5 (volume fraction) suspension of soil in either water (hereafter called pH H2O , ISO standard 10390).

Statistical analyses
We used linear regressions to test the effect of time since last harvest and former land use and the interaction between them on total MOTU diversity and the diversity of each fungal lifestyle while controlling for stand condition (G, Dg), deadwood volume (V CWD , V CWD.WD , V SWD ), climate (elevation, which is a good integrator of local variations in temperature and precipitation), landscape context (distance to the nearest forest edge), soil conditions (P, C org and pH H2O ) and harvesting intensity (G Harv and Dg Harv ) (Table 2).Stand characteristics and some soil variables may be related to harvesting or past land use so we analysed their correlations in Appendix Table 6 to ensure the absence of confounding effects and collinearity between predictors.
We fitted four linear models to predict 1 D for total diversity and for each fungal lifestyle: a null model (M 0 , with none of the predictors), an interaction model (M Inter , with only the main effects of time since last harvest and former land use and the interaction between them as predictors), a main model (M Main , same as M Inter without the interaction term) and the most parsimonious model (M Pars , which included the most explanatory variables among all control variables and variables of interest, Cf.Table 2).To compute this parsimonious model, we used a forward stepwise model selection algorithm (step function, Venables and Ripley 2002) to select only significant variables in the model.Finally, we compared these four models (M 0 , M Main , M Inter and M Pars ) based on their Akaike information criterion (AIC).If time since last harvest or former land use affect the biodiversity index after accounting for significant control variables, they will be selected in M pars .Comparing M 0 and M main or M inter allows us to test whether time since last harvest or former land use has a significant effect without the control variables.In this way, we ensure that the effect of harvesting or former land use is not mediated by the effect of other environmental predictors.For each model, spatial autocorrelation was checked using Moran's I index calculated on the residuals.
We studied the taxonomic composition of soil fungal communities by applying a canonical correspondence analysis (CCA, cca function, Legendre et al. 2012) to the MOTU abundance matrix, using the logarithm of the number of reads.We first calculated a CCA with only time since last harvest and former land use only (called CCA Main ), then a second CCA with the interaction between Last.Harv and Form.Use (called CCA Inter ) and a parsimonious CCA (called CCA Pars ) that included the most explanatory variables among the potential control variables, Last.Harvest, Form.Use and their interaction (Table 2).We selected the best subset of variables using permutations tests in a forward stepwise model selection algorithm (ordistep function, Blanchet et al. 2008).Significance of each CCA was assessed using permutation tests (anova.ccafunction, Legendre et al. 2011).For CCA Pars , variance partitioning was carried out to assess the unique and shared variance components of each group of selected predictors (varpart function).These analyses were performed using the vegan R package (Oksanen et al. 2020).MOTUs occurring only at one or two sites (n = 496 out of a total of 849 identified at the genus level) were discarded from this analysis.All the data used in this study are available in an online repository (Mollier et al. 2023).

Results
The MiSeq sequencing run of fungal ITS1 DNA yielded a total of 14,296,158 reads for the 72 samples (62 samples + 10 replicates).Eighty-two percent of these reads presented a quality score (Q score) equal or higher to 30 (i.e. less than 1 error per 1000 bases).Excluding controls, the processing and filtering of this dataset resulted in 1,773,967 reads that were consistently assigned to fungal species, with an average of 77,600 reads per sample.These reads were included in the subsequent analyses and clustered into 5933 molecular operational taxonomic units (MOTUs) assigned to 1097 distinct taxa.Among them, only 849 were identified at the genus level and were subsequently analyzed.Accumulation curves show that diversity coverage is good, with most samples reaching a plateau (Appendix Fig. 5), especially for 1 D, confirming that this index is less sensitive to sequencing depth than richness ( 0 7).Soil sample replicates from the 10 randomly selected sites showed that the average Sørensen-type dissimilarity within-site was 0.42 ± 0.12 (sd) and between-site dissimilarity was 0.66 ± 0.12 (Appendix Fig. 4).Twenty-eight percent of the MOTUs were detected in only one of the two soil sample replicates.
We found no evidence that time since last harvest or former land use were strongly correlated with the control variables.Time since last harvest was positively correlated with G and Dg and negatively correlated with V SWD , G harv and Dg harv , consistent with the expected effects of harvesting, but no correlation coefficient exceeded 0.6.In addition, there were no differences in edaphic and stand conditions between ancient and recent forests (Appendix Table 6).

MOTU diversity
On average, 88.3 MOTUs were found in each site.Ectomycorrhizae were the most diverse group in the community, followed by soil and litter saprotrophs (Table 3).
Total MOTU diversity ( 1 D) was not influenced by any environmental factor.Neither former land use nor harvesting intensity nor time since last harvest affected the diversity of any fungal lifestyle (Table 4).The diversity of ectomycorrhizae and soil saprotrophs was mainly influenced by soil conditions and elevation, while distance to the forest edge and harvesting intensity positively affected the diversity of wood saprotrophs, and stand structure influenced litter saprotrophs and pathotrophs (Fig. 2).Moran's I index showed no spatial autocorrelation of residuals for any of the models.

Taxonomic composition
Neither CCA Main nor CCA Inter were significant (p = 0.19 and p = 0.27, respectively), and the variable selection algorithm did not select harvesting intensity, time since last harvest or former land use in CCA Pars which only included elevation, pH, C org and P (p = 0.001).
The environmental variables selected in CCA Pars explained 11.4% of the total variability.Of this explained variability, soil conditions (pH, C org and P) accounted for 84.1%, elevation accounted for 7.6% and the shared effect of soil and elevation amounted to 8.3% according to the variance partitioning.
Only the first three axes were significant in CCA Pars .The first two axes correspond to a gradient of pH and elevation, respectively, while the last one reflects a gradient of organic carbon content (Fig. 3).This analysis confirmed the results observed for the diversity indices (Fig. 2), as ectomycorrhizae tended to be associated with high pH and high phosphorus content, while soil saprotrophs tended to be associated with high pH and elevation.In contrast, pathotrophic fungi tended to be associated with more acidic soils with higher soil organic carbon than ectomycorrhizae and soil saprotrophs.Ectomycorrhizae also tended to be associated with lower elevation than soil saprotrophs, a tendency which was not observed with diversity indices.88.3 (80.3, 95.7) 26.7 (23.5, 29.8) 5.5 (4.7,6.3) 21 (19.4, 22.6) 12.5 (11, 14) 5.8 (4.9, 6.6) 7.3 (6.5, 8.1)

Forest management and former land use do not affect soil fungal communities in uneven-aged high mountain forests
We found no evidence of direct or indirect effects of time since last harvest or former land use on the diversity of MOTUs, either on total diversity or on the diversity of main lifestyles.Indeed, the M Main and M Inter models were never better than the M 0 models .In addition, time since last harvest and former land use were never selected in M Pars .
Our results are consistent with Dove and Keeton (2015), who also found no effect of either individual tree selection or group selection harvesting on soil fungal communities in uneven-aged high forests in Vermont.In contrast, some studies in even-aged high forests managed by clearcutting have reported negative effects of harvesting on soil fungal communities, especially for symbiotrophic and saprotrophic fungi (Hartmann et al. 2012;Parladé et al. 2019).However, other authors found that a thinning operation in even-aged high forests had no effect, suggesting that fungal communities may be resilient to tree removal when harvesting is less intensive (Castaño et al. 2018).Thus, the effect of harvesting on soil fungal diversity may depend on the management type and harvesting intensity.
We initially hypothesised that recruitment limitation due to soil modification could alter community composition in recent forests, but we did not find any difference in soil properties between ancient and recent forests.This is probably because pastoralism has been the dominant agricultural activity in mountain areas since the early Middle Ages (Mouthon 2019).Therefore, all of our Table 4 Summary of the Akaike information criterion (AIC) for the null model (M 0 ), the model including the main effects of Form Use and Last.Harvest (M Main ), the model including the interaction between them (M Inter ) and the most parsimonious model including Form.Use, Last.Harvest and all the potential control variables as potential predictors (M Pars ).Predictors have a significant effect if the AIC of the model is lower than the AIC of M 0 (shown in bold).The P-value corresponds to the significance of the spatial autocorrelation according to a Moran's I-test performed on the residuals of the models.The expected Moran's I was −0.02 recent forest sites were located on former pastures.This type of former management had little effect on soil properties compared to former arable management, where fertilisation and ploughing would have enriched the soil in nutrients and modified the original forest soil structure, contributing to soil homogenisation (Koerner et al. 1997;Holmes and Matlack 2017;Abadie et al. 2018;Burst et al. 2020).Therefore, our results confirm a previous study we carried out in the same study area (Mollier et al. 2022), which showed that pastures had a weaker legacy effect on soil properties compared to cropland.We also hypothesised that dispersal limitation might have determined the fungal diversity and composition in recent forests.However, in our study area, forest recovery since the mid-nineteenth century has mostly been due to a steady expansion of the forests on the periphery of historic forest patches so that most recent forests are in strong continuity with ancient forests; this might have reduced the effects of dispersal limitation (Boeraeve et al. 2018).In addition, the dispersal capacity of fungi remains largely unknown, and some species are likely to produce large numbers of propagules that can disperse over long distances, and therefore, fungi may be less affected globally by forest continuity than are vascular plants (Nordén et al. 2014).

Fungal lifestyle
Our results contribute to a growing body of literature showing that the effects of past land-use and forest management on understorey communities are largely context dependent (Janssen et al. 2017;Depauw et al. 2019;Mollier et al. 2022).It would be interesting to investigate these effects under a wider range of historical and environmental conditions, in order to understand which factors control this context-dependent effect.For instance, it may be interesting to investigate the interaction between legacy effects of past land-use and forest management in lowland areas, where the effects of past and current human activities are expected to be stronger (dominance of former cropland and higher harvesting intensities).In particular, plantations and tillage may have a stronger effects on soil properties and soil fungal communities than the low-intensity management practised in the mountain forests of our study area (Bergès et al. 2017).

Soil fungal communities are mainly influenced by abiotic factors
The most parsimonious models included only abiotic factors.These results are consistent with the literature, as pH, soil moisture and nutrient levels are known to strongly influence fungal species composition (Taylor and Sinsabaugh 2015;Leclerc et al. 2023).Some studies have also suggested that soil temperature contributes to differences in fungal communities (Timling et al. 2014;Taylor and Sinsabaugh 2015), which may explain the effect of elevation found in our study.In addition, we showed that total phosphorus content positively affected the diversity of ectomycorrhizal fungi.This result contradicts Khalid et al. (2021), who showed that high phosphorus concentrations had a negative effect on the richness of symbiotrophic fungi; however, the phosphorus content gradient they sampled was 10 times larger than in our study.Fungal communities are known to be strongly structured by soil horizon and depth (Taylor and Sinsabaugh 2015), and we did not account for this variability by sampling only the first 10 cm of the topsoil.As the samples were located on stony leptosols, it was difficult to sample deeper layers under these conditions.In addition, we did not analysed the humus layer, which can be thick and harbour specific fungal communities that we were unable to characterise.It would be interesting to investigate the effect of harvesting and former land use on more evolved soils, in order to study fungal communities in deeper layers, and by taking into account humus-inhabiting fungi.

Conclusion
Our initial hypotheses were not confirmed, as we did not detect any effect of time since last harvest or former land use on fungal MOTU diversity and taxonomic composition.Our results suggest that (1) current abiotic factors have a stronger effect than past or present human activities in temperate uneven-aged mountain forests, as previously reported (Jansa et al. 2014;Janssen et al. 2018;Leclerc et al. 2023); (2) former land use, at least pastures, has limited legacy effects on soil and fungal communities in mountain areas, a finding that is likely to be widely applicable, as agropastoralism was the most widespread form of agriculture in the mountains of Europe; and (3) the silviculture, as applied in the uneven-aged high forests of these mountains, is lowintensity and compatible with the conservation of fungal diversity.In order to provide relevant management guidelines, further research is needed to determine the level at which management intensity may affect the soil fungal communities and to test the effects of other former land uses, such as former cropland or former meadows, which are more common in lowland areas.

Appendix
Fig. 4 Distribution of pairwise Sørrensen-type dissimilarity indices between soil replicates at each site and between sites Using the 10 randomly selected sites where soil sampling was replicated, we computed between-and within-site dissimilarity indices.
To assess within-site variability, we computed the alpha ( q D α ) and gamma ( q D γ ) diversity for q = 0, 1 and 2, for each pair of soil samples.We then calculated beta diversity ( q D β = q D γ / q D α ) and the Sørensen-type overlap as recommended by Alberdi and Gilbert (2019).The dissimilarity index was assessed by 1-C qN .The results of the two soil subsamples were then grouped together to calculate the composition of each site, and the same procedure was followed for each pair of sites to assess dissimilarity between sites.We used the nonparametric Kendall correlation coefficient to assess correlations between numerical variables, as some variables did not follow a normal distribution.We calculated a binomial regression and used the square root of a pseudo-R 2 based on the likelihood ratio with the signs of the regression coefficients to assess the sense of the correlations between numerical variables and former land use.The Kendall's test assessed the significance of the correlation for numerical variables, and the LR test assessed the significance of the correlations between forest ancientness and numerical variables.

Fig. 2 Fig. 3
Fig. 2 Effects of environmental variables (±IC 95% ) on 1 D (exponential of Shannon's diversity index) for each fungal lifestyle responding to environmental predictors.Predictors were scaled prior to calculation to provide standardised coefficients

Fig. 5 Fig. 6
Fig. 5 Accumulation curves obtained after 30 resampling iterations and for 10 sample sizes, for 0 D (species richness), 1 D (exponential of Shannon's diversity index) and 2 D (reciprocal of Simpson's diversity index).The plot "coverage" represents the Good's coverage index: 1-(number of singletons/total number of reads)

Table 1
Stratification of soil sampling by former land use and date of the last harvest Numbers without parentheses indicate the number of sites selected for sampling.Numbers in parentheses indicate the number of sites visited before selection

Table 2
Description of the predictors used in the models

Table 5
Pearson correlation coefficient between each order of diversity for total diversity and each fungal lifestyle 0 D species richness, 1 D exponential of Shannon's diversity index, 2 D reciprocal of Simpson's diversity index

Table 6
Kendall coefficient correlations for all the variables included in our study Indicates a significant correlation according to Kendall's test for numerical variables and according to a logistic regression for correlations between forest ancientness and numerical variables.Dist edge distance to the nearest forest edge, G basal area, Dg mean square diameter of living trees, C org total organic carbon, V CWD volume of coarse woody debris on the ground, V CWD.WD volume of well-decayed woody debris on the ground, V SWD volume of small woody debris on the ground, G harv basal area of harvested trees, Dg harv mean square diameter of harvested treesMollier et al.Annals of Forest Science  (2024)81:2 *

Table 7
Full names of fungal species in the dataset.Only the fungal MOTUs identified as species and present in more than 30% of the sites are shown