Introduction

Copper (Cu) compounds have traditionally been the means to combat phytopathogenic microbes, especially fungi, in crop plants1. In organic agriculture including organic vineyards, the use of Cu-based pesticides is currently the only chemical treatment allowed, although it is limited to a maximum of 6 kg Cu ha−1 per year in the EU (Commission Regulation, 2002). Such a strict regulation is explained by the adverse effects of Cu on soil organisms and its long-term persistence on surface horizons2. Another rational behind this is that Cu overuse may lead to the development of Cu resistance in pathogenic fungi3. Lastly, Cu residues on grapes may compromise wine quality4. In contrast, biocontrol agents, such as lactic acid bacteria, are not harmful to either the environment or health. Indeed, they are classified as “generally recognized as safe” (GRAS) by the Food and Drug Administration (FDA, USA) and have received the “qualified presumption of safety” (QPS) status from the European Food Safety Agency (EFSA)5.

The biocontrol potential of the lactic acid bacterium Lactobacillus plantarum has recently emerged as an important subject of study. The effectiveness of L. plantarum against several fungi has been emphasized in a large-scale screening study of over 7.000 lactic acid bacteria, which showed that most of the strains with antifungal properties belonged to L. plantarum6. Furthermore, many in vitro studies have demonstrated the antifungal properties of L. plantarum strains against plant parasites of the genera Botrytis3,7,8,9, Fusarium3,6,7,10,11, Aspergillus10,12, Rhizopus7,12, Alternaria, Phytophtora and Glomerella3, Sclerotium, Rhizoctonia and Sclerotinia7. Additionally, chili seeds infected by Colletotrichum gloeosporioides germinated well after treatment with L. plantarum13. In the field, a mix of L. plantarum and Bacillus amyloliquefaciens applied to durum wheat from heading to anthesis showed promising control of the fungal pathogens F. graminearum and F. culmorum14.

In vitro trials using different L. plantarum strains have shown a promising inhibition of plant pathogenic bacteria belonging to the genera Clavibacter15, Xanthomonas5,16,17, Erwinia16 and Pseudomonas8,17,18. In fact, the inhibitory activity of L. plantarum against P. syringae in planta has long been documented16, while that against Rhizobium radiobacter has only recently been reported19. Meanwhile, field sprays of L. plantarum on Chinese cabbage alleviated the severity of soft rot by Pectobacterium carotovorum subsp. carotovorum20 and contributed to a reduction of fire blight by Erwinia amylovora on apple and pear21. However, the aforementioned studies focused on the effects of L. plantarum on single organisms while its effects on the microbial community in a field remain to be investigated.

Next-generation sequencing (NGS) and quantitative PCR (qPCR) techniques facilitate rapid and cost-competitive mapping of complex microbiomes by tracking fastidious, unculturable or even unknown taxa and their roles22,23. In particular, through accurate, real-time enumeration of taxonomic or functional gene markers, qPCR has enabled microbiologists to quantify specific species or phylotypes out of an environmental “genetic soup”23. Nonetheless, qPCR assays rely exclusively on known genes and thus overlook taxa distinct from those already described24. This limitation is circumvented by NGS technologies and the approaches of either marker gene amplification (amplicon sequencing) or total (shotgun) sequencing of environmental DNA. Owing to the PCR amplification of specific molecular gene markers, amplicon sequencing can still be quite inefficient for drawing conclusions about the genus or species level25. Moreover, biases can be introduced due to horizontal gene transfer, interspecific gene copy number variation within a microbiome and underrepresentation25. Nevertheless, both NGS and qPCR have substantially contributed to explorations of the planet’s microbiome, and amplicon sequencing remains a valuable tool for comparative studies of microbial communities26,27.

From an anthropocentric viewpoint, there has been increasing interest in using the above technologies to decipher the microbial interactions that directly affect human health and resources (e.g. on crops and livestock). In this context, the positive role of probiotic bacteria has long been acknowledged and their potential is being investigated in some detail28. L. plantarum is a versatile lactic acid bacterium and probiotic, and is ubiquitous on plants. The bacterium has been isolated from various environments, such as fermented food products, the human mouth and grape must29, and hence constitutes a promising case study. L. plantarum has been examined for, among other things, its adequacy as a biocontrol agent against phytopathogens and, as mentioned previously, it is generally considered a good candidate for use in agriculture.

Even though available results of the antagonism of L. plantarum against phytopathogens in vitro are encouraging, harsh field conditions call for the careful design of field application trials16. For example, a lack of nutrients can be circumvented by repeated, surface spray applications to sustain high viable biocontrol numbers16. The scarcity of field studies makes it more difficult to evaluate this prediction, as well as the effect L. plantarum may have on a field’s native microbiome. The aim of this study was to investigate the dynamics of the epiphytic bacterial and fungal communities in two plots in the same vineyard, as the season progressed, also considering two different fungicide treatments applied: CuSO4 and a potential biocontrol strain of L. plantarum named MW-1. To authors’ knowledge, this study is the first evaluation of the dynamics of epiphytic bacterial and fungal community on grapevine treated leaves, sampling the whole growing season, from bud-break to post-harvest.

