Samples were collected in 2016, one week before the harvest period, from ungrafted grapevines of the MA and CS cultivars. Fourteen vineyards were sampled across four wine-producing provinces of Argentina: the South region (RN, Patagonia), the Centre-West region (SJ and MZA), and the North-West region (SAL).

In SJ, part of the Centre-West wine region, samples were obtained from two vineyards located in the Ullum Valley, 6 km apart. In MZA (also in the Centre-West wine region), samples were taken from vineyards in Agrelo (Luján de Cuyo), cultivated using a vertical shoot positioning system on loamy-clay soils. Both SJ and MZA have arid climates. In RN, located in the South wine region, samples were collected from Mainque (RN1, only MA samples) in the Upper Valley of Río Negro and Viedma (RN2, samples of both CS and MA) near the Atlantic Ocean. Both zones have semi-arid climates. Annual rainfall in RN1 ranges from 120 to 180 mm, with temperatures between 30 and 34 °C in summer and 10 to 14 °C in winter. RN2’s climate is influenced by its proximity to the Atlantic Ocean, where summer temperatures can reach 30 °C, and winter temperatures range from 2 to 12 °C. Samples were collected from two plots with loamy-sandy soils, situated 7 metres above sea level (masl).

In Salta Province, part of the North-West wine region, the climate is warm and dry. Samples from MA and CS vineyards were taken in the Calchaquí Valleys (southwest of the province) at three sites: Molinos (SAL1) and Cafayate (SAL2), while only MA samples were collected in Cachi (SAL3). The altitude of these sites ranged from 1700 masl at Molinos to 2500 masl at Cachi.

The study included nine vines sampled per plot within a 49 m ² area, with sampling points located at least 7 m from the plot edge. The sampling design comprised nine vines spaced 2.5 m apart within a 14 m ² quadrat. Samples were pooled to form three composite biological replicates. Samples were collected in sterile containers, transported on ice, and stored at −20 °C until analysis. A schematic of sample collection can be found in Figure S1 of Oyuela Aguilar et al. (2020). Further details regarding the locations, geographical characteristics, grapevine age, management type (organic or conventional), and inter-row management (bare ground or vegetation cover) can be found in our previous research (Oyuela Aguilar et al., 2021; Toscani et al., 2025).

Soil and rhizosphere samples were collected from each vineyard at a depth of 30 cm and 20–30 cm from the vine trunks. Sampling involved nine vines per plot, covering an area of 49 m ², with each plot located at least 7 m from the vineyard edge to avoid border effects. The first layer of surface soil was discarded prior to collecting root-associated samples. Soil was gently detached from the roots and collected using a sterile metal spoon; adherent soil was further dislodged by carefully rubbing the root surface with sterile scalpels. Samples were subsequently sieved (0.5 mm mesh) to remove residual roots and plant debris. DNA was extracted from 0.4 g of rhizosphere soil using the FastDNA Spin Kit for Soil (MP Biomedicals, LLC, Solon, OH, USA) according to the manufacturer’s instructions.

Leaves were washed with saline solution containing 0.01 % Tween 80 in 50 mL polypropylene tubes (2 or 3 samples). 3 or 4 leaves of each sample were maintained in a horizontal shaker for 1 h at 100 rpm. The resulting solution was centrifuged at 4000 × g, and the pellet was transferred to a 2 mL microtube for further analysis.

Grapes from different vineyards were hand-harvested and crushed under aseptic conditions in the laboratory. Only healthy grapes were used. A fraction of must derived from approximately 100 g of grapes was centrifuged at 16,000 rpm, and the pellet was washed twice with a 0.9 % saline solution. The pellet was transferred to a 2 mL microtube and stored.

Samples were stored at –20 °C until analysis. Nine samples from each location and cultivar were pooled to yield three biological replicates.

Soil samples were carefully collected from the root systems of grapevines in the vineyards. Adherent soil was gently removed into a sterile bag, using a scalpel when necessary. Detailed information on soil physicochemical characterisation is found in our previous publication (Toscani et al., 2025). In summary, SAL soils were distinguished from the others by their high sand content. SJ, MZA, and RN1 exhibited similar properties, whereas RN2 differed markedly from the other three regions, likely due to its proximity to the Atlantic Ocean coast (Toscani et al., 2025).

DNA extractions were conducted using the FastDNA Spin Kit for Soil (MP Biomedicals, LLC, Solon, OH, USA), following the manufacturer’s instructions. Extracted DNA was quantified with a Qubit 2.0 Fluorometer (Thermo Fisher Scientific), and its purity was assessed by 260/280 nm and 260/230 nm absorbance ratios using a NanoDrop spectrophotometer.

