Ten blueberry cultivars were selected from the blueberry germplasm collection of the Edmund Mach Foundation (FEM) (46.0744° N, 11.2334° E), Trentino-Alto Adige, Italy, and harvested during the 2023 growing season. To address the specific agronomic requirements of blueberry cultivation, a soil management strategy was applied to maintain the acidic pH necessary for optimal growth. For the cultivation, pellet elemental sulfur was broadcast at a rate of 100 kg/ha in November. Subsequently, a fertilization program was implemented from vegetative resumption until the beginning of harvest. Fertilization consisted of a broadcast application of a nutrient blend in the first week of each month from March to June, with the March application occasionally omitted depending on the phenological stage of the plants. The nutrient mixture was applied at a rate of 10 kg/ha (equivalent to 24 kg ammonium sulfate, 20 kg potassium sulfate, 20 kg magnesium sulfate, and 6 kg triple superphosphate per 1,000 m²), ensuring a balanced supply of macronutrients throughout key growth and fruiting phases.

For microbiome analysis, each cultivar was sampled in triplicate, collecting material from three individual plants. To control environmental variation, all samples were collected under stable weather conditions, even though sampling occurred on different days. Each cultivar was sampled in triplicate, with three individual plants per replicate. Statistical analyses accounted for this structure by including replicate as a random factor where appropriate. Four compartments were considered: leaves (30 samples), fruit surface (30 samples), rhizosphere (30 samples), and bulk soil (10 samples), resulting in a total of 100 samples. Fruits were harvested at stage 89 of the Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie Scale. Bulk soil was collected in triplicate for each plant at approximately 1 m from the base, after removing the topsoil, and pooled to form a composite sample. Rhizosphere samples were obtained by carefully exposing and cutting roots with sterile pruning shears, minimizing damage to the plants. All collections (leaves, fruits, rhizosphere, bulk soil) were performed using sterile tweezers, trowels, and shears, and samples were immediately placed in sterile plastic bags, stored on ice in the field, and subsequently kept at −20 °C for temporary storage. Samples were kept at −20 °C for up to 2 weeks before DNA extraction. Transport from the field to the laboratory was performed on ice to maintain low temperature and minimize microbial changes.

Samples were treated under sterile conditions. For the leaf and fruit samples, 150 mL of sterile washing solution (0.9% NaCl, 0.01% Tween 80) was added to the sample; then it was incubated on an orbital shaker at 400 rpm for 1 h (Behrens and Fischer, 2022). At the end of the incubation, the water was filtered using sterile filtration apparatus (2GPU02RE Stericup Quick Release Millipore Express PLUS 0.22 μm PES, 250 mL Merck Millipore) to recover as much biological material as possible. DNA was extracted from the filters using the DNeasy PowerWater kit QIAGEN according to the manufacturer’s protocol. A negative control was also analyzed and subjected to the same polymerase chain reaction (PCR) process as the other samples. No DNA amplification was observed in the negative control, indicating the absence of contamination during the extraction process.

Using sterile forceps and spatulas, we carefully scraped off adhering soil from the roots. 250 mg of rhizosphere and bulk soil were then used for DNA extraction using the DNeasy PowerSoil Pro Kit (QIAGEN) following the manufacturer’s protocol. A negative control was included and subjected to the same PCR amplification process as the samples. No DNA amplification was observed in the negative control, confirming the absence of contamination during the extraction process.

For fungi, the ITS1 region was amplified using ITS1F (5′-CTT GGT CAT TTA GAG GAA GTAA-3′) and ITS2 (5′-GCT GCG TTC TTC ATC GAT GC-3′) primers (White et al., 1990; Gardes and Bruns, 1993) with overhang Illumina adapters. The PCR were carried out with a total volume of 25 μL, containing 1 μL of each primer (10 μM of each one), 0.25 μL of FastStart High Fidelity PCR System (Roche), 18.75 μL of nuclease-free water (Sigma–Aldrich), 2.5 [10X] Buffer, 0.5 μL dNTP and 1 μL of DNA (5–10 ng/μL). PCR were performed by using the GeneAmp PCR System 9,700 (Thermo Fisher Scientific); the conditions were: initial denaturation at 95 °C for 3 min, 35 cycles of denaturation at 95 °C for 20 s, annealing at 50 °C for 45 s, extension at 72 °C for 90 s, followed by a final extension at 72 °C for 10 min.