Results

This study used a combination of amplicon library NGS for both 16S and ITS, and qPCR with universal and strain-specific primers. The outcome of the DNA sequencing is presented below, followed by the results of qPCR analyses on the vine leaves. This study evaluated the dynamics of the epiphytic bacterial and fungal community on grapevine leaves along the growing season, as they change accordingly to the different fungicide treatments applied.

DNA sequencing dataset description

The sequencing provided a dataset of 120 samples. After quality filtering, 14,028 amplicon sequence variants (ASV) were attained from 3,011,656 high-quality reads concerning just the 16S community on leaves. For the ITS sequencing of the leaves, 6,901 ASV were retained, coming from 4,854,508 reads. All the samples were sequenced with sufficient coverage to unravel the complexity of the microbial community harboured on the leaves, as seen by the rarefaction curves in Fig. S1. To minimise statistical issues that could arise from differences in terms of sequencing depth, the samples were normalised by randomly extracting 24,000 reads from each of the samples before any downstream analyses were performed.

Alpha diversity on leaves

This section focuses on two different parameters that describe the microbial community on leaves: phylogenetic diversity (PD) and evenness. Alpha diversity results are summarised in Fig. 1. Looking at the bacterial community on leaves, it was noted that PD was significantly higher in the Cu-treated samples compared to the biocontrol treated, as shown in Fig. 1a (p-value = 0.014), and that PD seemed largely unaffected by the sampling time during the season (p-value = 0.538), as shown in Fig. 1c. Looking at the evenness in Fig. 1a, it can be seen that the Cu samples had a greater evenness than the biocontrol treated samples (p-value = 0.0003). In contrast, no statistical effect on evenness was evident when looking at the collection time for 16S (p-value = 0.38), Fig. 1c.

Figure 1
figure 1

Boxplots of phylogenetic diversity (PD) (left) and evenness (right), grouped by treatments (above) and period (below) for bacteria and fungi. (a) PD (left) and evenness (right) of bacteria grouped by treatment, (b) PD (left) and evenness (right) of fungi grouped by treatment, (c) PD (left) and evenness (right) of bacteria grouped by period, (d) PD (left) and evenness (right) of fungi grouped by treatment.

In contrast, the fungal community seemed to be unaffected by the treatments in terms of PD and evenness (p > 0.05), as visualised in Fig. 1b. However sampling time produced a strong effect on PD (p-value  0.05) and evenness (p-value  0.05) on the ITS distribution, as revealed in Fig. 1d. In particular, the number of sequence variants in the community tended to increase over time from September to November, and then suddenly decreased from December to the post-harvest period in February and March.

Beta diversity on leaves

This section shows the differences between the samples using Bray-Curtis dissimilarity, looking at the distribution and relative abundances of the single feature between samples. The beta diversities of 16S and ITS sequence datasets, coloured by treatment or collection time, were visualised by PCoA plots and are shown in Fig. 2. With regard to the bacterial community in Fig. 2a, it can clearly be seen that the 16S distribution on leaves strongly clustered by treatment (p-value = 0.001), while no significant pattern (p-value > 0.05) was found in relation to the sampling collection time (Fig. 2c). The opposite trend was apparent for the fungal community. In fact, although ITS distribution does not appear to be affected by the treatment (p-value > 0.05) as displayed in Fig. 2b, it does strongly depend on the period of collection during the season, not only at month level (p-value = 0.001) as shown in Fig. 2d, but sometimes even at day level (p-value = 0.001), as presented in Fig. S2.

Figure 2
figure 2

PCoA plots displaying the beta diversity within the dataset for 16S (left column) and ITS (right column) with samples coloured by treatment (above) and period of collection (below). (a) Distribution of 16S data based on treatment, (b) distribution of ITS dataset grouped by treatment, (c) distribution of 16S dataset coloured by period, and (d) distribution of ITS coloured by period. All distance matrixes are calculated based on Bray-Curtis dissimilarity.

Taxonomical composition on leaves

The results below show the microbial composition in terms of taxa assigned to the different features found, highlighting the microbial representation on grapevine leaves. All the information regarding bacterial and fungal community are summarised in Fig. 3. Figure 3A shows the taxa bar plots of the bacterial population. In order, the dominant bacterial taxa inferred from the sequence features were Lactobacillaceae, Bacillus, Oxalobacteraceae, Enterobacteriaceae, Planococcaceae, Pseudomonas, Enterococcus, Sphingomonas and Staphylococcus. The main notable difference between the biocontrol-treated samples and the Cu-treated leaves was in the variation of Lactobacillaceae. This family includes the potential biocontrol agent and is an important indicator of its presence on leaves. Although members of Lactobacillaceae are generally present on the leaves, the difference between the biocontrol and Cu-treated samples was statistically significant (based on ANCOM). When the exact sequence variant assigned to this taxa was identified, an alignment was performed against the NCBI database using BLAST. The result of this confirmed the sequence belong to Lactobacillus plantarum. The presence of Lactobacillaceae on treated samples varied from 42% to 68% of the total bacterial community in the period from September to December. In the post-harvest months (February and March), the relative abundance decreased to between 0.1% and 3.8%, resembling the values found in the control/Cu-treated samples. This means that three months after the final spray (December) the relative abundance of Lactobacillaceae resembled the levels found on copper-treated leaves.

