The bacterial communities of Tuber aestivum: preliminary investigations in Molise region, Southern Italy

Truffles are colonized by a complex microbial community of bacteria, yeasts, and filamentous fungi, whose role has not yet been fully understood. The main purpose of the research was to characterize the bacterial communities associated with Tuber aestivum Vittad. fruiting bodies collected from natural truffle grounds in the Molise region (Southern Italy). Despite it is one of the Italian richest areas of truffles, little is known about truffles in Molise. Six ripe fruiting bodies of Tuber aestivum Vittad. and six soil samples were collected in July 2018 at Villa San Michele in the municipality of Vastogirardi, Molise region. Then, soil and truffle microbial communities were analyzed through 16S rRNA gene sequencing on the Illumina MiSeq platform and bioinformatics analyses. Consistently with previous studies, the main phyla retrieved in the investigated ascocarps were Proteobacteria and Actinobacteria, with the genus Bradyrhizobium particularly represented. Nevertheless, significant differences between soil and truffle microbiota and an unexpected heterogeneity across truffles were observed. It is likely that a specific recruitment of bacteria from soil to ascocarps occurs during the truffle formation and that local-scale factors play an important role in determining the structure of the investigated truffle microbial communities. Although further analyses (based on a larger soil and truffle sample size and aimed at defining in more detail microbial diversity, soil physical and chemical properties, microclimatic conditions, and vegetation) are required to better understand which are these factors and how they could influence the composition of truffle bacterial communities, this study represents the starting point for a deepened characterization of this economically important product.


Introduction
Most of terrestrial plant roots are colonized by mycorrhizal fungi which play a key role in soils by improving plant water and mineral nutrient uptake and receiving carbon compounds in return (Mello and Balestrini 2018). Among mycorrhizal symbioses, ectomycorrhizae (ECM) are mutualistic relationships between plants (both angiosperms and gymnosperms, as well as shrubs) and fungi (almost exclusively belonging to Basidiomycetes and Ascomycetes). A well-known example of ectomycorrhizal fungi is given by truffles, a polyphyletic group of hypogeous fungi whose fruiting bodies sequester their spores and develop underground, bringing benefits to the forest ecosystems and to the host plant (Pavić et al. 2013;Mäkelä et al. 2015;Vita et al. 2015;Benucci and Bonito 2016;Mello and Balestrini 2018). Different species of truffles are edible and highly sought after due to their organoleptic properties (Vita et al. 2015;Benucci and Bonito 2016). Among these, truffles belonging to the Tuber genus (including more than 180 species) are considered very precious and delicious food (Pacioni et al. 2014;Schmidberger and Schieberle 2017). Tuber magnatum Pico, the Italian white truffle, and Tuber melanosporum Vittad., the Périgord black truffle, represent the most valued species on the food market due to their excellent flavor. Tuber spp. are widely distributed across Europe, North America, South East Asia, and limited parts of Africa and South America (Barbieri et al. 2005;Vita et al. 2015;Benucci and Bonito 2016;Iotti et al. 2016;Ye et al. 2018).
Truffles are colonized by a complex microbial community of bacteria, yeasts, and guest filamentous fungi Vahdatzadeh et al. 2015). Bacteria can be found in inner and outer parts of truffle with densities ranging from a million to a billion cells per gram of fruiting bodies (dry weight) (Vahdatzadeh et al. 2015). They seem to be selected from the soil microbial communities during the early stage of truffle formation and could play an important role in the development, growth, and nutrition of fruiting bodies . Moreover, they contribute to truffle aroma through the production of volatile organic compounds (VOCs) Vahdatzadeh et al. 2015;Benucci and Bonito 2016;Ye et al. 2018). A combination of culture-dependent and independent methods has been used to extensively investigate truffle-associated bacteria Vahdatzadeh et al. 2015;Benucci and Bonito 2016) and showed that different Tuber species may harbor diverse bacterial communities, mainly constituted of Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria phyla (Vahdatzadeh et al. 2015). Nevertheless, recent studies highlighted that the genus Bradyrhizobium is particularly abundant in several Tuber species (Benucci and Bonito 2016).
Many of these investigations have been focused on T. borchii Vittad. and T. magnatum Pico (Sbrana et al. 2002;Barbieri et al. 2005;Barbieri et al. 2007;Pavić et al. 2013;Splivallo et al. 2015) whereas less is known on bacterial communities of T. aestivum Vittad. despite its wider geographic distribution (Büntgen et al. 2017). The main purpose of this research was to characterize the bacterial communities associated with Tuber aestivum fruiting bodies collected from truffle grounds in Molise region (Southern Italy), by using next generation sequencing techniques (NGS). Despite it has received less attention, Molise region represents one of the Italian richest areas of truffles.

