2.2. Design of the Field Survey
To assess the potential impact of urbanization on the soil fungal community in urban forests, we selected 20 deciduous forests belonging to the
Fagetum association [
37] along an urbanization gradient (
Figure 1). The forests studied ranged in size from 0.23 ha to 337.0 ha and differed in their historical development (
Table 1). Fifteen of them are surrounded by settlements and agricultural land and most of them are no longer connected to large continuous forests. Six of these forests are remnants (fragments) of former large continuous forests and six were planted after 1884, while three forests are part of large continuous forests (> 40 ha;
Table 1). The remaining five forests are situated in the rural surroundings and are part of large, continuous forests (> 76.2 ha;
Table 1). The most abundant tree species in these forests are European beech (
Fagus sylvatica), sycamore (
Acer pseudoplatanus) and European oak (
Quercus robur). The ground vegetation in the forests shows a high richness of vernal geophytes, including
Anemone nemorosa,
Ranunculus ficaria,
Polygonatum multiflorum and
Arum maculatum [
37]. Most forests were state owned and accessible to the public. Some forests are privately owned but managed by the forestry authorities.
In each forest, we chose an area dominated by European beech (80–90% of all tree individuals) and set up a study plot measuring 10 m x 10 m in its centre. The study plots had a minimum distance of 5 m to the nearest forest edge or to permanent paths to minimize possible edge effects. The forest management in the study plots (time since last thinning and management intensity) was similar among the forests investigated.
As a measure of degree of urbanization, we determined the percentage of sealed areas within a radius of 500 m from the centre of each study plot using satellite images from Google Earth [
38] and the pixel count function of Adobe Photoshop (version 10.0.1). Using the same method, we assessed the percentage of forest cover within a radius of 500 m from the centre of each study plot. The percentage of sealed areas around the study plots ranged from 1 to 69% in the forest examined, and the percentage of forest cover was between 1 and 92% (
Table 1).
2.4. Soil Sampling and Physiochemical Properties
In spring 2021, we collected five soil samples (one in each corner and one in the centre) in each sampling plot. We removed the litter layer and used a metal cylinder (diameter: 5.05 cm; soil volume: 100 cm3) to sample the soil to a depth of 5 cm. The five soil samples from a sampling plot were pooled and transported on ice to the laboratory where they were sieved (mesh size: 2 mm). This resulted in a total of 80 soil samples (4 sampling plots x 20 forests). One part of each soil sample was stored at –80 oC until DNA extraction, while the remaining part was stored at –20 oC for determination of physiochemical soil properties.
We determined soil moisture content (%) using the fresh to dry weight ratio and assessed soil pH in distilled water (1:2.5 soil:water) [
40]. We assessed the total soil organic matter content (SOM, %) as loss-on-ignition of oven dried soil at 750 °C for 16 h [
40]. Total soil nitrogen content (%) was determined using the standard method of Kjeldahl [
41]. We also determined the plant-available phosphorus content of soil (μg PO
43–/g) using the molybdenum blue method according to Sparks et al. [
42].
2.5. Soil Fungal Community
We extracted total soil genomic DNA from 0.5–0.6 g soil in triplicate using NucleoSpin Soil kit (Macherey-Nagel, Oensingen, Switzerland) according to the manufacturer’s instructions. The triplicate DNA extracts were combined into one sample and purified using NucleoSpin gDNA Clean-up kit (Macherey-Nagel, Oensingen, Switzerland). We quantified the concentration and purity of DNA using NanoDrop (NanoDrop Technologies Inc., Wilmington, USA), adjusted the sample to 5 ng/µl and stored it at –20 oC.
Fungal communities were analyzed by amplification and sequencing of the internal transcribed spacer 2 region (ITS2) using the primer pair ITS3/ITS4 [
43] added with the Illumina adapters. PCR reactions (25 μl) consisted of 5 μl template DNA, 12.5 μl Master Mix (HotStar Taq Master Mix kit; Qiagen, Switzerland), 2.5 μl Primer ITS3 (10 μM), 2.5 μl Primer ITS4 (10 μM) and 2.5 µl sterile water. The amplification was achieved in an Eppendorf Mastercycler Pro (Vaudaux-Eppendorf AG, Schönenbuch, Switzerland) under the following conditions: initial 15 min heat activation step at 94 °C, followed by 35 amplification cycles of denaturation at 94 °C for 40 s, annealing at 55 °C for 40 s and extension at 72 °C for 60 s, with a final extension step at 72 °C for 10 min. PCR reactions were done in triplicate.
PCR products were sent to Microsynth AG (Belgach, Switzerland) for sequencing and analysis. Briefly, PCR products were purified, quantified, and pooled at equimolar concentrations for sequencing on an Illumina Miseq using the 2 x 250 bp paired-end approach. Paired fungal sequences from Row Illumina Miseq were demultiplexed and merged using the USEARCH pipeline [
44]. The merged sequences were then quality-filtered and clustered into operational taxonomical units (OTUs) using the USEARCH pipeline and UPARSE [
45]. Singletons were removed prior to OTU determination at 97% sequence identity and chimeric representative sequences were removed with UCHIME [
46]. Finally, the original sequences were mapped to OTUs at a 97% identity threshold to obtain an OTU table. The taxonomy of each sequence was analyzed by Ribosomal Database Project Classifier [
47] against UNITE fungal database [
48] and the NCBI/GenBank [
49].
2.6. Data Analyses
All statistical analyses were performed in R [
50]. To avoid spatial pseudo-replication, we analyzed the data at the forest site level (n = 20) by using the combined data collected in the four sampling plots in each forest.
We used generalized linear mixed models (GLM) with Poisson-distributed errors to analyze the effects of degree of urbanization and forest size (log-transformed) on the number of plant, shrub and tree species. We applied analyses of covariance (ANCOVA) to assess the effects of degree of urbanization and forest size (log-transformed) as co-factor on ground vegetation cover (arcsine square root-transformed) and on soil moisture (%; arcsine square root-transformed), soil pH, SOM (%; Tukey-transformed), total soil nitrogen (%; Tukey-transformed) and plant-available phosphorus (square root-transformed).
To avoid bias due to differences in sequencing numbers between the soil samples, we rarefied the number of sequences of each sample to the lowest value for normalization using the procedure
rarefy in the
vegan package [
51]. All further analyses were conducted with rarefied OTU data. Alpha-diversity, Shannon-diversity and Pielou’s evenness were calculated using the
vegan package [
51]. Preliminary analyses revealed inter-correlations between several chemical and physical soil properties (soil pH vs. SOM:
rs = 0.72,
n = 20,
p < 0.001; soil pH vs. total soil nitrogen content:
rs = 0.51,
n = 20,
p = 0.022, and soil pH vs. plant-available phosphorus:
rs = 0.66,
n = 20,
p = 0.002). We therefore considered soil pH, soil moisture and the recorded forest vegetation characteristics as co-factors in the subsequent statistical analyses.
We applied generalized linear mixed models (GLM) with Poisson-distributed errors to analyze the effects of degree of urbanization, forest size, percentage of forest cover within a radius of 500 m from the centre of each study plot, vegetation characteristics, and soil properties on fungal OTU richness. Degree of urbanization was included as a fixed factor, forest size (log-transformed), percentage of forest cover within a radius of 500 m from the centre of each study plot (log-transformed), ground vegetation cover (Tukey-transformed) and plant, shrub and tree species richness (all log-transformed), as well as soil moisture and soil pH as cofactors in the GLM models. We used the same GLM models, but with gamma-distributed errors, to assess the effects of degree of urbanization, forest size, percentage of forest cover within a radius of 500 m from the centre of each study plot, vegetation characteristics, and soil properties on both the Shannon diversity and evenness of fungal OTUs.
To visualize differences in OTU composition, we plotted non-metric multidimensional scaling coordinates (NMDS) for first two dimensions based on Bray-Curtis dissimilarities matrices using metaMDS in the
vegan package [
51]. Permutational multivariate analysis of variance (PERMANOVA) was used to test whether the degree of urbanization affects the composition of fungal OTUs [
52]. The same co-factors as in the GLM models were included in the model. All PERMANOVA tests were based on 9999 permutations of the untransformed raw data, using the
adonis function in the
vegan package [
51]. Finally, individual OTU affinity with a given degree of urbanization was determined using indicator species analysis by the
multipatt function in the
indicspecies R package [
53], which tests the significance of the indicator species index through a permutations test with 9999 permutations.
We applied the same GLM models with gamma-distributed errors to assess the impact of degree of urbanization on the relative abundance of Basiodomycota, Chytridiomycota, Ascomycota and Morteriellomycota (all Tukey-transformed). PERMANOVA analyses, as described above, were conducted to test whether the degree of urbanization affects the composition and the relative abundance of fungal phyla. In addition, we used analysis of similarity (ANOSIM) in the
vegan package [
51] to examine differences in the composition of fungal phyla between the different degrees of urbanization. ANOSIM is a nonparametric permutation procedure that allows comparison of between-group and within-group dissimilarities [
54]. The procedure calculates
R statistics ranging from −1 to 1.
R = 0 indicates completely random grouping, while
R = 1 when all replicates within groups are more similar than all replicates between groups.
We also used GLM models with gamma-distributed errors as described above to analyze the impact of degree of urbanization on the relative abundance of saphotrophic (sqrt-transformed), symbiotrophic (Tukey-transformed) and pathotrophic (sqrt-transformed) fungi.