Differential Response of Wheat Rhizosphere Bacterial Community to Plant Variety and Fertilization

The taxonomic assemblage and functions of the plant bacterial community are strongly influenced by soil and host plant genotype. Crop breeding, especially after the massive use of nitrogen fertilizers which led to varieties responding better to nitrogen fertilization, has implicitly modified the ability of the plant root to recruit an effective bacterial community. Among the priorities for harnessing the plant bacterial community, plant genotype-by-microbiome interactions are stirring attention. Here, we analyzed the effect of plant variety and fertilization on the rhizosphere bacterial community. In particular, we clarified the presence in the bacterial community of a varietal effect of N and P fertilization treatment. 16S rRNA gene amplicon sequence analysis of rhizospheric soil, collected from four wheat varieties grown under four N-P fertilization regimes, and quantification of functional bacterial genes involved in the nitrogen cycle (nifH; amoA; nirK and nosZ) were performed. Results showed that variety played the most important role and that treatments did not affect either bacterial community diversity or bacterial phyla abundance. Variety-specific response of rhizosphere bacterial community was detected, both in relation to taxa (Nitrospira) and metabolic functions. In particular, the changes related to amino acid and aerobic metabolism and abundance of genes involved in the nitrogen cycle (amoA and nosZ), suggested that plant variety may lead to functional changes in the cycling of the plant-assimilable nitrogen.