The diversity of the fungal community in each sample was analysed by amplification of the eukaryotic ITS1 (Internal Transcribed Spacer) region, located between the 18S rRNA and 5.8S rRNA ribosomal genes, using the primers ITS1F (ACTTGGTCATTTAGAGGAAGTAA) and ITS2 (BGCTGCGTTCTTCATCGATGC) (White et al., 1990). A two-step PCR approach was employed for Illumina sequencing library preparation, as described by Gobbi et al. (2020). Sequencing was carried out on the Illumina MiSeq platform using the V2 500-cycle reagent kit.

The double PCR amplification was performed as follows. In the initial PCR, a final reaction volume of 25 μL included 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. The thermal cycling conditions were as follows: an initial denaturation at 95 °C for 2 minutes, followed by 33 cycles of 95 °C for 15 seconds, 55 °C for 15 seconds, and 68 °C for 40 seconds, with a final extension step at 68 °C for 4 minutes. The total DNA concentration was quantified using a Qubit 2.0 Fluorometer (Thermo Scientific). The second PCR added indices to the amplicons generated in the first PCR. This reaction, with a final volume of 28 μL, contained 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 the initial PCR amplicons. Both PCR steps were verified through 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 were uploaded to the SRA under the BioProject accession number PRJNA1282061.

Data analysis was performed using the 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 stringent quality and annotation thresholds to ensure high-fidelity analyses. During read processing, a Phred quality score cutoff of 20 was applied to trim low-quality ends, with a requirement that at least 80 % of nucleotides per read were correctly called. A minimum read length of 50 nucleotides was enforced. For OTU processing, reads were dereplicated using a prefix setting, with no upper limit on OTU length (maximum length set to 0) and a minimum OTU length of 50 nucleotides. Only sequences with a minimum abundance of 4 were retained during dereplication. Clustering was performed on both strands with a similarity threshold of 0.97.

OTU annotation was performed using the UNITE+INSD database (Abarenkov et al., 2024). The taxonomic assignment of ITS1 sequences was performed using BLAST against the UNITE database. First, a local BLAST database was created from the UNITE public release (version 21.04.2024) in FASTA format using the makeblastdb command with the nucleotide database type option. The database files were generated and stored in a specified directory. Then, a custom Python script was developed to automate the BLAST search and extract the best hits for each query sequence (Supplementary script). The script utilised the NcbiblastnCommandline function from the Biopython library to execute BLAST searches in tabular format (outfmt=6) with an e-value threshold of 1e-5. The output file containing the BLAST results was subsequently parsed to identify the best match for each query based on the lowest e-value. Only OTUs belonging to the Fungi and Stramenopila kingdoms (Phylum Oomycota) were used for subsequent analyses.

A weighted non-null normalisation by site and a filtering step, requiring a minimum of eight samples and a total abundance threshold of 2.4 (in logarithmic scale), were performed prior to analysis. Richness, Shannon and Inverse Simpson indices, principal coordinate analysis (PCoA), rarefaction plots, scatterplots, family barplots, and heatmaps were generated and analysed using the online Shaman platform.