For bacteria, the sub-region V3–V4 of the 16S ribosomal RNA gene was amplified using the 341F (5′-CCT ACG GGN GGC WGC AG-3′) and 805R (5′-GAC TAC NVG GGT WTC TAA TCC-3′) primers (Klindworth et al., 2013) with overhang Illumina adapters. The PCR were carried out with a total volume of 25 μL, containing 5 μL of each primer (1 μM), 12.5 μL MIX KAPA HiFi HotStart ReadyMix (Roche) and 2.5 μL of DNA (5–10 ng/μL).

PCR were performed by using the GeneAmp PCR System 9700 (Thermo Fisher Scientific); the conditions were: initial denaturation at 95 °C for 5 min, 35 cycles of denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, extension at 72 °C for 30 s, followed by a final extension at 72 °C for 5 min.

The obtained amplicons were purified using the CleanNGS (CleanNA), following the manufacturer’s instructions. Afterward, a second PCR was used to apply dual indices and Illumina sequencing adapters Nextera XT Index Primer (Illumina), by 7 cycles PCR (16S Sequencing Library Preparation, Illumina). The amplicon libraries were purified using CleanNGS (CleanNA), and the quality control was performed on a Tapestation 4150 platform (Agilent Technologies, Santa Clara, CA, United States). Finally, all barcoded libraries were pooled in an equimolar way and sequenced on an Illumina^® MiSeq (PE300) platform in pair ends (2 × 300 bp) (MiSeq Control Software 2.5.0.5 and Real-Time Analysis software 1.18.54.0). Library construction and Next-Generation Sequencing were performed at the sequencing platform at FEM. Negative controls were also sequenced to further exclude potential contaminants.

Physicochemical characterization of soil samples included texture, pH, salinity, contents of total NH₄⁺ and phosphorus (P). These characterizations were performed following standard methods (Siles et al., 2023).

The demultiplexed raw reads were analyzed using the MICCA (MICrobial Community Analysis) v1.7.2 bioinformatics pipeline (Albanese et al., 2015). Sequence pairs were merged to obtain consensus sequences. Primers were trimmed and sequences that did not contain these primers were discarded. Quality filtering was performed using an expected error rate of 0.75 for ITS and 16S, with the minimum read length set to 200 for ITS and 400 for 16S. Reads were clustered into amplicon sequence variants (ASVs) using the UNOISE protocol (Edgar, 2016). A total of 38,259 ASVs for bacteria and 4,658 ASVs for fungi were found. Finally, taxonomy was assigned using AMPtk (Amplicon ToolKit) v1.5.1 pipeline (Palmer et al., 2018) with the database UNITE v10.3 (Abarenkov et al., 2024) for the fungal component and Rdp_16s_v16 (Edgar, 2018; Wang and Cole, 2024) for the bacterial component. To ensure the absence of contaminants after sequencing, the R package “Decontam” version 1.26.0 (Davis et al., 2018) was used, with negative controls as reference. The analysis confirmed the absence of contamination in the dataset. To reduce the influence of extremely rare taxa and potential sequencing artifacts, we applied a prevalence filter prior to all compartment overlap analyses. Specifically, only ASVs present in at least 5% of samples within a given compartment were retained. This step ensured that downstream comparisons reflected biologically meaningful microbial taxa rather than low-abundance noise.

