Grapevine roots and soil samples were collected from three different vineyards in Emilia Romagna (Italy) across two timepoints. In particular, from each site, i.e., Bondeno (44.953 N/11.305 E, Ferrara), Finale Emilia (44.839 N/11.285 E, Modena), and Medolla (44.816 N/11.062 E, Modena; Figure 1), 15 plants and two bulk soils were retrieved both in June and November 2021 (immediately after the grape harvest), for a total of 90 root and 12 soil samples. Furthermore, for each root sample, the rhizospheric and endophytic microbial communities were both analyzed. The three vineyards considered in this study were characterized by different agronomic managements and biogeographical features. Specifically, both Finale Emilia and Medolla sites are located inside a protected designation of origin (PDO) area but differ in terms of the agricultural approach employed (chemical-based vs. organic, respectively) while the Bondeno site is found outside the PDO area and a traditional chemical-based farming approach is used. A schematic summary of samples distribution across the three sites and the two timepoints (June and November) is provided by Supplementary Table 1.

For the microbiome characterization, each plant root was investigated considering two different ecosystems, namely rhizospheric soil and root endophytic ecosystem, and samples of bulk soil were also analyzed, for a total of 102 samples (Supplementary Table 1). Grapevine thin lateral roots were collected after digging 10–20 cm under the plants. Bulk soil samples were collected near the area where plants were located, after removing the top centimeters of surface soil and the grass cover, if present. All samples were collected wearing sterile gloves, placed inside a sterile 50 ml Falcon tube and stored at −80°C until further processing.

In order to separate the two plant compartments (i.e., rhizosphere and endosphere), roots were thoroughly treated as previously described in D’Amico et al. (2018). In brief, approximately 3 cm of terminal roots portions, including tips, were dissected using sanitized scissors and tweezers to standardize the quantity of starting material. Then, root segments were placed in 15 ml Falcon tubes filled with 2.5 ml of modified PBS buffer (130 mM NaCl, 7 mM Na₂HPO₄, 3 mM NaH₂PO₄, pH 7.0, and 0.02% Silwet L-77) and left on a shaking platform at 180 rpm for 20 min to perform washing. After removing the roots, the washing buffer was centrifuged at 1.500 × g for 20 min and the resulting pellet was regarded as the rhizospheric soil. Roots were then re-washed under the same conditions and transferred to another 15 ml Falcon tube containing 2.5 ml of modified PBS buffer before undergoing 10 cycles of sonication as follows: 30-s pulses at 160 W with 30-s breaks in an ultrasonic bath (Branson 1800, Branson Ultrasonic Corporation, Danbury, CT, United States). After washing and sonication, roots were grinded by means of mortar and pestle in order to reach the root inner portions. This procedure led to a total of 180 samples (90 rhizospheres + 90 roots) that were analyzed together with the 12 bulk soil samples. All samples were kept frozen at −80°C until genomic DNA extraction.

Total genomic DNA was extracted from all the 192 samples, i.e., bulk soils (0.25 g), rhizospheres (approximately 0.25 g), and smashed roots, using the DNeasy PowerSoil Kit (QIAGEN, Hilden, Germany) following the manufacturer’s instructions with minor modifications: a FastPrep instrument (MP Biomedicals, Irvine, CA, United States) was used for the homogenization step with a cycle consisting of three 1-min steps at 5.5 movements per sec with 5-min incubations in ice between each run and, at the end of the protocol, DNA elution was preceded by a 5-min incubation in ice. Then, DNA was quantified by using NanoDrop ND1000 (NanoDrop Technologies, Wilmington, DE) and diluted in PCR grade water to the final concentration of 5 ng/μl before amplification. Five microliters of diluted DNA were used as template for the PCR reaction. PCR was performed in a final volume of 50 μl containing 25 ng of genomic DNA, 2X KAPA HiFi HotStart ReadyMix (Roche, Basel, Switzerland) and 200 nmol/L of 341F and 785R primers carrying Illumina overhang adapter sequences for amplification of the V3–V4 hypervariable regions of the 16S rRNA gene. Specifically, the thermal cycle consisted of initial denaturation at 95°C for 3 min followed by 25 cycles of denaturation (95°C for 30 s), annealing (55°C for 30 s), and extension (72°C for 30 s), with a final extension step at 72°C for 5 min (Turroni et al., 2016). PCR amplicons were cleaned up with Agencourt AMPure XP magnetic beads (Beckman Coulter, Brea, CA, United States). Indexed libraries were prepared by limited-cycle PCR using Nextera technology. Indexing was followed by a second clean-up step, as already described, and then libraries were quantified using Qubit 3.0 fluorimeter (Invitrogen), normalized to 4 nM and pooled. Before sequencing, the sample pool was denatured with 0.2 N NaOH and diluted to 4.5 pM with a 20% PhiX control. Sequencing was performed on Illumina MiSeq platform using a 2 × 250 bp paired end protocol, according to the manufacturer’s instructions (Illumina, San Diego, CA, United States). Sequencing reads were deposited in the ENA archive with the accession code PRJEB57815.

