Soil samples were collected from ungrafted grapevines of MA and CS cultivars one week prior to the harvest period in 2016. We sampled fourteen vineyards in four wine-producing regions of Argentina: the southern region (Río Negro Province, Patagonia), the western central region (San Juan and Mendoza Provinces), and the northwestern region (Salta Province) (Figure 1A). Information regarding the sampled sites is summarised in Table 1.

In the San Juan (SJ) Province (western central wine region), the samples came from two vineyards located in the Ullum Valley, 6 km from each other: Finca Norte (SJ1) and Finca Arriba (SJ2). In the Mendoza (MZA) Province (also western Central wine region), the samples came from vineyards with vines planted based on a vertical shoot positioning system in Agrelo (Luján de Cuyo); all the vineyards were on loamy-clay soils, and the climate of both San Juan and Mendoza is arid.

^{A) Map of Argentina indicating sampled sites. This map was created using OpenStreetMap (openstreetmap.org/copyright). The data is available under the Open Database License. B) PCA analysis of soil physicochemical characteristics. Samples from Río Negro (Green), San Juan (Blue), Mendoza (Purple) and Salta (Red, Orange, and Yellow). Two first dimensions are shown. Dots and squares represent Malbec and Cabernet-Sauvignon cultivars respectively.}

Soil samples were obtained by carefully removing soil from the roots sampled in the vineyards. It was rubbed off gently into a bag, using a scalpel when needed. Rhizosphere soil samples were sieved (pore size 0.5 mm) to eliminate the remaining plant debris. Soil samples were analysed at the Centro Nacional de Investigaciones Agropecuarias (CNIA-INTA) (Buenos Aires, Argentina). Samples were processed as previously described by Oyuela Aguilar et al. (2020). The pH was determined using the potentiometric method and a 1:2.5 distilled water ratio (ISO 10694-Soil quality determination of pH norms). A strong oxidising microscale mixture was used to assess organic carbon ( Walkley & Black, 1934), and the mass loss calcination method and the Kjeldahl method were used to assess organic material ( Ball, 1964) and organic nitrogen ( Bremner, 1960) respectively. Soil phosphorus content was determined using the Bray & Kurtz method ( Bray & Kurtz, 1945). Soil texture was analysed according to Kilmer et al. (1949). USDA soil texture diagram was used to determine the soil type in each vineyard, acquired lime, sand and clay content .

DNA extractions were performed using the FastDNA Spin Kit for Soil (MP Biomedicals, LLC, Solon, OH, USA) from 0.4 g of rhizosphere soil following the supplier’s instructions. The extracted DNA was quantified using a Qubit1 2.0 Fluorometer (Thermo Scientific). DNA purity at 260/280 nm and 260/230 nm absorbance ratios were assessed using a Nanodrop spectrophotometer.

The bacterial diversity was analysed by amplifying the hypervariable V3-V4 region of the 16S rRNA gene using 341F (CCTACGGGNGGCWGCAG) and 806R (GACTACHVGGGTATCTAATCC) primers. A 16S rRNA gene double PCR-step approach was used for Illumina sequencing library preparations ( Gobbi et al., 2020). Sequencing was performed on Illumina’s MiSeq platform using the V2 500 cycles reagent kit.

The double amplification was performed as follows. For the initial PCR, a final volume of 25 μL contained a mixture of 12 μL of AccuPrime™ SuperMix II (Thermo Scientific™), 0.5 μL of bovine serum albumin (BSA; final concentration, 0.025 mg/mL), 0.5 μL each of forward and reverse primers (10 μM stock), 1.5 μL of sterile water, and 5 μL of template DNA. This reaction mixture underwent pre-incubation at 95 °C for 2 min, followed by 33 cycles of 95 °C for 15 s, 55 °C for 15 s and 68 °C for 40 s, followed by a final extension at 68 °C for 4 min. Total DNA was quantified using a Qubit® 2.0 fluorometer (Thermo Scientific™).

The second PCR aimed to add indices to each amplicon from the first PCR. A final reaction volume of 28 μL was used, containing 12 μL of AccuPrime™ SuperMix II (Thermo Scientific™), 2 μL of primers with index sequences and Illumina P7 (CAAGCAGAAGACGGCATACGAGAT) and P5 (AATGATACGGCGACCACCGA) adapters, 7 μL of sterile water, and 5 μL of amplicons from the initial PCR. Cycling conditions included denaturation at 98 °C for 1 min, followed by 13 cycles of 98 °C for 10 s, 55 °C for 20 s and 68 °C for 40 s, followed by a final extension at 68 °C for 5 min. In both instances, PCR products were verified via 1.5 % agarose gel electrophoresis.