Figure 3
figure 3

Taxonomical bar plots for the bacterial and fungal community. (A) from left to right: taxa bar plots in which all the 16S samples are grouped by treatment, taxa bar plot representing t0 for both treatments and taxa bar plot in which samples are grouped by treatment during a specific period. Finally the legend representing the 18 most abundant bacteria. (B) from left to right: taxa bar plots in which the ITS samples are divided by treatment, taxa bar plot representing t0 for both treatments and taxa bar plot in which samples are grouped by treatment within a specific period. Finally the legend representing the 18 most abundant fungi.

The fungal communities, displayed in Fig. 3B, were relatively similar for the biocontrol and Cu-treated samples when looking at the entire growing season. The dominant taxa belonged to Pleosporaceae, Cladosporium, Alternaria, Capnodiales, Sporobolomyces and Aureobasidium pullulans. Although no differences appeared between the two treatments, another pattern was seen with respect to the collection time. In fact, during the growing season, there was a significant change in the taxonomical composition of the fungal communities. Several fungal outbreaks led to a variation in relative abundances during the period studied. For instance, Pleosporacea members tended to increase in relative abundance from September to December, reaching 55% of the total community, before suddenly decreasing to 9% in March, although they still appeared among the dominant taxa. Sporobolomyces, Sporidiobolus and Botryosphaeriaceae were more dominant in September, while Cladosporium and Alternaria were mostly represented in the post-harvest months. In light of this finding, it was of interest to examine the differences between the treatments just within the same sampling period. Interestingly Botryosphaeriaceae varied in relative abundance from 18% in the Cu-treated leaves to 0.67% to the biocontrol in September. Further analyses were performed using BLAST on the dominant sequence variant assigned to Botryosphaeriaceae and the read was assigned at species level to Diplodia seriata, a known pathogen that causes bot canker on grapevine. Conversely Leptosphaeriaceae detected in September ranged from 2.7% in Cu-treated leaves to 12% recorded in biocontrol-treated samples for the same period. In both cases this sudden change in relative abundance could be ascribed to the non-homogeneous distribution of these taxa between the two blocks at time zero, when no sprays were performed (Fig. 3) The known targets (in-vitro) of MW-1 Botrytis cinerea, Aspergillus niger, and Penicillium brevicompactum were undetected in both treatment and their respective higher taxonomical rank (family) were representing less than 1% of the total relative abundance.

ANCOM and Gneiss on leaves

To test which taxa change in a statistically relevant way, ANCOM was run using treatment and period as discriminants for these samples. The result was that for 16S between different treatments, the only taxa that changed in a statistically significant way was Lactobacillaceae, the family to which this study’s biocontrol agent belongs, with the main sequence variant assigned to Lactobacillus plantarum (Table S1). Only two statistically relevant differences were found when looking at the period regarding mitochondria and Ralstonia, but none were associated with grapevine and furthermore they appeared in low abundance (Table S2).

In the fungal community there were only two taxa that changed significantly between treatments. These taxa belong to the species Kondoa aeria and to the order of Filobasidiales (Table S3). They were in very low abundances and did not seem to be related to any known disease or have an impact on the plant itself. Instead, for this period a large amount of taxa were observed that changed during the season, 42 of which are classified at least at class level (Table S4). These findings provide information that can be used to understand the evolution of the fungal community as the season progresses. These 42 taxa included potential pathogens such as Alternaria, Cladosporium or members of Botryosphaeriaceae, with the main sequence assigned to Diplodia seriata using BLAST. As the period shaped the fungal community, an investigation was carried out on which taxa changed between different treatments during specific periods in the season. This ANCOM test returned 50 taxa, of which 46 were classified at order level (Table S5). Among the taxa that appeared to be differentially distributed between treatments, in a short period of time (such as September alone), there were Botryosphaeriaceae and Leptosphaeriaceae, which were previously highlighted when looking at the taxonomical composition in Fig. 3B.

Gneiss was then run on the bacterial community dataset to evaluate the impact of MW-1 on other bacteria. Interestingly, Lactobacillacea appeared to be sensitively different (p-value < 0.05) from the rest of the microbial community, based on the position occupied in the tree (Fig. S3), but the fact that this taxon branches out from the rest means that it does not interact with other species in the community within the same microbial niche. This is a further indication, which still requires biological confirmation, that the introduction of this potential biocontrol agent may not alter the overall bacterial composition normally detected on leaves sprayed with copper.