Study area and samplings
Six ripe fruiting bodies of Tuber aestivum and six soil samples were collected in July 2018 from natural truffle grounds at Villa San Michele in the municipality of Vastogirardi, Molise region (Southern Italy, Fig. 1a). The study area is located at an altitude of about 900 m above sea level (a.s.l.) with a vegetation composed by Turkey oak (Quercus cerris).
Tuber aestivum ascocarps were collected at a depth of approximately 10 cm in six different sites at a distance of some meters from each other. Truffles were dug out with the help of an expert person and a hunter/truffle dog by using a spade (Fig. 1b, c, d). Then, they were individually placed in sterile polypropylene containers. In addition, soil samples were also collected using sterile spoons and sterile polypropylene tubes under the fruiting bodies, at a depth of about 10-15 cm. Truffle and soil samples were then transported in a refrigerated container to the laboratory.
Tuber aestivum ascocarps were gently brushed with a sterile soft brush, rinsed with sterile distilled water and stored at − 20°C before proceeding with biomolecular investigations. Species identification and maturation stage assessment were performed on the basis of morphological traits by experienced mycologists Büntgen et al. 2017).

Biomolecular investigations DNA extraction
DNA was extracted from truffles (1t, 2t, 3t, 4t, 5t, 6t) and soil samples (1s, 2s, 3s, 4s, 5s, 6s) in order to assess bacterial diversity and taxonomic composition using Illumina amplicon sequencing of 16S rRNA genes (sample names with matching numbers indicate truffles and soils collected from the same site within the study area).
Tuber aestivum ascocarps (inclusive of peridium and gleba) were shredded and then pulverized by the addition of liquid nitrogen. Total genomic DNA was subsequently extracted from the fruiting bodies and soils by using the E.Z.N.A.® Plant DNA DS Mini Kit (Ye et al. 2018) and E.Z.N.A.® Soil DNA Kit (Omega Bio-tek, USA) respectively and following the manufacturer's instructions.

