Tanja Mrak
Nataša Šibanc
Philip Brailey-Jones
Ines Štraus1
Jožica Gričar
Hojka Kraigher- 1Department of Forest Physiology and Genetics, Slovenian Forestry Institute, Ljubljana, Slovenia
- 2Department of Agronomy, Biotechnical Faculty, University of Ljubljana, Ljubljana, Slovenia
- 3Department of Forest Yield and Silviculture, Slovenian Forestry Institute, Ljubljana, Slovenia
Temporal studies that would offer insight into resilience of ectomycorrhizal (ECM) communities in stress prone climates are scarce despite their role in tree nutrition and water supply. Our study characterized the vitality, community composition, diversity, and function of Quercus pubescens Willd. ECM fungi in the Sub-Mediterranean stress-prone environment for 2 consecutive years (June 2016–May 2018) and related the investigated measures to environmental parameters. ECM community was assessed for species actively associating with root tips and exploring the soil volume through the assessment of mycelial ingrowth into sand-filled mesh bags. The investigated period was characterized by a drier than average summer combined with wildfire in 2016 followed by another dry summer in 2017. The vital to non-vital ECM root tip ratio decreased below one in August 2016 and remained low until January 2018. This was ascribed to a series of stress events that occurred at the site including sequential droughts and wildfire. The most abundant ECM lineages on root tips were Tomentella and other Thelephoraceae, Sebacina, and Cenococcum while in mesh bags the most abundant were Tomentella, Sebacina, Pseudotomentella, Pyronemataceae, Inocybe, Cortinarius, Agaricales, and Boletales lineages. High intra-site variability was observed, with ECM communities directly associated with root tips and exploring the soil volume varying significantly among the plots. Community composition was stable over time, while species richness varied with mean air and soil temperature, relative air humidity, and solar radiation. The most abundant exploration type observed at this site was short distance, which was associated with precipitation along with long distance exploration type. The medium distance exploration type was temporally variable and responded to soil temperature and relative air humidity reflecting seasonality at the site. The presented results indicate complex relationships between environmental parameters, abiotic stress, and ECM fungi.
Introduction
Abiotic stresses are the major environmental factors that adversely impact tree growth and consequently forest productivity. They also play a major role in determining the geographic distribution of tree species (Harfouche et al., 2014). As mean global temperatures increase through climate change, the probability of extreme weather events causing heatwaves, droughts, and fires to occur simultaneously or sequentially is expected to rise in the future (Sutanto et al., 2020). This poses risks to the composition, structure, and biogeography of forests in many regions (Allen et al., 2010). Southern and Southeastern Europe are critical regions which could be drastically affected by these relative to other regions in Europe (Lehner et al., 2006; Sheffield and Wood, 2008; Gudmundsson and Seneviratne, 2016). The Slovenian Sub-Mediterranean is already a highly variable region in terms of climate, and this may become more unpredictable under climate change. Trees have evolved several physiological and anatomical mechanisms for coping with these abiotic stresses that may arise throughout the year. In recent years considerable progress has been made in the understanding of how the aboveground parts of trees respond to abiotic stresses. Although ECM fungi are an integral part of a tree, their symbiosis is often overlooked in studies of stress response, and specifically temporal studies that would offer insight into resilience of ECM communities in stress prone climates are scarce.
Trees invest 20–40% of their carbohydrate gain to mycorrhizal fungi in exchange for improved water and nutrient supply (Smith and Read, 2008). Through the ECM extramatical mycelium (i.e., hyphae of ECM fungi that extend beyond the fungal mantle of ECM root into surrounding soil), the symbiosis can improve water status under both ambient and drought conditions through an increased water-absorbing surface area, efficient conduction through mycelial strands, enhanced hydraulic conductivity at the soil-root interface, along with resulting hormonal and nutritional effects modifying stomatal regulation (Bréda et al., 2006). Some ECM fungi may also be beneficial to the wildfire response of host trees. Thelephora terrestris and Tomentella atramentaria have been shown to induce greater shoot weight in aspen following wildfires, while the same effect has been seen in Russula sp. associations in birch (Bent et al., 2011; Beidler et al., 2015).
Extramatrical mycelium is developed to a different extent across ECM fungal taxa, which can be classified into exploration types (Agerer, 2001). ECM fungi of the contact exploration type develop a smooth fungal mantle with only a few emanating hyphae. The short distance exploration type is characterized by a voluminous envelope of emanating hyphae from the mantle. The medium distance exploration type is characterized by the presence of undifferentiated or slightly differentiated rhizomorphs (aggregations of parallel-oriented hyphae), while number of emanating hyphae depends on the sub-type (fringe, mat, and smooth). Finally, the long-distance exploration type has differentiated rhizomorphs which spread throughout the soil up to several decimeters from the roots (Agerer, 2001). This morphological differentiation results in different relative importance of various species in nutrient and water uptake, and their effects on host carbon economy (Koide et al., 2005). It has been shown that the potential enzymatic activities of ECM root tips do not reflect fungal lineage or host species, but instead exploration type (Tedersoo et al., 2012b). Morphological differentiation also affects the ECM community composition observed, depending on whether root tips or soil-exploring mycelia are sampled (Koide et al., 2005; Kjøller, 2006). Fungi belonging to the long-distance exploration type are able to form common mycelial networks across tree hosts, which contribute to the connectivity among trees in forest soils and support the survival of seedlings by transporting water and nutrients from maternal trees to seedlings (Simard et al., 2012).
In the Slovenian Sub-Mediterranean, the adaptability of Quercus pubescens (Willd.) to ecosite characteristics (regions of similar soil characteristics) and inter-annual variability has been extensively studied in the past few years (Gričar et al., 2019, 2020; Vodnik et al., 2019). The region is characterized by two types of bedrock, flysch and limestone which meet in a sharp line along the Kraški rob (Karst edge). This is the geographic and climatic boundary where the limestone-based Karst Plateau region of Slovenia abruptly transitions to flysch terrain. The Karst region is defined by rugged terrain, a shortage of surface waters, and erosion by rainfall and bora wind in conjunction with at least 6000 years of intensive human transformation of local forest vegetation. These natural and anthropogenic disturbances have led to the deterioration of soils and the formation of a stony landscape (Zorn et al., 2015). However, after the second world war, social and demographic changes led to abandonment of pastures in this area and marked the beginning of spontaneous natural afforestation (Zorn et al., 2015) predominantly by Q. pubescens accompanied by Fraxinus ornus L. and Ostrya carpinifolia Scop. (Brus, 2012). Within this ecosite the ecological role of pubescent oak overrides its silvicultural importance due to its protective role against erosion (Brus, 2012). In contrast to the shallow soils sitting on limestone bedrock that are prone to quick water loss through deep percolation (Ferlan et al., 2016), soils sitting on flysch bedrock are deeper and have a greater water holding capacity (Vodnik et al., 2019). This is reflected in a greater height and circumference of Q. pubescens trees of the same age growing just a couple of meters away from each other on the two contrasting ecosites (Gričar et al., 2019). This region is also prone to regular wildfires due to the warm climate, highly permeable limestone bedrock, low precipitation, and high winds (Geršič et al., 2014). This region accounted for 78% of all wildfire affected area in Slovenia between 1995 and 2009. These wildfires are more common in the late Winter and Summer (Šturm and Podobnikar, 2017). Due to climate change even highly drought-adapted tree species will increasingly suffer from extended drought periods (Vodnik et al., 2019). The frequency of wildfires will also be expected to increase in this area through the association with lower than average quantity of precipitation and longer dry periods (Veble and Brečko Grubar, 2016). Drought stress may affect the structure of ECM communities or the abundance of individual taxa (Richard et al., 2011; Mrak et al., 2019) and wildfires have been reported either to cause reduction in diversity/richness of ECM fungi or have no impact, while decreased colonization of roots following wildfire by ECM fungi was confirmed by numerous studies (Taudiere et al., 2017). While Q. pubescens within this area is well-studied, the temporal dynamics of their associated ECM fungi are not. It is important to understand how fungi respond to inter- and intra-annual variation under the current climate in this region to be able to start assessing how they will respond under future climates.
In this study we took 11 cross-sectional measures to characterize the community of ECM fungi colonizing the root tips of Q. pubescens over a 2-year period within the highly variable Sub-Mediterranean climate. This was assessed through the vitality of ECM root tips, species diversity and community composition. The ECM mycelial network colonizing the bulk soil at the site was also assessed through the molecular identification of species colonizing mesh bags placed at the site to characterize the ECM species contributing to connectivity between host trees and the surrounding soil and trees within this Sub-Mediterranean ecosystem. Through this exploratory study we aimed to collect data on the ECM taxa that define the community to create a knowledge base from which we can further investigate the functional roles of such ECM taxa. Our time series assessments were made at the intersection of the two ecosites (flysch and limestone) in the region to maximize the variability present in this region. We selected three plots within these two contrasting ecosites only a few tens of meters away from each other. Two plots were located on the more heterogeneous limestone bedrock and one was located on more homogeneous flysch bedrock to incorporate bedrock with differing water holding capacities of soils in this area. Due to inter- and intra-annual variability of Sub-Mediterranean climate we hypothesized strong dynamics of parameters describing ECM community across this site in relation to environmental parameters. We expected that the extramatrical mycelium community composition will be comprised from a subset of taxa found on the root tips of Q. pubescens.
Materials and Methods
Study Area Description
The study was performed at Podgorski Kras (45°32′56.3″N, 13°54′36.1″E, 430 m a. s. l.), a karst region in SW Slovenia. The site was located at the edge of a plateau where the bedrock transitions from limestone to flysch. This area is characterized by a Sub-Mediterranean climate. The maximum precipitation occurs in the fall and late spring. The region tends to experience two distinct annual low-precipitation periods at the transition from winter to spring and the end of summer in July and August (Ogrin, 1996). ECM fungal communities inhabiting the root tips of Q. pubescens and the surrounding soil were assessed over a 2-year period, from June 2016 to May 2018. Long-term eddy covariance measurements in the area were established previously in 2008 and have shown Q. pubescens in the area as an important factor affecting soil carbon sequestration rates (Ferlan et al., 2011, 2016).
Three adjacent sampling plots of roughly equal size (15 x 18 m) were randomly selected in the study area, approximately 50 m apart from each other, one of them was on flysch (plot 1) and two on limestone (plot 2 and plot 3) (Supplementary Figure 1). Soil in the flysch plot was characterized as eutric cambisol- a soil with higher water retention capacity and depth of the soil profile > 0.6 m (Vodnik et al., 2019). Soils within the two plots on limestone were characterized as rendzic leptosol- shallow soil (0–0.5 m, but up to several meters in soil pockets) with up to 45% stoniness in the upper 0.5 m of topsoil and, consequently, very low water retention capacity (Vodnik et al., 2019). Both soils have a clay texture, with 7.8% of sand, 33.5% of silt, and 58.8% of clay (Vodnik et al., 2019). Soil chemical parameters for the three investigated plots are given in Supplementary Table 1. The mean annual soil water content in 2015 and 2016 was estimated as 45 m3 m−3 in both years for soils overlying flysch bedrock and 33 and 30 m3 m−3, respectively, for soil overlying limestone bedrock. Mean soil water content measured throughout the growing seasons (April–September) of 2015 and 2016 was 44 and 43 m3 m−3 for soil overlying flysch bedrock, respectively, and 32 and 24 m3 m−3 for soil overlying limestone bedrock, respectively (Vodnik et al., 2019). Q. pubescens trees within the three sampling plots were of similar age (55 ± 5 years), but the mean circumference at breast height differed significantly (30 ± 4.5 cm for flysch plot and 20 ± 1.5 cm for limestone plots (Gričar et al., 2019).
The 25-year average precipitation (1992–2017) and monthly precipitation throughout the sampling period (Supplementary Figure 2) was calculated from Slovenian Environment Agency (ARSO) weather data collected from the nearby Kozina weather station, while the 25-year mean temperature was averaged from Portorož, Postojna, and Godnje weather stations (ARSO, 2020). The 25-year (1992–2017) average precipitation for the area was 1,299 mm, and the average temperature was 11.8°C, with a mean maximum temperature of 27.9°C in July and a mean minimum temperature of −0.6°C in January (ARSO, 2020). Soil matrix potential and soil temperature (Supplementary Figure 3) were measured at the three plots using MPS-2 dielectric water potential sensors (Decagon devices Inc., USA) at 10-min intervals from May 2016 to June 2018. In addition to soil matrix potential and temperature, the mean precipitation, air temperature, relative air humidity, and solar radiation were calculated for the 30-day period before each sampling date (Supplementary Table 2).
2016 and 2017 were characterized by drier than average summers. There was also a high-spread low intensity wildfire that occurred at the site from 7 to 10 August 2016. The wildfire damaged the understory aboveground biomass and caused significant scorching of Q. pubescens crowns (Gričar et al., 2020). The fire effects on vegetation were not uniform. The fire caused little or almost no damage on flysch (plot 1). Trees had no visible heat damage of the stem bark or crown, and the herb and shrub layers were not affected by the fire. On limestone, the forest stand is very sparse, therefore certain parts of the sampled plots suffered little or no damage, whereas in other parts the understory aboveground biomass was significantly damaged. In the damaged part of the plot, the fire caused significant scorching of crowns of pubescent oak extending up to 2 m above the ground. In addition to this, the leaves up to a height of ca. 7 m dried and died but were mostly not burnt and remained on the trees. Small twigs were also damaged. Based on drone aerial photo estimation, severely burnt trees suffered more than 95% of leaf-loss. In these pubescent oaks, re-spouting of new shoots occurred in September 2016, after the first autumn rains (Gričar et al., 2020).
Sampling, Morpho-Anatomical, and Molecular Identification of Root-Tip ECM Fungi
Q. pubescens root tips were sampled at each plot every 2 months when the soil was not frozen providing altogether 11 cross-sectional measures (22 June 2016, 9 August 2016, 9 November 2016, 28 March 2017, 16 June 2017, 22 August 2017, 27 October 2017, 28 November 2017, 30 January 2018, 28 March 2018, 30 May 2018). Soil cores were collected under the Q. pubescens trees. For each sampling date five sub-samples were collected from each of the three sampling plots.
Cores were sampled by cutting into the topsoil (up to 20 cm depth) with a knife to extract soil, as the stony terrain prevented the use of a sampling probe with a definite volume. Each soil core was submerged into a measuring cylinder filled with tap water to obtain the exact volume of the sample. Subsequently, soil cores were left to soak in water overnight to enable easier separation of soil particles from roots. From each soil core Q. pubescens roots were isolated and divided under an Olympus SZH stereomicroscope (Olympus, Tokyo, Japan) into vital ECM root tips and non-vital ECM root tips. For vital ECM root tips all turgescent ECM tips were considered, while for non-vital ECM root tips all non-turgescent and dry ECM root tips were considered. Vital ECM root tips were separated into morphotypes according to Agerer (1991).
For each morphotype identified per soil core, mantle preparation was carried out, and anatomical characteristics determined under the Zeiss Axio Imager Z2 light microscope (Carl Zeiss Microscopy, Jena, Germany) using differential interference contrast following the protocols defined by Agerer (1991) and DEEMY1 (Agerer and Rambold, 2004-2020). ECM roots of each morphotype were scanned on an Epson scanner (Seiko Epson Corp., Suwa, Nagano, Japan) in a tray filled with water. This was also carried out separately for non-vital ECM and non-mycorrhizal roots. From each morphotype found per soil core DNA was extracted from one mycorrhizal system with DNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following manufacturer's instructions and sequenced using ITS1/ITS4 amplicons (Gardes and Bruns, 1993; Grebenc and Kraigher, 2007). Sequencing of pure DNA was performed at a commercial sequencing laboratory (Macrogen Inc., Seoul, South Korea). Forward and reverse sequences were assembled using Geneious Prime software 2020.0.4. and identified taxonomically using the UNITE (Version 8.2) database (Abarenkov et al., 2020). Representative sequences were deposited in GenBank under accession numbers MW027850-MW028132. Molecular identification was cross-checked with morpho-anatomical identification to ensure correct assignments. Exploration types were assigned to individual taxa based on Agerer (2001), DEEMY (Agerer and Rambold, 2004-2020), Suz et al. (2014), and our own image data and morphotype descriptions.
Scans were quantified using WinRhizo software (Regent Instruments Inc., Ville de Quebec, Canada) to obtain the number of root tips of each morphotype per sub-sample. Data were normalized to 250 mL sampling volume as this was the mean volume of collected samples. ECM fungi inhabiting tree root tips were found sparsely across sub-samples taken within each sampling plot and sampling date. Sub-samples from within each sampling plot at each date were therefore summed together to create aggregate samples that could be used as a replication of “plot,” “date,” and “year” for community composition analysis, while they were kept separate for other analyses to capture the full intra-site variation. ECM fungal normalized abundances were converted to relative abundances per sample prior to any further community analysis.
Sampling and Molecular Identification of Soil-Inhabiting Extramatrical ECM Fungi
To assess which ECM fungi were present and played an important role within the extramatrical mycelium across plots at the site, in-growth mesh bags (Wallander et al., 2001) were used. Extramatrical mycelium is comprised of emanating hyphae and rhizomorphs. Therefore, community composition of extramatrical mycelium includes only ECM fungi that form these structures and shows which ECM fungi form connections with surrounding soil and comprise to common mycelium networks among trees. Ingrowth mesh-bags were installed into the soil of the three plots over two periods during the experiment. The first set of bags were installed on the 21 June 2016 and removed on the 14 July 2017. The second set of bags were installed on the 9 November 2016 and removed on the 28 November 2017. Bags were left in the soil over these periods to allow adequate time for fungal mycelium to colonize them.
The triangular mesh bags (10 x 10 x 14.5 cm) were prepared from inox 304L woven metal cloth with pore openings of 46 μm diameter, 50 μm wire diameter, and filled with 25 g of silicate sand. The granulation of silicate sand varied over the two exposure periods. For the first exposure period this was 0.063–1.0 mm (Kema, Puconci, Slovenia) and for the second exposure period this was 0.71–1.25 mm (Aquasil Filtersand, Quarzwerke, Melk, Austria). The silicate sand used was first autoclaved, then washed with 2.5 M HCl to remove organic matter and subsequently with ultrapure water to remove any traces of the acid. The removal of organic matter was intended to reduce the possibility for ingrowth of decomposing saprotrophic fungi into the mesh bags by removing a potential nutrient source. The mesh bags were labeled with plastic-coated wire that was attached to one corner of the triangle. The 46-μm diameter of the pore openings in the mesh allows only hyphal ingrowth, excluding plant roots (and therefore directly root-associated fungi) from sampling. It is presumed that ECM fungi make up the majority of hyphae colonizing such a system with minor contributions from other saprotrophic hyphae based on previous comparisons of the isotopic signatures (δ15N or δ13C) of mesh-bag colonizing hyphae with those of known fruiting bodies known to belong to ECM and saprotrophic fungal species (Wallander et al., 2001, 2004; Bakker et al., 2009).
Mesh bags were inserted into soils of the same plots that were used for ECM root tip assessment. They were inserted into the soil at ~10 cm depth by lifting intact bulk soil with a spade and placing it back over the bags to minimize the disturbance of the soil layers. At each site, 18 mesh bags were inserted under the Q. pubescens trees, and for each of two sampling periods two mesh bags were used which were later pooled to obtain enough mycelium biomass. Once mesh bags were collected, they were kept frozen at −20°C until further processing. Sand with mycelium from inside of the mesh bags was mixed with distilled water from which mycelium floating on the water was picked out with tweezers and transferred into an Eppendorf tube. Afterwards, the collected mycelium was freeze-dried, pooled into one sample as previously described and homogenized using MillMix20 mixer mill (Domel, Železniki, Slovenia).
To quantify the total fungal community within the soil mycelium, DNA was extracted from up to 0.1 g of homogenized mycelium using DNeasy Plant Pro kit (Qiagen, Venlo, Netherlands), following the manufacturer's instructions. Fungal communities were quantified using Illumina MiSeq NGS sequencing of ITS2 region amplicons. To produce amplicon libraries for Illumina MiSeq NGS, the ITS2 fragment was first amplified by PCR using Q5® Hot Start High-Fidelity 2X Master Mix (New England BioLabs, Ipswich, Massachusetts, USA) and the primer pair ITS7f (Ihrmark et al., 2012) and ITS4r (White et al., 1990). Forward and reverse primers were modified to contain Illumina specific overhang adapter sequences. PCR was carried out in a 25 μl reaction volume with 2 μl of DNA template, 12.5 μl of Master Mix, and 1.25 μM of each primer. PCR conditions were 98°C for 30 s followed by 31 cycles at 98°C for 10 s, 57°C for 20 s, and 72°C for 20 s and final elongation at 72°C for 2 min on an Applied Biosystems Veriti Thermal Cycler (Thermo Fisher Scientific, Waltham, Massachusetts, USA). PCR products were purified using Agencourt AMPure XP magnetic beads (Beckman Coulter, Brea, California, USA). A secondary semi-nested indexing PCR was then used to attach Illumina sequencing adapters and multiplex indexes, using the Nextera XT Index Kit (Illumina, San Diego, California, USA) following Illumina's recommended protocols. Secondary PCR products were purified using Agencourt AMPure XP magnetic beads, before quantification using a PicoGreen dsDNA Assay Kit (Thermo Fisher, Waltham, Massachusetts, USA). Equimolar concentrations (4 nM) of successfully amplified samples were pooled, and the pooled library was quality checked with Bioanalyzer High Sensitivity DNA Analysis (Agilent, Santa Clara, California, USA). Sequencing was performed on the Illumina MiSeq (2x300 cycles, 5% PhiX).
Quality filtering of sequence reads and clustering into operational taxonomic units (OTUs) was performed as previously described by Mahnic et al. (2020). After quality filtering, we removed reads that were represented in an abundance of less than 5 sequences per OTU, or samples and OTUs where the taxonomy was not determined to at least the family level. We obtained a total of 227,508 sequences and 785 OTUs. Taxonomic classification of representative sequences was inferred using the UNITE version 8.2 ITS database (Abarenkov et al., 2020). Ecological guilds and trophic modes were also determined using the FUNGuild database (Nguyen et al., 2016). Two trophic modes assigned by FUNGuild were combined prior to analysis. These were pathotroph-saprotroph-symbiotroph and saprotroph-pathotroph-symbiotroph, as the assignments were composed of the same three trophic modes. FUNGuild assignments were used to assess the success of the mesh bags at capturing only ECM fungi, and subsequently to subset only the OTUs sequenced which belonged to the ECM functional guild for community analysis. Representative sequences of each OTU were deposited in GenBank under accession numbers MW039835-MW040018. Prior to further analysis the dataset was converted to relative abundances per sample.
Statistical Analysis
Statistical analyses were performed in R Version 3.6.3 (R Core Team, 2020) unless otherwise stated.
Fungal community composition was assessed through permutational multivariate analysis of variance (PERMANOVA) of Bray-Curtis dissimilarities between samples using the “vegan” R package (Oksanen et al., 2019). Where PERMANOVA analysis was found to be significant, post-hoc testing was conducted through pairwise PERMANOVA testing with Benjamini-Hochberg corrections for multiple comparisons. ECM root tip communities were compared between plots, years, sampling dates, and the interaction of plot and year. This was carried out for summed pooled samples. The interaction of plot and sampling date could not be assessed as each pooled plot x date sample was one replicate. The soil mycelium communities developed in mesh bags were compared between plots only, as the two sampling dates overlapped considerably in time.
The abundance of ECM mycelium OTUs assigned to each trophic mode and ECM root tip taxa/morphotypes assigned to each exploration type was compared among plots for trophic modes, and among sampling date, plot, year and plot x year interactions for explorations types using Kruskal-Wallis testing with Dunn's post-hoc testing and Benjamini-Hochberg corrections. For analysis of the effects of plot and sampling date on the ratio of vital to non-vital, and mycorrhizal to non-mycorrhizal roots Kruskal-Wallis was applied using Statistica Version 13.5.0.30 (TIBCO Software, Palo Alto, California, USA).
Alpha diversity metrics species richness (total number of species present), community evenness (inverse Simpsons index) and community dominance (reciprocal of the proportional abundance of the commonest/most abundant species; inverse Berger-Parker index) were estimated for all individual samples with the R package “vegan.” Kruskall-Wallis testing was used to compare alpha diversity of ECM and soil mycelium communities among plot, year, date, and the interactive variables plot x year and plot x date. Where Kruskall-Wallis testing was significant (P < 0.05), Dunn's post-hoc testing was conducted with Bonferroni corrections for multiple testing. Sample coverage, species richness and diversity were estimated across sampling units (plot and date) with the R package SpadeR using incidence-frequency data (Chao et al., 2015).
Both ECM root tip and soil mycelium fungal communities were visualized using non-metric dimensional scaling (NMDS). The relationship between environmental parameters and the ECM root tip fungal community composition was assessed using the “envfit” function of the “vegan” R package. This permutes linear correlations of environmental parameters against Bray-Curtis co-ordinate values for each NMDS axis to determine gradients of change associated with community differences. Envfit analysis was not performed for soil mycelium fungal communities as this was deemed inappropriate since mycelium were collected over a time and not as a snapshot which can be related to a specific environmental measure.
To determine the seasonal effect of environmental parameters on soil root tip species richness, diversity and presence across the experimental time period, general linear models (GLM) analysis was performed using the glm function in R with the family argument set to gaussian. Richness, evenness, dominance, vital-to-non-vital root tips ratio, and exploration type response variables were assessed against the environmental variables of mean temperature, mean soil temperature, mean soil matrix potential, sum of precipitation, relative air humidity, and sum of solar radiation. GLM models were originally conducted with all parameters included. Insignificant variables were stepwise removed from the original model using the step function in R to obtain a model with the lowest AIC value to be used as the final model. GLM analysis was performed on samples which were aggregated by calculating the mean values of the parameters across each plot x date interaction to match the environmental variables available. Response variables were all transformed according to the form of data (log transformations for alpha diversity measures, asin/sqrt transformations for all percentage variables. Due to the high proportion of zero values for the long-distance exploration type, proportions were also +1 transformed. Prior to analysis all response variables were checked for outliers using the boxplot method, resulting in two values being removed from the dataset for the abundance of long-distance exploration type. All model residuals when visually assessed through GLM models were randomly distributed around zero. Five samples were removed from the analysis due to missing environmental data.
Indicator species analysis was conducted to find OTUs which had a significant association with one or a combination of treatment variables. This form of analysis is preferred for multivariate data with large numbers of comparison factors, as multiple pairwise comparisons of OTU abundance are often inappropriate, lack statistical power and can result in multiple false positives without sufficient correction. Indicator species analysis was carried out using the “multipatt” function in the “indicspecies” R package (Cáceres and Legendre, 2009). P-values derived from the analysis were also subjected to Bonferroni corrections. The mean relative abundance OTUs with a significant association to one or more variables was visualized through heatmaps.
Results
Vitality of ECM Root Tips
There was a statistically significant effect of sampling date (H = 68.5, P < 0.0001), but not of plot (H = 5.05, P = 0.0802) on the ratio between vital and non-vital ECM root tips (Figure 1). The ratio of vital to non-vital root tips was above 1 only at the first sampling date in 2016 and the last three sampling dates (from January 2018 on).
Figure 1. Ratio of vital to non-vital ECM root tips of Q. pubescens in Slovenian Sub-Mediterranean in relation to sampling date. Data are not shown for sampling plots separately as there was no significant effect of sampling plot. Means ± stderrs are shown. Different letters mark significantly different ratios at P < 0.05 by Kruskal-Wallis.
ECM Root Tips Species Diversity and Community Composition
Using a combination of morphotyping and molecular analyses we observed altogether 93 ECM taxa-morphotypes. Of the 93 taxa-morphotypes found, 15 taxa-morphotypes belonged to Ascomycota and 77 to Basidiomycota, while one morphotype remained unidentified. Thirty-two genera were defined across 21 families. The most abundant ECM genera (>10% relative abundance) across the field site were Cenococcum (14.07% ± 2.46 SE), Tomentella (18.43% ± 2.32 SE), and Sebacina (17.86% ± 2.85 SE). At the family level, the most abundant lineages were Sebacinaceae (26.02% ± 4.44 SE) and Thelephoraceae (including Tomentella and other genera, 49.72% ± 4.33 SE). In plot 1 (flysch) the most abundant genera were Tomentella, Lactarius, and Sebacina. At plot 2 (limestone), the most abundant genera were Tomentella, Russula, and Cenococcum, along with the unidentified morphotype, that could be ascribed to Tuber based on morpho-anatomical characteristics. At plot 3 (limestone) the most abundant genera were Cenococcum, Russula, Sebacina, and Tomentella (Figure 2). All three plots shared only 16 taxa-morphotypes.
Figure 2. Taxa charts displaying the most abundant genera (>1%) associated with Q. pubescens root tips in the Slovenian Sub-Mediterranean per plot, sampling date and for the whole site (total).
The richest lineage was the Tomentella genus (containing 17 taxa-morphotypes), followed by the Russula genus (10 taxa-morphotypes). The highest number of taxa-morphotypes were detected at plot 1 (flysch). The lowest sample coverage was obtained for plot 3 (limestone) despite the same sampling effort (Table 1). The lowest estimated species richness per plot (Chao 2) was detected for plot 2 (limestone), while the highest was for plot 1 (flysch) (Table 1). Observed species richness, evenness and dominance per sampling unit were all significantly different among the three plots (Kruskall-Wallis richness: χ2 = 12.182, df = 2, P = 0.002; evenness = χ2 = 11.282, df = 2, P = 0.004; dominance χ2 = 9.329, df = 2, P = 0.009). Limestone plot 2 fungal community per sampling unit was significantly richer, more even, and less dominant compared to plot 3 (Table 1).
Table 1. Observed and estimated species richness and diversity (mean ± stderr) of ECM taxa of Q. pubescens root tips for three plots in Slovenian Sub-Mediterranean.
When samples were merged according to sampling dates, the estimated sample coverage per sampling date was lower, ranging from only 0.261 to 0.770 (Table 2). The observed number of taxa per sampling date ranged from 14 to 36 taxa—morphotypes, while the estimated species richness per sampling date (Chao 2) ranged from 20 to 156, showing strong seasonal variability. The lowest values of observed and estimated species richness were detected in August 2017 (Table 2). Only fungal richness per sampling unit was significantly different among sampling dates (Kruskall-Wallis χ2 = 38.743, df = 10, P < 0.001). This was not the case for evenness (Kruskall-Wallis χ2 = 13.921, df = 10, P = 0.177) or dominance (Kruskall-Wallis χ2 = 11.952, df = 10, P = 0.288). Over all three plots the highest richness per sampling unit was observed in May 2018, which was significantly richer compared to fungal community sampled between June 2017 and November 2017, while the lowest richness was observed in August 2017 and November 2017.
Table 2. Observed and estimated species richness and diversity (mean ± stderr) of ECM taxa of Q. pubescens root tips for 11 sampling dates in Slovenian Sub-Mediterranean.
Compositional changes in the root-tip ECM fungal community were assessed through PERMANOVA analysis of Bray-Curtis distances between pooled samples and visualized using non-metric multidimensional scaling (NMDS) of the same distances (Figure 3). PERMANOVA analysis revealed that intra-site variation was high, with the plot of origin being a significant determiner of ECM community composition (PERMANOVA: F = 4.911, df = 2, P = 0.001). The sampling year, date, and interaction between sampling year and plot had no observable impact on community composition (PERMANOVA year F = 1.329, df = 2, P = 0.082; date F = 0.749, df = 8, P = 0.749; plot x year F = 1.162, df = 4, P = 0.155). The NMDS displays a clear separation of plot 1 (flysch) from plot 2 (limestone), while most samples collected from plot 3 (limestone) were clustered between plot 1 and 2. Pairwise PERMANOVA post-hoc testing between plots confirmed the observed differences in ECM fungal community composition among all three of the plots as being significant. “Envfit” analysis was conducted to find environmental variables significantly correlated with the ECM root tip fungal community. These included maximum, minimum and mean values of air temperature, soil temperature and soil matrix potential for the months that samples were collected, along with the mean relative air humidity and the sum of solar radiation and precipitation over the same time period. None of the tested environmental variables were significantly correlated to ECM community composition (Supplementary Table 3).
Figure 3. NMDS ordination of ECM root tip community composition for the three plots in the Slovenian Sub-Mediterranean. Samples were pooled by plot and sampling date. Hulls delimit communities by the plot. Stress of the ordination is 0.196 for k = 2.
Indicator species analysis based on the relative abundances of root tip fungal species revealed six species with significant associations with one or two plots (Supplementary Figures 4, 5). Plot 3 (limestone) was characterized by increased abundances of Inocybe splendens and Tomentella bryophila, and plot 2 (limestone) by a morphotype of Boletales (an order belonging to the Agaricomycetes class) and an unidentified morphotype. Plot 1 (flysch) had no uniquely significantly associated ECM species, but was associated with Sebacina incrustans along with the plot 3 (limestone). This species was present in high relative abundances in plot 1 (flysch) and plot 3 (limestone) (0.26 ± 0.05 and 0.16 ± 0.04, respectively) and nearly completely absent from plot 2 (limestone, 0.02 ± 0.01). C. geophilum was also significantly associated with both plot 2 (limestone, 0.22 ± 0.03) and 3 (limestone, 0.16 ± 0.04) and occurred in very low abundance in plot 1 (flysch, 0.05 ± 0.02). Other indicator species relative abundances and analysis statistics are detailed in Supplementary Table 4.
Exploration Types
Across all root tip samples the short distance exploration type was the prevailing exploration type with 62.8% relative abundance. This was followed by contact (27.6%), medium distance (8.11%), and long distance (1.46%) exploration types. There was a statistically significant effect of plot on relative abundance of long distance (Kruskall-Wallis: χ2 = 6.6446, df = 2, P = 0.0361), contact (Kruskall-Wallis: χ2 = 8.7811, df = 2, P = 0.0124) and short distance exploration types (Kruskall-Wallis: χ2 = 12.066, df = 2, P = 0.0024, Figure 4). The prevalence of medium distance exploration types was similar across the three plots (Kruskal-Wallis: χ2 = 3.1924, df = 2, P = 0.2027). The contact exploration type was significantly more abundant within plot 2 (limestone) in comparison to plot 3 (limestone) and short distance exploration type was more abundant at plot 3 (limestone) compared to plot 2 (limestone) and plot 1 (flysch). Post-hoc testing of long-distance exploration prevalence among plots did not reveal any differences among plots, though it was generally more prevalent in plot 2 than the other plots. There was a significant effect of sampling date on the prevalence of ECM fungi belonging to the medium distance exploration type (Kruskall-Wallis χ2 = 19.76, df = 10, P = 0.0316), but no observed effect of sampling date on the prevalence of any of the other exploration types (Kruskall-Wallis Contact, χ2 = 7.7695, df = 10, P = 0.6513; Short distance, χ2 = 15.067, df = 10, P = 0.1296; Long distance, χ2 = 8.1626, df = 10, P = 0.613). The observed differences in medium exploration type prevalence was driven by a difference between abundance at the date 28/3/2018, and the dates 27/1/2017 and 22/8/2017 (Figure 5). There was no interactive effect of plot and date on the prevalence of exploration types (Kruskall-Wallis Contact, χ2 = 33.211, df = 32, P = 0.408; Short distance, χ2 = 37.126, df = 32, P = 0.2445; Medium distance, χ2 = 37.198, df = 32, P = 0.242; Long distance, χ2 = 29.462, df = 32, P = 0.5956).
Figure 4. Mean relative abundance (± stderr) of ECM exploration types of Q. pubescens at the three plots in the Slovenian Sub-Mediterranean. Significantly different (Kruskal-Wallis and Dunn's post-hoc test) mean relative abundances at P < 0.05 are marked by different letters.
Figure 5. Mean relative abundance (± stderr) of ECM exploration types of Q. pubescens over the 11 sample dates. Significantly different (Kruskal-Wallis and Dunn's post-hoc test) mean relative abundances at P < 0.05 are marked by different letters.
Seasonal Dynamics of ECM Root Tip Vitality, Species Richness and Exploration Types
The effect of environmental parameters on ECM root tip properties (including the percentage of vital root tips, alpha diversity parameters, and exploration types) was tested using linear modeling of these parameters against one another to unpick how seasonal shifts in the environment may contribute to the variation seen over time. The percentage of vital root tips did not show any association with the measured environmental properties. Species richness was found to be significantly associated with mean air and soil temperature, relative air humidity, and solar radiation. Community evenness and dominance were both associated with mean soil temperature. The percentage of ECM fungi belonging to the short distance exploration type was significantly associated with precipitation, medium distance exploration type was significantly associated with mean soil temperature and relative air humidity, whereas those belonging to the long exploration type was associated with precipitation. The strength and direction of association of each property is shown in Table 3.
Table 3. Results of linear regression analysis (through generalized linear modeling) between forward-selected environmental parameters and root tip ECM community parameters.
Soil-Inhabiting Extramatrical Mycelium Communities
Although the mycelium ingrowth mesh bags are supposed to allow only the growth of symbiotic fungi, FUNGuild assignments of trophic mode to OTUs showed that other trophic modes were present in addition to symbiotic ECM fungi detailed in the Supplementary Results and Supplementary Figure 6. Only fungi assigned to the ECM functional group were used for subsequent analyses.
One hundred eighty-four ECM fungal species belonging to 36 genera and 22 families were found within the mesh bags at the site. The mycelium of ECM fungi inside mesh bags was dominated by taxa belonging to the genera Tomentella, Sebacina, Pseudotomentella, which were present in relative abundances of 36.09% ± 4.41 stderr, 25.96% ± 4.44 stderr, and 13.62% ± 3.35 stderr across all samples and plots. Other abundant genera (>1% relative abundance) included Amanita, Xerocomus, Geopora, Trichophaea, Cortinarius, and Inocybe. Sixteen genera were present in relative abundances below 1%, accounting for only 2.38% of the community (Figure 6). Excluding families represented solely by one genus, abundant families (>1% relative abundance) included Sebacinaceae (represented primarily by the genus Sebacina with one OTU belonging to the genus Chaetospermum), Thelephoraceae (represented by the genera Tomentella and Pseudomentella), Boletaceae (represented by the genera Boletus, Xerocomellus, and Xerocomus) and Pyronemataceae (represented by the genera Genea, Geopora, Otidea, Tarzetta, and Trichophaea).
Figure 6. Taxa charts displaying the most abundant genera (>1%) found in soil-exploring ECM mycelium collected in mesh bags in the Slovenian Sub-Mediterranean per plot, sampling date and for the whole site (total).
High intra-site variation was observed, with the community of soil-exploring ECM within the extramatrical mycelium being found to be different between plots (PERMANOVA F = 6.8184, df = 2, P = 0.001, Figure 7). This explained 21.43% of the variation found between samples. Pairwise-PERMANOVA testing was conducted with Benjamini-Hochberg multiple comparison corrections, revealing that all three plots harbored divergent communities of ECM fungi (Figure 7). Based on the total number of OTUs found within each plot and the estimated coverage, Chao2 estimated species richness was found to be greater in plot 3 (limestone) than the other two plots (Table 4). The observed OTU richness, evenness, and dominance indices per sample were however not found to be significantly different among plots (Kruskall-Wallis OTU richness, χ2 = 2.470, df = 2, P = 0.2909; Evenness, χ2 = 0.2062, df = 2, P = 0.902; Dominance, χ2 = 0.9446, df = 2, P = 0.624, Table 4).
Figure 7. NMDS of ECM mycelium community developed in mesh bags at three Q. pubescens plots in the Slovenian Sub-Mediterranean. Hulls delimit communities by the plot and sampling date. Stress of the ordination is 0.19 for k = 2.
Table 4. No. of observed soil mycelium ECM OTUs, estimated sample coverage, estimated species richness, community dominance, and community evenness (mean ± stderr) for the two sampling dates and for the total community of ECM mycelia in mesh bags at three Q. pubescens plots in the Slovenian sub-Mediterranean.
Among the extramatrical ECM fungi, plot 3 (limestone) was characterized by greater abundance of OTUs belonging to Thelephoraceae, Tomentella pilosa, Tomentella bryophila, Pseudotomentella, Sebacina sp., and ascomycete Trichophaea woolhopeia. OTUs in plot 2 (limestone) belonged with greater abundance to Thelephoraceae, Tomentella punicea, Inocybe asterospora, and Amanita pantherina. Plot 1 (flysch) significant indicator OTUs were the most taxonomically diverse and belonged to Thelephoraceae, Tomentella coerulea and T. umbrinospora, Xerocomus subtomentosus, Inocybe luteifolia, Cortinarius sp., Sebacina sp., Lactarius quietus, and Geopora cervina. Again, plot 1 (flysch) and plot 3 (limestone) OTUs showed the highest similarity which corresponds to the results of ECM community composition. This is shown in Supplementary Figure 7, and other indicator species relative abundances and analysis statistics are detailed in Supplementary Table 5.
Across the root tip and ECM mycelial communities combined 50 genera of fungi were observed. Of these genera, 18 were shared between the two assessed compartments, while 18 were unique to the extramatrical mycelium, and 14 were unique to the root tips (Figure 8). Most mycelial fungal reads belonged to fungal genera which were shared between the two compartments, with those specific to the mycelial compartment comprising only 5.87% of all fungal reads. Based on normalized root tip counts, relative abundance of fungal genera found only within the root tip compartment was 29.93%.
Figure 8. Unique and shared ECM genera found on root tips and inside mesh bags. Abbreviations of exploration types were added in capital letters after genus name (C, contact; SD, short distance; MD, medium distance; LD, long distance; N, not known). Exploration types of genera in mesh bags were assigned by literature search (Agerer and Rambold, 2004-2020; Suz et al., 2014), where available.
Discussion
In this study we explored the temporal variability of the ECM fungal community inhabiting Q. pubescens roots and the surrounding soil across three plots with different soil conditions in the stress prone Sub-Mediterranean region of Slovenia, from June 2016 until the end of May 2018. Environmental conditions were highly variable over this period. The area received less than average precipitation in summer 2016 and 2017 and was affected by superficial wildfire in August 2016. Along with heat effects, wildfires impact carbon and nitrogen cycling and soil pH through the inputs of organic matter in the form of ash and charcoal into the soil system. Wildfires also result in the elimination of the soil organic layer and partial heat-sterilization of soil (Buscardo et al., 2010; Taudiere et al., 2017).
The ratio of vital to non-vital roots tips significantly varied over the course of the experiment. However, none of the measured environmental variables had an apparent effect on the vital ECM root tip abundance. This has also been reported in beech-ECM associations (Fagus sylvatica L.), with the ratio of vital to non-vital root tips instead displaying a significantly negative association with tree girdling, which decreases carbohydrate flow belowground (Pena et al., 2010). In this context, the decreased ratio of vital to non-vital ECM root tips in our study could hypothetically be associated with the sequence of stress events at the site which resulted in reduced carbon flow belowground. The vital to non-vital ratio remained low from after drought and fire sequence in 2016 onwards and did not exhibit an improvement despite the normal amount of precipitation in autumn 2016 and spring 2017. Spring 2017 was again followed with less than average amount of precipitation in summer 2017. The consistently low ratio of vital to non-vital root tips below one over this period denotes a low investment into new mycorrhizal associations. Only from the beginning of 2018 did the vital to old ECM root tip ratio increase beyond a value of one, indicating that the majority of mycorrhizal root tips were vital and that trees were re-investing considerably in the ECM association. It appears that the trees required a prolonged time period to begin reinvesting into the ECM symbiosis after this sequence of extreme events. The previous assessment of tree growth and physiology at this site over the same period observed trees with scorched crowns from the August 2016 wildfire re-sprouted in September of 2016, potentially re-using carbohydrates that would otherwise have supported ECM over the next year. Tree growth was also depressed over the subsequent growing season due to unfavorable growing conditions—the site received 43% less precipitation than the long-term average (Gričar et al., 2020). This growth depression could also explain the lack of investment in ECM fungal associations. Tree responses to environmental stressors are complex and are not yet fully understood. Integrating studies of tree physiology and ECM associations may help better elucidate both tree and ECM responses in future studies.
The most abundant taxa found associated with Q. pubescens root tips at our study site belonged to the Tomentella genus and other Thelephoraceae genera, Sebacina, and Cenococcum. Species belonging to Thelephoraceae, Sebacinaceae and the Inocybe genus were also previously reported as dominant components of wooded meadows in Estonia (Tedersoo et al., 2006). Encroaching or fragmented karst grasslands in Slovenia are similar to Estonian wooded meadows in terms of high biodiversity of vascular plants and may therefore be a somewhat analogous system (Kaligarič and Ivajnšič, 2014). Fungi belonging to the family Thelephoraceae and the genus Cenococcum are characterized by melanized cell walls (Fernandez and Koide, 2013). This feature protects fungi from environmental stress, such as drought and was reported favorable against desiccation of fungi for Estonian wooded meadows with shallow soil and reduced organic matter content (Tedersoo et al., 2006). Protection from drought would be also beneficial in our investigated area where summer droughts are common. Cenococcum sp. was significantly more abundant at plot 2 (limestone) compared to both other plots and had high relative abundance already at the first sampling date before summer droughts in 2016 and 2017. Sebacina is a very widespread and host-rich fungal genus, that also associates with orchids. In temperate deciduous forests, Sebacinales ectomycorrhizae seem to co-dominate with those of Russulales and Thelephorales (Oberwinkler et al., 2013). Sebacinales or Sebacina s. str. were reported from a wide spectrum of oak species (Oberwinkler et al., 2013). Across the three plots assessed in this study several Sebacina ECM fungi were detected, but only the relative abundance of Sebacina incrustans was significantly different among plots, with the highest occurrence in plot 1 (flysch). Thelephoraceae were also reported as dominant ECM in Q. douglasii in Mediterranean California along with Inocybe (Smith et al., 2007).
The overall community composition of root-tip inhabiting ECM fungi was not affected by the sampling date, which is consistent with findings of Smith et al. (2007) and Richard et al. (2011) who also investigated ECM community composition of Quercus in the Mediterranean. The sampling date did however have a significant effect on fungal richness in our case, but not in the case of Smith et al. (2007) and Richard et al. (2011). The species richness of our site was significantly negatively associated with mean air temperature and solar radiation, and positively with relative air humidity, while community evenness and dominance were both negatively associated with mean soil temperature. Soil temperature has been previously shown to affect fungal diversity across wide latitudinal gradients (Shi et al., 2014), but we observed this happening temporally within this study. Globally, both plant and animal diversity are determined to some degree by temperature and precipitation, which also extends to ECM fungi (Hawkins et al., 2003; Tedersoo et al., 2012a). Our results indicate that even in smaller local setting seasonal variations in these properties can also drive intra-site community diversity over time.
The high prevalence of ECM fungi belonging to the short distance exploration type, followed by those of the contact exploration type indicated that Q. pubescens trees in our study area relied mainly on close contact with neighboring substrates in topsoil. Short-range (contact and short distance) exploration types mainly take up N as nitrate and ammonium, and due to their biochemistry contribute to humus build-up in soil, leading to high C sequestration in soils (Rosinger et al., 2018). Fungi belonging to contact and short distance exploration types have generally broad environmental ranges and are not as C demanding as those associated with medium distance and long-distance exploration types (Rosinger et al., 2018). The abundance of ECM fungi belonging to the short distance type was found to be negatively associated with precipitation, whereas of the long-distance type was found to be positively associated with precipitation. This is unexpected, as the ability of long-distance exploring fungi to explore the surrounding soil volume and increase water acquisition under drier scenarios (Rosinger et al., 2018) would lead us to expect the opposite relationship. Despite this relation to precipitation, long-distance exploration type taxa were present at low abundances across all investigated plots at the site. As carbon resources are scarce, at least within the limestone plots, where lower water retention capacity of soil result in stomatal limitation of photosynthesis (Vodnik et al., 2019), investments into long distance exploration types may be too costly. Alternatively, long distance exploration types might be more abundant in deeper soil. Taxa belonging to the medium distance exploration type were the only group found to vary in abundance significantly across the timescale of the experiment. This was related to the mean soil temperature and relative air humidity, from which we suggest a strong seasonal dynamic in the association of taxa belonging to this exploration type with the host tree. Our results indicate that environmental parameters may induce functional shifts on a temporal scale and are supported by the finding of another European study that showed 36% variation in exploration types could be ascribed to environmental variables (soil pH, mean annual temperature, and host taxon) (Rosinger et al., 2018).
Ingrowth mesh bags are usually used to quantify biomass of extramatrical mycelium (Wallander et al., 2001). In our case, they were used to identify taxa that contribute to the connectivity in investigated plots. Assessing ECM communities through NGS resulted in the detection of non-ECM fungal taxa in the mesh bags. Although the acid-washed sand used to fill the bags is proposed to prevent ingrowth of saprotrophs, there was still a significant proportion of them detected. The presence of saprotrophs could be related to the presence of dead organic matter (e.g., extramatrical mycelium) that accumulates in mesh bags over time or due to the presence of hyphal exudates. This itself is an important finding as it demonstrates that mesh bags are not as exclusive in recruiting ECM fungi as previously thought through isotopic labeling studies (Wallander et al., 2001, 2004; Bakker et al., 2009). This calls into question the validity of the assumptions of this experimental approach. The mesh bags did however successfully select for taxa that were exploring the soil environment rather than directly associated with Q. pubescens roots. The majority of taxa found in the mycelium mesh bags belonged to medium distance or short distance exploration types, though some contact types were detected (Figure 8). The mycelial network in mesh bags was dominated by fungi shared between this compartment and the root tips with only 5.9% of fungal reads in this compartment found to be unique to the mycelial compartment. Presence of unique taxa within mycelium compartment could be explained by fungal spores that are shed from sporocarps into the topsoil (Koide et al., 2005) or alternatively the mycelium from other ECM hosts in the surrounding (O. carpinifolia and Pinus nigra). Some genera that were unique to the mycelium community are less known as ECM fungi and we could not find their exploration type in the literature (Figure 8). However, the low percentage of fungal reads belonging to unique taxa in mesh bags demonstrates the nested nature of the mycelial community, as it was made up primarily of fungi which are found in Q. pubescens root tips. This shows the relative importance of Q. pubescens as a keystone ECM host in structuring the common mycelial network at this site.
Plots within the site were distinguished by extramatrical mycelium of different representatives of the same taxonomical groups of fungi, indicating redundancy inside higher taxonomical groups. The extramatrical mycelium community composition was characterized by a high abundance of unidentified Thelephoraceae, different Tomentella taxa, Inocybe spp., and Sebacina spp. with specificities for each plot. These taxa - except for Inocybe- were also the most abundant taxa found in the root tips of our plots. This demonstrates a robustness and complementarity of assessing both root tips and the wider mycelial community to fully assess the importance of ECM taxa on a site. In another study incorporating mesh bags carried out within Danish beech forests, Tomentella taxa were found to occur with equal frequencies in both mesh bags and neighboring roots tips, whereas boletoid species were more frequent as mycelia and Cortinarius species were absent from mesh bags (Kjøller, 2006) despite belonging to medium distance fringe exploration type (Agerer, 2001). The absence of Cortinarius from mesh bags was explained by the fact that Cortinarius species may preferentially grow their mycelia into organic substrates rather than into the pure sand in the mesh bags (Kjøller, 2006). In contrast, our mesh bags contained Cortinarius sequences even with the apparent removal of organic material.
Conclusion
Our study characterized the temporal variability in vitality of ECM root tips, species richness and the fungal consortia found within the Slovenian Sub-Mediterranean climate. We assessed both the fungi actively associating with the roots of Q. pubescens and those exploring the surrounding soil matrix both taxonomically and through functional characteristics. This approach demonstrated that different lineages and functional groups of fungi may be important in this stress-prone environment to maintain the nutritional status of Q. pubescens and the connectivity of the wider mycelial network. It demonstrated the association of environmental parameters to temporally variable ECM parameters such as species richness and percentage of selected exploration types but not to variability in ratio of vital to non-vital ECM root tips that could be tentatively related to decreased flux of carbohydrates belowground due to sequential stress events. We found out that the community composition of ECM root tips is temporally stable, but highly variable among the investigated plots, indicating relative stress resilience of the fungal communities at the investigated site. We demonstrated high importance of Q. pubescens at this site for structuring the extramatrical mycelium community. Our study and others at this site demonstrate the complex interplay of environmental parameters, abiotic stressors, tree vitality and ECM fungi within the Slovenian Sub-Mediterranean, which require further investigation to assess how this region will be affected under future climate change scenarios.
Data Availability Statement
The sequences presented in this study can be found in GenBank database: https://www.ncbi.nlm.nih.gov/genbank/, accession numbers MW027850-MW028132 and https://www.ncbi.nlm.nih.gov/genbank/, accession numbers MW039835-MW040018.
Author Contributions
JG and HK: funding acquisition. JG, HK, and IŠ: experiment conceptualization. IŠ, TM, and NŠ: data acquisition. NŠ, PB-J, and TM: data analysis. TM: original draft preparation. All authors: writing review and editing.
Funding
This study was funded by the Slovenian Research Agency projects J4-7203 (Short- and long-term responses of oak in Sub-Mediterranean region to extreme weather events using a tree-anatomical analysis and ecophysiological measurements) and J4-9297 (The agreement and synchrony between woody carbon sequestration and eddy covariance estimates of net ecosystem productivity in a heterogeneous open woodland ecosystem) and Slovenian Research Agency Research Programme P4-0107 (Forest biology, ecology, and technology).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
Authors would like to thank Dr. Mitja Ferlan (Slovenian Forestry Institute) for field measurements of soil temperature and soil matrix potential and Dr. Aleksander Mahnic and Dr. Maja Rupnik (National laboratory of health, environment and food, Slovenia) for NGS sequencing service and bioinformatic analyses. Melita Hrenko and Barbara Štupar (Slovenian Forestry Institute) are acknowledged for technical assistance during lab work.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2021.599946/full#supplementary-material
Footnotes
1. ^DEEMY – An Information System for Characterization and Determination of Ectomycorrhizae (2004–2007). Available online at: www.deemy.de
References
Abarenkov, K., Zirk, A., Piirmann, T., Pöhönen, R., Ivanov, F., Nilsson, R. H., et al. (2020). UNITE QIIME release for Fungi. doi: 10.15156/BIO/786385
Agerer, R. (2001). Exploration types of ectomycorrhizae Mycorrhiza 11, 107–114. doi: 10.1007/s005720100108
Agerer, R., and Rambold, G. (2004–2020). DEEMY – An Information System for Characterization and Determination of Ectomycorrhizae. Available online at: www.deemy.de
Allen, C. D., Macalady, A. K., Chenchouni, H., Bachelet, D., McDowell, N., Vennetier, M., et al. (2010). A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. For. Ecol. Manage. 259, 660–684 doi: 10.1016/j.foreco.2009.09.001
ARSO (2020). Available online at: http://meteo.arso.gov.si/met/sl/archive/
Bakker, M. R., Jolicoeur, E., Trichet, P., Augusto, L., Plassard, C., Guinberteau, J., et al. (2009). Adaptation of fine roots to annual fertilization and irrigation in a 13-year-old Pinus pinaster stand. Tree Physiol. 29, 229–238. doi: 10.1093/treephys/tpn020
Beidler, K. V., Taylor, B. N., Strand, A. E., Cooper, E. R., Schönholz, M., and Pritchard, S. G. (2015). Changes in root architecture under elevated concentrations of CO2 and nitrogen reflect alternate soil exploration strategies. New Phytol. 205, 1153–1163. doi: 10.1111/nph.13123
Bent, E., Kiekel, P., Brenton, R., and Taylor, D. L. (2011). Root-associated ectomycorrhizal fungi shared by various boreal forest seedlings naturally regenerating after a fire in interior Alaska and correlation of different fungi with host growth responses. Appl. Environ. Microbiol. 77, 3351–3359. doi: 10.1128/aem.02575-10
Bréda, N., Huc, R., Granier, A., and Dreyer, E. (2006). Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences. Ann. For. Sci. 63, 625–644. doi: 10.1051/forest:2006042
Brus, R. (2012). Drevesne vrste na Slovenskem (Tree Species in Slovenia), 2nd Edn. Ljubljana: Mladinska knjiga & Robert Brus, 406.
Buscardo, E., Rodríguez-Echeverría, S., Martín, M. P., De Angelis, P., Pereira, J. S., and Freitas, H. (2010). Impact of wildfire return interval on the ectomycorrhizal resistant propagules communities of a Mediterranean open forest. Fungal Biol. 114, 628–636. doi: 10.1016/j.funbio.2010.05.004
Cáceres, M. D., and Legendre, P. (2009). Associations between species and groups of sites: indices and statistical inference. Ecology 90, 3566–3574. doi: 10.1890/08-1823.1
Chao, A., Ma, K. H., Hsieh, T. C., and Chiu, C. H. (2015). Online Program SpadeR (Species-richness Prediction and Diversity Estimation in R). Program and User's Guide. Available online at: http://chao.stat.nthu.edu.tw/wordpress/software_download/
Ferlan, M., Alberti, G., Eler, K., Batič, F., Peressotti, A., Miglietta, F., et al. (2011). Comparing carbon fluxes between different stages of secondary succession of a karst grassland. Agric. Ecosyst. Environ. 140, 199–207. doi: 10.1016/j.agee.2010.12.003
Ferlan, M., Eler, K., Simončič, P., Batič, F., and Vodnik, D. (2016). Carbon and water flux patterns of a drought-prone mid-succession ecosystem developed on abandoned karst grassland. Agric. Ecosyst. Environ. 220, 152–163. doi: 10.1016/j.agee.2016.01.020
Fernandez, C. W., and Koide, R. T. (2013). The function of melanin in the ectomycorrhizal fungus Cenococcum geophilum under water stress. Fungal Ecol. 6, 479–486. doi: 10.1016/j.funeco.2013.08.004
Gardes, M., and Bruns, T. D. (1993). ITS primers with enhanced specificity for basidiomycetes–application to the identification of mycorrhizae and rusts. Mol. Ecol. 2, 113–118.
Geršič, M., Repe, B., Blatnik, M., Brečko Grubar, V., Kovač, B., Pozvek, N., et al. (2014). Geografija in rastlinska sukcesija - izbrani primeri iz slovenskih pokrajin vol 23. Georitem. Ljubljana: ZaloŽba ZRC.
Grebenc, T., and Kraigher, H. (2007). Changes in the community of ectomycorrhizal fungi and increased fine root number under adult beech trees chronically fumigated with double ambient ozone concentration Plant Biol. 9, 279–287. doi: 10.1055/s-2006-924489
Gričar, J., Hafner, P., Lavrič, M., Ferlan, M., Ogrinc, N., Krajnc, B., et al. (2020). Post-fire effects on development of leaves and secondary vascular tissues in Quercus pubescens. Tree Physiol. 40, 796–809. doi: 10.1093/treephys/tpaa030
Gričar, J., Zavadlav, S., Jyske, T., Lavrič, M., Laakso, T., Hafner, P., et al. (2019). Effect of soil water availability on intra-annual xylem and phloem formation and non-structural carbohydrate pools in stem of Quercus pubescens. Tree Physiol. 39, 222–233. doi: 10.1093/treephys/tpy101
Gudmundsson, L., and Seneviratne, S. (2016). Anthropogenic climate change affects meteorological drought risk in Europe. Environ. Res. Lett. 11:044005. doi: 10.1088/1748-9326/11/4/044005
Harfouche, A., Meilan, R., and Altman, A. (2014). Molecular and physiological responses to abiotic stress in forest trees and their relevance to tree improvement. Tree Physiol. 34, 1181–1198. doi: 10.1093/treephys/tpu012
Hawkins, B. A., Field, R., Cornell, H. V., Currie, D. J., Guégan, J., Kaufman, D. M., et al. (2003). Energy, water and broad-scale geographic patterns of species richness. Ecology 84, 3105–3117. doi: 10.1890/03-8006
Ihrmark, K., Bödeker, I. T. M., Cruz-Martinez, K., Friberg, H., Kubartova, A., Schenck, J., et al. (2012). New primers to amplify the fungal ITS2 region – evaluation by 454-sequencing of artificial and natural communities. FEMS Microbiol. Ecol. 82, 666–677. doi: 10.1111/j.1574-6941.2012.01437.x
Kaligarič, M., and Ivajnšič, D. (2014). Vanishing landscape of the “classic” Karst: changed landscape identity and projections for the future. Landsc. Urban Plan. 132, 148–158. doi: 10.1016/j.landurbplan.2014.09.004
Kjøller, R. (2006). Disproportionate abundance between ectomycorrhizal root tips and their associated mycelia. FEMS Microbiol. Ecol. 58, 214–224. doi: 10.1111/j.1574-6941.2006.00166.x
Koide, R. T., Xu, B., and Sharda, J. (2005). Contrasting below-ground views of an ectomycorrhizal fungal community. New Phytol. 166, 251–262. doi: 10.1111/j.1469-8137.2004.01313.x
Lehner, B., Döll, P., Alcamo, J., Henrichs, T., and Kaspar, F. (2006). Estimating the impact of global change on flood and drought risks in europe: a continental, integrated analysis. Clim. Change 75, 273–299. doi: 10.1007/s10584-006-6338-4
Mahnic, A., Breskvar, M., Dzeroski, S., Skok, P., Pintar, S., and Rupnik, M. (2020). Distinct types of gut microbiota dysbiosis in hospitalized gastroenterological patients are disease non-related and characterized with the predominance of either enterobacteriaceae or enterococcus. Front. Microbiol. 11:120. doi: 10.3389/fmicb.2020.00120
Mrak, T., Štraus, I., Grebenc, T., Gričar, J., Hoshika, Y., Carriero, G., et al. (2019). Different belowground responses to elevated ozone and soil water deficit in three European oak species (Quercus ilex, Q. pubescens and Q. robur). Sci. Total Environ. 651, 1310–1320. doi: 10.1016/j.scitotenv.2018.09.246
Nguyen, N. H., Song, Z., Bates, S. T., Branco, S., Tedersoo, L., Menke, J., et al. (2016). FUNGuild: an open annotation tool for parsing fungal community datasets by ecological guild. Fungal Ecol. 20, 241–248. doi: 10.1016/j.funeco.2015.06.006
Oberwinkler, F., Riess, K., Bauer, R., Selosse, M.-A., Weiß, M., Garnica, S., et al. (2013). Enigmatic sebacinales. Mycol. Prog. 12, 1–27. doi: 10.1007/s11557-012-0880-4
Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2019). vegan: Community Ecology Package. R package version 2.5–6. Available online at: https://CRAN.R-project.org/package=vegan
Pena, R., Offermann, C., Simon, J., Naumann, P. S., Geßler, A., Holst, J., et al. (2010). Girdling affects ectomycorrhizal fungal (EMF) diversity and reveals functional differences in EMF community composition in a beech forest. Appl. Environ. Microbiol. 76, 1831–1841. doi: 10.1128/AEM.01703-09
Richard, F., Roy, M., Shahin, O., Sthultz, C., Duchemin, M., Joffre, R., et al. (2011). Ectomycorrhizal communities in a Mediterranean forest ecosystem dominated by Quercus ilex: seasonal dynamics and response to drought in the surface organic horizon. Ann. For. Sci. 68, 57–68. doi: 10.1007/s13595-010-0007-5
Rosinger, C., Sandén, H., Matthews, B., Mayer, M., and Godbold, D. L. (2018). Patterns in ectomycorrhizal diversity, community composition, and exploration types in European beech, pine, and spruce forests. Forests 9:445. doi: 10.3390/f9080445
Sheffield, J., and Wood, E. F. (2008). Projected changes in drought occurrence under future global warming from multi-model, multi-scenario, IPCC AR4 simulations. Clim. Dyn. 31, 79–105. doi: 10.1007/s00382-007-0340-z
Shi, L., Mortimer, P. E., Ferry Slik, J. W., Zou, X.-M., Xu, J., Feng, W.-T., et al. (2014). Variation in forest soil fungal diversity along a latitudinal gradient. Fungal Divers. 64, 305–315. doi: 10.1007/s13225-013-0270-5
Simard, S. W., Beiler, K. J., Bingham, M. A., Deslippe, J. R., Philip, L. J., and Teste, F. P. (2012). Mycorrhizal networks: mechanisms, ecology and modelling. Fungal Biol. Rev. 26, 39–60. doi: 10.1016/j.fbr.2012.01.001
Smith, M. E., Douhan, G. W., and Rizzo, D. M. (2007). Ectomycorrhizal community structure in a xeric Quercus woodland based on rDNA sequence analysis of sporocarps and pooled roots. New Phytol. 174, 847–863. doi: 10.1111/j.1469-8137.2007.02040.x
Smith, S. E., and Read, D. (2008). “Chapter 8 - Growth and carbon allocation of ectomycorrhizal symbionts,” in Mycorrhizal Symbiosis (Third Edition), eds S. Smith, D. Read (London: Academic Press), 295–320. doi: 10.1016/B978-012370526-6.50010-6
Šturm, T., and Podobnikar, T (2017). A probability model for long-term forest fire occurrence in the Karst forest management area of Slovenia. Int. J. Wildland Fire 26, 399–412. doi: 10.1071/WF15192
Sutanto, S. J., Vitolo, C., Di Napoli, C., D'Andrea, M., and Van Lanen, H. A. J. (2020). Heatwaves, droughts, and fires: exploring compound and cascading dry hazards at the pan-European scale. Environ. Int. 134:105276. doi: 10.1016/j.envint.2019.105276
Suz, L. M., Barsoum, N., Benham, S., Dietrich, H. P., Fetzer, K. D., Fischer, R., et al. (2014). Environmental drivers of ectomycorrhizal communities in Europe's temperate oak forests. Mol. Ecol. 23, 5628–5644. doi: 10.1111/mec.12947
Taudiere, A., Richard, F., and Carcaillet, C. (2017). Review on fire effects on ectomycorrhizal symbiosis, an unachieved work for a scalding topic. For. Ecol. Manage. 391, 446–457. doi: 10.1016/j.foreco.2017.02.043
Tedersoo, L., Bahram, M., Toots, M., Diédhiou, A. G., Henkel, T. W., Kjøller, R., et al. (2012a). Towards global patterns in the diversity and community structure of ectomycorrhizal fungi. Mol. Ecol. 21, 4160–4170. doi: 10.1111/j.1365-294X.2012.05602.x
Tedersoo, L., Naadel, T., Bahram, M., Pritsch, K., Buegger, F., Leal, M., et al. (2012b). Enzymatic activities and stable isotope patterns of ectomycorrhizal fungi in relation to phylogeny and exploration types in an afrotropical rain forest. New Phytol, 195, 832–843. doi: 10.1111/j.1469-8137.2012.04217.x
Tedersoo, L., Suvi, T., Larsson, E., and Kõljalg, U. (2006). Diversity and community structure of ectomycorrhizal fungi in a wooded meadow. Mycol. Res. 110, 734–748. doi: 10.1016/j.mycres.2006.04.007
Veble, D., and Brečko Grubar, V. (2016). Pogostost in obseg požarov v naravi na Krasu in v slovenski Istri / The frequency and extent of wildfires on the Kras and in Slovenian Istria. Geogr. Vestnik 88, 9–20. doi: 10.3986/GV88101
Vodnik, D., Gričar, J., Lavrič, M., Ferlan, M., Hafner, P., and Eler, K. (2019). Anatomical and physiological adjustments of pubescent oak (Quercus pubescens Willd.) from two adjacent sub-Mediterranean ecosites. Environ. Exp. Bot. 165, 208–218. doi: 10.1016/j.envexpbot.2019.06.010
Wallander, H., Göransson, H., and Rosengren, U. (2004). Production, standing biomass and natural abundance of 15N and 13C in ectomycorrhizal mycelia collected at different soil depths in two forest types. Oecologia 139, 89–97. doi: 10.1007/s00442-003-1477-z
Wallander, H., Nilsson, L. O., Hagerberg, D., and Bååth, E. (2001). Estimation of the biomass and seasonal growth of external mycelium of ectomycorrhizal fungi in the field. New Phytol. 151, 753–760. doi: 10.1046/j.0028-646x.2001.00199.x
White, T., Bruns, T., Lee, S., Taylor, J., Innis, M., Gelfand, D., et al. (1990). Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. 31, 315–322.
Keywords: pubescent oak, abiotic stress, drought, wildfire, extramatrical mycelium, mesh bags, exploration types, ECM fungi
Citation: Mrak T, Šibanc N, Brailey-Jones P, Štraus I, Gričar J and Kraigher H (2021) Extramatrical Mycelium and Ectomycorrhizal Community Composition of Quercus pubescens in a Sub-Mediterranean Stress-Prone Environment. Front. For. Glob. Change 4:599946. doi: 10.3389/ffgc.2021.599946
Received: 28 August 2020; Accepted: 11 January 2021;
Published: 11 February 2021.
Edited by:
Weile Chen, Zhejiang University, ChinaReviewed by:
Roger Tai Koide, Brigham Young University, United StatesJian Huang, Northwest A and F University, China
Copyright © 2021 Mrak, Šibanc, Brailey-Jones, Štraus, Gričar and Kraigher. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Tanja Mrak, tanja.mrak@gozdis.si