Quantifying fungal and MW-1 biomass on leaves

The results of the quantitative evaluation by qPCR of the microbial community on the samples are summarised in Figs. 4 and S4. The detection of MW-1 was determined by strain-specific qPCR combined with a high-resolution melting curve. Furthermore, the total relative abundance of fungi during the period and between the treatments was compared using universal primers for fungi (with the same primer set used for the generation of the sequencing reads to reduce biases). MW-1 was detected in all the treated samples after the first spray (26 September) and up to 2 February, one month after harvest. The cell numbers showed a range spanning from 103 to 105 cells/leaf. The maximum amount of Lactobacillus plantarum MW-1 cells was detected after the spray on 15 November, while on two dates the number of cells detected was below the detection limit. The first sample where MW-1 was not detected correspond to our t0 (15th of September), when no spray were previously done. The second sample in which MW-1 was not detected correspond to the 9th of March, two months after harvest. This gives an indication about the persistence of MW-1 on leaves.

Figure 4
figure 4

(a) Positive correlation between index-ratios IR1 and IR2; letters identify the collection date as reported in Table 1. (b) Normalised genomes count for MW-1 (light blue triangles), estimated fungal genomes on community treated with biocontrol agent (dark blue squares) and fungal community genomes count under Cu treatment (black diamonds); curves are expressed as log10 of the genome count. The red dashed line indicates the harvest.

The total fungal abundance could be estimated at cell level number without introducing significant biases due to an uneven ploidy variation between the different taxa detected. However, since the community composition between the treatments appeared to be very similar, it was assumed to have a negligible difference in terms of genome distribution. For this reason, a relative quantification was performed between the two fungal communities grouped by treatments. With this in mind, the results obtained showed that the fungal communities between the two treatments followed a similar trend, with no profound variations between the different treatments. It was decided to estimate the number of fungal genomes, assuming a fungal-genome average size of 27.5 Mb. Accordingly, the number of genomes ranging from 102 to 104 genomes/leaves during this period was calculated (Fig. 4b). After normalising the raw-qPCR signals using the gDNA amount, an ANOVA was run and return a p-value of 0.56. This confirmed that the variation between the two fungal communities coming from the two treatments was not significant. Finally a correlation analysis was performed between the number of cells of MW-1 detected and the corresponding fungal genome amount for each sample. In this case, the analyses returned a coefficient of −0.39, which could be interpreted as a moderate negative correlation between fungi and MW-1.

To further demonstrate the link between initial gDNA amount and qPCR signal, an additional correlation analysis was performed. First two index ratios (IR) were calculated: IR1 between input gDNA in Cu versus biocontrol-treated leaves, and IR2 by dividing qPCR genome-copies detected on Cu versus biocontrol-treated leaves. Finally, the correlation between IR1 and IR2 during the period was determined. The correlation coefficient was 0.81, confirming that the two IRs were highly positively correlated, as shown in Fig. 4a. A t-test on the two IR series (the assumption of equal variance was tested with a F-test that resulted in a p-value of 0.165) confirmed that the difference was not significant (p-value = 0.699). This means that the tendency of the two IRs did not differ significantly between the treatments. In conclusion, this indicated that not only did the two fungal communities not differ significantly in composition between treatments, but they were also quantitatively comparable.

Discussion

Understanding the structure and dynamics of microbial communities in the phyllosphere is important because of their effects on plant protection and plant growth. Since the structure and dynamics of the microbiome are affected by environmental factors such as geography30, grape cultivar31 and climate32, this study focused on different treatments applied in two blocks of the same vineyard. Both treatments, Cu and MW-1, were applied to control fungal diseases. The treatments were applied throughout the growing season and were followed by the periodic collection of samples. The number of sprays performed was decided based on the visible symptom evaluation on the plants within the two blocks during the season. This is important to note since the impacts of any chemical and biological treatments depend on the dosage and frequency of applications, as well as other biotic and abiotic factors32.

This study revealed the development of the grapevine leaf microbiome over the growing season and displayed the impact of different types of vineyard management practices using NGS technology combined with qPCR. The combined use of these approaches allowed quantitative and qualitative information to be recovered from the compositional dataset32.

In the present study, the biocontrol-treated phyllosphere-bacterial community had a reduced phylogenetic diversity and evenness compared with Cu-treated leaves, seemingly due to the presence of highly dominant taxa (i.e. MW-1, which appears as Lactobacillaceae, being dominant in almost all the treated samples). These highly abundant taxa probably compact the remaining populations, reducing evenness and phylogenetic diversity, and excluding some of the rare species that are likely to be missed at a fixed sequencing depth. In fact, phylogenetic diversity positively correlates with sequencing depth, as shown by33. The beta diversity of the bacteria displayed significant differences since both communities, by treatment, scattered in two distinct directions along axis 1, which explained 28.5% of the variation.