Illumina reads were demultiplexed using bcl2fastq V.2.17.1.14 (Illumina). Adapters were trimmed with Trim Galore v0.4 (https://github.com/FelixKrueger/TrimGalore.git) running cutadapt v1.8.3 ( Martin, 2011) and primer sequences were deleted from the 5’ ends of each read using the custom script (https://github.com/padbr/asat/blob/master/strip_degen_primer.py). Raw sequencing data associated with this work was uploaded to the SRA under the BioProject accession number PRJNA1087980.

Data analysis was performed using Shaman platform (https://shaman.pasteur.fr/), following the pipeline described by Volant and collaborators ( Volant et al., 2020). Read and operational taxonomic unit (OTU) processing parameters were set with strict quality and annotation thresholds to ensure high-fidelity analyses. For read processing, a Phred quality score cutoff of 20 was applied to trim low-quality ends, with a minimum requirement of 80 % correctly called nucleotides per read and a minimum read length of 50 nucleotides. In the OTU processing stage, reads were dereplicated with a prefix setting, no upper limit on OTU length (maximum OTU length set to 0) and a minimum OTU length of 50 nucleotides. Only sequences with a minimum abundance of 4 were retained at dereplication, and clustering was performed on both strands with a 0.97 similarity threshold. For OTU annotation, both strands were used, and identity thresholds were progressively more stringent across taxonomic levels: a minimum of 0.75 for Kingdom, 0.75–0.785 for Phylum, 0.785–0.82 for Class, 0.82–0.865 for Order, 0.865–0.945 for Family, 0.945–0.98 for Genus, and 0.98–100 for Species annotations. These parameters ensured accurate read trimming, clustering and taxonomic assignment throughout the analysis.

Taxonomic annotation of OTUs was performed using SILVA database ( Pruesse et al., 2007) version 138. A weighted non-null normalisation by site and a filtering step for a threshold of a minimal number of 8 samples and 2.4 on the total abundance (in log) was perform prior to data analysis. Richness, Shannon and Inverse Simpson indexes, principal coordinate analysis (PCoA), rarefaction plots, scatterplots, families’ barplot and heatmaps were analysed and plotted using the on-line Shaman platform.

Principal component analysis (PCA), clustering analysis using Euclidean similarity index and Canonical Correspondence Analysis (CCA) were performed using Past software ( Hammer et al., 2001). Data were normalised by Standard score ((X-µ)/σ) before analysis ( Han et al., 2012).

Linear Discriminant Analysis Effect Size (LEfSe) ( Chang et al., 2022) was performed to identify overrepresented OTUs in soil samples. Alpha values for the factorial Kruskal-Walli’s test among classes and for the pairwise Wilcoxon test between subclasses used were 0.05. Taxonomic units above a logarithmic LDA score of 4 were considered as overrepresented and are shown in the graph. LEfSe was run within the Galaxy Metabiome environment (http://mbac.gmu.edu/mbac_wp).

PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) ( Douglas et al., 2020)–run within the Galaxy Metabiome environment–was used to predict metabolic pathways in 16sRNA NGS sequencing data. The presence of differential pathways was inferred using LEfSe.

Microbial networks were constructed using the SparCC algorithm ( Friedman & Alm, 2012) in the Integrated Network Analysis Pipeline 2.0 (iNAP 2.0. https://inap.denglab.org.cn) ( Peng et al., 2024). Network graphs were generated with Cytoscape v3.10.3 ( Shannon et al., 2003). The SparCC algorithm was employed to infer microbial associations from compositional data. Pseudo p-values were calculated through 100 permutations, with 20 inference iterations averaged per calculation. To reduce the influence of strongly correlated pairs, 10 exclusion iterations were applied, using a correlation strength exclusion threshold of 0.7. This approach ensured robust estimation of microbial associations. A SparCC correlation matrix was generated and filtered based on a threshold of 0.6. Correlation coefficients with absolute values below this threshold were excluded. Additionally, the SparCC pseudo p-value matrix, derived from the 100 permutations, was filtered at a p-value threshold of 0.05 to retain statistically significant correlations. Network modularisation was performed using the Greedy modularity optimisation method, as implemented in the Galaxy tool (Version 1.0.0). The network matrix, pre-filtered after applying correlation and p-value cutoffs, was analysed to identify modular structures and key hub nodes. The modularity analysis also incorporated Z-P scores to classify module hubs and assess their significance within the microbial association network.

The figures were created using the Inkscape (Inkscape Project, 2020. https://inkscape.org) and GIMP (The GIMP Development Team, 2019. https://www.gimp.org) software packages.

Soil physicochemical properties were determined for the fourteen sampled sites (Table S1). Like in any edaphic study on soils, and as a starting point for visualising any variation in the soil composition of different sites, we performed a principal component analysis (PCA). This two-dimensional PCA plot (Figure 1B) clearly distinguished the samples, showing a clear difference between locations. Component 1 differentiates SAL soils from the rest of the sampled sites, while Component 2 discriminates RN2. This analysis showed a clear grouping of SAL and RN2, while samples from the western region (SJ and MZA) are dispersed in quadrants 1 and 4. In contrast to RN2, RN1 is positioned between SJ and MZA, showing a clear distinction between the Patagonian samples.

After performing the 16S rRNA sequencing, we obtained a total of 903,462 amplicons, with a remaining 668,603 after dereplication. Followed by the sequence removal of singletons and chimeras, 14,469 sequences remained for classification into OTUs, resulting in a total of 2,945. Further description of the sequenced samples can be found in Table S2. Although three composite samples were collected for each vine, one sample from the SAL3MA set could not be successfully sequenced (Table 1). Additionally, because the rarefaction curve for sample run.MO.A31_S31_L001.fastq - which belonged to one sequence from the SJ2CS - did not reach saturation, it was excluded from the results (Figure S1), leaving a total of 40 sequences for the final analysis. Samples were normalised per site prior to the statistical analysis.

Three estimators of α-diversity were used: species (OTUs) richness, Shannon diversity index and the inverse Simpson index. As there were no observed differences between MA and CS indexes (Figure S2), we decided to pool samples in terms of diversity according to the sites.

To evaluate the effect of location on rhizosphere microbial composition, diversity indexes were plotted accordingly (Figure 2). Although the species richness comparison using the non-parametric Kruskal-Walli’s test showed no significant differences in bacterial communities, pairwise-comparisons through a post hoc Mann-Whitney Bonferroni adjusted test showed differences for RN1 with respect to SAL1 and RN2 (Figure 2A). While global comparison given by the Shannon index of Kruskal-Walli’s test showed no significant differences, pair to pair comparisons showed statistically significant differences between RN1 with respect to SJ1, SJ2, MZA, SAL1 and RN2 (Figure 2B). Similar results were obtained for the inverse Simpson index (Figure 2C).

^{Boxplots (A), (B) and (C) represent the prokaryotic species α-diversity indexes. The statistical} analysis ^{was performed using a Kruskal Walli’s test. The identified pairwise-comparisons results were obtained through a post hoc Mann–Whitney Bonferroni adjusted test.}

As previously mentioned, 2,945 OTUs were obtained after sequencing. All these OTUs were annotated using SILVA ( Pruesse et al., 2007) at the Kingdom level, of which 2,934 (99.63 %) were assigned to the phylum level, 2,864 (97.25 %) to the class level, 2,584 (87.74 %) to the order level, 2,113 (71.54 %) to the family level and 1,273 (43.23 %) to the genus level. Only 514 (17.45 %) of the OTUs could be assigned to the species level.

^{PCoA summarizes the differences between rhizospheric bacterial 16S rRNA gene sequence analysed. The two first axis are shown. Permanova test was performed for the analysis of variance using distance matrices. p-value = 0.001.}

Network analysis can identify disproportionately influential species, or keystones, within communities that contribute to the system's resilience to external disturbances. A SparCC-based network was therefore constructed to explore the interactions between the OTUs in the analysed samples. This network was compared to random networks generated using Greedy modularity optimisation and Short random walks. In the comparison of the empirical network with both randomised networks, the values of the Clustering Coefficient, Average Path Distance, Geodesic Efficiency, Centralisation of Stress Centrality, Modularity, Transitivity, Hierarchy and Lubness validated the network generated with our data (Table S4). The empirical network has a significantly higher clustering coefficient than the random networks and larger average distances between nodes. Nodes in the empirical network show high stress centrality, implying that certain species are crucial for connectivity and may act as “hotspots” within the network. The higher transitivity observed in the empirical network suggests a greater tendency for a node's neighbours to also be connected to each other, which is typical in organised networks and less common in random networks.

^{A. Colours} show ^{modules that comprise the OTUs’ network. B. Colour code indicates connector, module, network and peripheral hubs}.

The results of various metagenomics studies have supported the connection between the grape microbiome and factors such as vineyard location, climatic conditions, and other vineyard-related variables ( Mezzasalma et al., 2018; Oyuela Aguilar et al., 2020; Zarraonaindia et al., 2015). To assess the possible existence of patterns linking the physicochemical characteristics of the soil and the families of bacteria found, a canonical correspondence analysis (CCA) was performed (Figure 6). Quadrant 1 is characterised by high sand contents, with a high number of bacteria from families such as Pseudomonadaceae, Bacillaceae and Paenibacillaceae. The opposite quadrant (Quadrant 3) is characterised by soils rich in lime, clay, sodium, potassium and carbonate and in bacteria such as Rhizobiaceae, Sphingomonadaceae and Xanthomonadaceae, among others.

This therefore indicates a clear positive correlation between the families of Gram-positive bacteria Bacillaceae and Paenibacillaceae and soils with high sand content (Figure S6A and S6B). Consistent with the results of the CCA, these Gram-positive bacterial families negatively correlated with Xanthomonadaceae (Figure S6C and S6D) and Steroidobacteraceae (Figure S6E and S6F). Thus, the results regarding the Xanthomonadaceae and Steroidobacteraceae families are related to the amount of clay in the soil (Figure S6G and S6H), reinforcing the hypothesis that a soil’s physicochemical characteristics condition its microbiota.

^{CCA analysis of the most abundant bacterial families and their relationship with the physicochemical features of soil. Only families representing more than 0.5% of total counts are shown. Dots and squares represent Malbec and Cabernet-Sauvignon cultivars respectively. Perm (n = 999) Trace = 0.6385, p = 0.001.}

This study provides valuable insights into the prokaryotic diversity of Argentine vineyard soils, with a particular focus on the rhizosphere microbiota of two grapevine cultivars, MA and CS, across four key wine-producing provinces: Mendoza, San Juan, Río Negro, and Salta. The physicochemical characteristics of the soils, including the sand-to-clay ratio and concentrations of sodium, potassium and carbonates, play a pivotal role in defining regional groups and shaping microbial diversity. Taxonomic assignment and LEfSe analysis identified a unique taxonomic signature for each site, highlighting the prevalence of specific microbial families associated with particular regions. These patterns likely reflect variations in environmental factors, soil properties and vineyard management practices. While the microbial communities associated with the different grapevine regions displayed distinct compositions, a subset of beneficial bacteria–including Bacillus, Paenibacillus, Pseudomonas, and Rhizobium–was consistently identified across all studied locations. These bacteria, recognised as plant growth-promoting rhizobacteria (PGPR), may enhance nutrient availability, confer pathogen resistance and promote overall vine health. The findings underscore the potential of microbial communities to support grapevine health and contribute to wine quality. This comprehensive analysis sheds light on the intricate relationship between soil characteristics, microbial diversity and the concept of terroir in Argentine vineyards. Understanding these interactions is essential for optimising vineyard management practices and ultimately influencing the organoleptic properties of wines produced in different regions.

We would like to thank the wineries in Argentina for allowing us to collect samples from their vineyards. We would also like to thank the “MicroWine” team for their support and input. MO carried out her Marie-Curie ITN fellowship at the IBBM- (Instituto de Biotecnología y Biología Molecular) CONICET-UNLP. AMT is a researcher at UNLP; CR and RRW are fellows of the Research Career of CONICET; MFDP and MP are members of the Research Career of CONICET; and LS is a Head Professor of the Universidad Nacional de Quilmes and a member of the Research Career of CIC-PBA.

This study was funded by the European Commission Horizon 2020 programme in the Marie Skłodowska-Curie Innovative-Training-Network “MicroWine” (grant number 643063) and the Agencia Nacional de Promoción Científica y Tecnológica PICT-2017-2833.