Canonical Correspondence Analysis (CCA) was conducted using the PAST software package. (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 employed to identify overrepresented OTUs in the analysed samples. An alpha value of 0.05 was applied for the factorial Kruskal-Wallis test among classes and the pairwise Wilcoxon test between subclasses. Taxonomic units with a logarithmic LDA score greater than 4 were considered overrepresented and visualised in the resulting figure. LEfSe analyses were conducted using the Galaxy Metabiome environment (http://mbac.gmu.edu/mbac_wp).

To assign ecological functions to the fungal OTUs identified, we employed the tool FungalTraits (Tanunchai et al., 2023). This database was used to assign lifestyle-related functional traits at the genus and species hypothesis level. This resource integrates expert-curated information from FUNGuild and FunFun (Nguyen et al., 2016; Krivonos et al., 2023), and allows for trait-based assignments across a wide range of fungal taxa.

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

Fungal 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), run on the Galaxy platform (v1.0.0). Network graphs were generated in Cytoscape v3.10.3 (Shannon et al., 2003). SparCC was utilised to infer microbial associations from compositional data. A SparCC correlation matrix was computed using FastSpar (Watts et al., 2019), averaging 20 inference iterations. An initial filter was applied, excluding OTUs present in less than 50 % of the samples. To stabilise estimates, we used 10 exclusion iterations with a correlation-strength exclusion threshold of 0.7. Two-sided SparCC pseudo p-values were estimated with 100 permutations. Edges were retained when |r| ≥ 0.6 and p ≤ 0.05. Networks were generated with the “Generate networks from SparCC” Galaxy tool (version 1.0.0) as a single-mode adjacency network and exported to Cytoscape. Network modularisation used the greedy modularity optimisation algorithm, and module roles were classified using Z–P scores.

The integrated soil prokaryotic and fungal network was generated by normalising the respective OTU count tables using Scaling with Ranked Subsampling (SRS) prior to merging the tables ( Beule & Karlovsky, 2020). Prokaryotic OTUs were obtained by 16S rRNA gene amplicon sequencing and were reported in our previous publication (Toscani et al., 2025). The SRS method for normalising species count data minimises subsampling error while preserving the original community structure. The SparCC algorithm was employed to infer microbial associations from compositional data. Pseudo p-values were estimated using 100 permutations, and a correlation strength threshold of 0.6 was applied.

Following ITS sequencing, 9,857,782 amplicon sequences were obtained, which were reduced to 1,398,814 after dereplication (Table S1). After removing singletons and chimeras, 98,259 sequences remained for classification into Operational Taxonomic Units (OTUs), yielding 4669 OTUs (Table S2). Since ITS metabarcoding can detect a wide range of eukaryotes, we excluded OTUs assigned to taxa outside the phylum Fungi and the phylum Stramenopila, order Oomycota, leaving 4422 OTUs for further analysis (Table S3). Extended details on the sequenced samples are provided in Table S1. Although three replicates were collected per vine, 11 samples could not be successfully sequenced (Table 1). Furthermore, as the rarefaction curve did not reach saturation, six of the sequenced samples were excluded from further analysis (Figure S1), resulting in 39 rhizosphere, 34 leaf, and 38 must sequences. Samples were normalised per site prior to statistical analysis.

To gain insight into organ-specific diversity, three estimators were used: OTU richness, Shannon diversity index, and the inverse Simpson index. Since no significant differences were observed between the MA and CS indices (Figure S2–4), we opted to pool samples by site for subsequent analyses.

To assess differences in fungal composition among plant organs, diversity indices were plotted accordingly (Figure 1A–C). A Mann–Whitney test with Bonferroni adjustment identified significant differences and revealed that the rhizosphere was the most diverse sample. Must samples exhibited the lowest diversity, whereas leaf samples showed intermediate values. The rhizosphere was the most diverse compartment, leaf samples showed intermediate diversity, and must samples showed the lowest values. PCoA analysis showed that Axis 1 significantly separated rhizosphere samples, whereas must and leaf samples largely overlapped (Figure 1D). Greater fungal diversity in the soil, both bulk and that associated with the rhizosphere, has been previously reported by other authors (Swift et al., 2021). Environmental factors, such as solar UV radiation, to which aerial organs are more exposed, are key drivers of shifts in fungal communities worldwide, exerting their greatest influence on the presence of free-living fungi (Egidi et al., 2023). The selective pressure imposed by these abiotic factors may account for the lower diversity observed in more exposed plant organs, as only fungi tolerant of these conditions can proliferate in these tissues. The richness observed in the rhizosphere may be due to routine application of fungicides, which imposes stress on fungal growth on leaves and fruits, and to the fact that the rhizosphere is a richer, physicochemically less hostile environment than exposed organs. Additionally, rhizosphere microbes are fuelled by sustained inputs of labile carbon released by grapevine roots, as well as by carbon derived from the rapid turnover and decomposition of fine roots (Darriaut et al., 2024). In parallel, arbuscular mycorrhizal fungi channel substantial host-derived carbon into the soil through their hyphal networks, which also turn over rapidly, thereby supporting a diverse microbial food web and contributing to soil organic matter formation (Kakouridis et al., 2024). Collectively, these factors may contribute to the greater diversity that characterises the rhizosphere.

Violin plots illustrating: A. Richness, B. Shannon, and C. Inverse Simpson indices for fungal diversity. Statistical analysis was performed using the Kruskal–Wallis test, with pairwise comparisons identified through a post hoc Mann–Whitney test adjusted using the Bonferroni correction. D. PCoA summarising the differences in ITS rRNA gene sequences across various organs. The first two axes are displayed. A PERMANOVA test was used to assess variance based on distance matrices (p-value = 0.001). E. Venn diagram depicting fungal OTUs identified in Root (2716 total), Leaf (1405 total), and Must (1534 total) samples. F. Distribution of the fungal families identified in each organ. Families representing more than 1 % are shown individually. G. Cladogram plotted from LEfSe analysis showing the phylogenetic relationships of taxa significantly enriched across rhizosphere, leaf, and must samples. Each circle represents a taxonomic level, and colours indicate the group in which the taxon is significantly overrepresented. The coloured fan shapes represent all possible taxa to which the OTU.

In contrast to the pattern observed for rhizospheric prokaryotic diversity, where RN1 was the least diverse environment (Toscani et al., 2025), the analysis of rhizospheric fungal diversity revealed that SAL2 and SAL3 exhibited the lowest richness values (Figure S5A–C). Notably, only SJ1 showed significantly higher diversity, as indicated by both the Shannon and inverse Simpson indices. Furthermore, PCoA revealed that the fungal community at SAL3 was clearly distinct from those at the other sampling sites (Figure S5D). These findings suggest that fungal and prokaryotic communities in vineyard rhizospheres may be shaped by different environmental or edaphic factors, resulting in contrasting diversity patterns and site-specific community structures.

To further investigate fungal diversity in the rhizosphere of the sampled vines, OTUs were taxonomically annotated and integrated with physicochemical data and previously reported information on rhizospheric prokaryotic diversity (Toscani et al., 2025) (Table S5). Lophiostomataceae, Nectriaceae, Pleosporaceae, and Cladosporiaceae were the most frequently detected families, accounting for approximately 50 % of the annotated OTUs. Lophiostomataceae and Nectriaceae predominated in SJ1, SJ2, MZA, RN1, and RN2, whereas Saccotheciaceae and Dothioraceae were enriched in SAL (Figure 2A). Pleosporaceae, Cladosporiaceae, and Metschnikowiaceae were dominant in SAL2 and SAL3, whereas Lophiostomataceae was not detected. In contrast, SAL1 showed higher relative abundances of Lophiostomataceae and Cladosporiaceae. These families – Saccotheciaceae, Dothioraceae, Pleosporaceae, Cladosporiaceae, and Metschnikowiaceae – were associated with sandy soils, higher humidity, and increased precipitation, as observed in SAL (Figure 2C), supporting the fungal taxonomic distribution patterns reported herein.

To refine the characterisation of such fungal communities, a LEfSe analysis was conducted (Figure 2B). This analysis identified differential OTUs mainly associated with saprotrophic or pathogenic lifestyles (Table S6), including members of the Nectriaceae family in MZA and SJ1, where an OTU affiliated with the genus Fusarium, a genus of filamentous fungi that contains many agronomically important plant pathogens, mycotoxin producers, and opportunistic human pathogens (Ma et al., 2013), was overrepresented.

To gain insight into fungal interactions in the grapevine rhizosphere, we constructed a SparCC network that comprised 11 modules, with Modules 1 and 2 containing the majority of OTUs (Figure 3A). The topology is dominated by a densely connected core (Module 1, red), which concentrates most edges and contains several common vineyard families ( e.g., Didymellaceae, Pleosporaceae, Cladosporiaceae, and Nectriaceae). Surrounding this core, multiple satellite modules (notably the cyan, orange, and purple clusters) form smaller subgraphs with comparatively fewer between-module links. Positive correlations predominate within modules, whereas negative edges are rarer and tend to bridge modules, indicating niche partitioning among groups that otherwise exhibit strong intramodular co-selection.

A. Distribution of the fungal families identified in rhizosphere samples. Families representing more than 1 % are shown individually. B. Taxonomic distribution of differential OTUs determined by LEfSe analysis. C. CCA plot of the most abundant fungal families and their relationship with the physicochemical features of soil. Perm (n = 999) Trace = 1.745, p = 0.012.

When analysing the diversity of fungi present in grapevine leaves, it was observed that samples from RN1 exhibited high values for richness, Shannon index, and inverse Simpson index, whereas SAL1 and SAL2 showed the lowest values (Figure 4A–C). Additionally, PCoA revealed differences between Cuyo (SJ and MZA) and Patagonian (RN1) grapevine leaf samples (Figure S6A). Due to the low number of sequenced samples (Table 1), SAL3 was not included in the analysis.

Violin plots depict: A. Species richness, B. Shannon diversity, and C. Inverse Simpson indices for fungal communities associated with Leaf samples. Statistical comparisons were conducted using the Kruskal–Wallis test, followed by pairwise Mann–Whitney tests with Bonferroni correction. D. Distribution of the fungal families identified in Leaf samples. Families representing more than 1 % are shown individually. E. CCA analysis of fungal families and their relationship with the climatological features. Perm (n = 999) Trace = 0.7496, p = 0.001.

Similarly to leaf samples, must from RN1 showed the highest values for species richness, Shannon, and inverse Simpson index, indicating a more diverse fungal community (Figure 6A–C). This consistency in the patterns observed across both aerial plant samples may be related to the fact that RN1 is the only organic vineyard sampled. This distinctive production approach, which relies on alternative pest control methods, may account for the higher fungal diversity observed in RN1 compared to the other vineyards analysed. Furthermore, PCoA showed differences between must samples from SJ and RN (Figure S6B). In contrast to the pattern observed for leaves, must from MZA differed from that of SJ and clustered more closely with the RN samples.

A. Species richness, B. Shannon diversity, and C. Inverse Simpson indices for fungal communities associated with must samples. Statistical comparisons were performed using the Kruskal–Wallis test, followed by pairwise Mann–Whitney tests with Bonferroni correction. D. Taxonomic distribution of fungal families detected in must samples. Families contributing more than 1 % to the overall community are shown individually. E. Taxonomic distribution of differential OTUs determined by LEfSe.