In contrast to the bacterial communities, the fungal communities were mainly affected by the sampling period, during which PD and evenness showed similar trends for both the biocontrol and Cu-treated leaves. This is in line with Singh et al.31 who showed that the shift of the microbial community in samples collected at different times was detected for fungi only. The fact that phylogenetic diversity increased, together with evenness, until December suggested that the complexity of the fungal community increased due to a higher number of similarly abundant taxa. In December, the fall in PD and evenness was probably due to the presence of one or a few taxa with a high relative abundance that replaced rare species following the same mechanism highlighted for the bacterial community. This hypothesis is based on the fact that the sequencing depth was comparable between samples. The beta diversity of the fungal communities showed that the two treatments had an insignificant effect on the mycobiome when the whole period was considered. This result is consistent with the study of 32 where a chemical and a biocontrol agent against downy mildew were tested in different locations and climatic conditions. They showed that the fungal community was impacted mostly by geography and was resilient to the treatments tested. Instead, it is clear from the present study that the fungal community was strongly affected by the sampling time, suggesting an evolution of the harboured microbiota due to seasonal changes such as temperature, solar exposition, rainfall and other abiotic fractions. The fact that microbial seasonal shift only impacts the fungal community and not the bacterial community is consistent with the work of 31 in which two sets of samples from the same vineyards were analysed and the impact of different sampling times on the microbial community was measured.

The differences in beta diversity reflected the change in taxonomical composition in both bacterial and fungal communities. In particular, the bacterial communities differed significantly only for the relative abundance of Lactobacillaceae, as shown in Table S1, obtained by applying ANCOM. This taxon appeared in the biocontrol-treated leaves in relative abundance up to 40% and contained reads from MW-1. After blasting the exact sequence variant against NCBI, the analyses could be refined to species level, identifying Lactobacillus plantarum as the most likely assignment. The present potential biocontrol strain, which here appeared as the dominant sequence variant, was not the only representative of Lactobacillaceae. The presence of other sequence variants assigned to Lactobacillaceae was also retrieved on copper-treated samples. This is consistent with the existing literature in which a community of LAB is found living on the phyllosphere naturally34. Furthermore, the strain-specific qPCR applied to recognise and quantify MW-1 produced a precise signal peak only on the biocontrol-treated samples. This allowed the conclusion to be drawn that the different amount detected in relative abundance for the Lactobacillaceae family was due to the presence of reads belonging to MW-1. Finally the bacterial community, which appeared to be relatively stable throughout the season, reflected the composition already reported by31 and35. The dominance of Proteobacteria and Firmicutes highlighted in their paper, was consistent with the phyla-distribution shown in the present study. In our study though, the most dominant phyla is Firmicutes instead of Proteobacteria. This difference is probably due to the different primers in use. Furthermore31, analysed two different sets of samples, one collected in spring and one during the harvest period, which displayed similar bacterial community composition with a few differences in relative abundance, e.g. the taxa Cyanobacteria. In addition to35, we retrieved also Actinobacteria and on leaves. At genus and species level the population is also reporting some similarity with35 where Sphingomonas, Pseudomonas are the most represented taxa which also are reported to influence plant productivity and health36,37,38. Furthermore, our study confirm the presence on leaves of Agrobacterium, Staphylococcus and Burkholderia which were commonly reported as endophyte bacteria but likely to be found on the phyllosphere as well35,39. Other dominant bacterial families are represented by Enterobacteriaceae, Oxalobacteraceae and Enterococcaceae confirmed by40. Among these, Enterobacteriaceae was very abundant at the t0 within the copper treated block while it was barely detected in the biocontrol treated; instead, t0 within the biocontrol treated block was dominated Pseudomonas and Sphingomonas. Finally the family Lactobacillaceae which contains reads of MW-1, is reported increasing also on leaves treated with Cu during the season. This is also reported in40 even though the period was different but still corresponding to the same phenological state of our plants.

