454 Pyrosequencing Analysis of Fungal Assemblages from Geographically Distant , Disparate Soils Reveals Spatial Patterning and a Core Mycobiome

Identifying a soil core microbiome is crucial to appreciate the established microbial consortium, which is not usually subjected to change and, hence, possibly resistant/resilient to disturbances and a varying soil context. Fungi are a major part of soil biodiversity, yet the mechanisms driving their large-scale ecological ranges and distribution are poorly understood. The degree of fungal community overlap among 16 soil samples from distinct ecosystems and distant geographic localities (truffle grounds, a Mediterranean agro-silvo-pastoral system, serpentine substrates and a contaminated industrial area) was assessed by examining the distribution of fungal ITS1 and ITS2 sequences in a dataset of 454 libraries. ITS1 and ITS2 sequences were assigned to 1,660 and 1,393 Operational OPEN ACCESS Diversity 2013, 5 74 Taxonomic Units (OTUs; as defined by 97% sequence similarity), respectively. Fungal beta-diversity was found to be spatially autocorrelated. At the level of individual OTUs, eight ITS1 and seven ITS2 OTUs were found in all soil sample groups. These ubiquitous taxa comprised generalist fungi with oligotrophic and chitinolytic abilities, suggesting that a stable core of fungi across the complex soil fungal assemblages is either endowed with the capacity of sustained development in the nutrient-poor soil conditions or with the ability to exploit organic resources (such as chitin) universally distributed in soils.