The ITS dataset included 4,289 unique ASVs, with a total of 605,540 raw, non-rarefied reads, which decreased to 501,813 reads after rarefaction. The 16S dataset comprised 38,152 unique ASVs, with 4,959,760 raw, non-rarefied reads, reduced to 3,204,000 reads after rarefaction. The number of raw reads for each sample is reported in Supplementary Table 1. Rarefaction was applied to standardize sequencing depth across samples, which enabled reliable comparisons of microbial and fungal community diversity. The rarefaction depth was chosen based on a trade-off between retaining samples and capturing true diversity, which balanced coverage and diversity loss. Samples from all compartments were rarefied to 9,000 reads for bacteria and 10,000 for fungi. Two alpha diversity indices were computed: Chao1 and Shannon diversity (Kim et al., 2017; Willis, 2019; Finn, 2024), using the “estimate_richness” function from the R package Phyloseq v1.46.0 (McMurdie and Holmes, 2013). Chao1 and the Shannon indices were used to estimate richness and diversity, respectively. The Wilcoxon-Mann–Whitney test was used to assess significance of difference in alpha diversity among groups with Holm–Bonferroni correction for multiple testing. For beta diversity analysis, the Bray–Curtis dissimilarity matrix was used to assess the difference between the plant tissues, which were visualized with a principal coordinate analysis (PCoA). Statistical significance was assessed using permutational multivariate analysis of variance (PERMANOVA) with pairwiseAdonis v0.4 R package (Martinez-Villegas et al., 2024). The parameters considered in this analysis were difference in tissues, and difference in plant cultivar compared to the reference set of pairwise distances between samples from the same plant cultivar and tissues.

For differential abundance analysis and Spearman’s correlations between alpha diversity and edaphic factors, the Microeco v1.13.0 R package (Liu et al., 2021) was used, applying Random Forest classification along with a zero-inflated log-normal (ZILN) model-based differential test, particularly suitable for identifying differentially abundant taxa in sparse datasets (Paulson et al., 2013). Tissue overlap was calculated by converting the ASVs table into a presence/absence matrix, where ASVs were considered present if their abundance was greater than zero. ASVs were defined as present in a tissue if they were detected in at least one sample and present in at least 9 out of 10 cultivars. For visualization, the UpSetR v1.4.0 package was used (Conway et al., 2017). Distance-based redundancy analysis (dbRDA) was performed with Microeco v1.13.0 R package; this constrained ordination method related community dissimilarity matrices to environmental variables, enabled the identification of factors shaping community composition (Yao et al., 2023), moreover we employed envfit correlation from “vegan” to fit the environmental variables onto the PcOA ordination (Kuusisto et al., 2025). The environmental variables used in envfit were NH₄⁺, P, pH, salinity and texture.

We first explored microbial alpha and beta diversity to assess how communities differ in richness across plant compartments. After rarefaction, we computed the Chao1 estimator of species richness and the Shannon entropy index to characterize microbial richness and diversity in each compartment. Richness was highest in bulk soil and the rhizosphere both for fungi (Figure 1A) and bacteria (Figure 1D), with a marked decline on the leaf and fruit surfaces. Shannon diversity index, both for fungi (Figure 1B) and bacteria (Figure 1E), shows that bulk soil and rhizosphere exhibit significantly higher diversity compared to the leaf and fruit epiphytes (Supplementary Table 2).

To explore how soil properties influence belowground microbial communities, we measured soil physico-chemical parameters (Supplementary Figure 1 and Supplementary Table 3). Most variables differed significantly except texture. Pairwise correlations revealed a negative correlation between pH and both salinity and ammonium, and a positive correlation between texture and phosphorus (Supplementary Figure 3). The influence of environmental factors on alpha diversity was assessed using Spearman’s correlation analysis (Figures 1C–F). The only significant correlation was observed in the bacterial community of the rhizosphere, where Chao1 was positively correlated with pH (ρ = 0.517, p = 0.034), whereas no significant correlations were found for Shannon. In bulk soil no statistically significant correlations were detected after p-value adjustment, although some weak to moderate non-significant trends were present.

PCoA of beta diversity revealed two main clusters separating aboveground from belowground communities for both fungi (Figure 2A) and bacteria (Figure 2B). PERMANOVA also indicated significant differences among aboveground and belowground parts of the plants, both for fungi (p-value = 0.001, R² = 0.3963) and bacteria (p-value = 0.001, R² = 0.3352). Additionally, pairwise PERMANOVA confirmed that microbial communities differed significantly between plant compartments. For both fungi and bacteria, aboveground compartments (fruit and leaf) were distinct from belowground compartments (bulk soil and rhizosphere), whereas differences between bulk soil and rhizosphere were smaller (Supplementary Table 5), highlighting the close connection between rhizosphere and surrounding soil.

