Grafted Falanghina grapevine plants (Vitis vinifera L. cultivar Falanghina) were sampled in September 2021 (just before the harvest season) at the farming co-operative La Guardiense (Benevento, Sothern Italy), which is one of the biggest in Italy and a representative of the wine-growing region of Sannio. Specifically, two different rootstocks were selected: V. berlandieri × V. rupestris PAULSEN 1103 (1103P) and V. berlandieri × V. riparia KOBER 5BB (5BB) ( Table S1 in the Electronic Supplementary Material ). We investigated seven vineyards subjected to the same pedo-climatic conditions, similar age, training system (double Guyot), spacing, and management ( Figure 1 ). Rhizosphere soils were taken from seven 1103P and five 5BB grapevine plants, randomly sampled from three rows of the same vineyard field. Specifically, to sample the rhizospheric soil, three soil cores were taken near the stem of each plant (up to 50 cm). The roots trapped in the soil cores collected near the stem were considered for rhizosphere soil sampling. The three rhizosphere sub-samples were then combined to obtain a homogeneous sample at each plant location. Then, from the three homogeneous samples, a composite one (deriving from the randomly selected plants) was collected. This type of sampling allows us to account for biological variability and reduce sampling bias as much as possible. Overall, 12 rhizosphere soil samples were collected and transported immediately to the Department of Science and Technology (University of Sannio) while maintaining the cold chain. For each sample, 5 g was stored at -80°C for metagenomic analyses.

Additionally, for each vineyard, a composite soil sample was collected from the same three soil cores used for rhizophere soil samples. Specifically, 200 g of soil was stored at -20°C until the agronomic characterization of the soil.

The physical and chemical properties of the soil were determined at each vineyard (seven) following official analysis methods of the soils (Ministero per le Politiche Agricole e Forestali, 1999): texture, pH, cation exchange capacity (CSC), organic carbon (OC), electrical conductivity (EC), and calcium (Ca), magnesium (Mg), potassium (K), phosphorous (P), and nitrogen (N) content. The soil texture was determined according to Method II (Ministero per le Politiche Agricole e Forestali, 1999). The pH reaction was obtained following the DM 13/09/99 G.U. no. 248 of 21/10/99 Method III.1. The cation exchange capacity (CEC) was obtained according to Method XIII.1 using a molar of ammonium acetate. The exchange bases (calcium, magnesium, potassium, and sodium) removed from the exchange sites with barium chloride solution buffered to pH 8.2 were determined by flame atomization atomic absorption spectrophotometry (FAAS), according to Method XIII.5. The electrical conductivity was determined following Method IV.1 (Ministero per le Politiche Agricole e Forestali, 1999). Total N and P were determined by XIV.3 and XV.1 methods, respectively.

Total genomic DNA was extracted using E.Z.N.A.^® Soil DNA Kit (Omega Bio-tek, Inc.) according to the manufacturer’s instructions. Nucleic acid quantification was performed using the Qubit 3.0 fluorometer (Thermo Fisher Scientific, Carlsbad, CA, USA) with Qubit DNA High Sensitivity Assay Kit (Thermo Fisher Scientific). 16S rRNA genes were amplified through PCR using primers for hypervariable regions V2, V3, V4, V6–7, V8, and V9 by using the “Ion 16S Metagenomics Kit” (ThermoFisher Scientific, 5781 Van Allen Way Carlsbad, CA, USA). The libraries were prepared using the Ion Plus Fragment Library Kit (ThermoFisher Scientific) and then sequenced by using Ion 510™ & Ion 520™ & Ion 530™ Kit – Chef on the Ion Torrent ChefTM and the Ion Torrent™ Genestudio S5™ Plus System. We produced around 500,000 reads per sample, on average, by loading 16 samples on 520 ChipTM. The generated unaligned binary data files (Binary Alignment Map, BAM) have been analyzed using Ion Reporter™ Software 5.20 platform with the workflow of Metagenomics 16S w1.1 v5.14. The Ion Reporter workflow is a custom workflow that uses both MicroSEQ™ ID 16S Reference Library v2013.1 and Greengenes v13.5 as reference libraries. All reads shorter than 150 (in base pairs) were filtered out after trimming the primers. A minimum alignment coverage of 90% and a read abundance filter of 10 were set; the threshold value for percentage identity was 99%. 16S rRNA sequence reads were finally clustered into operational taxonomic units (OTUs) using the closed-reference method. The Ion Reporter metagenomics workflow also provides primer information, classification information, percent ID, and mapping information. Concerning the workflow of Metagenomics 16S w1.1 v5.14, a detailed overview is reported in https://assets.thermofisher.com/TFS-Assets/LSG/Vector-Information/ion-reporter-16s-metagenomics-algorithms-whitepaper.pdf.