Introduction
The plant bacterial community, and especially the rhizospheric bacterial community, has been clearly related to crop yield and resistance to biotic and abiotic stresses [1][2][3][4][5][6]. Domestication and breeding have strongly influenced plant bacterial community [7,8], resulting in some cases in a reduced ability of the plant to attract an adequate rhizosphere bacterial community [9]. Since the diversity of the plant bacterial community is contributing to the resilience and sustainability of agroecosystems to multiple stressors, including climate change, increased diversity of the crop-associated bacterial community can contribute to conjugating sustainable yield with a reduced footprint of agriculture on the environment. Indeed, in general terms, plant genotype is a relevant factor affecting diversity and functions of the plant bacterial community and influencing soil bacterial community [3,[10][11][12]. Consequently, there is a growing focus on considering the recruited plant bacterial community for future crop breeding programs [3]. Among the priorities for harnessing the plant bacterial community, plant genotype × microbiome × environment × management interactions are stirring attention [13]. In fact, not only the species, but also the individuality of genotypes, viz. germplasm or varieties, show peculiarities in the associated bacterial community stresses [5,8,14].

Physicochemical Characterization of the Soil
Results of the physicochemical analysis of the soil samples are reported in Table 1. For each plot, the analysis was performed on bulk soil samples. A Principal Component Analysis (PCA) was performed on soil physicochemical parameters and the location of sample collection (i.e., latitude and longitude of the plots). No evident clustering was detected, either with respect to physicochemical parameters and geographical coordinates ( Figure 1).

Diversity of the Soil and Rhizosphere Bacterial Community
A total of 3,180,595 16S rRNA reads were obtained, and 2,131,486 (67% of total reads) passed quality filtering (Table S1). After the clustering step, a total of 12,795 amplicon sequence variants (ASVs) were obtained (Dataset S1). To Bacteria, 10,271 ASVs were assigned, 59 ASVs to Archaea, and 2,465 were unassigned (NA). ASVs assigned to Archaea and NA were removed from the following analyses. Rarefaction curves obtained from ASVs reached plateaus for all samples (Figure S1), indicating a satisfactory survey of the bacterial diversity (Good's coverage Table S2), which allowed to estimate alpha diversity indices (Table S3). Bulk soil and rhizosphere samples did not differ in alpha diversity nor in relation to both  Table 1.

Diversity of the Soil and Rhizosphere Bacterial Community
A total of 3,180,595 16S rRNA reads were obtained, and 2,131,486 (67% of total reads) passed quality filtering (Table S1). After the clustering step, a total of 12,795 amplicon sequence variants (ASVs) were obtained (Dataset S1). To Bacteria, 10,271 ASVs were assigned, 59 ASVs to Archaea, and 2,465 were unassigned (NA). ASVs assigned to Archaea and NA were removed from the following analyses. Rarefaction curves obtained from ASVs reached plateaus for all samples ( Figure S1), indicating a satisfactory survey of the bacterial diversity (Good's coverage Table S2), which allowed to estimate alpha diversity indices (Table S3). Bulk soil and rhizosphere samples did not differ in alpha diversity nor in relation to both variety and fertilization (Permanova p > 0.05). Varieties differed in relation to the Shannon index, but not in relation to the Simpson index ( Figure 2, see specific pairwise comparison for p-values), indicating that species richness and species abundance distribution (equitability) changed among varieties, while the dominance of bacterial groups did not. Concerning the taxonomic diversity among samples, alpha and beta diversities of soil and rhizosphere samples were not separated (Figure 3), while plant varieties were, with evident grouping for Andriolo and Bologna, and less evident, but significant for Verna and Sieve (Table S4). Concerning the taxonomic diversity among samples, alpha and beta diversities of soil and rhizosphere samples were not separated (Figure 3), while plant varieties were, with evident grouping for Andriolo and Bologna, and less evident, but significant for Verna and Sieve (Table S4)  Fertilization did not affect the taxonomic diversity of the samples. Statical analy was performed using a permutational multivariate analysis of variance (PERMANOV on centered log-ratio-transformed samples (Table S4). Figure 4 reports the relative abundances of taxonomies at the phylum level for t top 20 ASVs detected (which represent relative abundances from 0.1% to 0.08%). Tax nomic abundances in relation to variety and fertilization are shown in Figure S2. The m represented phyla were Acidobacteriota, Actinobacteriota, Chloroflexi, Proteobacter and Verrucomicrobia. Differences among varieties were present in relation to phyla abu dance ( Figure 5). In fact, differential abundance analyses showed remarkable differenc in Actinobacteriota, Chloroflexi, Desulfobacterota, Entotheonellaeota, Firmicutes, Myx coccota, Nitrospirota, and the candidate division RCP2-54, in particular for Verna varie as confirmed by a Permanova on the PCA to the centered log-ratio ( Figure S3). When co sidering single ASVs, 37 of them were found with statistically different abundance in lation to variety (Table S5). These ASVs included members of the orders Burkholderial Rhizobiales (Beijerinckiaceae, Hyphomicrobiaceae, Xanthobacteraceae), Blastocatellal Pyrinomonadales, as well as of Actinobacteriota. Fertilization did not affect the taxonomic diversity of the samples. Statical analysis was performed using a permutational multivariate analysis of variance (PERMANOVA) on centered log-ratio-transformed samples (Table S4). Figure 4 reports the relative abundances of taxonomies at the phylum level for the top 20 ASVs detected (which represent relative abundances from 0.1% to 0.08%). Taxonomic abundances in relation to variety and fertilization are shown in Figure S2. The most represented phyla were Acidobacteriota, Actinobacteriota, Chloroflexi, Proteobacteria, and Verrucomicrobia. Differences among varieties were present in relation to phyla abundance ( Figure 5). In fact, differential abundance analyses showed remarkable differences in Actinobacteriota, Chloroflexi, Desulfobacterota, Entotheonellaeota, Firmicutes, Myxococcota, Nitrospirota, and the candidate division RCP2-54, in particular for Verna variety, as confirmed by a Permanova on the PCA to the centered log-ratio ( Figure S3). When considering single ASVs, 37 of them were found with statistically different abundance in relation to variety (Table S5). These ASVs included members of the orders Burkholderiales, Rhizobiales (Beijerinckiaceae, Hyphomicrobiaceae, Xanthobacteraceae), Blastocatellales, Pyrinomonadales, as well as of Actinobacteriota.

Prediction of Potential Functions and Metabolic Pathways
Functional profiles of the bacterial community were inferred. All the identified metabolic pathways present in the taxa retrieved in the bacterial community are listed in Table S6. We identified a total of 422 pathways, all distributed in all samples. However, pathways' abundance showed a clustering according to varieties (Figure 6a). Using a SIMPER test, we then selected the top three pathways that mostly differentiate the four varieties (in terms of amount of variance explained). Six different pathways-(VALSYN-PWY, ILEUSYN-PWY, BRANCHED-CHAIN-AA-SYN-PWY, PWY-3001, PWY-3781, and PWY-7111) were selected, most of them were related to the metabolism of amino acids and aerobic respiration. In particular, the pathways of L-valine biosynthesis (VALSYN-PWY), L-isoleucine biosynthesis (ILEUSYN-PWY and PWY-3001) and aerobic respiration (PWY-3781) were found. Interestingly, this subset of pathways also suggested a clustering related to, in particular, Verna being grouped apart from the other three varieties (Figure 6b). Pathways related to the nitrogen cycle showed only a minor contribution to the overall differences (i.e., the denitrification pathway accounted to 0.05% of total variance).

Quantitative PCR Analysis of the Nitrogen Cycle Genes
In order to overcome the limitation of PICRUSt-based pathway predictions (which are based on overall pathways and on phylogenetic reconstruction) and more clearly define the effect of plant variety and fertilization on the bacterial community fraction involved in the nitrogen cycle, bacterial genes related to nitrogen fixation (from N 2 to NH 3 , nifH), nitrification (conversion of ammonia NH 3 to hydroxylamine, amoA), and denitrification (nitrite to nitric oxide, nirK and nitrous oxide to N 2 , nosZ) were detected, and their abundance was estimated by quantitative PCR on rhizosphere samples (Figure 7). While for nifH no differences among varieties or fertilizations were detected, amoA, nirK, and nosZ were different with respect to varieties (Scott-Knott test p < 0.05, Table S7). In particular, Andriolo and Bologna showed higher abundance of nirK and amoA genes than Sieve and Verna did. Verna also had lower nosZ values than the other three varieties did, suggesting that the denitrification pathway could be less effective on the rhizosphere of this variety.
ways' abundance showed a clustering according to varieties (Figure 6a). Using a SIMPER test, we then selected the top three pathways that mostly differentiate the four varieties (in terms of amount of variance explained). Six different pathways-(VALSYN-PWY, ILE-USYN-PWY, BRANCHED-CHAIN-AA-SYN-PWY, PWY-3001, PWY-3781, and PWY-7111) were selected, most of them were related to the metabolism of amino acids and aerobic respiration. In particular, the pathways of L-valine biosynthesis (VALSYN-PWY), Lisoleucine biosynthesis (ILEUSYN-PWY and PWY-3001) and aerobic respiration (PWY-3781) were found. Interestingly, this subset of pathways also suggested a clustering related to, in particular, Verna being grouped apart from the other three varieties ( Figure  6b). Pathways related to the nitrogen cycle showed only a minor contribution to the overall differences (i.e., the denitrification pathway accounted to 0.05% of total variance).

Quantitative PCR Analysis of the Nitrogen Cycle Genes
In order to overcome the limitation of PICRUSt-based pathway predictions (which are based on overall pathways and on phylogenetic reconstruction) and more clearly define the effect of plant variety and fertilization on the bacterial community fraction involved in the nitrogen cycle, bacterial genes related to nitrogen fixation (from N2 to NH3, nifH), nitrification (conversion of ammonia NH3to hydroxylamine, amoA), and denitrification (nitrite to nitric oxide, nirK and nitrous oxide to N2, nosZ) were detected, and their abundance was estimated by quantitative PCR on rhizosphere samples (Figure 7). While for nifH no differences among varieties or fertilizations were detected, amoA, nirK, and nosZ were different with respect to varieties (Scott-Knott test p < 0.05, Table S7). In particular, Andriolo and Bologna showed higher abundance of nirK and amoA genes than Sieve and Verna did. Verna also had lower nosZ values than the other three varieties did, suggesting that the denitrification pathway could be less effective on the rhizosphere of this variety.

Discussion
The taxonomic assemblage and functions of the plant bacterial community are strongly influenced by soil and host plant genotype [5,7,8]. Crop breeding, especially after the massive use of nitrogen fertilizers which led to varieties responding better to nitrogen fertilization [36], has implicitly decreased the pressure toward plant genotypes highly ef-

Discussion
The taxonomic assemblage and functions of the plant bacterial community are strongly influenced by soil and host plant genotype [5,7,8]. Crop breeding, especially after the massive use of nitrogen fertilizers which led to varieties responding better to nitrogen fertilization [36], has implicitly decreased the pressure toward plant genotypes highly efficient in the use of native soil nutritional resources, then indirectly reducing the ability of the plant root to recruit an effective bacterial community [7,8,37]. Tough physicochemical analysis of soil did not detect evident groupings according to wheat varieties (apart from Andriolo), and the analysis of the bacterial community did, clearly differentiating the four varieties under study (Andriolo being the most heterogeneous, Figure 3). In fact, the overall body of our results showed that variety played the most important role and that fertilization did not affect either bacterial community diversity or bacterial phyla abundance. Interestingly, rhizosphere bacterial community and bulk soil bacterial community were highly similar, suggesting that under the plant density and the depth we sampled, the root apparatus could extend its effects at a certain distance from the plant stem (30 cm). This bacterial community was in agreement with the extensive root apparatus of wheat, which allows affecting the soil at a distance from the plant [38]. Consequently, we may consider "bona fide" our bulk soil samples as very close to a low-bound rhizosphere soil. It would be interesting to evaluate the extent in space around the plant of this assimilation of soil on the rhizosphere bacterial community, which can have relevance for the setting-up and interpretation of results from intercropping trials [39]. The three old varieties and the modern Bologna variety showed a differential abundance of several microbial taxa. Among them, Acidobacteriota, Bacteroidota, Chloroflexi, Firmicutes, and Nitrospirota were particularly relevant. Previous studies on wheat rhizosphere showed that N fertilization decreased the abundance of Acidobacteria and increased Bacteroidetes [20]. In our work, a decrease of Acidobacteriota was shown by Andriolo variety in the plots with increased P-fertilization ( Figure S2), which may suggest that this variety may respond differently to P fertilization, in turn affecting its rhizosphere bacterial community. However, both Andriolo and Bologna showed levels of abundance of members of phyla Acidobacteriota and Bacteroidota higher than that of Sieve and Verna varieties. Since both Sieve and Verna derive from breeds with "Est Mottin 72", we could not a priori exclude that the lower levels of Acidobacteriota and Bacteroidota in their rhizosphere may reflect this partially common genetic basis, highlighting a contribution of plant genotype to the interaction with specific bacterial phyla. If confirmed in following experiments, these bacterial phyla could be good candidates for quantitative trait loci (QTL) analysis identifying in breeds from the old Verna variety and the modern Bologna. The high abundance of members of Acidobacteriota also agrees with previous reports (see, for instance, [33,40]). Moreover, Bacteroidota, Chloroflexi, and Firmicutes phyla were recognized as important in differentiating wheat species (T. durum Desf. and T. aestivum L.) [33] and soil-varieties combinations [31], suggesting the presence of loci in the wheat genome which can be identified for interaction with the rhizosphere microbiota. To date, few loci/genes in a relatively small number of plant species have been identified as related to the recruitment of microbiota [41], and no studies have been performed on wheat. Interestingly, ASVs affiliated to the uncultured phylum RCP2-54 were found more abundantly in Sieve and mainly in Verna, compared to Bologna and Andriolo. The functional role of this phylum is still unknown; however, previous works showed its presence in the rhizosphere and in differentiating tomato and maize varieties under water deficit conditions [42], again suggesting a possible genetic basis in the plant for its differential recruitment. At the level of ASVs it was relevant to find members of Blastocatellales and Pyrinomonadales. A recent paper investigating soil suppressiveness to Fusarium culmorum in wheat, found Blastocatellales in the network associated with pathogen suppression [43], and Pyrinomonadales were reported to vary in abundance in wheat rhizosphere in relation to soil nitrate level [44]. This evidence may suggest that the four varieties in our study may be able to differentially recruit a fraction of the soil bacterial community able to improve plant health and contributing differently to plant use of nitrogen. Bologna variety had the highest levels of both Blastocatellales and Pyrinomonadales, although high variation between plots was present ( Figure S4). The role of these bacterial groups on plant tolerance to pathogenic fungi and nitrogen use may deserve further attention.
Under the light of possible functional differences among the bacterial community recruited by the four plant varieties, PICRUSt2 results ( Figure 6) mirrored in part the taxonomic profiling in varieties differentiation. In particular, Verna grouped separately from the other three varieties, similarly to the pattern shown in the nMDS analysis of ASVs ( Figure 3). Interestingly, the metabolic pathways which mostly differentiated Verna from the other varieties were related to aerobic respiration and amino acid biosynthesis. We may speculate that root exudates from the four varieties differed in the amino acid composition, hence differentially promoting the growth of bacteria harboring specific biosynthesis pathways. Indeed, root exudates are rich in amino acids, and previous work showed that amino acid content strongly varies among different plant varieties. Concerning the aerobic respiration pathway, Verna showed a more abundant representation than the other varieties did, though with differences among fertilization regimes. This may suggest that the other varieties could hold a higher fraction of anaerobic bacteria/anaerobic metabolism. Verna rhizosphere bacterial community contained a higher amount of Actinobacteriota and Firmicutes, phyla including several aerobic bacterial groups thriving in soil, as Streptomyces and Bacillus.
Concerning possible differences in nitrogen cycling by the rhizosphere bacterial community (which is partially related to anaerobic metabolism), variation in the abundance of members of Nitrospirota was found in Verna and Sieve varieties compared to that in Bologna and Andriolo varieties. Nitrospirota includes Nitrospira, a well-known genus of nitrite oxidizers present in agricultural soil and in the plant rhizosphere [45]. Interestingly, this result agrees with qPCR estimation of functional gene abundance. In fact, the analysis of genes involved in the nitrogen cycle showed that amoA (nitrification) and nosZ (denitrification) were differentially abundant with respect to varieties, suggesting that variety may lead to functional changes in the cycling of the plant-assimilable nitrogen. However, we could not a priori exclude fertilization may partially contribute to this difference. Experiments with bulk soil samples taken at higher distances from the plant may be helpful to clarify how much of the differences resulted from the sole fertilization. Indeed, it would highly relevant from the perspective of low-input precision agriculture to clarify the extent and impact on plant productivity of such variety-fertilization interaction in the rhizosphere bacterial community, to exploit the potentialities of variety-specific bacterial community. Multiple seasons and long-term field experiments in different soil and climatic conditions would be needed to fully address this question.
Soil tillage was carried out to a depth of 0.4 m with a moldboard plow in October 2019, followed by a tandem disk harrow (0.1 m depth) to break clods. The organic fertilizer "Endurance N8" (Unimer s.p.a., Milano, Italy), containing organic nitrogen from animal protein source, and the organic soft rock phosphate "GAFSA 27" (Panfertil s.p.a, Ravenna, Italy) were used as N and P 2 O 5 sources, respectively. On 27 November, 2019, both fertilizers were broadcasted on the soil surface, immediately followed by spike tooth harrow (0.05 m depth) to incorporate the fertilized and prepare the false seedbed. Weeds were uprooted to a depth of 0.05 m using a spring tine harrow on 7 and 17 December 2019. Winter wheat seeds were sown 90 kg ha −1 on 17 December 2019 with a row spacing of 0.13 m. Grain harvesting was performed at crop commercial maturity (grain moisture lower than 13%) on 15 July 2020.

Sampling and Samples Treatment
Sampling of bulk soil (thereafter indicated as soil) and rhizosphere soil was performed in May 2020. For each plot, five individual plants at flowering stage were selected, separated at least 5 m from each other. After digging up the entire plant (10 cm radius around the plant, and down to 20 cm soil depth) and shaking off the plant roots, the soil remaining attached to the roots was detached with a sterile scalpel and used as rhizosphere soil. Bulk soil was collected using a stainless steel sterile corer to sample the top 10 cm of soil in patches of bare soil between plants (ca. 30 cm apart from the plant) and sieved at 2 mm. Bulk and rhizosphere soil samples for each plot were combined to create composite bulk and rhizosphere soil samples. A total of 32 combined samples from 16 plots were analyzed ( Table 2). Samples were stored at −80 • C until DNA extraction. Table 2. Samples analyzed. The code, wheat variety, the fertilization of the plot, and their geographical position is reported. Numbers in the sample code refer to the plots. "S" and "R" in the sample code indicate soil and rhizosphere fraction, respectively.

Soil Chemical and Physical Analyses
Soil pH was determined in a 1:5 soil-to-distilled water solution. Soil total nitrogen (STN) was determined by applying the dry combustion method on 5.0 mg soil samples, using a CHNS analyzer (CHN-S Flash E1112, Thermo-Finnigan LLC, San Jose, CA, USA) [48,49]. Soil samples (5.0 mg) were pretreated with chloride acid solution (10%) within silver cups to volatize the carbonate carbon and then analyzed by means CHNS to determine soil organic carbon (SOC) [50,51]. Ammonium nitrogen and nitrate nitrogen were determined by spectrophotometric determination (Lambda 20 spectrophotometer; PerkinElmer, Waltham, MA, USA) after extraction with the calcium chloride (CaCl 2 ) procedure by Houba et al., (1995), which has the advantage of extraction uniformity for the considered N forms [52]. Aliquots of soil extract solution were treated according to the nitrate copper-cadmium reduction method [53], and the resulting solutions were spectrophotometrically analyzed at 540 nm after Griess reaction [54] for the determination of nitrate nitrogen. Further, aliquots of soil extract solution were treated following the Nessler method [55] and then absorbance analyzed at 420 nm for the determination of ammonium nitrogen. According to [56], soil samples were extracted with 0.5 M NaHCO 3 at pH 8.5 to solubilize inorganic P forms which were spectrophotometrically analyzed at 882 nm after reaction with ammonium molybdate to determine the soil assimilable phosphorous (AP). Exchangeable cations (calcium, iron, magnesium, potassium, selenium, sodium, and zinc) were extracted from the soil using a 1 N NH 4 OAc extracting solution pH 7.0 and then determined using an ICP-OES analyzer (ICAP™ 7400 ICP-OES Analyzer, Thermo Scientific, Waltham, MA, USA). The soil cation exchange capacity (CEC) was determined using the BaCl2-triethanolamine method as proposed by Mehlich (1939) [57].

eDNA Extraction and 16S rRNA Gene Amplicon Sequencing
Environmental DNA (eDNA) was extracted from 500 mg of soil using a DNeasy PowerSoil Pro (Qiagen). From the extracted DNA, the bacterial V4 region of 16S rRNA genes was amplified with primers 515F (5 -GTGCCAGCMGCCGCGGTAA-3 ) and 806R (5 -GACTACHVGGGTATCTAATCC-3 ) [58] in a 25 µL total volume with KAPA HiFi HotStart ReadyMix, 1 µM of each primer, and 10 ng of template DNA with 25 cycles with the following temperature profile: 30 s at 95 • C, 30 s at 55 • C, and 30 s at 72 • C. PCR products were sequenced in a single run using Illumina MiSeq technology with pair-end sequencing strategy and a MiSeq Reagent Kit v3 (Illumina, USA). Library preparation (Nextera XT, Illumina, San Diego, CA, USA) and demultiplexing were performed following Illumina's standard pipeline as previously reported [10].

Quantitative PCR of Genes Involved in the Nitrogen Cycle
Quantitative PCR was performed on eDNA using primers for the amplification of nifH, amoA, nosZ, nirK genes (primer sequences and citations are reported in Table S8). Total bacterial load was estimated on the 16S rRNA gene. Reactions were performed with 5 ng of template DNA in 10 µL total volume on a QuantStudio7 Flexi apparatus (Applied Biosystems) with Sybr Green technology (Maxima SYBR Green qPCR Master Mix, Thermofisher). Reactions were performed in triplicates. Relative abundance of nifH, amoA, nosZ, nirK genes was expressed in relation to 16SrRNA, as previously reported [59,60].

Bioinformatic and Statistical Analyses
Illumina reads were trimmed using the "Trim Galore!" tool on Galaxy server (https: //usegalaxy.org/, last accessed on 21 December 2021). Paired-end sequences were clustered into Amplicon Sequence Variants (ASVs) following the DADA2 pipeline (version 1.16) [61]. After filtering the sequences and removing the chimeras, the taxonomy assignment was carried out comparing our data against the SILVA NR99rel138 standard database of bacteria [62] using "DECIPHER" R package (version 2.18.1) [63] as implementation of DADA2 (SSU version 138 available at: http://www2.decipher.codes/Downloads.html, last accessed on 21 December 2021). Annotated ASVs count tables were processed in Phyloseq package in R environment version 4.0.5 [64].
The analysis of microbial communities was performed through the "Phyloseq" R package (version 1.34.0) [64]. For alpha diversity analysis, the Shannon and Simpson indices were calculated and plotted using the function "diversity()" within "microbiome" R package (version 1.12.0) [65]. Good's coverage and Evenness indices were calculated through the R functions "goods()" and "evenness()", respectively, within the "microbiome" R package (version 1.12.0). A Wilcoxon test for multiple comparison of averages was performed on alpha diversity indices using the "ggviolin()" and "stat_compare_means()" within the "ggpubr" R package (version 0.4.0). For ordination plots of phyloseq objects, a multivariate analysis based on Bray-Curtis distance and NMDS ordination generated using the "ordinate" function was performed, and plots were generated using the "plot_ordination()" function within phyloseq package. Rarefaction curves were generated using the "ggplot2" (version 3.3.3) [66] and "ranacapa" (version 0.1.0) [67] R packages using the "ggrare()" function on the phyloseq object. The "ggplot2" R package (version 3.3.3) was used to generate relative abundance plots.
Statistical analyses for quantitative PCR data were performed in R environment (version 4.0.5) calculating the ratio between the threshold cycles of the target genes and 16S gene. Graphics of the ratio averages were performed using the "ggboxplot()" function within the "ggplot2 R package (version 3.3.3). A Scott-Knott test was performed on averages using the functions "with()" and "SK()" within the "ScottKnott" R package (version 1.2-7) [68].

Prediction of Functional Abundances
The software PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States, https://github.com/picrust/picrust2, last accessed on 21 December 2021) was utilized to investigate and predict the functional potential for each determined sequence in the microbial community. PICRUSt2 pipeline returned as output the metabolic pathway abundances, and functional gene profiling by comparison with MetaCyc and Kyoto Encyclopedia of Genes and Genomes databases [69]. Ranking of the most relevant pathways for varieties differentiation was conducted by running a SIMPER test using the function "simper()" with the "vegan" R package (version 2.5-7). For each contrast between varieties, the top three differentially expressed pathways were selected, and analyses were carried out using the function "metaMDS()" in the "vegan" R package (version 2.5-7) to perform nonmetric multidimensional scaling (nMDS) and the function "heatmap()" in the "stats" R package (version 4.0.5) to graphically represent results.

Links to Deposited Data
The sequences dataset was deposited in the SRA database under the BioProject PRJNA720503 (last accessed on 21 December 2021).