The analysis of the fungal community throughout the season, offered interesting insight into the possible correlations between fungal outbreaks and seasonal variations. When looking at the whole period, there was a clear succession of dominant fungal taxa with Sporobolomyces in September, Pleosporaceae in December and Cladosporium, Alternaria and Aureobasidium pullulans in February and March. All of these taxa are listed as being commonly recovered on grapevine plants in the review of41. However, focusing only on the differences between treatments, that appeared only at a specific period of the season, can hide the presence of other taxa which could be sensitive to the treatments but are detected only in a short amount of time. For instance, in September, just before spraying started, there were differences in the dominant taxa between the two blocks which had been treated in two different ways. The block subject to biocontrol sprays showed a reduced relative abundance of Botryosphaeriaceae compared to the Cu-treated block; the dominant sequence variant was assigned to Diplodia seriata, a known grapevine pathogen that can cause bot canker42,43. Conversely, within the family Botryosphaeriaceae, there is also the genus Botryosphaeria, whose anamorph stage has been confused with Diplodia44,45. This genus appears not to be strongly affected by Cu46 but from this study its presence decrease, in relative abundance, within few weeks. Botryosphaeriaceae includes known grape-pathogen and it has been reported mainly in wood samples as it is connected with the incidence of several grape trunk diseases (GTD)47. This fungal taxa has rarely been reported on leaves; in fact only Pinto et al. in 2014 reported a low relative abundance associate to this taxa (0.23%). In this case instead, the relative abundance is much higher but only for a short range of time. Another observation, which can be made focusing on the post-harvest period, reveal that to a rise in relative abundance of Aureobasidium pullulans correspond a decrease in Pleosporaceae; this family include the genus Alternaria which was reported to be negatively affected by the presence of Aureobasidium pullulans by40. However, other reads were assigned to Alternaria sp. but here the negative correlation is not consistent. These outbreaks support the hypothesis that fungal community is more dynamic during time compared with bacterial. At last, a 3% of relative abundance was attributed to Erysiphaceae, a family that contains phytopathogens such as powdery-mildew agent Erysiphe necator (syn. Uncinula necator) and appeared only in September, before and right after the first spray in both treatments. However, targeted quantitative methods, such as qPCR, are required to confirm or deny these variations.

Finally, the outcome of the NGS approach was combined with a targeted strain-specific qPCR, to detect and quantify MW-1, and with qPCR applied on the fungal communities for the two treatments. This combined approach has previously been reported32 as an efficient way of drawing conclusions about microbial ecology when biocontrol agents are involved. The number of cells detected for MW-1 was consistent with the existing literature in which biocontrol agents on the phyllosphere are studied and quantified, and could eventually be used as an indicator of the effectiveness of the treatment32. The detection limit in this study was 102 cells/g of leaves, a frequent limit of detection when qPCR is used48. The number of cells detected during the season varied between 103 and 105. A sudden tenfold decrease was expected 1–6 days after the spray, as previously reported in other field studies such as34 using L. plantarum PM411 on apple, kiwi, strawberry and pear plants. Comparing culture-dependent results, based on isolation of LAB in biocontrol studies such as20, our findings are also confirmed. In fact, Tsuda and colleagues reported 105–104 CFU/g detected right after foliar spraying, followed by fast decrease in the next 7 days. When investigating the amount of fungi, the number of cells was not estimated due to the biases that could have been introduced by different ploidy between species. However, an average amount of fungal genomes was estimated based on the information reported in49, with the average size for ascomycetes and basidiomycetes estimated at 13 and 42 Mb respectively. Accordingly, since members of both phyla were retrieved, 27.5 Mb was used as the average size for genomes. However, the results could be slightly different if these parameters were changed, according to other studies on fungal genomes such as50, which reported different genome sizes. All the analyses established that the two fungal communities were quantitatively similar. Finally, when comparing the two normalised qPCR signals, an ANOVA test confirmed that there were no significant differences between treatments.

Based on these results, leaves treated with L. plantarum MW-1 showed minor differences at specific periods of the year in terms of taxonomical composition, and there were no significant differences in terms of fungal load when compared with Cu-treated leaves. In both treatments, the microbes that were considered targets for MW-1 display a comparable low relative abundance, which might be due to either a comparable antagonistic effect of Cu and MW-1 or, in a more cautious hypothesis, to the low disease pressure occurred in that specific year. In light of this, further evidences are required before considering Lactobacillus plantarum MW-1 as a potential biocontrol alternative or supplement to the use of Cu on grapevine leaves against Botrytis cinerea, Aspergillus niger and Penicillium brevicompactum. Notwithstanding this and to conclude, the methodology applied in this study allowed to underline the dynamics of the microbial community on grapevine leaves, throughout the growing season. This raise awareness on the best timing to apply treatments to the plants, before the development of symptoms, within the framework of a more sustainable precision agricultural practices.

Material and Methods

Materials

Leaves were collected periodically from September 2016 to March 2017 from a commercial vineyard in ine Western Cape (South Africa). Grape cultivar was Shiraz, grafted with US8-7 rootstock. Two plots (each of them of about two hectars) in the same field, separated by 5 rows of grapevine, were managed in two different ways: one block was sprayed with Cu in a conventional vineyard management system (6 kg Cu ha−1), while the other was sprayed with a Lactobacillus plantarum MW-1 as a potential biocontrol agent. A third plot was initially sampled as no-treated control but in few weeks, after the appearance of symptoms of downy mildew and Botrytis cinerea, the winery decided to spray this plot with copper sulfate following the same scheme of the Cu-treated block. For this reason the control-plot was excluded from the analyses. The Lactobacillus plantarum MW-1 strain has previously been isolated from grapes in the same wine region (data not shown). In laboratory, previous tests highlighted an effect of growth inhibition against Botrytis cinerea, Aspergillus niger and Penicillium brevicompactum. The volume of MW-1 sprayed vary from 350 L/Ha in September to 1000 L/Ha at the end of October and remained stable til December. During the studied period, 12 different samplings were carried out; leaves were collected in five spots for each of the treatment areas within the vineyards. For each spot and each treatment there were three independent biological replicates. Each sample consisted of five leaves, picked from mixed position along the shoot, that were placed in a 50 ml Falcon tube. The tubes were frozen at −20 °C and shipped to the laboratory in Denmark. Sample collection followed the spray calendar for the biocontrol and copper treatment. After the grape harvest in January, there was a follow-up with two extra samplings in February and March. All the plants from which leaves had been collected were marked and did not change throughout the experimental plan. A complete calendar of the sampling is given in Table 1. A visual evaluation was carried out between sprays to detect the onset of symptoms and in both blocks the leaves were healthy.