We then investigated the relative abundance of the main microbial groups. We identified a total of 9 fungal phyla (Supplementary Table 6) and 10 bacterial phyla (Supplementary Table 7), together with 14 genera for both fungi and bacteria (Figures 3A,B). At the phylum level, fungi (Supplementary Figure 2A) showed a dominance of Ascomycota across all compartments except for the rhizosphere. Following in decreasing order, we identified Basidiomycota, Mortierellomycota, and Rozellomycota. For bacteria (Supplementary Figure 2B), Proteobacteria emerged as the dominant phylum, followed by Actinobacteria in all the compartments; whereas both the rhizosphere and bulk soil displayed a higher relative abundance of Acidobacteria. This pattern already highlights a compositional difference among the different plant-associated compartments. At the genus level (Supplementary Table 7), compositional differences became even more pronounced for fungi such as Mortierella, which was particularly abundant belowground; Serendipita and Pezoloma, typical of the rhizosphere; whereas Ascochyta and Epicoccum were characteristic of aboveground compartments, revealing a clearly different community structure. A similar pattern was observed for bacteria (Supplementary Table 8): genera such as Aciditerrimonas and Actinomadura were more prevalent in the belowground compartments, whereas Massilia and Hymenobacter were more abundant in the aboveground compartment. Notably, a significant proportion of the fungal community, exceeding 25%, remained unclassified across all plant tissues analyzed. This result contrasts sharply with observations in bacterial communities, where the abundance of unclassified taxa was considerably lower. At the genus level, bacterial communities showed a greater degree of diversity, with the “Other” category which represents the most abundant group across all compartments and including all taxa not within the top 15 most abundant.

To more precisely define the taxa specifically associated with each compartment, or predictive biomarkers taxa, we trained a random forest model that predicts the compartment using taxonomic profiles at the genus level as input (Figures 3C,D). We identified several predictive biomarkers taxa that exhibit compartment-specific associations. For fungi, the bulk soil was characterized by the presence of Linneamannia, Talaromyces and Mortierella, while we found in the rhizosphere genera such as Trichoderma and Sporothrix as predictive biomarkers taxa. The leaf surface exhibited a distinct set of predictive biomarkers taxa, including Epicoccum, Ascochyta, Cladosporium, Aureobasidium, and Peniophora as well as Vishniacozyma. In contrast, the fruit epiphytic community was enriched in Taphrina, Alternaria and Curvibasidium. Similarly, bacterial predictive biomarkers taxa varied across compartments. In bulk soil, we identified Parachlamydia while the rhizosphere harbored Aciditerrimonas, Acidibrevibacterium, Acidibacter, Mycobacterium, Iamia. The leaf surface was dominated by Methylobacterium, Sphingomonas and Massilia. The fruit surface microbiota was enriched in Pseudomonas, alongside Pantoea, Micrococcus, and Acinetobacter.

To assess the influence of edaphic factors on microbial taxa across different soil niches we analyzed Spearman’s correlation coefficients (ρ) between microbial genera (Supplementary Tables 10, 11) (in both rhizosphere and bulk soil) and physico-chemical parameters (NH₄⁺, P, pH, salinity and texture). We found that pH was the primary influencer of the abundance of fungal genera (Figure 4A). For some genera relative abundance increased with pH (as for example Plectosphaerella) while in other cases (e.g., Cladophialophora) the correlation was negative. Additionally, salinity and phosphorus were identified as influential factors, both showing negative correlations (respectively with Acidomyces and Pyrenochaetopsis), whereas different concentrations of NH₄⁺ and texture did not show any effect. In the rhizosphere, several fungal taxa were notably influenced by the edaphic factors analyzed (Figure 4C). For example, in the rhizosphere Oidiodendron and Serendipita were influenced by P concentration, pH and salinity; specifically, Oidiodendron showed a negative correlation with pH and a positive correlation with salinity, while Serendipita was negatively correlated with P and positively correlated with salinity. Interestingly, we found some genera that were influenced by different factors depending on the compartment analyzed. For example Acidomyces was negatively correlated with salinity in bulk soil, while in the rhizosphere was negatively correlated with NH₄⁺ and texture. The correlation with salinity was not significant.