Introduction
Uncovering a core microbiome (the suite of members shared among microbial consortia from similar habitats) is critical to understanding the stable, consistent components across complex microbial assemblages [1].
The hunt for the soil core microbiome is connected with explaining the patterns of beta-diversity (differences in community composition among sites) and, hence, distribution of microorganisms across geographically separated sites.Soils represent a huge reservoir of biodiversity (one gram of soil being estimated to contain more than 10 10 prokaryotes [2] and approximately 1,000 Gbp of microbial genome sequences [3]).Therefore, only high-throughput and high-resolution detection methods, such as those offered by next-generation sequencing (NGS) techniques, can adequately assist in the task of extensively and intensively investigating patterns of distribution of microbial communities in this environment.A number of studies taking advantage of such methods have revealed distinct biogeographical patterns in soil bacterial diversity [4][5][6][7][8].Most of the current knowledge on microbial diversity patterns in soil, therefore, concerns bacteria, whereas the fungal communities (mycobiomes) have received little attention.This is a serious omission, given that fungi comprise a major portion of the biodiversity and biomass in soils and play crucial roles in maintaining soil processes, which ultimately affect the functioning of terrestrial ecosystems [9][10][11][12][13].
While the early debate on the forces shaping microbial biogeography has been dominated by the comparison between macro-and micro-organisms, the focus is now shifting to identifying the mechanisms determining the biogeographical patterns observed [7].Lifestyle and dispersal properties can be critical determinants of spatial patterns [14,15].Given their unique biological and ecological features, fungi may therefore deviate from the relationships documented for prokaryotes.For instance, while bacterial diversity and biomass are mainly correlated with soil pH at different spatial scales [16][17][18][19][20], the latter environmental factor exerts a weaker influence on fungal community composition [21].Although the hyphae of some fungal species are subjected to a high turnover [22], hyphae are generally more persistent in soil than bacterial cells [23,24].The fungal:bacterial ratio of soil communities may change in response to physical disturbance, soil nutrient availability and moisture [25,26].Such key differences between fungi and bacteria have been shown to lead to distinct diversity and distribution patterns.For instance, differences in body size between bacteria and fungi translate into different local abundance and diversity patterns along physico-chemical gradients in soil [27].In their comparison of bacterial and fungal community composition in rhizospheric and bulk soil, Hovatter and colleagues [28] found that, in contrast to the mechanisms observed for bacteria, almost none of the variability in fungal communities was explained by geographic coordinates or soil characteristics.Indeed, patterns in fungal diversity profiles are often difficult to discern, due to high variability overall and low consistency between replicate samples, e.g., [29].
However, none of the studies focusing on spatial patterns in fungal diversity carried out so far have taken advantage of the potential of NGS techniques, relying instead on community diversity profiling tools (e.g., DGGE, ARISA, t-RFLP).Although NGS technologies are not exempt from a number of potential sources of bias [30], they allow to investigate microbial diversity at an unprecedented level of resolution.By contrast, while allowing to gain information on how the dominant members of microbial communities differ in composition across landscapes, diversity profiling techniques suffer from inherently low levels of taxonomic discrimination [20].Taxonomic resolution has a critical influence on the interpretation of both biogeographic patterns and the processes driving such patterns.Indeed, patterns observed at a given level of taxonomic discrimination might not be found at lower or higher levels [5,7,14].Studies on the diversity of soil fungal communities carried out so far by means of NGS have either focused on a phylogenetically restricted set of taxa (e.g., arbuscular mycorrhizal fungi [31][32][33][34]) or have been performed on a restricted (local or regional) spatial scale (e.g., [21,[35][36][37][38][39][40][41][42][43][44]).To bridge the gap between local and wide-scale assessments, additional studies are required, encompassing greater spatial sampling across multiple soil biomes.
To this aim, in this study, we compare the diversity of soil fungal communities (assessed by 454 pyrosequencing of both ITS1 and ITS2 sequences) at distinct localities in Italy and France, separated by up to 760 km.At both the local and the broader (trans-regional) scale, the studied samples represented a diverse array of soil and environmental characteristics.Our specific goals were: (1) to assess patterns of variation of soil fungal communities on a wide spatial scale and (2) to address the occurrence of a set of fungal OTUs (-core‖ mycobiome) shared by the disparate soils examined at the latter scale.

Comparison of the Mycobiomes in the Different Soil Samples
We examined the distribution of fungal ITS1 and ITS2 sequences in a dataset of 454 libraries assembled from a variety of soil types and distinct geographic localities ( [39,41,43] and this study; Table 1).Consistently with the findings by other authors [39,43], overall, a higher number of ITS1 than ITS2 sequences were obtained (32,673 and 21,667 reads, respectively).Following quality trimming, denoising and chimera removal, 16,758 ITS1 and 11,234 ITS2 sequences were subjected to downstream analyses (Supplementary Material 1).To account for the different sequencing depth per sample, ITS1 and ITS2 datasets were rarefied to even sequencing depth (254 and 158 ITS1 and ITS2 sequences per sample, respectively).The resulting sequences were assigned to 1,660 ITS1 and 1,393 ITS2 fungal Operational Taxonomic Units (OTUs), defined using a 97% sequence identity threshold (depending on the soil sample, 63-385 and 59-338 ITS1 and ITS2 OTUs, respectively; Supplementary Material 1).Rarefaction analysis for all samples and datasets indicated undersampling of actual richness (Figure 1).When considering the distribution of OTUs across all 16 sampled plots (Supplementary Material 2a), more than 73% and 80% of all ITS1 and ITS2 OTUs, respectively, were not detected in more than one soil sample and only approx.1% of all OTUs was observed in ≥50% of assemblages for both ITS1 and ITS2.Similarly, all sites revealed a distribution pattern featuring a few widely dispersed OTUs and many more confined OTUs (at each geographic site, 68.6-89.9%ITS1 and 81.6-91.4% ITS2 OTUs were retrieved in just one sampled plot, and 0.8-20.2%ITS1 and 1.1-18.1% ITS2 OTUs were retrieved in all plots) (Supplementary Material 2b).
Principal Components Analysis (PCA) was carried out to compare OTU distribution among the different soil samples (only OTUs occurring in either all five soil sample groups-WTG, BTG, SAR, SER and CIS-or all samples within a group, were included in the analyses).In the case of ITS1 OTUs (Figure 2), the first ordination axis (58.1% of total variance) separated fungal communities in the black truffle grounds from those in the other sample groups (mainly due to 65 OTUs occurring exclusively in the black truffle soil samples; Supplementary Material 3), whereas the second axis (13.6% of the total variation) distinguished the white truffle grounds (that featured seven exclusive OTUs) from the other sample groups.Similarly for ITS2, the black truffle ground mycobiomes segregated along the first axis (48.6% tot.var.; 44 exclusive OTUs), while the second axis (20.8% tot.var.) distinguished the white truffle ground samples (14 exclusive OTUs).In particular, both the ITS1 and ITS2 OTUs segregating the black truffle ground assemblages from the other mycobiomes included a high proportion of the OTUs assigned to ectomycorrhizal (ECM) taxa (43.3 and 54.5% of total ECM ITS1 and ITS2 OTUs included in the analysis, respectively; the ECM OTUs exclusively found in the BTG samples being ITS1 OTUs 108, 113, 132, 220, 285, 323, 415, 482, 501, 528, 542, 700, 884 (assigned to Inocybe, Russula, Scleroderma, Tomentella, Cenococcum, Pheangium, Trichophaea and Tuber spp.) and ITS2 OTUs 40, 239, 254, 328, 378, 631 (assigned to Hymenogaster, Sarcodon, Scleroderma, Geopora and Trichophaea spp.)).However, when the overall degree of similarity between pairs of soils was examined (Jaccard indices), each mycobiome generally exhibited higher similarity to the other mycobiomes from the same sample group than to the mycobiomes from the other groups.This result was constant over different sequencing depths.Figure 3 illustrates data obtained with the rarified datasets of 150 sequences (150 sequences per sample for both ITS1 and ITS2).The same pattern was observed when singleton OTUs were excluded from the analysis (data not shown).
Correlations between community and geographic distances were therefore tested by significance by means of Mantel tests.For both, the ITS1 and ITS2, Mantel correlation between community and geographic distances was significant (correlation coefficients being 0.43 and 0.28, with associated probabilities of 0.0002 and 0.008 for ITS1 and ITS2, respectively).Similarly, highly significant correlations were found for log-transformed distances (p = 0.0001 for both ITS1 and ITS2).Significant correlations were also obtained, with dissimilarity matrices derived from the other rarefied datasets, as well by excluding singleton OTUs from the analyses (data not shown).
In spite of their overall dissimilarity, the mycobiomes from different sample groups shared some features.The phylum-level taxonomic assignment of ITS1 and ITS2 OTUs in each fungal community is reported in Figure 4.In most soil samples, Ascomycota were represented by the largest number of both ITS1 and ITS2 OTUs (accounting for 36.7-92.9% of total ITS1 and 38.5-92.9% of total ITS2 OTUs).However, the mycobiomes of the soil samples with aboveground vegetation dominated by ectomycorrhizal plants (BTG1-2, WTG1-3 and SAR4-5) featured considerable proportions of OTUs assigned to Basidiomycota (accounting for 21.6-56.7%and 12.1-51.3% of total ITS1 and total ITS2 OTUs, respectively).Overall, Basidiomycota OTU numbers were significantly higher in the former soil samples than in the other soil samples (Mann Whitney U test, p = 0.001 for both ITS1 and ITS2).
Comparison between the ITS1 and ITS2 datasets was performed following rarefaction to the same sequencing depth for both datasets (150 sequences per sample).Consistently, with its higher fungal specificity, the ITS1F/ITS2 primer pair yielded a lower number of non-fungal sequences than the ITS3/ITS4 primer pair (six and 269 sequences, respectively).The ITS1F/ITS2 primer pair is reported to be biased towards amplification of basidiomycetes [47].However, such a bias was not consistently observed for the samples we amplified, as, for instance, no significant difference between the numbers of basidiomycetous sequences obtained with the ITS1F/ITS2 and the ITS3/ITS4 primers was observed for samples WTG1,3, BTG1, SAR2,5, SER2 and CIS2 (Supplementary Material 4).

Fungal OTUs Shared by the Different Sample Groups
Given that higher beta-diversity was found between rather than within the five sample groups, we examined patterns of OTU occurrence across the five groups.
Eight ITS1 and seven ITS2 OTUs were found in all groups (Tables 2 and 3, Figure 5).Among such core OTUs, ITS1 OTU 274 and ITS2 OTU 185 were retrieved from all but one analyzed plots.The other core OTUs were found in all sample groups, but not in each of their constitutive plots.The less widely distributed shared OTUs were ITS1 OTU 426 and 472, as well as ITS2 OTU 394, occurring in an average of 56% and 49% of ITS1 and ITS2 assemblages, respectively (Figure 5).The best BLAST hits for each of these core OTUs (Supplementary Material 5) were sequences identified in a diverse array of soils, continents and regions.
Although read numbers were not considered in this study, such ITS core sequences were present multiple times in each dataset, from 50 to 712 read numbers for the ITS1 core OTUs and from 38 to 475 for the ITS2 ones (Tables 2 and 3).

Soil Mycobiomes Are Spatially Structured Over a Broad Scale
In this study, comparison of the mycobiomes of disparate soils across Italy and France (separated by 0.1 to 760 km) revealed non-random distributions at the examined scales, demonstrating the occurrence of spatial patterning.
Spatial variation in species assemblages can be explained by either contemporary environmental conditions (present-day attributes of the environment) or historical contingencies (past events related to origin, dispersal and extinctions of species) [5,8].According to the first scenario, high beta-diversity (i.e., large differences in community composition among sites) is caused by contemporary interactions among organisms, as well as with their physical and biotic environments [8], and community composition is regulated by niche partitioning and environmental conditions in a habitat patch (so-called -species sorting‖, [8]).By contrast, the second hypothesis postulates a dominant role for -dispersal limitation‖ of taxa within habitats, and community composition is assumed to be regulated by spatial proximity to other populations [8].The relative influence of historical versus environmental factors seems to be related to the scale of sampling [5,48].
In our study, one of the main environmental factors separating most of the examined soil samples was the nature of the plant cover.Plants provide both the resources for obligate root-associated pathogenic and mutualistic symbionts and the organic carbon required for the functioning of the decomposer subsystem, which is responsive to the nature of organic matter entering the soil [49].There are, therefore, direct correlative links between plant and soil fungal community composition (e.g., [50][51][52].In our study, several sampling plots (BTG, WTG and SAR4-5 plots) were dominated by ectomycorrhizal plants, which recruit specific guilds of symbiotic fungal species.Thus, species sorting is suggested by the clustering and segregation from the other mycobiomes of the fungal consortia in both the black and white truffle grounds (BTG and WTG).Higher numbers of OTUs assigned to Basidiomycota (including several ectomycorrhizal fungi) were found in the truffle ground mycobiomes, as well as in the mycobiomes of the other two samples collected under an ECM-type vegetation (SAR5 and SAR4, the Sardinian soil samples collected in an oak forest and wooded pasture, respectively), than in the remaining fungal consortia.However, the fungal assemblages of the SAR5 and SAR4 samples were more similar to the mycobiomes of the other Sardinian soil samples, despite strong differences not only in plant cover, but also in management practices, as in the case of the tilled vineyard, impacted by pesticides and fertilization.This observation indicates that beta-diversity was, at least to some degree in this case, decoupled from environmental variability.Indeed, fungal beta-diversity and spatial distance were found to be significantly correlated.Similarly, in a regional-scale study of Australian desert ascomycete assemblages, Green and colleagues [53] found that geographic distance was a better predictor for community turnover than was habitat (as defined by soil and vegetation type).Arbuscular mycorrhizal fungal community composition in field soils was also shown to be the product of both dispersal and environmental variables [54].In different studies, a range of physico-chemical parameters was found to structure the mycobiomes of soils sampled along a 2,370 km latitudinal gradient in the southern maritime Antarctic [55], in alpine soils [38], as well as field soils in Texas [56] and Japan [57].By contrast, Hovatter et al. [28] did not find evidence for either species sorting or dispersal limitation to be acting on a regional scale, on fungal community composition in the rhizosphere of Lobelia siphilitica plants or adjacent bulk soil.
As previously mentioned, while differing at the level of OTU composition, the mycobiomes of the three localities with vegetation dominated by ECM plants (BTG, WTG and SAR4-5) shared higher numbers of OTUs assigned to Basidiomycota relative to the other fungal consortia.A likely explanation to this observation is that most if not all ectomycorrhizal species guilds are indeed dominated by Basidiomycota taxa.This observation indicates, as already reported by other authors (e.g., [14]), that differentiation among microbial communities can decrease when taxonomic resolution decreases.The level of microbial taxonomic resolution also influences the interpretation of the processes driving biogeographic patterns [5,7,14].In particular, the dominant effect of edaphic control may be manifest only when assessing sequences at a coarse taxonomic level, and caution is therefore needed in interpreting the strengths of environment-space relationships when communities are assessed at low levels of taxonomic resolution [20].

Core Soil Fungi Likely Exhibit Nutritional Versatility
Biogeographic signature patterns may be life history-dependent, depending on traits related to dispersal and colonization mode/efficiency and generally lifestyle [14].At the level of individual OTUs, we found both fungi with a restricted distribution in the studied localities and seemingly soil generalists (the -core‖ OTUs).Some of the latter were ubiquitous in all samples within a group, as well as across groups, while others exhibited different frequencies among the different sample groups, suggesting either patchy local distributions or insufficient sampling of taxa that are present in low abundance.Microbial endemism is difficult to assess, because many rare taxa will not be detected in a sample, even when they are present in the location, falsely suggesting a more restricted distribution than is actually the case [7].Widespread distribution is, on the contrary, an obvious feature, although difficult to fully quantify.Indeed, unequal sampling and undersampling, as well as differences in DNA extraction efficiency and sequence yield for the different soil samples could underestimate the actual number of OTUs shared by the soil mycobiomes.For instance, Qin et al. [58] showed that a 3X sequencing depth revealed a 25% larger core than did 1X coverage.The distribution of the -core‖ fungi we identified most likely extends beyond the sites investigated.The best BLAST hits for the core OTUs were indeed sequences from numerous and a diverse array of soils and regions (Supplementary Material 5).
Such -core‖ OTUs encompassed very few unidentified fungi, such as ITS1 OTU 426, which, based on BLAST results, could only be ascribed to Zygomycota.This is at odds with the general belief that fungal diversity is largely unknown.Furthermore, in some fungal groups, the high number of unidentifiable OTUs from environmental sequencing projects likely represent species that have been described, but for which there are no reference sequences in GenBank [59,60].However, based on the estimate of 1.5 million extant fungal species by Hawksworth [61], the odds that an unidentified OTU represents one of the currently described species for which there are no data in GenBank are about eighteen-to-one or higher, based on the diversity estimates of O'Brien et al. [62,63].Since soil is one of the most bio-diverse environments on earth, it is likely that fungi, such as ITS1 OTU 426, represent truly undescribed, possibly high-ranking taxa, such as other fungal clades detected from environmental DNA sequences, e.g., [64,65].
Other core fungal OTUs could only be identified at the genus level.For most of them, the reason was the unsuitability of the ITS for identifying species in some fungal genera [66] or ambiguous BLAST matches within a given genus.Species in these genera (Fusarium, Alternaria, Cladosporium, Mortierella, Tetracladium; Tables 2 and 3) are amongst the most commonly isolated soil fungi [67].Although the possibility that the ITS1 and ITS2 shared OTUs found in our work could represent uncultured species within the latter genera cannot be ruled out, it is likely that they encompass species (such as Fusarium oxysporum or F. solani) found amongst the commonest soil fungi, based on culture-dependent methods.
Strikingly, the remaining core OTUs, identified at the species level (Mortierella alpina, Bionectria ochroleuca, Paraphoma chrysanthemicola), are easily culturable fungi [67].Likewise, several culturable bacteria were found among the dominant 16S sequences in a large survey (approx.139,000 sequences) in four randomly chosen soils from distinctly different sites in South and North America, separated by up to 9,000 km [68].The ease of cultivation on standard nutritional media has been related to a copiotrophic strategy (as opposed to oligotrophy), since traditional culturing methods are likely to select for microorganisms that can grow rapidly in high resource environments [69].However, microbial strains may switch from copiotrophic to oligotrophic strategies, depending on the culture conditions and stage in the lifecycle [70,71], and intermediate/secondary strategies may well exist along the copiotroph-oligotroph continuum.Most of the core fungal species we identified seem indeed to combine the capacity to grow well in simple nutrient-rich media, with either oligotrophic abilities or the capacity of utilizing refractory compounds, such as chitin, an abundant N source in soils.Relatively high growth rates on silica gel containing no added carbon or even in distilled water on carefully cleaned glassware (-oligocarbonotrophy‖, [72]) have, for instance, been reported for Gliocladium roseum (=Bionectria ochroleuca), as well as several Fusarium, Penicillium, Cladosporium and Phoma (Paraphoma = Phoma) species [73][74][75][76].Very good chitin decomposition capacity (often with relatively little variation in this capacity between isolates) and/or mycoparasitic behavior are well known for Mortierella alpina and Bionectria ochroleuca [67].Multiple substrate utilization, together with the ability to scavenge low levels of nutrients, would confer an obvious ecological advantage in soil.Bulk (non-rhizospheric), sub-litter soil is indeed regarded as an environment poor in available organic carbon, since much organic material entering it has already been partially exploited and, as a consequence, readily available substrates have been removed.Organic materials in this environment are therefore composed of refractory compounds, persisting for long periods [77].
It is generally assumed that, due to such predominant low nutrient conditions, microbes mostly live in the soil as dormant propagules and that the soil microbiome is mainly composed of survivors, their numbers being maintained via periods of transient activity (when fresh organic matter enters the ecosystem), concluded by the restoration of the dormant condition [77].However, both oligocarbonotrophy and the capacity of utilizing refractory compounds, such as chitin, would permit continuous growth in soil, allowing the microorganisms endowed with such abilities to utilize new nutrient sources more effectively than microorganisms relying largely on spores for survival.The former biota would fit a concept of autochthonous species (S-selected, capable of sustained development in soil and little affected by fluctuations in resource availability), as opposed to zymogenous (R-selected) species that burst into sporadic activity when presented by a suitable resource, then return to a quiescent state after sporulation or the formation of resting structures [77].Collectively, the ecological traits outlined above are compatible with both the likely truly resident nature and the widespread occurrence of the -core‖ identified fungi.Dead and alive mycelia and insects comprise a significant proportion of soil biomass everywhere, and interestingly, chitinolytic bacteria assigned to Chitinophaga have also been found to be prevalent in dissimilar soils [68].The widespread occurrence of the -core fungi‖ could therefore be related to the universal distribution and abundance of organic resources suitable for exploitation by them.

Study Sites
Published and unpublished fungal ITS1 and ITS2 sequences from 16 soil samples were examined.These samples had been collected either in environmentally different plots at a single geographic site or in environmentally similar soils at different geographic sites (Table 1).
Two geographic sites encompassed white (Tuber magnatum) and black (T.melanosporum) truffle grounds.Three plots were sampled from a white truffle ground (WTG) in Italy.Two of the plots corresponded to productive areas, where fruit bodies were collected; one plot to a nearby non-productive area [78][79][80].The sequences obtained by Mello et al. [39] from inside and outside -brûlé‖ (burned) soils in a French black truffle ground (BTG) were also analyzed.
Sequences from a Mediterranean site located in the island of Sardinia (Italy) were also studied [43,81,82].Five plots, each corresponding to a specific vegetation type, were sampled at this site.Vegetation types were a tilled vineyard (SAR1), a covered vineyard (SAR2), a managed meadow (SAR3), a wooded pasture (SAR4) and a cork-oak (Quercus suber) semi-natural forest (SAR5).
Four sites located in the Western Alps in Italy (one plot per site each), characterized by serpentine substrates (SER), were sampled [41].Two corresponded to disused asbestos mines, and the two other were pristine ophiolitic sites with serpentine rocks and fibrous asbestos outcrops.
Finally, two contaminated industrial soils (CIS), located on the alluvial deposits of a river bank in the Western Po Plain (Italy), were sampled.Due to the past industrial activities, these soils featured a mixed pollution from organic (mainly nonylphenols) and inorganic (metals) contaminants [83].

Sampling, DNA Extraction, PCR and Pyrosequencing
Sequences were obtained from several soil cores samples (5 cm diameter and 20 cm depth) in each plot (two cores per plot).Samples were independently packed in ice upon collection and transported to the laboratory.Soil cores were then sieved (2 mm) to remove fine roots and large organic debris and stored at −80 °C before use for DNA extraction.
Genomic DNA was extracted from two 500 mg aliquots (three for the SER samples) of each soil sample using the FastDNA SPIN for Soil Kit (Q-Biogene, Rome, Italy), following a modified protocol described by Luis et al. [84].The two DNA samples extracted per soil core were independently amplified using two primer sets fused with the 454 pyrosequencing adapters A (GCCTCCCTCGCGCCATCAG) and B (GCCTTGCCAGCCCGCTCAG).These primer pairs target the non-coding nuclear rDNA ITS region, which has been chosen as the universal DNA barcode marker for fungi [65].The first primer pair (ITS1F (5'-(A)CTTGGTCATTTAGAGGAAGTAA-3')/ITS2 (5'-(B)GCTGCGTTCTTCATCGATGC-3')) is fungal-specific for the ITS1 region and amplifies a fragment of c. 400 bp.The second set (ITS3 (5'-(A)GCATCGATGAAGAACGCAGC-3')/ITS4 (5'-(B)TCCTCCGCTTATTGATATGC-3')) is eukaryote-specific for the ITS2 region and amplifies a fragment of c. 350 bp.The choice of a universal primer set (ITS3-ITS4) for the ITS2 spacer was due to the fact that fungal-specific primers targeting this ITS subregion are only available for some fungal phyla (Ascomycota and Basidiomycota) [85].Polymerase chain reaction (PCR) mixes contained 17.1 μL of sterile water, 2.5 μL 10X of reaction buffer (Sigma), 2.5 μL of each deoxyribonucleotide triphosphate (dNTP 2.0 μM), 0.5 μL of each primer (10μM), 0.4 μL of DNA polymerase (High Fidelity Taq, Roche) and 2 μL of DNA template in a final volume of 25 μL.The DNA was amplified using a T3000 thermal cycler (Biometra, Göttingen, DE).The following program was used: initial denaturation at 94 °C for 3 min, followed by 35 cycles of denaturation at 94 °C for 45 s, annealing at 50 °C for 45 s, extension at 72 °C for 1 min and a final extension at 72 °C for 7 min.The PCR products obtained with the two primer pairs were purified with the Agencourt ® AMPure ® Kit (Beckman Coulter, CA, USA).The quality of these samples was assessed through (i) gel electrophoresis of 5 μl subsamples on 1.5% agarose gel, (ii) evaluation of the AD260/280 ratio calculated using the ND-1000 Spectrophotometer NanoDrop ® (Thermo Scientific, Wilmington, DE) and (iii) analysis with the Experion™ System (Bio-Rad, Hercules, CA, USA), using a DNA1K Chip.At this stage, the different amplicons obtained from each studied plot were pooled before sequencing.
Pyrosequencing was performed by BMR Genomics s.r.l.(Padua, Italy) using 15 out of the 16 available lanes in the 454 Genome Sequencer FLX System (454 Life Science Branford, CT, USA), by means of the GS FLX Standard Reagent Series Kit.The obtained number of sequences differs among samples, since in some cases, the 1/16th line has been used to sequence other genomic regions [43].For the second white truffle ground soil sample (WTG2), problems during pyrosequencing occurred, leading to an absence of reads for the ITS2 subregion.
single isolate [90,91].In phylogenetically diverse bacteria, the rRNA operon copy number per genome can fluctuate in response to resource availability in the environment [92].Taken together, these observations caution against using read numbers as an estimate of actual abundance in nature when dealing with environmental samples [93,94].We therefore decided not to take read numbers into consideration.Hence, data from the QIIME OTU table files were reduced to binary data (OTUs were counted as present or absent) in Excel.
Beta-diversity was estimated using Jaccard indices for pairwise comparisons-determined as the ratio of the number of OTUs shared and the total number of OTUs in both samples.Jaccard indices were computed for each of the four rarefied datasets.The fungal assemblages retrieved in the different soil samples were also compared by means of Principal Components Analysis (PCA), with -standardized‖ option, carried out using the SYNTAX 2000 package.Correlations with the original variables were also analyzed with the -Mixed (Rohlf) biplot‖ option.Since replication is imperative [95], only OTUs occurring in either all five soil sample groups (WTG, BTG, SAR, SER and CIS) or all samples within a group were included in the analyses.The Jaccard indices were also used to calculate Mantel statistics (correlation coefficients) between community and geographic distance matrices.Community distance D was computed as the equivalent of the inverse of similarity as expressed by the Jaccard index J (D = 1 -J).Geographic distance matrices were calculated using both straight line distances and using log-transformed distances [28].Mantel analyses (simple Mantel tests with 10,000 iterations) were performed in ZT [96].
Comparisons between the ITS1 and ITS2 datasets were performed following rarefaction to 150 sequences per sample for both regions.Since the ITS1F/ITS2 primer pair is reported to be biased towards amplification of basidiomycetes [47], differences in the proportions of basidiomycetous sequences in the ITS1 and ITS2 datasets were tested for significance by means of the chi-square test (degrees of freedom = 1) performed with the XLSTAT software (Addinsoft; [97]).

Conclusions
To our knowledge, this study represents the largest attempt so far to comprehensively identify patterns of soil fungal beta-diversity at multiple scales over long distances, taking advantage of the high resolution of NGS techniques, at a relatively narrow taxonomic breadth.By means of comparative analysis of pyrosequencing data, we have found evidence of spatial distance acting on the assembly of soil fungal communities.We have also uncovered a core mycobiome shared by unlike soil types.The ecological features of the generalist fungal species we identified are consistent with the potential for continuous, active growth in the oligotrophic soil environment.Technical advances in NGS, as well as more extensive sampling, will allow to unearth further core taxa.The analysis of soil fungal metagenomes and metatranscriptomes (e.g., [40,98]) will contribute to deciphering fungal strategies in soil, increasing our understanding of microbial ecology in this environment.

Figure 1 .
Figure 1.Accumulation curves describing the observed number of ITS1 and ITS2 operational taxonomic units (OTUs) as a function of the sequencing effort per soil plot (rarefaction index values were calculated at intervals of 10 sequences).WTG: white truffle grounds; BTG: black truffle grounds; SAR: Sardinian soils; SER: serpentinitic substrate samples; CIS: contaminated industrial soils.

Figure 2 .
Figure 2. Principal Component Analysis of fungal assemblages in WTG (WTG 1-3), BTG (BTG 1-2), SAR (SAR 1-5), CIS (CIS 1-2) and SER (SER 1-4) soil samples, based on presence/absence data of ITS1 or ITS2 OTUs occurring in either all five soil sample groups or all samples within a group.Percentage variance values accounted for by the two first ordination axes are reported along each axis.

Figure 3 .
Figure 3. Pairwise Jaccard similarity indices between fungal communities from the same sample group or the other groups.For both ITS1 and ITS2, Jaccard indices were computed following rarefaction to 150 sequences per sample for both ITS1 and ITS2, to compare each fungal community with each of the other fungal communities from either the same sample group (blue columns) or the other groups (red columns).Mean values (columns) and standard deviations (bars) are reported for each comparison.Couples of columns marked with an asterisk indicate a significant difference between indices computed for intra-group pairwise comparisons and indices computed for inter-group pairwise comparisons (p < 0.05; Kruskal-Wallis test).

Figure 4 .
Figure 4. Proportion of the ITS1 and ITS2 fungal OTUs assigned to different fungal phyla (according to the NCBI Taxonomy) within each sample (following rarefaction to even sequencing depth: 254 and 158 ITS1 and ITS2 sequences per sample, respectively).WTG: white truffle grounds; BTG: black truffle grounds; SAR: Sardinian soils; SER: serpentinitic substrate samples; CIS: contaminated industrial soils.

Figure 5 .
Figure 5. Heat maps of the OTUs occurring in all localities showing the frequency of each OTU in each sample group.Each column represents a different sample group; each row represents a different OTU; the color of the cells represents the frequency of that OTU within the sample group (percentage of the group samples the OTU was retrieved from).

Table 1 .
Origin of the soil samples examined.

Table 2 .
ITS1 OTUs occurring in at least three soil sample groups (WTG, BTG, SAR, SER and CIS).Occurrence in the examined soil samples and sample groups, total number of clustered reads and taxonomic assignment are reported for each OTU.OTUs are listed according to decreasing frequency (no.soils and samples they were retrieved from).

Table 3 .
ITS2 OTUs occurring in at least three soil sample groups (WTG, BTG, SAR, SER and CIS).Occurrence in the examined soil samples and sample groups, total number of clustered reads and taxonomic assignment are reported for each OTU.OTUs are listed according to decreasing frequency (no.soils and samples they were retrieved from).