Table 1 Sampling dates in this study; second column report the dates of MW-1 sprays, third column the code used during data-processing and analyses; fourth column identify the phenological state of the vines at the time of the samples collection; fifth column contains notes about the sampling moment.

Sample preparation and DNA extraction

The tubes containing the leaves were thawed at room temperature and 20 ml of a washing solution was subsequently added (Singh et al.31) and placed on a rocket inverter, applying one hour of shaking rotation. After this, the leaves were discarded using sterile pincettes. The washing solution was centrifuged for 15 minutes at 6000 rpm to create a pellet. Following removal of the washing solution, the pellet was resuspended in 978 ul of phosphate buffer solution and DNA was extracted according to the manufacturer’s protocol using FAST DNA Spin Kit for Soil (MP-Biochemical, CA). After the extraction, the DNA quality and concentration were measured with Nanodrop (Thermo Fisher Scientific) and Qubit 2.0 fluorometer (Thermo Scientific).

Library preparation for sequencing

Amplicon library preparation for 16S bacterial gene was performed as described by51, with minor modifications hereby reported. To reduce the amount of plastidial DNA from the grapevine leaves, specific mPNA and pPNA were used with the same sequences as those recommended by52. Each reaction contained 12 µL of AccuPrime SuperMix II (Thermo Scientific), 0.5 µL of forward and reverse primer from a 10 µM stock, 0.625 ul of pPNA and mPNA to a final concentration of 0.25 mg/mL, 1.5 µL of sterile water, and 5 µL of template. The reaction mixture was pre-incubated at 95 °C for 2 min, followed by 33 cycles of 95 °C for 15 sec, 75 °C for 10 sec, 55 °C for 15 sec and 68 °C for 40 sec. A further extension was performed at 68 °C for 10 min.

The fungal community was sequenced using the same double-step PCR approach for library preparation described in53 and 54, but the primers in use were ITS1 (TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-GAACCWGCGGARGGATCA) and ITS2 (GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-GCTGCGTTCTTCATCGATGC) with adapters for Illumina MiSeq Sequencing. Each reaction for the first PCR on ITS contained 12 µL of AccuPrime SuperMix II (Thermo Scientific), 0.5 µL of forward and reverse primer from a 10 µM stock, 0.5 µL of bovine serum albumin (BSA) to a final concentration of 0.025 mg/mL, 1.5 µL of sterile water, and 5 µL of template. PCR cycles were the same as above, but without the step at 75 °C for 10 sec. The number of cycles was set to 40. The second PCR of fragment purification using beads and Qubit quantification were common for 16S and ITS and followed the protocol reported in51. The final pooling was proportional to the length of the fragment in order to sequence an equimolar amount of gene fragments for all the samples. Thus by sequencing all the samples within the same sequencing run, biases due to run variations were minimised. Sequencing was performed using an in-house Illumina MiSeq instrument and 2 × 250 paired-end reads with V2 Chemistry.

Strain-specific primer for MW-1

In order to find a region that was unique for MW-1, part of the genome was uploaded and scanned using PHASTER (PHAge Search Tool Enhanced Research)55 to detect phage and prophage genomes within the bacterial genome. Following this, the uniqueness was tested of an amplicon generated from one primer that targets the genome and another targeting the phage sequence against the NCBI database and a private genome database repository of the strain provider. Thus a fragment of 300 bp was obtained that did not show any hit in NCBI and only one full hit in the private genome repository. To amplify this fragment, the primers MWP3F (CATCCCAACCGCTAACAA) and MWP3R (CGCAGAAAAGGTAGCAAA) were used. These primers were tested further against other strains of L. plantarum from the authors’ collection, showing amplification only toward MW-1. To create a standard curve, the DNA was extracted from a pure culture of MW-1 and serial dilutions were performed from 10−2 to 10−7. Knowing the length of the fragment of 300 bp, the qPCR signal and the fact that the fragment appears only once in a genome, the number of cells of L. plantarum MW-1 was estimated for each reaction. This information was used to estimate the number of cells in leaf extract samples.

qPCR