16S rRNA gene amplicon library preparation and sequencing
Partial 16S rRNA gene sequences were amplified from extracted DNA using primer pair Probio_Uni (5′-CCTACG GGRSGCAGCAG-3′) and Probio_Rev (5′-ATTACC GCGGCTGCT-3′), which target the V3 region of the 16S rRNA gene sequence (Milani et al. 2013;Di Luccia et al. 2018). Illumina adapter overhang nucleotide sequences were added to the partial 16S rRNA gene-specific amplicons, which were further processed employing the 16S Metagenomic Sequencing Library Preparation Protocol (Part #15044223 Rev. B-Illumina). Amplifications were carried out using a Verity Thermocycler (Applied Biosystems). The integrity of the PCR amplicons was analyzed by electrophoresis on a 2200 TapeStation Instrument (Agilent Technologies, USA). DNA products obtained following PCR-mediated amplification of the 16S rRNA gene sequences were purified by a magnetic purification step involving the Agencourt AMPure XP DNA purification beads (Beckman Coulter Genomics GmbH, Bernried, Germany) in order to remove primer dimers. DNA concentration of the amplified sequence library was determined by a fluorimetric Qubit quantification system (Life Technologies, USA). Amplicons were diluted to a concentration of 4 nM, and 5 μL quantities of each diluted DNA amplicon sample were mixed to prepare the pooled final library.

Bioinformatics analysis
Following sequencing, the .fastq files were processed using a custom script based on the QIIME software suite (Caporaso et al. 2010). Paired-end read pairs were assembled to reconstruct the complete Probio_Uni/Pro-bio_Rev amplicons.
Quality control retained those sequences with a length between 140 and 400 bp and mean sequence quality score > 20 while sequences with homopolymers > 7 bp and mismatched primers were omitted.
In order to calculate downstream diversity measures (alpha and beta diversity indices, Unifrac analysis), 16S rRNA Operational Taxonomic Units (OTUs) were defined at 100% sequence homology using DADA2 (Callahan et al. 2016); OTUs not encompassing at least 2 sequences of the same sample were removed. Notably, this approach allows highly distinctive taxonomic classification at single nucleotide accuracy (Callahan et al. 2016). All reads were classified to the lowest possible taxonomic rank using QIIME2 (Caporaso et al. 2010;Bokulich et al. 2018) and a reference dataset from the SILVA database (Quast et al. 2013).
Biodiversity within a given sample (alpha-diversity) was calculated with Shannon and Chao1 indices. Similarities between samples (beta-diversity) were calculated by weighted UniFrac (Lozupone and Knight 2005). The range of similarities is calculated between values 0 and 1. PCoA representations of beta-diversity were performed using QIIME2 (Caporaso et al. 2010;Bokulich et al. 2018).
To identify the core taxa retrieved in the soil and truffle habitats, microbial genera were analyzed for their presence or absence in each sample and the related Venn diagram drawn by Venny 2.1 (Oliveros 2007(Oliveros -2015. For differential abundance analysis between soil and truffle microbial communities at phylum level, a negative-binomial-based approach in tandem with paired Wald test, as available in DESeq2 version 1.24.0 (Love et al. 2014) in R environment (R core team 2019), was performed.
In order to structure soil, truffle, or soil-truffle microbiomes, co-occurrence patterns were determined by applying rcorr function from Hmisc package in R (Harrell Jr et al. 2019) and calculating correlations among species abundances by Pearson method. Correlation coefficient significance was also assessed by using the same R function.
Networks were composed by selecting highly significant correlations (P < 0.01) between species abundances. Specifically, networks were derived from abundances of species identified in the soil or in the truffle habitat and also from a dataset composed by the species identified in the soil and truffle habitats together. The analysis and visualization of the networks and the related statistical analysis were performed by Cytoscape (Shannon et al. 2003) and R.
Network structure was assessed by measuring the node degree of each node that is the number of partners each node has. Node degree is related to sparsity, a property that relies on the number of connections observed in the network (Busiello et al. 2017). In a "complete network," each node is connected with all other nodes; thus, the average number of node degree is equal to the number of nodes minus 1. A lower mean node degree is a sign of a sparser network. The larger is the network the higher is the sparsity (Busiello et al. 2017). The density, as the number of observed edges respect all possible ones, shows how dense are the connections per node. In biological networks, density has been demonstrated to be lower than 0.1 (Leclerc 2008). This observation indicates that the structure of a biological network fits with a general sparse connected graph which gives evolutionary advantages in terms of resilience to possible network damages (Leclerc 2008;Pavlopoulos et al. 2011).
To have a stronger evidence about network structure, the centrality of the microbial populations present in the combined habitat network was measured in order to define the relative importance of soil species with respect to the truffle ones. To achieve this goal, node degree distributions of soil-truffle network were analyzed by considering first the totality of the nodes present, then only the nodes related to species identified in the soil habitat, and finally only the species identified in the truffle habitat. The frequency distribution of node degree provides information about the structure of a network (sparsity), and it is often used to compare the nature of networks (Liu et al. 2011;Nacher and Akutsu 2013;Suweis et al. 2013;Grilli et al. 2017). In order to achieve this goal, Welch's unequal variance t test (Welch 1947) on node degree distributions was applied.
Nodes relevant to the structure of the microbiomes (keystone nodes) were identified according to their node degree and betweenness centrality that is a measure of the shortest paths that pass through the node (Freeman 1977).
To test the robustness of the network and its ability to be resilient to species loss, an approach based on the progressive removal of species (nodes) according to their relevance in the network structure (Iyer et al. 2013;Ruiz et al. 2017), both considering node degree and betweenness centrality, was applied. The area under the curve (AUC) was calculated according to Gini's formula.

Results
NGS results allowed to obtain detailed information about the composition of microbial communities in soil samples and associated with summer black truffle (Tuber aestivum) ascocarps collected from six different truffle grounds in the Molise region (Fig. 1).
Actually, even though it seems that truffle could host peculiar bacterial communities including taxa absent in soil, those genera retrieved only in the truffle group were not found in all the analyzed six ascocarps and were often present at very low relative abundances.
Overall, important differences were found between soil and truffle groups even though a significant heterogeneity of ascocarp microbiota has also been detected. Principal Coordinate Analysis (PCoA) based on weighted UniFrac index revealed that the truffle microbial communities varied from those of the soil (Fig. 2). The rarefaction curves obtained by using Shannon and Chao1 indices for each sample highlighted a greater microbial diversity in soil samples compared to the Tuber aestivum fruiting bodies and differences among samples of the truffle group in terms of α-diversity (Additional file 2). Differential abundance analysis at phylum level revealed the presence of 14 categories with significantly different relative abundances between soil and truffle (Fig. 3). Among these, the phyla Tectomicrobia, Nitrospirae, Fibrobacteres, Planctomycetes, Gemmatimonadetes, Chloroflexi, and Acidobacteria were more abundant in soil whereas Proteobacteria, Saccharibacteria, Firmicutes, Cyanobacteria, Fusobacteria, Tenericutes, and other unclassified Fig. 2 Principal Coordinate Analysis (PCoA). Plot was generated using weighted UniFrac distance matrix. Soil (black spheres) and truffle (gray spheres) samples are shown members of the Bacteria domain were found at significantly higher relative abundances in truffle (Fig. 3).
When analyzing the bacterial community composition on the whole, the main phyla retrieved in the analyzed Tuber aestivum fruiting bodies were Proteobacteria, Actinobacteria, Bacteroidetes, and Firmicutes (Fig. 4). Proteobacteria ranked first in all truffle samples, with a relative abundance between 57.6 and 95.9%. The phylum Actinobacteria showed percentages ranging from 2.2 to 25.1%, with Nocardioides being one of the main bacterial genera present in the ascocarps. Bacteroidetes were found in one of the truffle samples with a relative abundance of 25.9%, while in the other fruiting bodies they ranged between 0.8 and 4.7%. Firmicutes occurred in percentages ranging from 0.2 to 12.1%.
Among Proteobacteria, the main families were represented by Bradyrhizobiaceae, Rhizobiaceae, Comamonadaceae, and Hyphomicrobiaceae, with the genera Bradyrhizobium, Rhizobium, and Devosia present in all the analyzed summer black truffles at relatively high percentages (up to 40.9%, 6.5%, and 14.5%, respectively). The family Comamonadaceae showed variable values ranging from 1.9 to 73.8%, with the genera Acidovorax and Polaromonas poorly represented in most of fruiting bodies but with a relative abundance of 67.6% and 4.0% in sample 3t, respectively (Fig. 4).
In order to have an insight on the relevance and impact of cross-habitat associations between species, three microbial networks were generated for soil and truffle communities both separately and taken together, based on correlation between species abundances (Fig. 5). The topology was investigated in order to identify the most significative features of each network and keystone species relevant for their structures and organizations.
The network obtained for the soil habitat (Fig. 5a) was characterized by 747 nodes, each representing a species, and 8326 connections among them. The nodes were organized in a large connected component (clique) made of 99.19% (741) of the nodes and three dyads (0.81%). The mean number of partners for each node of this network (mean node degree) was equal to 22.29 whereas network density was 0.0149.
The network obtained for truffle habitat (Fig. 5b) contained 683 nodes and 28,057 connections, and it was characterized by a mean node degree of 82.16 and a Fig. 3 Differentially abundant phyla in soil (on the left) and truffle (on the right) microbial communities network density of 0.0601. The largest component of the network, the clique, was made of 649 nodes which represented the 95.02% of the total nodes. A component of 28 nodes (4.01%) was also present, and there were three dyads (0.97%).
The network inferred by soil and truffle microbial communities combined together (Fig. 5c) was characterized by 1439 nodes organized in a large connected component (clique) comprising 99.72% of the nodes and two dyads (0.38%). This network was also characterized by 53,721 connections, of which 8326 (15.51%) were between species identified in soil habitat, 28,057 (52.27%) were between species from truffle habitat, and 17,299 (32.22%) were across the habitats (a partner from soil habitat, the other from the truffle habitat). The average node degree of the soil-truffle network was 74.59, and the density was 0.026. A high number of significant cross-habitat correlations was found, suggesting that the species from the two habitats interacted or were similarly influenced by the environmental conditions. Subsequently, the centrality of the microbial populations present in the combined habitat network was measured in order to define the relative importance of soil species with respect to the truffle ones. As shown in Fig. 5, panel d, the number of network connections between soil species was strongly and significantly lower with respect to the connections found between truffle species (Welch t test p value < 10 −14 ) within the soiltruffle network. It is likely that those microorganisms recruited in the truffle from the soil must have a closer and stricter linkage to optimize the colonization of that specific ecological niche.
Based on betweenness centrality scores, the top five keystone OTUs in soil communities included unclassified members of the Sedimentibacter, Nitrospira, Gemmatirosa, and Dyadobacter genera, and Acidiferrobacteraceae family (Fig. 6a). For truffles, Kosakonia cowanii, Massilia sp. UMI-21, unclassified members belonging to the genera Piscinibacter and Blastococcus, and to the Polyangiaceae family represented the more relevant bacterial groups (Fig.  6b). On the other hand, the higher stability observed in the network showing soil and truffle microbial communities together seemed to rely on the presence and abundance of OTUs (such as Pseudomonas brassicacearum Fig. 4 Soil (1s,2s,3s,4s,5s,6s) and truffle (1t, 2t, 3t, 4t, 5t, 6t) microbial composition. Relative abundance at the phylum (a) and genus (b) levels. U.m. refers to unclassified microorganisms subsp. brassicacearum, Paenibacillus polymyxa, beta proteobacterium DC2b-18, and unclassified members of the Sedimentibacter genus and Saccharibacteria phylum) defining both the environments, demonstrating a more complex inter-relationship of microbial taxa (Fig. 6c) which were more closely related to one or to the other habitat. In addition, this higher stability was also shown by the analysis of the influence of species loss on the network connectivity through the measure of the "attack robustness" (Fig. 7). This analysis clearly showed that the size of the largest component of soil microbial community network decreased faster compared to truffle or soil-truffle communities whose response to robustness attack differed only of a light skew. Among all, the soil-truffle robustness attack curve showed the highest AUC values (Fig. 7). This suggested also a higher structural similarity between the truffle and soil-truffle microbial communities that is coherent with a specific recruitment of the truffle microbiota from the soil.

Discussion and conclusions
Several methods (cultivation-dependent and molecular approaches) have been employed to reveal microbial community composition and responses to environmental changes in various environments and in different contexts (Bucci et al. 2011(Bucci et al. , 2014(Bucci et al. , 2015a(Bucci et al. , 2015bCrescenzo et al. 2017;Di Luccia et al. 2018;Petrella et al. 2018;Pietrangelo et al. 2018). In fact, an understanding of the temporal and spatial structures, functions, interactions, and population dynamics of microbial communities is critical for many aspects of life, including scientific discovery, biotechnological development, sustainable agriculture, environmental protection, and human health ).
In the present study, we used NGS technology to investigate microbial communities associated with summer black truffle ascocarps collected in Molise region (Southern Italy), one of the richest Italian areas of this product. Despite its economically important value, T. aestivum from Molise has received less attention from a scientific point of view compared to truffles from other Italian regions. Thus, the research aimed at filling the gap in the current knowledge by analyzing bacterial communities associated with fruiting bodies as an initial step for a further and deepened characterization. Each node represents a species, and node color is related to the habitat the species belongs to (blue for soil habitat and yellow for truffle one). Edges (lines) between nodes represent a significant correlation between the abundances of the species connected by. Edge width is proportional to the significance of correlation between species. Networks are plotted by organic layout in Cytoscape. d Species identified in the truffle habitat have a higher connectivity (measured by node degree centrality) than species identified in the soil habitat (respectively yellow and blue boxplots). The network derived from abundances combined from species identified in soil and truffle habitats has a connectivity lower respect the only truffle component and higher respect the only soil component (striped boxplot) The main bacterial phyla retrieved in the Molise truffle were the same found in the fruiting bodies of Tuber spp. of different geographic origin (Vahdatzadeh et al. 2015;Benucci and Bonito 2016;Ye et al. 2018). In fact, the analyzed communities were dominated by Proteobacteria and Actinobacteria, with the genus Bradyrhizobium particularly represented (Gryndler et al. 2013;Benucci and Bonito 2016;Ye et al. 2018).
As expected, summer black truffle microbiota was significantly different from that of soil which showed a higher α-diversity although a high number of shared taxa at genus level. Some bacterial genera, such as Bradyrhizobium and Devosia, were detected with relatively higher abundance values in most of the fruiting bodies compared to the soil, demonstrating that most likely a recruitment of bacteria from soil to ascocarps occurs during the truffle formation, in agreement with previous researches (Antony-Babu et al. 2013;Splivallo et al. 2015).
The presence of genera comprising nitrogen-fixing bacteria (such as Bradyrhizobium and Rhizobium), as already reported for other Tuber species (Le Tacon et al. 2016), is relevant: in fact, the ability to modify nutrient availability during their biological cycle could be of particular importance during the development of fruiting bodies which need nutrients in order to complete the maturation process independently of the host plant (Barbieri et al. 2010).
The network analysis performed to elucidate the interactions between microbial taxa showed a higher robustness of the system when microbial communities were examined as a whole rather than individually. To this greater stability contributed bacterial taxa more strictly related both to one and to the other habitat, indicating Fig. 6 Keystone species obtained from the network analysis. Nodes represent species that are shown agreeing to their centrality metrics (triangles indicate the top 5 keystone species based on node degree and betweenness centrality). Each node of a plot represents a species in the related network (a for soil, b for truffle, and c for soil-truffle). Nodes are colored according to the network they belong to (blue for soil, yellow for truffle, and red for soil-truffle). Keystone species names are showed for some of the most relevant nodes of each network. In the soil-truffle network the apex beside of species name indicates the habitat in which the species was detected: (s) for soil and (t) for truffle Fig. 7 Attack robustness strategy to measure network resilience. Node removal was performed in agreement with their centrality, considering both node degree (a) and betweenness centrality (b) or by a random selection (c). The robustness was measured as percentage of nodes remaining in the main connected component of the network (clique). Here we plotted the results of this analysis performed on each of the three networks as a line (blue for soil, yellow for truffle, and red for soil-truffle). On the X-axis are reported the percentages of nodes removed, while on the Y-axis the percentage of the nodes remaining in the clique. AUC indicates the area under the curve calculated according to Gini's formula complex connections among species that could somehow determine and maintain the equilibrium of this peculiar ecosystem.
Nevertheless, a surprising and remarkable heterogeneity across truffle samples, in terms of microbial community composition and relative abundance of the main taxa, was also observed. Since the study area is small and the host plants belong to the same species, we suppose that local-specific factors could play an important role in determining the structure of the investigated truffle microbial communities. Further analyses are, thus, required to better understand which these factors are and how they could influence the composition of microbial communities. In our opinion, the full comprehension of the role of bacteria in mycorrhization process and their contribution to the development of specific traits in truffles as well as of the factors that drive the establishment of specific microbial communities in the ascocarps, could have a significant impact on truffle industry, mainly at regional scale: for the detection and identification of specific quality marks in the foodstuff (product promotion) and also for the improvement of cultivation techniques in artificial grounds.
Additional file 1. Venn diagram of bacterial communities with shared and unique genera between soil and truffle samples Additional file 2. Rarefaction curves of soil (black lines) and truffle (gray lines) samples. Alpha-diversity plots obtained by using the Shannon (left) and the Chao1 (right) indices