Soil chemical analyses were presented as mean ± standard deviation; moreover, coefficient of variation (CV) was used to quantify soil feature variations and evaluate the homogeneity of soils.

Microbial community analysis was performed using R v. 3.1.0 (R Core team, 2021) with general dependency on the following packages: reshape2 (v.1.4.4), ggplot2 (v.3.4.2), plyr (v.1.8.8) (Wickham, 2007; Wickham, 2009; Wickham, 2011) Hmisc (v.5.1-0) (Harrell and Dupont, 2019), and phyloseq (v.1.42.0) (McMurdie and Holmes, 2013). The 16s data processed into OTUs were imported into R using phyloseq to perform exploratory analyses, examining the taxonomic composition and diversity of the samples, and assessing the dissimilarity between groups (5BB and 1103P rootstocks). Specifically, the input data were imported as phyloseq object, which is the most commonly used data structure for microbiome data in R environment. For the creation of the phyloseq object, the OTU abundance and OTU taxonomy tsv files were downloaded from Ion Reporter and subsequently converted into csv files, and both were arranged with the additional csv file, including metadata, into a combined object.

In particular, microbiomeMarker, an R/Bioconductor package (Cao et al., 2022), for microbiome marker identification and visualization was used to process both taxonomic and functional metagenomic data and perform differential analysis (DA).

Liner discriminant analysis effect size (LEfSe) was conducted for metagenomic biomarker discovery after total sum scaling (TSS) data normalization. This was chosen since it addresses the challenge of finding taxa or pathways that consistently explain the differences between two or more microbial communities (in our case study, two different rootstock-associated microbial communities), which is a central problem to the study of metagenomics (Segata et al., 2011). LEfSe was based on wilcoxon_cutoff = 0.05, kw_cutoff = 0.05, and lda_cutoff = 2.

Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) (Langille et al., 2013) was applied to infer functional categories associated with taxonomic composition (Ijoma et al., 2021). The generated OTU table was also imported into the PICRUSt package, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to predict the functional potential of bacterial communities. The output from the KEGG database containing the predicted function was further analyzed: a single-factor statistical comparison was performed by Mann–Whitney U -test to evaluate the functional differences between groups. Then, genes involved in plant growth promotion (PGP) and encoding key enzymes in nutrient cycling were accordingly identified in the resulting profiles using their KEGG orthologs (Sashidhar and Podile, 2010; Kaiser et al., 2016; Marasco et al., 2018) and plotted using ggplo2 package.

The soil information were prior retained in the information system managed by the Agriculture Research Council of Italy (Costantini et al., 2014). The available soil subsystems geodatabase (1:250,000) was considered. The sampled vineyard soils are grafted onto flood plains and terraces and were classified as Calciustepts and Haplustepts according to the Keys to Soil Taxonomy (Soil Survey Staff, 2015). All soil samples were sandy clay loam and silty clay loam ( Table 1 ). Based on our analyses, the soil pH was 8.02, on average, typical of alkaline soils, with a low CSC (<10). The organic carbon content ranged from low to normal, with an average value of 11.28 g/kg, and was comparable to data reported for European Mediterranean cropland semiarid soils (Tóth et al., 2013; Grilli et al., 2021). The total N mean value was 0.76 g/kg. To exclude biases due to different soil features, which play a key role in shaping the microbiome layout, the variability of soil characteristics was inspected by coefficient of variation (CV, %), which is a useful statistical parameter for comparing the degree of variation of the samples. The CV data ( Table 1 ) showed that the soil samples were quite homogenous as spread around, and the mean is overall less than 20% (Selmy et al., 2022). However, the statistical analysis demonstrated a high variability of EC with a coefficient of variation of 62.5%. The soil pH and CEC levels showed minimum variability in comparison to other soil properties.

A total of 4,293,111 amplicon sequences from 12 samples were collected ( Table S2 in the Electronic Supplementary Material ). The bacterial community dataset included 883,864 quality-filtered reads, with 18,988 to 130,421 reads per sample ( Table S2 ). We identified 458 taxa based on the conventional criterion of 99% sequence similarity. Alpha-diversity measures ( Figure 2 ) indicate that the two rootstocks’ rhizosphere microbiota are significantly different (p < 0.05). Overall, the bacterial alpha-diversity (observed OTUs) significantly increased in the 1103P rhizosphere microbiota. Interestingly, the 5BB soil microbiome showed the lowest within-group variability, while greater dispersion was observed for the rhizosphere microbiome of 1103P. The Shannon diversity measure ( Figure S1 in the Electronic Supplementary Material ) indicates no significant difference between groups.