Raw sequences were analyzed using a pipeline which combines PANDAseq 2.11 (Masella et al., 2012) and QIIME2 2021.8.0 (Bolyen et al., 2019) for all 192 samples. High-quality reads (min/max length = 350/550 bp) were retained thanks to the “fastq filter” function of the Usearch 11.0.667 algorithm (Edgar, 2010) and then binned into amplicon sequence variants (ASVs) using DADA2 2021.8.0 (Callahan et al., 2016). Samples with less than 1,000 high-quality reads were discarded and not used for subsequent analyses. The VSEARCH algorithm 2021.8.0. (Rognes et al., 2016) and the SILVA database (December 2017 release; Quast et al., 2012) were employed for taxonomic classification. All unassigned and eukaryotic sequences were discarded. ASVs table was then rarefied to retain a number of 1,138 sequences per sample. The QIIME2 feature table rarefy plugin was used to perform rarefaction. Statistical analyses were carried out using the R software (R Core Team; www.r-project.org—last access: April 2022), v. 4.2.0, implemented with the packages “Made4” 1.72.0 (Culhane et al., 2005), “vegan” 2.6–4 (https://cran.r project.org/web/packages/vegan/index.html—last access: October 2022), “pairwiseAdonis” 0.4 (Martinez Arbizu, 2020), and STAT 0.1.0 (https://cran.r-project.org/web/packages/STAT/index.html—last access: April 2019). Beta diversity based on unweighted UniFrac distances was computed, and the separation of data in the Principal Coordinates Analysis (PCoA) was assessed with a permutation test with pseudo-F ratios (function “adonis” in the vegan package and function pairwiseAdonis in the homonymous package). Kruskal–Wallis test was used to assess significant differences in alpha diversity distributions between groups. Bacterial genera with the largest contribution to the ordination space were detected by the function envfit of the R package vegan on the genus relative abundances. p values, when necessary, were corrected for multiple testing by means of the Benjamini-Hochberg method, with a false discovery rate (FDR) ≤ 0.05 considered to be statistically significant. Linear discriminant analysis (LDA) effect size (LEfSe, Segata et al., 2011), aimed at identifying discriminant rhizospheric taxa between vineyards located inside and outside the Lambrusco DOC PDO viticultural area, regardless of the agricultural practices employed, was performed on genus-level relative abundance tables, retaining only taxa with LDA score threshold of ±2 (on a log10 scale) and value of p ≤ 0.2. The online Galaxy Version interface (https://huttenhower.sph.harvard.edu/galaxy/, last accessed in October 2022) was used to run LEfSe. All taxa identified by LEfSe, thus significantly enriched either inside or outside PDO sites, were then tested for their putative ability to support plant growth by the presence of some well-known plant growth-promoting (PGP) genes. For this purpose, starting from the QIIME2 genera level taxonomic assignment, an oligotyping procedure (Eren et al., 2014) was implemented to detect the species belonging to the genera previously identified by LEfSe through the “Minimum Entropy Decomposition” (MED) module and the global earth microbiomes (GEM) catalogue (November 2020 release; Nayfach et al., 2021). For each genus, the command line was “decompose <ASVs representative sequences fasta.file> − g-M 1-V5.” The-M integer defines the minimum substantive abundance of an oligotype, and the-V integer defines the maximum variation allowed in each node. The node representative sequence of each oligotype was used for species profiling with the QIIME2 feature-classifier plug-in (Bolyen et al., 2019), selecting the VSEARCH algorithm 2021.8.0 (Rognes et al., 2016) and the GEM database (Nayfach et al., 2021). Then, the aminoacidic sequence of some well-characterized PGP proteins, obtained from the reference sequence of the NCBI protein database (https://www.ncbi.nlm.nih.gov; accessed from the 1st to the 31st October 2022), was recovered and blasted against the non-redundant protein sequences NCBI database selecting as target organisms for our queries the bacterial taxa identified by the oligotyping procedure. Unclassified members of a specific taxon were considered when it was impossible to assign the ASVs at the species level, or when the number of oligotypes not assigned at the species level but assigned at higher taxonomic levels overcame the number of species-level matches. Alignments were filtered according to a query coverage of at least 40% and an alignment percentage of identity of at least 20%. The PGP functions selected for our analysis were: nitrogen fixation, phosphorous solubilization, iron chelation, production of the phytohormone indole-3-acetic acid (IAA), and production of the enzyme 1-aminocyclopropane-1-carboxylic acid (ACC) deaminase. For each of these functions, we selected some marker genes from previous scientific literature resulting in 15 protein sequences recovered and blasted. Specifically, the chosen marker genes were NifB, NifE, NifH, NifN, NifV, and NifU (i.e., nitrogen fixation genes) for nitrogen fixation, the alkaline phosphatase phoA, and the glucose dehydrogenase GDH for phosphorous solubilization, three markers of two relevant bacterial siderophores for iron chelation (namely EntF/EntS for enterobactin and FslA for rhizoferrin), three genes directly involved in IAA synthesis (i.e., ipdC, aro10, and aldH) and the AcdS gene encoding the enzyme ACC deaminase (see Supplementary Table 3 for genes accession and version numbers). Moreover, the presence of the same marker genes was verified across the entire rhizospheric microbiome by means of Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2 v. 2.5.0) analysis (Douglas et al., 2020). Notably, during the process of ASVs sequences matching to the KEGG database (Wixon and Kell, 2000, queried on January 10, 2023), two out of 15 reference proteins (i.e., FslA and aro10) were not found in the KEGG database (Supplementary Table 4).

Bacterial co-abundance groups (CAGs) were determined as formerly described by Schnorr et al. (2014). In brief, the Kendall correlation test was used to evaluate the associations among bacterial genera, which were visualized using hierarchical Ward clustering with a Spearman correlation distance metric and used to define CAGs at the genus level. The significant associations observed were controlled for multiple testing with the q value method (FDR ≤ 0.05; Dabney et al., 2013). Permutational multivariate ANOVA (PERMANOVA; Anderson, 2001) was employed to verify whether the CAGs were significantly different from one another. The Wiggum plot network analysis was carried out using cytoscape software v. 3.9.1 (http://www.cytoscape.org/, last accessed in November 2022) as previously described (Claesson et al., 2012).

A total of 90 V. vinifera Cultivar Lambrusco roots samples and 12 bulk soil samples were taken from three different viticultural farms in June and November 2021 in Emilia Romagna, Italy (Figure 1). In particular, from each vineyard (located in Bondeno, Finale Emilia and Medolla) 30 roots (15 in June and 15 in November) and four bulk soils (two in June and two in November) were retrieved, resulting in 102 samples. Among those, all the 90 roots were treated as previously described in order to separate the rhizospheric from the endophytic compartment, leading to a total of 180 V. vinifera samples and 12 bulk soil samples. The selected farms were characterized by different designation of origin and by different agricultural practices: (i) Bondeno (non-PDO area, conventional farming), (ii) Finale Emilia (PDO area, conventional farming), and (iii) Medolla (PDO area, organic farming). For the three sites and the two timepoints, microbiome compositional structure was investigated by NGS sequencing of the 16S rRNA gene (V3–V4 hypervariable regions), resulting in ≃1,5 M high-quality reads, with an average of 9,581 ± 2,329 reads per sample (mean ± SD), which were binned in 31,264 ASVs (samples with less than 1,000 high-quality reads were not analyzed).

Firstly, the bulk soil microbiome was characterized by a significantly higher degree of biodiversity with respect to both rhizospheric and endophytic compartments (p ≤ 0.05, Kruskal–Wallis test, Supplementary Figure 1). When we sought for differences among farms, we only observed a gradual increase of the soil biodiversity from Bondeno, to Finale Emilia and Medolla, with a trajectory that mirrored the path from non-PDO to PDO area and from conventional to organic management. However, these differences are only appreciable at soil level, with the rhizosphere and root compartments from the different farms showing comparable levels of biodiversity (Figure 2).

Beta-diversity analysis revealed a clear pattern toward segregation of the rhizospheric microbial communities according to the sampling location, but not to the sampling season, as shown by the unweighted UniFrac distances (permutation test with pseudo F-ratio, p ≤ 0.001; Figure 3A). Interestingly, the same PCoA indicates that a similar trend can be observed also for the bulk soils, as if the differences detected into the rhizosphere compartment mirrored differences in the soil. In order to identify those bacterial genera most contributing to the separation of the rhizospheric samples in the PCoA, the relative abundances of such taxa were superimposed in the unweighted UniFrac beta diversity plot (Figure 3B). Our results indicate that some bacterial genera are more represented in a particular farm regardless of the season. Specifically, the genus Pirellula, Micromonospora, and Nocardioides are the most characteristic of the Bondeno farm, while Pseudomonas, Flavobacterium, Acinetobacter, Pir4 lineage, and Planctomyces can be associated with Finale Emilia samples and, finally, Skermanella, Gaiella, Solirubrobacter, and Rubrobacter are the most distinguishing of the Medolla farm. Conversely, it is noteworthy to point out that the only ASVs detected across the entire rhizospheric cohort were assigned to uncultured members of the Planctomycetaceae and to uncultured members of the Tepidisphaeraceae. For these taxa, coefficients of variations were 0.4 (mean ± SD % rel. ab., 3.6 ± 1.5) and 0.6 (mean ± SD % rel. ab., 1.7 ± 1.0), respectively, meaning that these taxa were present in the rhizospheres at comparable levels, independently of site and season, constituting a sort of core bacterial group for V. vinifera cultivar Lambrusco. Interestingly, the root samples show no significant structural differences across different sites and seasons (permutation test with pseudo F-ratio, p > 0.05) and a sharp segregation appears in the PCoA only when comparing the endophytic cluster with the entire set of the bulk soil samples (permutation test with pseudo F-ratio, p ≤ 0.001; Supplementary Figure 2).

The LEfSe was finally used to identify rhizospheric bacterial genera that discriminated PDO-associated from non-PDO microbiomes, regardless of farming site, season and type of management (Figure 4). In particular, genera associated with Lambrusco DOC PDO area were Bacillus, Pseudarthrobacter, unclassified members of the order Gaiellales, Planctomyces, Skermanella, Pir4 lineage, Microlunatus, and Paenibacillus. On the other hand, genera less representing the PDO area were Nocardioides, Micromonospora, and Pirellula, unclassified members of the family Gemmatimonadaceae and of the order Acidimicrobiales, Mycobacterium, Legionella, and Chthoniobacter. Notably, such taxa were also generally more represented into the correspondent bulk soil microbiome, with the exception of Pseudarthrobacter and Pir4 lineage for what concerns the PDO area and Nocardioides, Legionella, and Chthoniobacter for what concerns the non-PDO area (Supplementary Table 2).

To identify specificities of the rhizospheric microbiome structure associated with the Lambrusco DOC PDO region and seasonality, we established co-abundance associations of genera and then clustered correlated bacterial taxa into four co-abundance groups (CAGs), describing microbiome configurations across the entire dataset (Supplementary Figure 3). The dominant (i.e., the most abundant) genera in these CAGs were Pirellula (red), Nocardioides (blue), Pseudomonas (pink), and Bacillus (green). The network-establishing CAGs relationships are named Wiggum plots, where genera abundances are represented as a circle proportional to the genus normalized over-abundance (Figure 5). The microbiome variation from Bondeno to Finale Emilia and Medolla through the two different seasons was accompanied by distinctive CAGs dominance, and most relevantly by abundances of the Pirellula and Nocardioides CAGs (Bondeno), the Pseudomonas CAG (Finale Emilia) and the Bacillus CAG (Medolla). When we sought for shared network topological features among Finale Emilia and Medolla microbiome structure, distinctive of the PDO region and not included in the control site (Bondeno), we found that nodes corresponding to Bacillus and Rhizobium were over-abundant during spring, whereas Pseudarthrobacter and Microlunatus nodes were over-abundant in the fall season. When combined, such results underline a sort of seasonal dynamic, very peculiar to the Lambrusco DOC PDO region independently of the type of management. Remarkably, most of these taxa constitute a subgroup of the species previously identified by LEfSe.

Plant growth-promoting microorganisms regulate plant physiological reactions and foster plant growth with several mechanisms. Here, we sought for some of these functions within the reference genomes of the taxa revealed by LEfSe. Specifically, we first used oligotyping (Eren et al., 2014) to identify the bacterial species (or higher taxonomic levels in some cases, as explained above) nested by the ASVs sequences belonging to the genera identified by LEfSe. In particular, ASVs sequences coding for Bacillus, Pseudarthrobacter, Planctomyces, Paenibacillus, Microlunatus, Skermanella, Pir4 lineage, and uncultured members of Gaiellales (i.e., the PDO-related taxa identified by LEfSe) along with ASVs sequences coding for Chthoniobacter, Nocardioides, Micromonospora, Pirellula, Legionella, Mycobacterium, and uncultured members of Acidimicrobiales and Gemmatimonadaceae (i.e., the non-PDO taxa identified by LEfSe) were processed using the “Minimum Entropy Decomposition” (MED) module and the global earth microbiomes (GEM) catalogue (Nayfach et al., 2021). We found that Bacillus korlensis, Bacillus mediterraneensis, Bacillus tuaregi, Bacillus niacini, Bacillus jeotgali, Bacillus lonarensis, and Bacillus litoralis, unclassified species of the genus Bacillus, species of the genus Pseudarthrobacter, species of the genus Planctomyces, unclassified species of the order Gaiellales, Azospirillum brasilense, Azospirillum thiophilum, species of the Pirellulales order, Microlunatus phosphovorus, Paenibacillus castaneae, Paenibacillus harenae, Paenibacillus ferrarius, Paenibacillus beijingensis, and Paenibacillus uliginis, and unclassified species of the genus Paenibacillus, were the taxa characterizing the PDO-area. Notably, when applying oligotyping and GEM database, the ASVs sequences previously assigned to Skermanella were assigned to A. brasilense and A. thiophilum. We chose to retain both Skermanella and Azospirillum genomes for the following analysis, also because of the high level of overlapping found between the 16S rRNA sequences of these two taxa (Zhu et al., 2014). On the other hand, Nocardioides massiliensis, Nocardioides allogilvus, Nocardioides exalbidus, Nocardioides halotolerans, and Nocardioides szechwanensis, unclassified species of the genus Nocardioides, Micromonospora cremea, Micromonospora marina, Micromonospora nigra, Micromonospora sediminis, Gemmatirosa kalamazoonesis, Pirellula staleyi, Legionella fallonii, Legionella saoudiensis, Mycolicibacterium moriokaense, and Mycolicibacterium sphagni, unclassified species of the order Acidimicrobiales and Chthoniobacter flavus were the taxa most distinguishing the non-PDO area. Interestingly, the oligotyping procedure and the GEM database identified Mycolicibacterium species nested in the ASVs belonging to the genus Mycobacterium. In this regard, a recent comprehensive phylogenomic study by Gupta et al. (2018) revealed that Mycolicibacterium can be actually regarded as a distinct clade previously classified as Mycobacterium and now forming a novel microbial genus. Then, in order to investigate the presence of potential PGP traits related to all these microorganisms, the NCBI reference genomes of all of these taxa were scanned for genes associated with nitrogen fixation (essential for plant growth), phosphorous solubilization (important for plant P uptake), siderophore production (for growth in iron-limiting conditions), indole-3-acetic acid (IAA) phytohormonal secretion (beneficial to increase water and nutrient absorption), and 1-aminocyclopropane-1-carboxylic acid (ACC) deaminase (for ethylene precursor degradation and regulation of plant stress response). Among the taxa identified by LEfSe as PDO-characterizing, species belonging to the genera Bacillus, Skermanella/Azospirillum, Paenibacillus, and unclassified species of the order Pirellulales contained most of the PGP features, whereas species from the genera Pseudarthrobacter, Planctomyces, and Microlunatus, together with unclassified members of the order Gaiellales, contained only one or two out of the five investigated PGP traits (Figure 6A). Conversely, if we look at the microbial taxa related to the non-PDO area (Figure 6B), species of the genus Nocardioides (i.e., unclassified Nocardioides and N. exalbidus) are the only ones in which more than two PGP traits out of five have been detected. Furthermore, two important PGP functions, namely nitrogen fixation and IAA production, have been scarcely observed in non-PDO related taxa (with the first only detected in species of unclassified Nocardioides while the latter entirely absent in non-PDO related taxa). PICRUSt2 confirmed most of the findings (Douglas et al., 2020), with some exceptions, above all for what concerns siderophore production and ACC deaminase production (Supplementary Table 5). This can be attributed at least in part to the fact that two markers used in our analysis and necessary for predicting the functionalities are absent in the databases provided with PICRUSTt2.