All PCR reactions were prepared using UV sterilised equipment and negative controls (NTC) were run alongside the samples. The qPCR with specific primers for MW-1 and the fungal community was carried out on a CFX Connect Real-Time PCR Detection System (Bio-Rad).

The primers used for the bacterial agent were MWP3F and MWP3R, while for the total fungal community the same primers ITS1 and ITS2 were used as during library preparation. To reduce biases the adapters for Illumina MiSeq Sequencing were also maintained in qPCR. Single qPCR reactions contained 4 µL of 5x HOT FIREPol EvaGreen qPCR Supermix (Solis BioDyne, Tartu, Estonia), 0.4 µL of forward and reverse primers (10 µM), 2 µL of bovine serum albumin (BSA) to a final concentration of 0.1 mg/mL, 12.2 µL of PCR grade sterile water, and 1 µL of template DNA. Since DNA concentration in some of the samples was very low (around 0.5 ng/ul), normalisation by diluting was avoided before qPCR, but the qPCR signal have been normalized using the total input DNA, previously measured by Qubit. The qPCR cycling conditions included initial denaturation at 95 °C for 12 min, followed by 40 cycles of denaturation at 95 °C for 15 sec, annealing at 56 °C for 30 sec, and extension at 72 °C for 30 sec; a final extension was performed at 72 °C for 3 min. Quantification parameters showed an efficiency of E = 87.7% and R2 of 0.999. The same qPCR cycle conditions were applied for both bacterial and fungal qPCR. All qPCR reactions were followed by a dissociation curve in which temperature was increased from 72 °C to 95 °C, rising by 0.5 °C in each cycle, and fluorescence measured after each increment. Since the initial concentration of DNA was not homogeneous and in a few samples was below 1 ng/ul, it was decided to normalise the qPCR result using the initial concentration of total DNA, as reported in56.

Bioinformatics

Sequencing data were analysed and visualised using QIIME 2 v. 2018.1157 and the same pipeline described in51. After demultiplex, reads were processed using DADA2 denoise.paired plugin58. For both 16S and ITS data a multiple-sequence alignment was performed with MAFFT59 and the following phylogenetic tree was built using FastTree60. Diversity analyses were performed using q2-diversity plugin (https://github.com/qiime2/q2-diversity). Alpha-diversity was measured in terms of Phylogenetic Diversity (PD)61 while evenness was reported using Pielou score62. Beta Diversity was measured using PERMANOVA and the PCoA plots were visualized using Emperor63. Taxonomic assignments were done using qiime feature-classifier classify-sklearn with a pre-trained Naïve-Bayes classifier64 with Greengenes v_13.865 for 16S and UNITE v7.2 for ITS66. All the high-quality reads assigned to mitochondria and chloroplast were removed before the analyses of diversity. The raw data of this study are available on SRA (Study Accession Number PRJEB30496).

Statistics

Statistical evaluation of these results was performed separately for qPCR and the sequencing dataset. All statistical evaluations and visualisations regarding the qPCR dataset were performed on Microsoft Office Excel 2010, while the NGS dataset was analysed through QIIME2 v. 2018.11. The qPCR signals obtained allowed the calculation of the number of cells of MW-1 detected on the sprayed leaves throughout the period, while for the fungal community it was possible to perform a relative quantification between the two different treatments of Cu and biocontrol. For both dataset the signal obtained was normalised using the total DNA as input56.

To relate the differences in qPCR signal on the fungal community to the presence/absence of the biocontrol treatment, the qPCR signal had to be correlated with the initial DNA amount before checking the differences and their statistical relevance. To do this, two different index ratios (IR) were created: IR1 dividing the raw qPCR signal from Cu samples from that from biocontrol treated samples, and IR2 dividing the initial amount of DNA of Cu-treated with that of biocontrol-treated leaves. Two series of IRs were then obtained to test for equal variance with a F-test. This was intended to establish whether the two IRs, deriving from two different datasets, had a comparable variance within their samples. The F-test for variances, with a p-value = 0.16, proved that the two variances were equal and then a t-test was run on two series assuming equal variances between the IRs on the two treatments, Cu and biocontrol. A correlation analysis was also performed between the initial DNA amount and the resulting qPCR signal. Finally, correlation analyses and a t-test were repeated on the samples treated with the bacterial agent to compare the variation between the total fungal community on biocontrol-treated leaves and the number of cells detected of the biocontrol bacteria. Sequencing data, after QIIME 2 pipeline processing, were statistically evaluated using the Kruskal-Wallis test for alpha diversity. The resulting p-value of the alpha diversity comparison was based on two parameters (phylogenetic diversity and evenness) calculated between the different series analysed. Beta-diversity analyses were evaluated using PERMANOVA with 999 permutations. Finally a statistical evaluation of differentially abundant features was performed based on analysis of composition of microbiomes (ANCOM)67. This test is based on the assumption that few features change in a statistical way between the samples and hence is very conservative. To individuate microbial niches, Gneiss68 was also applied on the bacterial dataset, using treatment as the determining variable.