Irrespective of the cultivar, Pseudomonadota (synonym Proteobacteria) was the dominant phylum ( Figure 3 ), followed by Actinomycetota (synonym Actinobacteria) > Bacteroidota (synonym Bacteroidetes) > Thermoproteota (synonym Crenarchaeota). Pseudomonadota ranged from 69.19% (V9 sample) to 93.34% (V12 sample). Actinomycetota ranged from 2.56% (V12 sample) to 26.71% (V9 sample). At the class level (Figure 4), the most abundant taxa in all samples were Gammaproteobacteria (64.92% and 57.37%, on average, in the 5BB and 1103P rhizosphere samples, respectively) > Actinobacteria (12.70% and 21.08%, on average, in the 5BB and 1103P rhizosphere samples, respectively) > Alphaproteobacteria (12.37% and 14.62%, on average, in the 5BB and 1103P rhizosphere samples, respectively).

At the family level ( Figure S2 in the Electronic Supplementary Material ), Pseudomonadaceae was the most abundant taxa (with a relative abundance ranging from 48.61% to 88.07%) in all samples. The second most representative family were Rhodospirillaceae (V1, V10, V4, and V5 samples), Nitrososphaeraceae (V11 and V12 samples), Flavobacteriaceae (V2 and V3 samples), Nocardioidaceae (V6, V7, and V8 samples), and Micrococcaceae (V9 samples).

In order to obtain further information on the recruitment of the microbiome that thrives in the rhizosphere of the two different rootstocks, a DA was performed.

The LEfS robustly identified rhizosphere microbial taxa that are statistically different among the two investigated rootstocks. The LEfSe detected 26 differentially abundant taxa in the rhizospheres ( Figure 5 ), which discriminated the bacterial communities between the different root genotypes.

Actinomycetota was the major phyla that contributed to differentiate the rhizosphere microbial communities associated with the different rootstock types ( Figure 5 ). In particular, at the genus level, g:Nocardioides contributed to the relative high dissimilarity between 5BB and 1103P (enriched in the latter) with a LDA score above 4.5. Among Actinomycetota, the genus g:Mycobacterium, g:Agromyces, g:Streptomyces, g:Modestobacter, and g:Arthrobacter also determined the dissimilarities among rootstocks in the rhizosphere fractions, being enriched in 1103P (with a LDA score ranging from 3 to 3.6). We also found different microbial biomarkers belonging to Alphaproteobacteria class, such as the genus g:Sphingomonas, g:Hyphomicrobium, and g:Lacibacterium. Finally, among Gammaproteobacteria, the 1103P rhizosphere microbiome revealed that the genus g:Lysobacter was also enriched with respect to 5BB.

prediction of functional content of microbiome from the rhizosphere of both rootstocks was obtained by processing the OTU table on PICRUSt using the KEGG database. We found 2,179 molecular functions in terms of functional KEGG orthologs (KO) in our samples. However, we focused on the pathways involved in carbon, nitrogen, and phosphorus cycling and PGP activities ( Figure 6 ) and found 40 KO to be involved. A total of 15 KO were associated to C cycle, 12 to N cycle, three to P cycle, and four to S cycle. In addition, we also found five KO to be associated to PGP activities. Comparing the soil microbiomes associated with the two rootstocks by Mann–Whitney U-test (p < 0.1), we found a significant effect of plant genotype.

Most key enzyme-encoding genes involved in the cycling of N were more abundant in the 5BB rootstock rhizosphere soil (Mann–Whitney test, P < 0.1; Figure 6 ). Specifically, these were involved in dissimilatory nitrate reduction (KO00362) and Anammox (KO10535) pathways and were 1.7- and 7. 8-fold higher in the functional profile associated to 5BB.

As regards the C cycle, chitinase genes (involved in the chitin breakdown pathway) showed a 1.4-fold higher abundance in the 1103P-associated rhizobiome compared to 5BB. In addition, 1103P also had a 2.9-fold higher abundance of carbon monoxide dehydrogenase (K00192). Interestingly, we found that the 1103P rhizosphere was also significantly enriched in microbes with a PGP functionality compared to 5BB. Specifically, for potential PGP traits involved in plant drought tolerance (K01575) (Vigani et al., 2019), a 3.8-fold higher abundance in 1103P was observed. Genes that encode enzymes linked to either nitrilase route in indole acetic acid (IAA) biosynthesis (K01501) and 1-aminocyclopropane-1-carboxylate (ACC) deaminase activity (K01505) were enriched in 1103P (2.3- and 2. 2-fold, respectively).

However, considering the PGP key enzyme-encoding genes for biostimulation, IAA production (K01426), and auxin production (K00817), no rootstock effect was identified. In addition, no differences of functionality related to the S cycle were observed.