Regarding bacteria, they were primarily influenced by salinity and soil texture, with additional contribution from pH and phosphorus levels (Figure 4B), while NH₄⁺ did not exhibit any significant effect. In the rhizosphere, pH emerged as the most influential factor shaping bacterial community composition, showing both positive and negative correlations with different taxa. Phosphorus and soil texture also influenced the taxonomic composition, while salinity appeared to have a relatively minor effect. As observed in bulk soil, NH₄⁺ showed no significant influence on bacterial taxa in the rhizosphere (Figure 4D). In the rhizosphere, Gp 2 was correlated negatively with phosphorus concentration while it was positively correlated with both salinity and texture. Iaima showed a negative correlation with pH. Also for bacteria we found that some taxa, such as Acidibacter, responded differently depending on the compartment analyzed. In particular, in the rhizosphere Acidibacter was negatively correlated only with pH while in bulk soil was positively correlated only with soil texture.

The distance-based redundancy analysis (dbRDA) of beta diversity revealed distinct differences in fungal community composition between rhizosphere and bulk soil samples and how the soil physico-chemical parameters significantly explained the variation in microbial community structure (Figures 5A,B). Both for bacteria and fungi dbRDA analyses showed a clear separation between rhizosphere (R) and bulk soil (B) samples, accounting for 47.3% of the variance in fungal communities and 58.4% in bacterial communities. This separation was primarily driven by edaphic factors such as salinity, ammonium (NH₄⁺) concentration, and pH. Salinity and NH₄⁺ showed a positive association with rhizosphere samples, whereas pH was negatively correlated with them. The regression of environmental variables with ordination axes (function “envfit,” See methods) revealed significant associations between environmental variables and microbial community structure in the dbRDA ordination space (Supplementary Table 11). For fungal communities, pH, available phosphorus, ammonium (NH₄⁺), and salinity were significantly correlated with community composition, while soil texture was not significant. pH aligned primarily with dbRDA1, while phosphorus was oriented along dbRDA2, indicating distinct environmental gradients influencing fungal distribution. For bacterial communities, pH and ammonium were the strongest predictors. Salinity and phosphorus were also significant but explained less variance. Moreover, we assessed the influence of plant genotype on microbial community composition, which revealed significant differences for both fungi and bacteria. Genotype explained 42.6% of the variance in bacterial communities and 40.4% of the variance in fungal communities (Supplementary Table 13).

After identifying the different predictive biomarkers taxa, we further analyzed the distribution of ASVs to measure the overlap between aboveground and belowground compartments, for both bacterial and fungal communities (Figures 6A,B) in order to assess whether there is a flow of taxa between the belowground and aboveground compartments. We found only 9 ASVs and 12 ASVs shared among all compartments for fungi and bacteria, respectively. In the belowground compartments, we identified 235 ASVs exclusive to the fungal rhizosphere and 3,977 ASVs for bacteria, 65 ASVs exclusive to fungal bulk soil and 3,061 ASVs for bacteria, and 158 fungal ASVs and 1,258 bacterial ASVs shared between rhizosphere and bulk soil. In contrast, in the aboveground compartments, 72 fungal ASVs and 241 bacterial ASVs were unique to the leaf surface, 14 fungal and 43 bacterial ASVs were exclusive to the fruit surface, and 74 fungal ASVs and 192 bacterial ASVs were shared between these two compartments.

To better understand the distribution of ASVs shared across all compartments (9 for fungi and 12 for bacteria), we analyzed their relative abundances for fungi (Supplementary Table 14) and bacteria (Supplementary Table 15). Fungal ASVs (Supplementary Figure 4A) exhibited distinct abundance patterns between aboveground (green) and belowground (brown) compartments. Typically, ASVs showed asymmetrical abundance distributions. For instance, Penicillium was notably enriched in the rhizosphere but less abundant aboveground. In contrast, Fusarium predominated in bulk soil, showing reduced abundance aboveground. Conversely, genera such as Epicoccum and Cladosporium consistently showed higher abundance in aboveground compartments. A similar trend was observed for bacterial communities (Supplementary Figure 4B). For instance, Methylobacterium and Sphingomonas were more abundant in the epiphytic communities of aboveground compartments, whereas Bacillus and Microvirga exhibited higher relative abundance in belowground habitats and were relatively scarce in the phyllosphere.